Solving Stokes equation with NewtonSolver

jegor.kravchenko · November 29, 2023, 5:24pm

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)

nate · November 29, 2023, 6:48pm

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"

jegor.kravchenko · November 29, 2023, 6:53pm

jegor.kravchenko:

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()

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

dokken · November 29, 2023, 7:00pm

There are several issues with your implementation

The sign of the second term in the variational formulation. You have integrated both terms by part, so it should have the opposite sign
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)

dokken · November 29, 2023, 7:02pm

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

jegor.kravchenko · November 29, 2023, 7:16pm

Thank you so much for fixing my issues again

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.

dokken · November 29, 2023, 7:20pm

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).

Topic		Replies	Views
RuntimeError: Newton solver did not converge because maximum number of iterations reached Errors	1	2124	April 20, 2022
Krylov solver's option max_it and a few questions General	25	601	November 26, 2023
Newton Solver Does not converge in mixed formulation Errors error	6	291	December 20, 2023
Looking for documentation on nonlinear solver dolfinx	4	262	October 23, 2023
How to use solution of linear problem as initial guess of non linear problem General	7	309	October 18, 2023

Solving Stokes equation with NewtonSolver

Related Topics