Setting Vel upper limit?

wadewei

Hello everybody,

I am trying to refine one of my modelled glaciers and start from previous thickness/temperature fields. However, although solving for stressbalance returns healthy results, it would blow up after several transient steps with velocities at some vertices rapidly increasing to 1E7 m/yr.

I am trying to zoom in and take a closer look, but plotmodel gets very slow since I have almost three million elements (anything I can do to speed it up?). I wonder if it's possible to set an upper limit to the velocity so that I can run it for a while and see if the surface would smooth itself to a point that the velocity would drop to a reasonable value to allow for a smooth run.

Thank you!

Best regards,
Wade

mathieumorlighem

Hi Wade,

unfortunately there is no limiters implemented. If you would like, you can try to do it yourself in src/c/analyses/StressbalanceAnalysis.cpp (one of the InputUpdateFromSolution). That being said, there is always a reason for the speed to blow up, either a problem with the model set up, the boundary conditions, or the time stepping that may be too large. Alternatively, you can constrain (spcvx and spcvy) the vertices that blow up, but I think it is worth investigating the issue further.
If you have an idea of the area that blows up (using find(md......Vx>1e7)), you can "extract" a model around that area and use this smaller model to plot everything.

I hope this helps
Best
Mathieu

wadewei

mathieumorlighem
Hi Mathieu,

Thank you for you reply. I am still exploring ways to make it work. However I noticed something interesting and would like to hear your thoughts. Both cases are continental ice sheet only with nothing to do with the ocean. spcvx,vy,vz are set to 0 on lateral boundaries and NaN on the base.

In a previously working model setup, md.mask.ice_levelset was set to -1 in ALMOST everywhere, although no ice should exist on a large portion of the model domain. These places have thickness==masstransport.min_thickness. I wonder if this is potentially problematic. Should I set ice_levelset to -1 strictly on ice-covered areas? I will try it and compare the results right away but would like your confirmation, does ice_levelset actually evolve every once in a while as specified by levelset.reinit_frequency?
I realized my higher-resolution model, which borrowed almost everything from the said working model through interpolation functions provided and yet still refuse to run, has mask.ice_levelset=-1 everywhere. In a desperate attempt, I changed mask.ice_levelset to 1 in places where thickness==masstransport.min_thickness. Miraculously the model would solve now without blowing up. could it be because some places with thickness equal to masstransport.min_thickness=1 that was causing nonconvergence is now ignored?

When you mention blow up during velocity calculation is largely due to model set up and boundary conditions, is there a universal rule to follow, and what would you inspect/modify when it blows up? Once spcvx,y,z along with some other stuff are set by invoking md=SetIceSheetBC(md), is there anything in particular that usually cause nonconvergence and needs special attention? Something that resembles an emergency checklist for pilots would be very helpful.

Best regards,
Wade

mathieumorlighem

Hi Wade,

the location where your level set is 0 is absolutely critical, because it determines where water (or air) pressure is applied at the ice front. Level set = 0 is basically your ice front, and regions where level set >0 are simply ignored in the calculation. I think your approach #2 is more appropriate since you should not include elements where there is no ice in your model.

For your second point, we tried to work hard on ismodelselfconsistent which checks every field in your model, and catches most problems but not all. It is virtually impossible to catch everything, but one thing I strongly recommend is to plot your boundaries conditions with plotmodel(md,'data','BC'); before you run to make sure that things look good.

Best
Mathieu

wadewei

mathieumorlighem

the location where your level set is 0 is absolutely critical, because it determines where water (or air) pressure is applied at the ice front. Level set = 0 is basically your ice front

Hi Mathieu,

I compared two model runs with identical setup, except for the size of (sub-)domain with levelset=-1. They are all land-ice only, but was not specified with any vertices with levelset=0, although ice did not fully cover the entire domain. The first one has levelset=-1 at vertices with ice thickness>min_thickness only, while the second one has levelset=-1 in areas slightly larger than the former.
After some transient time steps, nowhere previously without ice became ice-covered in the first model. However, in the second setup, ice did extend into areas not previously occupied by ice (thickness=min_thickness & levelset=-1).

This puzzled me a lot. I did peek into FemModel::ResetLevelset and FemModel:: DistanceToFieldValue but wasn't able to fully understand the mechanism at a first glance. However, looking at this reminder, I wonder if the reason why my model would not reinitialize or evolve mask.ice_levelset is that I did not specify any vertices with levelset=0?

Best,
Wade

mathieumorlighem

Hi Wade,

you do not need (and should not) have vertices with level set = 0 explicitly, ISSM takes care of figuring out where in the mesh the boundary is. ISSM will not change your level set if you have md.transient.ismovingfront = 0 (I recommend not allowing moving boundaries unless you know what you are doing). If you have level set = -1 everywhere, then after the first step ISSM will change the ice thickness to min_thickness automatically.
Does that answer your question?
Mathieu