Hello,
I am using ISSM Version 4.24 with Python 3.11.11 on Ubuntu 24.04. I am encountering an error when I try to run a transient simulation of a 1D (flowline) free surface ice sheet model with periodic boundary conditions between inflow and terminus.
To implement the BCs I based myself on the procedure used for the IsmipF.py experiment (the BCs are found in its respective runme.py script). The error occurs when the simulation is running its first timestep. The exact error message is this:
call computational core:
iteration 1/5 time [yr]: 1.00 (time step: 1.00)
computing new velocity
computing slope...
extruding SurfaceSlopeX from base...
computing slope
computing basal mass balance
maximum number of nonlinear iterations (100) exceeded
computing basal mass balance
computing mass transport
call free surface computational core
[0] ??? Error using ==> ./classes/Loads/Penpair.cpp:187
[0] PenaltyCreateKMatrix error message: analysis 1539 (FreeSurfaceBaseAnalysis) not supported yet
-----------------------------
AI says that it appears to be a conflict between my periodic boundary conditions and the free surface configuration in the mass transport module.
This is my setup_periodic_boundary method:
def setup_periodic_boundary(self, md):
"""Configure periodic boundary conditions between inflow and terminus"""
# Find vertices at the inflow (minimum x) and terminus (maximum x)
inflow_vertices = np.where(np.isclose(md.mesh.x, self.x_params['start']))[0]
terminus_vertices = np.where(np.isclose(md.mesh.x, self.x_params['end']))[0]
# Ensure we have matching numbers of vertices to pair
min_vertices = min(len(inflow_vertices), len(terminus_vertices))
if min_vertices == 0:
print("Warning: Could not find boundary vertices for periodic conditions")
return md
inflow_sorted = inflow_vertices
terminus_sorted = terminus_vertices
# Pair the vertices (add 1 because ISSM uses 1-based indexing)
paired_vertices = np.vstack((inflow_sorted[:min_vertices] + 1,
terminus_sorted[:min_vertices] + 1)).T
# Set the vertex pairing for both stress balance and mass transport(?)
md.stressbalance.vertex_pairing = paired_vertices
md.masstransport.vertex_pairing = paired_vertices
print(f"Set up periodic boundary conditions with {min_vertices} paired vertices")
return md
This is my setup_mass_balance module:
def setup_mass_balance(self, md):
"""Configure mass balance settings"""
# Initialize surface mass balance
md.smb.mass_balance = np.zeros((md.mesh.numberofvertices))
# Set up mass transport configuration
md.masstransport.spcthickness = np.nan * np.ones((md.mesh.numberofvertices,))
md.masstransport.min_thickness = self.mass_transport_params['min_thickness']
md.masstransport.hydrostatic_adjustment = self.mass_transport_params['hydrostatic_adjustment']
md.masstransport.isfreesurface = 1
md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
md.basalforcings.geothermalflux = np.zeros((md.mesh.numberofvertices))
return md
In the runme.py file for IsmipF.py the vertex pairing is set like I have set it in my BCs module:
if ParamFile=='IsmipF.py':
# if we are dealing with IsmipF the solution is in masstransport
md.masstransport.vertex_pairing = md.stressbalance.vertex_pairing
However, the main difference between the example and my model are that
in my 1D model I explicitly set md.masstransport.isfreesurface = 1 and md.masstransport.hydrostatic_adjustment = 'Incremental', which triggers the FreeSurfaceBaseAnalysis.
NOTE: I've verified that the model runs successfully with periodic boundary conditions when mass transport is disabled (md.transient.ismasstransport = 0), which confirms the periodic boundary implementation works properly for stress balance.
I have already tried some workarounds:
Using a different hydrostatic adjustment method: 'Absolute'
This yielded the same error as previously
[0] ??? Error using ==> ./classes/Loads/Penpair.cpp:187
[0] PenaltyCreateKMatrix error message: analysis 1539 (FreeSurfaceBaseAnalysis) not supported yet
Setting isfreesurface = 0 and handling surface elevation changes differently
This yielded the following error:
solver residue too high!: norm(KU-F)/norm(F)=3.61814e-06 > 1e-06 (md.settings.solver_residue_threshold)
WARNING: Solver failed, Trying Recovery Mode
solver residue too high!: norm(KU-F)/norm(F)=3.61814e-06 > 1e-06 (md.settings.solver_residue_threshold)
[0] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[0] Solverx error message: Recovery solver failed...
--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpiexec detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:
Process name: [[36558,1],0]
Exit code: 1
--------------------------------------------------------------------------
Since none of my solutions worked, I was wondering if if there's a recommended approach for handling free surface evolution with periodic boundaries in a 1D flowline model with these settings?
If you need more information please let me know and I'll happily provide it. Thanks for your time and guidance 🙂
Cheers,
Ana
