Hi Mathieu,
thanks a lot for your help 🙂
I'm not sure if I understand your reply entirely: with Neumann BC and
md.mask.ice_levelset(md.mesh.vertexonboundary==1)=-1;
I'm in the second case above where the blow-up occurs.
I've tried
md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
md.mask.ocean_levelset(md.mesh.vertexonboundary==1) = -1;
though with Neumann BC. Maybe this is what you meant? The runtime for was much faster (6-7 min) compared to previously where the ocean levelset was 1 everywhere (20 min). It's still not as fast as with Dirichlet boundary conditions, but at least it's closer.
Depending on how you have the Dirichlet conditions implemented I guess you are solving for much fewer degrees of freedom, resulting in smaller linear systems and faster runtime. If 15% of the vertices are on the boundary, ignoring those for a direct solver with O(N3) complexity would reduce the runtime to 0.853 = 61%. Here we have an iterative solver of course and sparse linear systems but also nonlinearities in the governing equations so this "analysis" doesn't carry over directly, but it's at least an explainable guess.
Anyway, I'd be happy with the 7min/5years runtime, plus I tested running it for 100 years and it didn't blow-up either - so a huge improvement to what I had previously. Is the change to the mask what you had in mind or did I misunderstand?
Best regards,
Nicole