Currently doing a MISMIP Full-Stokes run using the Python interface and running into errors, here is the part of the script which configures the run.
if org.perform('GlenFS'):
# If using a 2D mesh (e.g. 500m case), use the steady state model
if modelnum == 5 or modelnum == 6:
md = org.loadmodel('Transient_Steadystate')
else:
# or extruded model
md = org.loadmodel('Transient_steadystate3')
# or load from GlenHO to reuse geometry
# md = org.loadmodel('GlenHO')
md.initialization.vx = md.results.TransientSolution[-1].Vx
md.initialization.vy = md.results.TransientSolution[-1].Vy
md.initialization.vel = md.results.TransientSolution[-1].Vel
md.geometry.thickness = md.results.TransientSolution[-1].Thickness
md.geometry.base = md.results.TransientSolution[-1].Base
md.geometry.surface = md.results.TransientSolution[-1].Surface
md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
md = md.extrude( 10, 1.1)
md = setflowequation(md, 'HO', 'all')
md.transient.isthermal = 0
md.transient.issmb = 0
md.transient.ismasstransport = 0
md.initialization.temperature[:] = 273.0
md = setflowequation(md, 'FS', 'all')
md.stressbalance.shelf_dampening = 1
md.masstransport.isfreesurface = 1
md.transient.isgroundingline = 0 # or 1, depending on your test
md.groundingline.migration='Contact' # or 'SoftMigration'
md.timestepping.time_step = 0.00001
md.timestepping.final_time = 0.0001
md.settings.output_frequency = 1
# md = md.SetInput(md, 'Bed', md.geometry.bed)
md.toolkits = toolkits.addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
md.flowequation.fe_FS = 'TaylorHood'
md.stressbalance.maxiter = 20
md.stressbalance.restol = 0.5
md.stressbalance.reltol = 0.001
# If you know surface and thickness:
#md.geometry.bed = md.geometry.surface - md.geometry.thickness
#md.geometry.surface = InterpFromGridToMesh(x, y, usrf, md.mesh.x, md.mesh.y, 0);
#md.geometry.bed = InterpFromGridToMesh(x, y, topg, md.mesh.x, md.mesh.y, 0);
# Or if you already interpolated base (from NetCDF):
#; % Sometimes 'base' is used internally
md.stressbalance.abstol = float('nan')
md.verbose.convergence = False
md.cluster = cluster
md.miscellaneous.name = 'MISMIP_' + modelname + '_GFS'
print(md.mesh)
solutiontype = 'tr' # 'tr' = transient, 'sb' = stress‑balance, etc.
md = solve(md, solutiontype,'loadonly', loadonly, 'lock', lock, 'runtimename', False)
#org.savemodel(md)
and the output I get when running it is
call computational core:
iteration 1/10 time [yr]: 0.00 (time step: 0.00)
computing new velocity
computing slope
Updating constraints and active domain of analysis L2ProjectionBaseAnalysis for time: 1e-05
Allocating matrices (Kff stiffness matrix size: 54534 x 54534)
Assembling matrices
Create nodal constraints
Dirichlet lifting applied to load vector
Solving matrix system
Merging solution vector from fset to gset
Updating inputs from solution for L2ProjectionBaseAnalysis
Updating constraints and active domain of analysis L2ProjectionBaseAnalysis for time: 1e-05
Allocating matrices (Kff stiffness matrix size: 54534 x 54534)
Assembling matrices
Create nodal constraints
Dirichlet lifting applied to load vector
Solving matrix system
Merging solution vector from fset to gset
Updating inputs from solution for L2ProjectionBaseAnalysis
[31] ??? Error using ==> ./classes/Elements/Element.cpp:1476
[31] GetInputListOnVertices error message: input Bed not found in element
any idea on what the issue is?
