Converting a tutorial to Python

EnnaD

Hi !
I am trying to run the Subglacial Channel Formation from a Single Moulin (SHaKTI) tutorial with Python. As I could not find Python input files for it I tried translating the matlab ones. I believe I was successful for steps 1 and 2 but I don't think I am setting-up the moulin properly. Or maybe I am doing a bunch of different things poorly !
I wanted to attach my .py translation of it here but don't see how I can do that so I'll paste it in a second post. If someone could help me identify where I made mistakes that would be very helpful.

Many thanks !
Anne

EnnaD

import sys                                                                                                                                            
import numpy as np
sys.path.append('blabla/0_CODES/trunk/bin')
sys.path.append('blabla/0_CODES/trunk/lib')
from issmversion import issmversion
from model import *
from plotmodel import *
from triangle import *
from parameterize import *
from setmask import *
from hydrologyshakti import *
from transient import *
from solve import *

model = model()

--STEP 1
print("  Step 1: Mesh")
md=triangle(model,'./outline.exp',20);
plotmodel(md,"mesh", 'data','mesh')

--STEP 2
print("  Step 2: Parameterization")
md=setmask(md,'','');
md=parameterize(md,'moulin.par')
--HYDROLOGY SPECIFIC PARAMETERIZATION:
md.hydrology=hydrologyshakti()
--Define initial water head such that water pressure is 50% of ice overburden P
md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base
--Initial subglacial gap height of 0.01m everywhere
md.hydrology.gap_height = 0.01*np.ones((md.mesh.numberofelements,1))
--Typical bed bump bump spacing (2m)
md.hydrology.bump_spacing = 2*np.ones((md.mesh.numberofelements,1))
--Typical bed bump height (0.1m) 
md.hydrology.bump_height = 0.1*np.ones((md.mesh.numberofelements,1))
--Define distributed englacial input to the subglacial system (m/yr)
md.hydrology.englacial_input = 0.0*np.ones((md.mesh.numberofvertices,1))
--Initial Reynolds number (start at Re=1000 everywhere)
md.hydrology.reynolds= 1000*np.ones((md.mesh.numberofelements,1))
--Set up atmospheric pressure Type I boundary condition at left edge of domain (outflow, i.e. h=zb at x=xmin) 
md.hydrology.spchead = float('NaN') * np.ones((md.mesh.numberofvertices,1))
pos = np.nonzero(np.logical_and(md.mesh.vertexonboundary, (md.mesh.x==min(md.mesh.x))))[0]
for i in range(len(pos)):
    md.hydrology.spchead[i]=md.geometry.base[i]
plotmodel(md,"init", 'data', md.geometry.base, 'title', 'Bed elevation in m',
                     'data', md.geometry.surface, 'title', 'Surface elevation in m',
                     'data', md.hydrology.head, 'title', 'Head in m',
                     'data', md.hydrology.gap_height, 'title', 'Gap height in m')

--STEP 3
print("  Step 3: Solve!")
md.transient.deactivateall()
md.transient.ishydrology=1
--Specify that you want to run the model on your current computer
--Change the number of processors according to your machine
md.cluster=generic('np',1);
--Time stepping
md.timestepping.time_step=3600./md.constants.yts
md.timestepping.final_time=30.0/365.0
num_ts = int(md.timestepping.final_time/md.timestepping.time_step)
print("final time ?", md.timestepping.final_time)
print("time stepping ?", md.timestepping.time_step)
print("number of ts ?", num_ts)
--Add one moulin with steady input at x=500, y=500 
time = []
curr_time = 0.0
for i in range(num_ts):
    time.append(curr_time)
    curr_time += md.timestepping.time_step
md.hydrology.moulin_input = np.zeros((md.mesh.numberofvertices+1, num_ts))
md.hydrology.moulin_input[-1,:] = time
pos = min(np.sqrt((md.mesh.x - 500.)**2+(md.mesh.y - 500.)**2))
md.hydrology.moulin_input[pos,:] = 4
--Specify no-flux Type 2 boundary conditions on all edges except for outflow
md.hydrology.neumannflux=np.zeros((md.mesh.numberofelements+1, num_ts))
md.hydrology.neumannflux[-1,:] = time
md.verbose.solution=1
md=solve(md,'Transient')
plotmodel(md,"finalSol", 'data', md.results.TransientSolution[-1].EffectivePressure, 'title', 'Pressure effective in Pa',
                     'data', md.results.TransientSolution[-1].HydrologyHead, 'title', 'Head in m',
                     'data', md.results.TransientSolution[-1].HydrologyBasalFlux, 'title', 'Basal water flux in m2s-1',
                     'data', md.results.TransientSolution[-1].HydrologyGapHeight, 'title', 'Gap height in m ')

bdefleurian

Hi Anne.
It would help if you could post your error message here too. That would help pinpoint the issue.

So a few issues with your file, I am writing those as I read so they may not be related to your problem specificaly.

You should avoid to append paths in the test, do that in you profile.
You can remove model = model() and have md = triangle(model(), ... instead
We use .py for the parameters file in python, There may be issues in this file but you have not provided .
All the vectors in python are vectors, so you don't need the ,1 in their definition, np.ones((md.mesh.numberofelements,1)) should be np.ones((md.mesh.numberofelements))
You don't need to iterate to fix the spc, you can have md.hydrology.spchead[pos]=md.geometry.base[pos]
curr_time can be written : curr_time = np.arange(md.timestepping.start_time, md.timestepping.final_time + md.timestepping.time_step, md.timestepping.time_step)

EnnaD pos = min(np.sqrt((md.mesh.x - 500.)2+(md.mesh.y - 500.)2))

is not giving you a position just the minimum distance to the (500, 500) point, you should use something like that :pos = np.where((np.sqrt((md.mesh.x - 500.)**2+(md.mesh.y - 500.)**2))==np.nanmin(np.sqrt((md.mesh.x - 500.)**2+(md.mesh.y - 500.)**2)))
I expect that this last one is your issue.
Also we are interested into the translation when you will be done. If you could send the files my way I would have them on our repository for others to use.

EnnaD

And the moulin.par:

import numpy as np                                                                                                                                    
from frictionshakti import *
from verbose import *
from SetIceSheetBC import *
#Start defining model parameters here

# Set up bed topography and ice geometry for a tilted 500m thick slab
md.geometry.base = .02*md.mesh.x
md.geometry.bed = md.geometry.base
md.geometry.surface = .02*md.mesh.x + 500 
md.geometry.thickness = md.geometry.surface - md.geometry.bed

# Define ice sliding velocity (m/yr) 
md.initialization.vx = 10e-6*md.constants.yts*np.ones((md.mesh.numberofvertices,1))
md.initialization.vy = np.zeros((md.mesh.numberofvertices,1))
md.initialization.pressure=np.zeros((md.mesh.numberofvertices,1))

# Materials
# Ice flow law parameter (note that the standard parameter A=B^(-3))
md.materials.rheology_B= (5e-25)**(-1/3)*np.ones((md.mesh.numberofvertices,1))
md.initialization.temperature=(273)*np.ones((md.mesh.numberofvertices,1))
md.materials.rheology_n=3.*np.ones((md.mesh.numberofelements,1))

#Calving
md.calving.calvingrate=np.zeros((md.mesh.numberofvertices,1))

# Friction law and coefficient
md.friction=frictionshakti(md)
md.friction.coefficient=20.*np.ones((md.mesh.numberofvertices,1))

#Numerical parameters
#md.stressbalance.viscosity_overshoot=0.0;
md.masstransport.stabilization=1.
md.thermal.stabilization=1.
md.verbose=verbose(0)
md.settings.waitonlock=30
md.stressbalance.restol=0.05
md.steadystate.reltol=0.05
md.stressbalance.reltol=0.05
md.stressbalance.abstol=float('NaN')
md.timestepping.time_step=1.
md.timestepping.final_time=3.

#GIA: 
md.gia.lithosphere_thickness=100.*np.ones((md.mesh.numberofvertices,1)) # in km
md.gia.mantle_viscosity=1.0*10**21*np.ones((md.mesh.numberofvertices,1)) # in Pa.s
md.materials.lithosphere_shear_modulus=6.7*10**10                   # in Pa
md.materials.lithosphere_density=3.32                              # in g/cm^-3
md.materials.mantle_shear_modulus=1.45*10**11                       # in Pa
md.materials.mantle_density=3.34                                   # in g/cm^-3

#Boundary conditions:
md=SetIceSheetBC(md)

#Change name so that no test have the same name
#A=dbstack
#if (length(A)>2):
#    md.miscellaneous.name=A(3).file(1:end-2)

Thanks

EnnaD

Thanks a lot for your input and help.
So to clarify, the simulation was not crashing, it was just producing garbage results, if I compare to what is provided on the website here.
I am unsure of what you mean by

bdefleurian We use .py for the parameters file in python, There may be issues in this file but you have not provided

As I thought my third post was providing just that ?

With the modifications you suggested, the results seem better in that I now clearly have an input flux in the center of the domain. However the head and water flux seem off by several order of magnitudes...
I wanted to attach a picture of what I get, hoping it would help, but I can't seem to be able to do that so I don't really know how to describe the solution. I am running for 30 days (at least I think I am)

Many thanks for your help again
Anne

EnnaD

Ok only way I could think of was host my img somewhere ... hope you can see it and maybe think of what might be going wrong ? Thanks again.
https://ibb.co/mqr1ykD

bdefleurian

Hei Anne.

Sorry I overlooked your par file, you must have uploaded it while I was answering. So on this one,

as in the runfile you should fix the shape of the vectors (without the ,1)
You messed up the conversion of the exponents, 10^-6 in matlab translates to 1.0e-6 in python
Try to keep the floats as floats and not integers (for example the temperature should be 273. (with the dot) rather than 273
1.0*10**21 can be reduced to 1.e21, same for 6.7*10**10 and 1.45*10**11

The change in velocity might make a big enough difference that you would get the proper results.

EnnaD

Hi, thanks a lot for your input.
I have reviewed my code and fixed these matlab to python conversions -at least I believe I have rectified all of them, but my simulation keeps producing weird results. I don't see a clear channel formation, and the head and water flux are still order of magnitues off from what I see reported on the tutorial webpage. Is there a possibility that the tutorial and paper parametrization vary ?
I am attaching my conversion per your request, maybe you will be able to spot either a mistake or a bug ! At this time I cannot run in matlab unfortunately so it's difficult for me to debug efficiently.
Thanks !
Anne

EnnaD

This is the runme.py

import sys 
import numpy as np
from issmversion import issmversion
from model import *
from plotmodel import *
#from mesh2d import *
from triangle import *
from parameterize import *
from setmask import *
from hydrologyshakti import *
from transient import *
from solve import *


# STEP 1
print("  Step 1: Mesh")
md=triangle(model(),'./outline.exp',20);
## plot
plotmodel(md,"mesh", 'data','mesh')

# STEP 2
print("  Step 2: Parameterization")
md=setmask(md,'','');
md=parameterize(md,'moulin.par.py')
##  HYDROLOGY SPECIFIC PARAMETERIZATION:
md.hydrology=hydrolyshakti()                                                                                                                            
## Define initial water head such that water pressure is 50% 
## of ice overburden P
md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base
## Initial subglacial gap height of 0.01m everywhere
md.hydrology.gap_height = 0.01*np.ones((md.mesh.numberofelements))
## Typical bed bump bump spacing (2m)
md.hydrology.bump_spacing = 2*np.ones((md.mesh.numberofelements))
## Typical bed bump height (0.1m) 
md.hydrology.bump_height = 0.1*np.ones((md.mesh.numberofelements))
## Define distributed englacial input to the subglacial system (m/yr)
## Change the value 0.0 to add distributed input 
md.hydrology.englacial_input = 0.0*np.ones((md.mesh.numberofvertices))
## Initial Reynolds number (start at Re=1000 everywhere)
md.hydrology.reynolds= 1000*np.ones((md.mesh.numberofelements))
## Set up atmospheric pressure Type I boundary condition at left edge
## of domain (outflow, i.e. h=zb at x=xmin) 
md.hydrology.spchead = float('NaN') * np.ones((md.mesh.numberofvertices))
pos = np.nonzero(np.logical_and(md.mesh.vertexonboundary, (md.mesh.x==min(md.mesh.x))))[0]
md.hydrology.spchead[pos]=md.geometry.base[pos]
## plot
plotmodel(md,"init", 'data', md.geometry.base, 'title', 'Bed elevation in m',
                     'data', md.geometry.surface, 'title', 'Surface elevation in m',
                     'data', md.hydrology.head, 'title', 'Head in m',
                     'data', md.hydrology.gap_height, 'title', 'Gap height in m')

#STEP 3
print("  Step 3: Solve!")
md.transient.deactivateall()
md.transient.ishydrology=1
## Specify that you want to run the model on your current computer
## Change the number of processors according to your machine
md.cluster=generic('np',1);
## time stepping
md.timestepping.time_step=3600./md.constants.yts
md.timestepping.final_time=30.0/365.0
time   = np.arange(md.timestepping.start_time, md.timestepping.final_time + md.timestepping.time_step, md.timestepping.time_step)
num_ts = len(time)
print("final time ?", md.timestepping.final_time)
print("time stepping ?", md.timestepping.time_step)
print("number of ts ?", num_ts)
## Add one moulin with steady input at x=500, y=500 
md.hydrology.moulin_input = np.zeros((md.mesh.numberofvertices+1, num_ts))
md.hydrology.moulin_input[-1,:] = time
pos = np.where((np.sqrt((md.mesh.x - 500.)**2+(md.mesh.y - 500.)**2))==np.nanmin(np.sqrt((md.mesh.x - 500.)**2+(md.mesh.y - 500.)**2)))
md.hydrology.moulin_input[pos,:] = 4
## BC
md.hydrology.neumannflux=np.zeros((md.mesh.numberofelements+1, num_ts))
md.hydrology.neumannflux[-1,:] = time
## verbose
md.verbose.solution=1
## solve
md=solve(md,'Transient')
## plot
plotmodel(md,"finalSol", 'data', md.results.TransientSolution[-1].EffectivePressure, 'title', 'Pressure effective in Pa',
                     'data', md.results.TransientSolution[-1].HydrologyHead, 'title', 'Head in m',
                     'data', md.results.TransientSolution[-1].HydrologyBasalFlux, 'title', 'Basal water flux in m2s-1',
                     'data', md.results.TransientSolution[-1].HydrologyGapHeight, 'title', 'Gap height in m ')

This is the moulin.par.py

mport numpy as np                                                                                                                                      
from frictionshakti import *
from verbose import *
from SetIceSheetBC import *

#Start defining model parameters here
# Set up bed topography and ice geometry for a tilted 500m thick slab
md.geometry.base = .02*md.mesh.x
md.geometry.bed = md.geometry.base
md.geometry.surface = .02*md.mesh.x + 500.
md.geometry.thickness = md.geometry.surface - md.geometry.bed

# Define ice sliding velocity (m/yr) 
md.initialization.vx = 1.0e-6*md.constants.yts*np.ones((md.mesh.numberofvertices))
md.initialization.vy = np.zeros((md.mesh.numberofvertices))
md.initialization.pressure=np.zeros((md.mesh.numberofvertices))
md.initialization.pressure=np.zeros((md.mesh.numberofvertices))

# Materials
# Ice flow law parameter (note that the standard parameter A=B^(-3))
# Not sure about this one in matlab it is (5e-25)^(-1/3)
md.materials.rheology_B= (5.e-25)**(-1/3)*np.ones((md.mesh.numberofvertices))
md.initialization.temperature=(273.)*np.ones((md.mesh.numberofvertices))
md.materials.rheology_n=3.*np.ones((md.mesh.numberofelements))

#Calving
md.calving.calvingrate=np.zeros((md.mesh.numberofvertices))

# Friction law and coefficient
md.friction=frictionshakti(md)
md.friction.coefficient=20.*np.ones((md.mesh.numberofvertices))


#Numerical parameters
#md.stressbalance.viscosity_overshoot=0.0;
md.masstransport.stabilization=1.
md.thermal.stabilization=1.
md.verbose=verbose(0)
md.settings.waitonlock=30
md.stressbalance.restol=0.05
md.steadystate.reltol=0.05
md.stressbalance.reltol=0.05
md.stressbalance.abstol=float('NaN')
md.timestepping.time_step=1.
md.timestepping.final_time=3.

#GIA: 
md.gia.lithosphere_thickness=100.*np.ones((md.mesh.numberofvertices)) # in km
md.gia.mantle_viscosity=1.e+21*np.ones((md.mesh.numberofvertices)) # in Pa.s
md.materials.lithosphere_shear_modulus=6.7*1.0e+10                   # in Pa
md.materials.lithosphere_density=3.32                              # in g/cm^-3
md.materials.mantle_shear_modulus=1.45*1.0e+11                       # in Pa
md.materials.mantle_density=3.34                                   # in g/cm^-3

#Boundary conditions:
md=SetIceSheetBC(md)

#Change name so that no test have the same name
#A=dbstack
#if (length(A)>2):
#    md.miscellaneous.name=A(3).file(1:end-2)

bdefleurian

Hi Anne
I just came back from holidays sorry for the long wait on that matter. I am able to run the model you pasted above and get the results that are expected. The only issue that I noted is with your call to hydrologyshakti but that should throw an error if it does not work so not sure if it is an error or just a copy past issue.
md.hydrology=hydrolyshakti() should be md.hydrology=hydrologyshakti()
Also minor points but you should not need to import sys and issmversion.
Basile