Hi Ben,

On 06/06/2018 09:06 PM, Ben Melosh wrote:

I am trying to model Stokes flow at very high viscosities (~10^20), from my results it looks like I cannot achieve this and am better suited using viscosities of ~10^5 or so. I know there is this older post [1] which states: "A word of warning though: SfePy does not have a dedicated Navier-Stokes solver (the Newton (+ and direct) solver use that is in the navier_stokes.py example is rather naive, and works for small problems with only, with viscosity "high enough")"

Has there been a different solver implemented since this was posted? If not, what do you recommend to work with such large viscosities?

I have only some general recommendations:

• having a matrix with very different magnitudes of items in different blocks is problematic - you might need to use some scaling, or choose different physical units for your input.

• You can try using the FieldSplit preconditioner from PETSc for the Stokes problem, if it is so large that it cannot be solved using a direct solver (Umfpack, or, in the git version, also MUMPS). Some of the relevant options of the solver are mentioned in the docstring of examples/multi_physics/biot_parallel_interactive.py, see PETSc docs for more, and examples/multi_physics/biot_short_syntax.py shows how to pass any additional options to the PETSc solvers from a problem description file.

Does it help?

Cheers, r.

Thank you

[1] "Problem with the mesh while trying to solve the Navier Stokes equations with sfepy"

---

SfePy mailing list -- sfepy@python.org To unsubscribe send an email to sfepy-leave@python.org https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/

On 11/6/18 1:11 AM, Ben Melosh wrote:

Thanks for the help Robert,

I will look into these different options. I am wondering what specifically is problematic about using very different or very large magnitudes of items in different blocks? Is it a problem with processing power? I am hoping for specifics so I can better understand what is going on in the model and improve it.

It is a problem of machine precision - when using the 64 bit float values, 1 + 1e-17 is still just 1, so having a block with values of the order of magnitude 1e20 and another or the order 1e3, the second block is essentially zeros from the global matrix point of view.

r.