Solving Stokes equation with NewtonSolver

dolfinx
1

Hello everyone,
I have trouble solving the Stokes equation with the NewtonSolver. As far as I understand, the Newton solver should converge in one iteration, since the Stokes equations are linear. I do plan on extending to a non linear model in the future though, so that’s why I try to use the Newton solver in the first place. I tried to use the same settings for the krylov solver as in this post from the forum Incompressible Navier-Stokes with particular scheme & mixed elements FEniCSx, but with no success.

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

I also printed solver.A[:,:] and solver.b[:], which I assume are the stiffness matrix and the residual in my case. The matrix A seems to be fine, but the vector b contains some pretty high numbers like 1e+15 and they grow with each Newton iteration. I’m not sure, if I messed something up with my mixed finite elements or whether my settings for the krylov solver are not right for my problem.

Thank you in advance.

Here is my code:

import numpy as np
import ufl
from ufl import inner, grad, dx, div, VectorElement, FiniteElement, TestFunctions
import dolfinx
from dolfinx import default_scalar_type
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc

# %% mesh and function spaces
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (1.0, 1.0)],
                          (3,3), mesh.CellType.quadrilateral)

# combined function space, taylor hood elements with 2nd order for momentum
velocity_element = VectorElement("Lagrange", domain.ufl_cell(), degree=2)
pressure_element = FiniteElement("Lagrange", domain.ufl_cell(), degree=1)
space = fem.functionspace(domain, velocity_element*pressure_element)

v, q = TestFunctions(space)

# w should hold the solution
w = fem.Function(space)
u,p = ufl.split(w)

# %% boundary conditions
U = space.sub(0).collapse()[0]
Q = space.sub(1).collapse()[0]

def left(x):
    return np.isclose(x[0], 0.0)
def right(x):
    return np.isclose(x[0], 1.0)
def top(x):
    return np.isclose(x[1], 1.0)
def bottom(x):
    return np.isclose(x[1], 0.0)
def corner(x):
    return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))

bc = [] # list of Dirichlet boundary conditions

# velocities for cavity flow
noslip = fem.Constant(domain, (0.0, 0.0))
dofs = fem.locate_dofs_geometrical(U, left)
bc.append(fem.dirichletbc(noslip, dofs, U))
dofs = fem.locate_dofs_geometrical(U, right)
bc.append(fem.dirichletbc(noslip, dofs, U))
dofs = fem.locate_dofs_geometrical(U, bottom)
bc.append(fem.dirichletbc(noslip, dofs, U))
dofs = fem.locate_dofs_geometrical(U, top)
bc.append(fem.dirichletbc(fem.Constant(domain, (1.0, 0.0)), dofs, U))

# fix pressure in one corner
dofs = fem.locate_dofs_geometrical(Q, corner)
bc.append(fem.dirichletbc(fem.Constant(domain, 0.0), dofs ,Q))

# %% solve the problem
# residual of weak form of Stokes
F = inner(grad(u), grad(v))*dx + p*div(v)*dx + div(u)*q*dx
problem = NonlinearProblem(F, w, bc)
solver = NewtonSolver(domain.comm, problem)

solver.rtol = 1e-6
solver.report = True
solver.error_on_nonconvergence = False
solver.max_it = 1

info = solver.solve(w)
print(info)
2

What happens if you use an LU solver instead of CG preconditioned with GAMG? E.g.,

opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
1 Like
3

Newton doesn’t converge with these settings as well. Is there somewhere documentation on what all the possible settings for PETSc.Options are?

4

There are several issues with your implementation

  1. The sign of the second term in the variational formulation. You have integrated both terms by part, so it should have the opposite sign
  2. Setting of boundary conditions. You are setting them on a collapsed sub space, not the mixed space.
    Here is a functional version of the code:
import numpy as np
import ufl
from ufl import inner, grad, dx, div, VectorElement, FiniteElement, TestFunctions
import dolfinx
from dolfinx import default_scalar_type
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc

# %% mesh and function spaces
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (1.0, 1.0)],
                               (15, 15), mesh.CellType.quadrilateral)

# combined function space, taylor hood elements with 2nd order for momentum
velocity_element = VectorElement("Lagrange", domain.ufl_cell(), degree=2)
pressure_element = FiniteElement("Lagrange", domain.ufl_cell(), degree=1)
space = fem.functionspace(domain, velocity_element*pressure_element)

v, q = TestFunctions(space)

# w should hold the solution
w = fem.Function(space)
u, p = ufl.split(w)

# %% boundary conditions
U, _ = space.sub(0).collapse()
Q = space.sub(1)


def left(x):
    return np.isclose(x[0], 0.0)


def right(x):
    return np.isclose(x[0], 1.0)


def top(x):
    return np.isclose(x[1], 1.0)


def bottom(x):
    return np.isclose(x[1], 0.0)


def corner(x):
    return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))


bc = []  # list of Dirichlet boundary conditions

# velocities for cavity flow
noslip = dolfinx.fem.Function(U)
noslip.x.array[:] = 0
fdim = domain.topology.dim-1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, left_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, right_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, bottom_facets)
bc.append(fem.dirichletbc(noslip, dofs, space.sub(0)))
top_facets = mesh.locate_entities_boundary(domain, fdim, top)
dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, top_facets)
slip = dolfinx.fem.Function(U)
slip.interpolate(lambda x: (np.ones(x.shape[1]), np.zeros(x.shape[1])))
bc.append(fem.dirichletbc(slip, dofs, space.sub(0)))

# fix pressure in one corner
corner_vertex = mesh.locate_entities_boundary(domain, 0, corner)
dofs = fem.locate_dofs_topological(Q, 0,  corner_vertex)
bc.append(fem.dirichletbc(fem.Constant(domain, 0.0), dofs, Q))

# %% solve the problem
# residual of weak form of Stokes

F = inner(grad(u), grad(v))*dx - p*div(v)*dx - div(u)*q*dx
problem = NonlinearProblem(F, w, bc)
solver = NewtonSolver(domain.comm, problem)

solver.rtol = 1e-6
solver.report = True
solver.error_on_nonconvergence = False
solver.max_it = 1
info = solver.solve(w)
print(info)

with dolfinx.io.VTXWriter(domain.comm, "u.bp", [w.sub(0).collapse()], engine="BP4") as vtx:
    vtx.write(0.0)

image

2 Likes
5

Documentation of PETSc options can for instance be found at: KSP: Linear System Solvers — PETSc v3.20.1-276-g06d7d4ca38 documentation

1 Like
6

Thank you so much for fixing my issues again :smiley:

I have two questions about your implementation. The first one relates to this line

U, _ = space.sub(0).collapse()
Q = space.sub(1)

Is collapse only necessary when the subspace is based on a VectorElement?

dofs = fem.locate_dofs_topological((space.sub(0), U), fdim, left_facets)

Why do I pass a (space.sub(0),U) as a functions space base? For the pressure you specified only Q.

7

The reason for using collapse is related to the second question really.
In DOLFINx, there are some limitations of when one can use Constants as input to dirichlet bc.
One of them is when using a mixed space where you want to enforce a bc on a sub space where you are using a vector element. For this reason, one has to use a Function.

Since you preferably want to use a function from the sub space, we collapse the sub space, and interpolate the boundary condition into the collapsed sub space. This space has no connection to the parent space, and therefore we pass both the collapsed space and the sub space into the fem.locate_dofs_topological, which then returns a tuple, mapping bc dofs from U to space.sub(0).

2 Likes