Hey all
I am performing a full-Stokes run of the MISMIP configuration and running into issues. The run crashes after the third time step after failing to invert the stiffness matrix.
analysis for time: 3e-05
Renumbering degrees of freedom
Get solution from inputs
Reduce vector from g to f set
Input updates from constant
Updating inputs from solution for StressbalanceAnalysis
Allocating matrices (Kff stiffness matrix size: 12725643 x 12725643)
Assembling matrices
[
[16] ??? Error using ==> ./classes/Elements/PentaRef.cpp:260
[16] GetJacobianDeterminant error message: negative jacobian determinant!
Any thoughts on how this could be fixed? Here is the run config:
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 = 1 # 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)
