HO transient runs blow up

wadewei

Hi everyone,

I have a transient run that always blew up after some steps with large calculated velocities. It's been run with version 4.16 precompiled binary for quite a long time before crashing.The topography near the location where abnormally large velocity appears has been smoothed slightly and it only seemed to postpone the problem rather than solving it.
I've tried to switch to SIA or FS instead of HO and it doesn't seem to work either. Any suggestion is welcomed

Here's a zoomed-in plot of the area in question and some settings that may be relevant.
https://pkueducn-my.sharepoint.com/:i:/g/personal/wadew_pku_edu_cn/EbBiqp68D6xMikDubM9gCEABK5-Q8fI4wco2PyhxtQRW5g?e=aAdki7
https://pkueducn-my.sharepoint.com/:i:/g/personal/wadew_pku_edu_cn/EVMcaV0TIR9NljT3mqUtm58BkJypLfk5-Kd33wpfbrBH-Q?e=J5KgwA

md.initialization.pressure=md.materials.rho_icemd.constants.gmd.geometry.thickness;
md.friction.coefficient=80*ones(md.mesh.numberofvertices,1);
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=ones(md.mesh.numberofelements,1);
md.stressbalance.restol=0.005;
md.stressbalance.reltol=0.05;
md.stressbalance.abstol=NaN;
md.stressbalance.maxiter=40;

---log----
solver residue: norm(KU-F)/norm(F)=5.16767e-07
mechanical equilibrium convergence criterion 229.749 > 0.5 %
Convergence criterion: norm(du)/norm(u) 22.4449 > 5 %
Convergence criterion: max(du) 2.04841e+07
number of unstable constraints: 0

solver residue: norm(KU-F)/norm(F)=2.77637e-07
mechanical equilibrium convergence criterion 85.0726 > 0.5 %
Convergence criterion: norm(du)/norm(u) 14.4541 > 5 %
Convergence criterion: max(du) 1.61737e+07
number of unstable constraints: 0

solver residue: norm(KU-F)/norm(F)=1.69396e-07
mechanical equilibrium convergence criterion 39.738 > 0.5 %
Convergence criterion: norm(du)/norm(u) 9.41795 > 5 %
Convergence criterion: max(du) 1.20681e+07
number of unstable constraints: 0

solver residue: norm(KU-F)/norm(F)=4.5098e-07
mechanical equilibrium convergence criterion 24.8383 > 0.5 %
Convergence criterion: norm(du)/norm(u) 6.184 > 5 %
Convergence criterion: max(du) 8.67112e+06
number of unstable constraints: 0

solver residue: norm(KU-F)/norm(F)=4.39168e-06
solver residue too high!: norm(KU-F)/norm(F)=4.39168e-06 => Trying recovery solver
WARNING: Solver failed, Trying Recovery Mode

solver residue: norm(KU-F)/norm(F)=2.54363e-07
mechanical equilibrium convergence criterion 19.7379 > 0.5 %
Convergence criterion: norm(du)/norm(u) 4.08125 < 5 %
Convergence criterion: max(du) 6.07546e+06
number of unstable constraints: 0 solver residue: norm(KU-F)/norm(F)=5.16767e-07
mechanical equilibrium convergence criterion 229.749 > 0.5 %
Convergence criterion: norm(du)/norm(u) 22.4449 > 5 %
Convergence criterion: max(du) 2.04841e+07
number of unstable constraints: 0

solver residue: norm(KU-F)/norm(F)=4.39168e-06
solver residue too high!: norm(KU-F)/norm(F)=4.39168e-06 => Trying recovery solver
WARNING: Solver failed, Trying Recovery Mode

[0] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[0] Solverx error message: Recovery solver failed...

[22] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[22] Solverx error message: Recovery solver failed...

[3] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[3] Solverx error message: Recovery solver failed...

mathieumorlighem

Hi Wade,

most of the time when the solver does not converge (or even fails, like here) it is a problem of boundary condition or model set up. You should run the stress balance with md.stressbalance.maxiter = 2; and take a look at the results.
I hope this helps
Best
Mathieu

wadewei

mathieumorlighem
Good point! I realized md.basalforcings.geothermalflux was set 1000 times larger than realistic value (0.04-0.08 W/m2). Also corrected md.thermal.spctemperature form something like repmat(surfacetemperature,numoflayers) to md.thermal.spctemperature(surfacevertices)=surfacetemperature and md.thermal.spctemperature(elsewhere)=NaN.
Looks like velocity calculation converges now. However, thermal calculation refuses to converge. Is there something in particular that may be suspicious given this error, or do you have a recommended troubleshooting path?

Thanks,
Wade

iteration 2/80 time [yr]: 162.2 (time step: 0.125)
updating boundary conditions...
computing temperatures

solver residue: norm(KU-F)/norm(F)=2.83844e-16
number of unstable constraints: 3634

solver residue: norm(KU-F)/norm(F)=2.83258e-16
number of unstable constraints: 5281

solver residue: norm(KU-F)/norm(F)=2.85347e-16
number of unstable constraints: 5108

solver residue: norm(KU-F)/norm(F)=2.83767e-16
number of unstable constraints: 2792

solver residue: norm(KU-F)/norm(F)=2.82612e-16
number of unstable constraints: 391

solver residue: norm(KU-F)/norm(F)=2.83158e-16
number of unstable constraints: 102

solver residue: norm(KU-F)/norm(F)=2.82362e-16
number of unstable constraints: 84

solver residue: norm(KU-F)/norm(F)=2.82724e-16
number of unstable constraints: 84

solver residue: norm(KU-F)/norm(F)=2.82211e-16
number of unstable constraints: 83

solver residue: norm(KU-F)/norm(F)=2.84148e-16
number of unstable constraints: 90
maximum number of iterations (10) exceeded
computing melting

mathieumorlighem

With the penalty method, it is not "unexpected" that the steady state solver would not converge. You may want to try md.thermal.isenthalpy = 1; you may get better convergence. If not, you can just use md.thermal.penalty_lock = 10; or something like that (it will freeze penalties if they keep switching after 10 iterations).
Best
Mathieu

wadewei

mathieumorlighem
md.thermal.penalty_lock = 10 did not seem to work. Setting isenthalpy=1 does converge the thermal model although I had to set below parameters to zeros per nightly run scripts. Do they actually matter for land glaciers?
md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);

Best,
Wade

mathieumorlighem

ok the enthalpy formulation is more accurate and do not worry about these initial conditions, they don't matter for what you are trying to do.
Mathieu

wadewei

mathieumorlighem
Good to hear. One more thing before I let you go (hopefully). Is there a 4.17 binary avavilable that works with CentOS 7? 4.16 works with stock glibc 2.17 but the newer binary seem to be compiled with a higher version one and it complained about /lib64/libm.so.6: version `GLIBC_2.27' not found (required by /opt/ISSM/4.17/lib/IssmConfig_matlab.mexa64).
I am interested to see what kind of performance improvement is available and yet still remember previous compilation attempt that failed with mysterious and seemingly random petsc segmentation faults.

Wade

justinquinn

wadewei This is very good to know as we have a limited number of Linux configurations on which we can test our precompiled packages. I will work on including our version of libm in the package (we do so already with libgfortran) and will let you know here when it is available.

justinquinn

wadewei Which minor release version of CentOS are you running? If updating the OS or GLIBC is an option, I would go with that. I don't know how far back we can support varied Linux builds, but I will bring it up in our team meeting this week. Alternatively, you might try compiling from source.

wadewei

justinquinn I am currently using CentOS 7.5/7.6. If I recall correctly, glibc version remains the same throughout CentOS minor releases (e.g. 7.0-7.x) so updating to 7.7 will not help and involves updating InfiniBand and RAID controller drivers which may not be feasible. With CentOS 8 released only last year, 7 will remain mainstream for another couple years I think.
I will consider compiling from source if 4.17 does provide significant improvement but let me know if you guys could come up with something for redhat users :-)

Thanks,
Wade

justinquinn

wadewei Hi Wade, we're not really in a position to spend the resources building out the infrastructure to support a large number of Linux variants/versions. I understand that CentOS is not that old, but we, for example, only support back through Ubuntu Linux 18 (Ubuntu 16 uses an older version of glibc as well). Really, these precompiled packages are intended for users who have little experience with the command line or compiling from source. If we get a large number of requests for backward compatibility on the precompiled packages we may need to reassess.