I've been trying to run some small test simulations of an ice shelf using ParStokes, with an eye towards scaling up to looking at the long-term evolution of Larsen C. Unfortunately, I don't seem to be able to get the parallel solver to converge to anything remotely sensible. Identical runs with the standard Navier-Stokes solver work fine, as do similar runs with grounded geometry (both ParStokes and Navier-Stokes). As far as I can tell, it's just the combination of ParStokes and ice shelf type boundary conditions which breaks down.

So, two questions:

Is there anything in particular that has to be done to get ParStokes to play nice with Neumann boundary conditions? As far as I can tell from the code, it should "just work", but is there something I've missed?

Has anybody successfully used ParStokes in an ice shelf setup? Is there any special configuration required?

Dear Martin,
if your "normal" Stokes solver produced a solution, then - in principle - the ParStokes should do so, too.
Could you please indicate which part of ParStokes refuses to converge: the outer (GCR) system or one of the block-preconditioning sub-problems?
Do you have grounding-line dynamics switched on? Perhaps there is an incompatibility in these procedures, which, to my knowledge, haven't been tested with the block pre-conditioner.

Initially, using the actual ice shelf geometry, I thought that the GCR was failing to converge. However, I've been trying to reduce things to a simpler test case, and when I use a simple square geometry, I get (fairly slow) convergence of the GCR, but to a nonsensical result.

I haven't been using any grounding line dynamics - the domain is fully floating, with velocities specified along the upstream and lateral boundaries, and pressure boundary conditions (using the Buoyancy USF) on the downstream and lower boundaries.

I've attached a test case you can try for yourself. The "normal" Stokes solution is generated by seq.sif and the ParStokes by par.sif.