Details about the Full-Stokes solver

nicare

Hi,
over the last few weeks I've been trying to understand how to use ISSM, in particular how to use the Full-Stokes model. I understand that FS is tricky, so I've given it some time to test around with different discretizations and settings. Unfortunately, I don't quite trust the results I've been getting with FS, details below. We had a project meeting today, and we agreed we'd really want to use FS but that it doesn't seem to work thus far. My background is in model order reduction, so it's an obvious choice because of its cost, especially long term.

I guess what I'm looking for at this point is more details about how the FS model is solved internally, and if there are any settings that might help with stability. For instance, what's the difference between MINI and MINI-condensed? I've tried Taylor Hood, but the solver wouldn't converge regardless of mesh size. I've been going through the documentation on the website, but I couldn't find much in regard to computational details outside the use of ISSM.

A big question that came up today was mass conservation in the stressbalance equations. If I understand correctly, for SSA and HO the pressure is computed from the velocity such that mass is conserved. For the Stokes problem this is only weakly enforced unless the finite element scheme is chosen appropriately e.g. such that the pressure space contains P0 (e.g. https://link.springer.com/article/10.1007/s10915-011-9549-4 ). For the FE schemes implemented in ISSM, is mass conservation guaranteed at the element level? I don't think Taylor-Hood (P2-P1) does it, but maybe the MINI schemes? Have you tried using P0 or P1+P0 as pressure space?

I'm happy to discuss per email if there's a specific contact person on the ISSM side.

Best regards,
Nicole

nicare

The reason I don't trust the FS solver yet is the following: I've been using the Greenland outline and parameterization from the tutorial with different discretizations (code below) to get an idea for how much solutions differ between models and meshes. While SSA and HO show good consistency in regards to velocity and pressure, the FS results are somewhat different:

The FS velocity maximum is twice as high as with the other models.
the general behavior of the FS velocity is different, it is by far not as smooth, and it has a lot of local fluctuations. I've not been able to get rid of them via mesh refinement, and I'm wondering if they are caused by numerical instability.
While the FS pressure agrees with SSA and HO overall, the SSA and HO pressure is non-negative but the FS pressure has some huge local negative values (order -3e+7). In particular, the HO pressure on the surface is a constant zero, whereas the FS surface pressure ranges from -1.7252e+07 to 1.9075e+07 with mean -2.0137e+03.

Velocity comparison (at surface for HO and FS):

Pressure comparison (at base for HO and FS):

Mesh: bamg(md,'hmax',25000,'hmin', 100,'gradation',1.5,'field',vel,'err', 2);

FS pressure at surface:

Matlab script:

stressbalacegrids.txt

9kB

mathieumorlighem

I will send you an email with some documentation that we wrote a long time ago. Happy to hop on Zoom if you have further questions.
Mathieu

mathieumorlighem

While I am here...
There is no Dirichlet Condition applied to the pressure. It should reach 0 at the surface but we can't force it. Yes, negative pressures are not realistic, it may be that the model is too coarse, but it is hard to say. FS Greenland-wide is difficult, I am not surprised it does not work out of the box.
It is normal if the speed of HO is higher than SSA, and FS higher than HO. As the ice is non-Newtonian, the viscosity decreases with the effective strain rate. Since FS considers "more" components in the strain rate tensor, the effective strain rate is higher, and the viscosity lower, hence a higher speed. Inversion are generally model dependent. You need more friction for HO, and even more friction for FS to fit the surface speed as we have more internal deformation.
I hope this makes sense!
Mathieu