Iterative method (Krylov solver and preconditioner) for stationary N-S equation with high Reynolds number

Question

I am trying to solve stationary N-S equations with Re=5*10^5 and the equations were solved successfully with a mesh consist of 1million triangles using the direct solver. However, it costs quite a lot computer memory (about 20GB). Thus, I tried various Krylov solvers and preconditioners but they never converge. My code is like

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from numpy import linalg as la
from numpy.random import rand

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Inlet(SubDomain):
    def inside(self, x, on_boundary):
        tol=1e-14
        return on_boundary and near(x[0], -0.5,tol)

class Outlet(SubDomain):
    def inside(self, x, on_boundary):
        tol=1e-14
        return on_boundary and near(x[0], 1.25,tol)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        tol=1e-14
        return on_boundary and (near(x[1], 0.0,tol))# and between(x[0], (-0.5, 0.0+tol)))

class Top(SubDomain):
    def inside(self, x, on_boundary):
        tol=1e-14
        return on_boundary and near(x[1], 0.2,tol)

class Plate(SubDomain):
    def inside(self, x, on_boundary):
        tol=1e-14
        return on_boundary and (near(x[1], 0.0,tol) and between(x[0], (0.0-tol, 1.25)))

# Initialize sub-domain instances
inlet = Inlet()
top = Top()
outlet = Outlet()
bottom = Bottom()
plate = Plate()

# Create mesh
mesh=Mesh("plate_500000.xml")
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
inlet.mark(boundaries, 1)
top.mark(boundaries, 2)
outlet.mark(boundaries, 3)
bottom.mark(boundaries, 4)
plate.mark(boundaries, 5)
# Define new measures associated with the interior domains and
# exterior boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = MixedElement([P2,P1])
V = FunctionSpace(mesh, TH)

# Create functions
(v, q) = TestFunctions(V)
dw = TrialFunction(V)
w = Function(V)
(u, p) = split(w)

#Create parameters
nu=Constant(1.0/(5.0*10**5))
R_noslip=Constant((0.0,0.0))
R_inlet=Constant((1.0,0.0))
R_symm=Constant(0.0)

#Boundary connditions
BC_inlet=DirichletBC(V.sub(0), R_inlet, boundaries, 1)
BC_noslip=DirichletBC(V.sub(0), R_noslip, boundaries, 5)
BC_top=DirichletBC(V.sub(0).sub(1), R_symm, boundaries, 2)
BC_bottom=DirichletBC(V.sub(0).sub(1), R_symm, boundaries, 4)
bcs=[BC_inlet,BC_noslip,BC_bottom,BC_top]

# Define variational form
n=FacetNormal(mesh)
F=(nu*inner(grad(u),grad(v))+inner(grad(u)*u,v)-div(v)*p+q*div(u))*dx-nu*inner(grad(u)*n,v)*ds(2)-nu*inner(grad(u)*n,v)*ds(4)+dot(v,n)*p*ds(2)+dot(v,n)*p*ds(4)

solve(FF == 0, w, bcs,solver_parameters={'symmetric':True,'newton_solver':
    {'linear_solver':'gmres','preconditioner':'hypre_amg',
    'krylov_solver':{'absolute_tolerance':1e-9,'relative_tolerance':1e-7,'maximum_iterations':5000,
    'nonzero_initial_guess':False,'error_on_nonconvergence':True,'monitor_convergence':True,'report':False,'divergence_limit':10000.0}}})
#solve(F == 0, w, bcs,solver_parameters={'newton_solver':{'linear_solver':'mumps',"relative_tolerance": 1e-10,"relaxation_parameter": 1.0}})

The mesh can be downloaded by the link https://www.dropbox.com/s/2pbgelnt73wrlox/plate_500000.xml?dl=0

Comments

I found a package called FENaPack which implements preconditioners for Navier-Stokes problem using FEniCS and PETSc packages. It may help solve N-S equations using Krylov solver.

Answer 1

It seems that you solve the N.S. equations in a coupled manner. I do not know for sure, but I believe, the system matrices resulting from this problem are typically not symmetric. However, you do set some setting to symmetric (even though the gmres solver should be suitable for non-symmetric problems). Does it help if you change this to 'symmetric': False?

Comments

Still not converge. The preconditioned residual norm is converged but the solution is still not converged. I am not sure if the method or the preconditioner is suitable for this nonlinear problem. I tried the same solver setting for a nonlinear Poisson equation and it worked. I was wondering if the problem is caused by the Mixed Function Space.

---

I just remembered that I once tried to solve a similar problem (i.e. stationary coupled Stokes) and also could not solve the problem using Krylov solvers (note that the Stokes problem does not require the Newton iterations) . I finally solved it by following this fenics demo. The documentation of this problem, says:

"The Stokes equations as formulated above result in a system of linear equations that is not positive-definite. Standard iterative linear solvers typically fail to converge for such systems. Some care must therefore be taken in preconditioning the systems of equations. Moreover, not all of the linear algebra backends support this. We therefore start by checking that either "PETSc" or "Tpetra" (from Trilinos) is available. We also try to pick MINRES Krylov subspace method which is suitable for symmetric indefinite problems. If not available, costly QMR method is choosen."

I can imagine that you face a similar problem for the full Navier-Stokes equations., but you'd have to check.