Convergence issue with parallel on MG Multigrid Solver

[HanyMKhalil]

Dear Romain, I ran a simple 2D model on parallel 8 cores using default MG solver, yet the model runs fine for half of the time period then it couldn't converge any more (reach iterations limit of 500), I tried change solver parameters but did not work, I attached you the simple code (I added txt to upload it) for a quick look and also SLURM job and output file. note: the model runs fine with mumps. any advice? thanks so much.

scenario3-Ref-10-1-8core-defaultmg.py.txt

scenario3-Ref-10-1-8core-defaultmg.slm.txt

5016706.output.txt

[rbeucher]

Hi @HanyMKhalil

There are many reasons why this could happen. That could be due to some Viscosity contrasts due to Temperature / pressure discontinuities, scaling can also be a pb. That could also be the model being unable to pick a solution between equivalent options... Solver can go back and forth between activating one fault or another... etc.

If you say that Mumps is working properly, that might be because of numerical instabilities. I am mentioning this because some of my 3D models have shown that recently. I had some good result using the penalty method with mg:

Just add

Model.solver.set_penalty(1.0)

Before init_model

Give it a try and let me know

R

[HanyMKhalil]

Thanks Romain, ok I will try this

[rbeucher]

If it does not work we can have a look at the details of your model.

[HanyMKhalil]

Hi Romain, it works very well actually without any problems. So I just want to know what exactly the penalty did and what does it mean to have penalty 1? like how this affects the computations and results?

Thanks a lot really

[rbeucher]

Hi @HanyMKhalil

So I think we need to work on the documentation about that penalty method. To make it short... The penalty method helps to relax the constraint on incompressibility. If the model boundary conditions are such that the incompressibility might be too constraining... if you have an inflow for example.. The penalty can help to find a solution.

@julesghub @jmansour any comments?

R

[HanyMKhalil]

Hi @rbeucher you think is there another way to get around this without penalty for example smooth viscosity so that no big contrast between the viscosities happens, or another way? in my models I put minimum and maximum viscosity values. @julesghub @jmansour ?

[HanyMKhalil]

@jmansour Hi John as we discussed just a reminder on how to calculate the divergence between the inflow and outflow so that to check which penalty parameter is the best. I attached you my simple code, it does work and converge without penalty for first couple of time steps then it failed, but with penalty 1 it is ok. I also though of making the air layer comprissible but in that case I won't be able to use penalty. I also though of adjusting the CFL parameter?

mgfull25 2.p.txt

[jmansour]

Hi @HanyMKhalil

I'd be interested to know if adding an air layer restores solver convergence without using the penalty method.

@rbeucher, in @HanyMKhalil's model he hasn't set a top condition. Does this do the same as Underworld would, resulting in the open condition? Also, does LecodeIsostasy result in a Neumann or Dirichlet type condition?

Also, how should Hany obtain the velocity field? I assume he'll construct the divergence as we do in UW:

grad_v = velocityField.fn_gradient
div_v = grad_v[0] + grad_v[3]

https://underworld2.readthedocs.io/en/v2.8.1b/build/underworld.mesh.html#underworld.mesh.MeshVariable.fn_gradient

You can then visualise div_v to see how divergent your flow fields are and whether your penalty is sufficient. Apparently 10% of your maximum viscosity contrast is a reasonable starting place.

You might also like to see what effect switching to a constant viscosity has on your solver convergence. I don't think CFL should matter.

[rbeucher]

@jmansour , yes no condition on the top means it is open.

@HanyMKhalil You can get the velocityField from Model.velocityField

R

[HanyMKhalil]

Hi @jmansour Yes, I calculated the divergence and it shows values really close to zero, attached in the picture one test, I will try with different penalty magnitudes and see which one will give me the closest to zero. But I am happy that it helps my 2D and 3D models converge smoothly. maybe I will try to calculate Peclet number from the models however I do not know exactly how to do that, any clues?