Good choice of solver parameters on parallel Stokes mixed formulation resolution?

Linear Algebra
1

Hello,
I encounter some problems in scaling of transport problems to large mesh: how original ! :slight_smile:
Like a lot of threads, my codes actually work for small meshes but not for larger ones.

So, I ran in parallel (with about 20 processors available on my machine). The problem is that some of my codes don’t converge anymore and I wonder if this is due to the capacity of the machine or to my choices of solver (iterative, tolerance etc) and preconditioning. I’m not yet familiar with this world (the solver world), that’s why I need your help.

For example, in the case of a Stokes mixed formulation, with an undefined matrix, I know that it is possible to use the minres or tfqmr solver with an amg preconditioner. However, I get an error message like "PETSc reason DIVERGED_INDEFINITE_MAT with a residual norm ||r||=0.00000e0 for the solver minres. The tfqmr code is spinning his wheels.

Please note that I tested the direct mumps solver (without preconditioning) in parallel (notifying or not the number of threads to 1 with the command OMP_NUM_THREADS=1 as some threads in the forum recommend it) and this time I get an error message of the form “DIVERGED_PC_FAILED” with a residual norm ||r||=0.00000e0.

Thank you in advance for your help

Here is the code (i did not put the mesh files)

from dolfin import *
import sys

mesh = Mesh()
with XDMFFile(sys.argv[1]) as infile:
infile.read(mesh)
File("Dolfin_circle_mesh.pvd").write(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile(sys.argv[2]) as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("Dolfin_circle_facets.pvd").write(mf)

#mixed space for the stokes equation
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# creates measure for the inflow on the left of my geometrie
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
ds = Measure("ds")(subdomain_data=boundaries)
ds_G = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=100001)

# constant and geometrical entities
n = FacetNormal(mesh)
f = Constant((0, 0))
mu1 = 0.001
mu = Constant(mu1)

# boundaries conditions on interior obstacle (bc_os) and extern square boundaries (outflow + symmetry condition)
bc_os = DirichletBC(W.sub(0), (0.0,0.0), mf, 100000)
bc_noflow = DirichletBC(W.sub(0).sub(1), 0.0, mf, 100004)
bc_2noflow = DirichletBC(W.sub(0).sub(1), 0.0, mf, 100002)
bcp_outflow = DirichletBC(W.sub(1), 0.0, mf, 100003)
bc = [bc_2noflow,bc_noflow, bc_os,bcp_outflow]

# weak formulation of the stokes problem with a left inflow and right outflow
a1 = mu*inner(grad(u),grad(v))*dx\
- div(v)*p*dx - q*div(u)*dx
L1 = -dot(1*v,n)*ds_G

# iterative solver parameter
prm = parameters['krylov_solver']
prm['absolute_tolerance'] = 1E-10
prm['relative_tolerance'] = 1E-6
prm['maximum_iterations'] = 1000


# resolution of the same problem with MINRES or TFQMR and amg
w1 = Function(W)
solve(a1 == L1, w1, bc, solver_parameters={'linear_solver':'minres','preconditioner': 'amg'})
(u1, p1) = w1.split()

w2 = Function(W)
solve(a1 == L1, w2, bc, solver_parameters={'linear_solver':'tfqmr','preconditioner': 'amg'})
(u2, p2) = w2.split()
2

Looks like you’re missing an appropriate preconditioner for the symmetric positive indefinite problem.

Consider the dolfin demo and dolfinx demo.

I’d strongly advise against using the syntactic sugar of solve(a == L, ...) when you need fine control over the linear algebra solution method.

2 Likes
3

Hi Nate,

Thank you really much for your answer : that’s working now !

I always thought that the two syntax were equivalent.

I have another question : Is there any interest in running in parallel with “assemble”? I see that it is quite slow in my case and I would like to speed up the process.

Thank you

4

Scalable performance is contingent on a well configured system with appropriate MPI bindings, in addition to careful consideration of an iterative solver’s design and a corresponding preconditioner.

The performance tests provide implementations of problems and corresponding solvers which are scalable. This is typically a good place to begin benchmarking your build and system.

Regarding performance of FEniCS reported in the literature, consider here and here.

1 Like