Constructing a Nonlinear Problem on DOLFIN for PETSc

Hello everyone,

I am using DOLFIN to solve a nonlinear problem with PETSc.SNES instead of NonlinearVariationalSolver to have more control over settings. I am constrained to use FEniCS (not FEniCSx) because my problem requires the real element (an element with only one degree of freedom across the entire domain).

Since there is nothing like assemble_matrix_block in DOLFIN, I need to permutate the PETSc matrix and vector in the constructor I pass to SNES (for applying fieldsplit in KSP). So I tryed to create a problem class by myself.

I wrote a basic version to solve NS to check if it works. The domain is a rectangle. I set a rightward inflow at the left inlet with a parabolic velocity profile. No-slip boundary conditions are applied on the top and bottom walls, and a pressure of zero is set on the right to approximate a do-nothing boundary condition. While the solution converges, there are some suspicious behaviors(the output in paraview is shown below):

  • The boundary conditions are flipped by a negative.
  • The pressure space remains zero.
  • Despite being a nonlinear problem, it converges in a single SNES iteration. Even reducing viscosity only increases the number of KSP iterations, but it always halts after one SNES iteration.

I think something is wrong in the problem class. Are there any demos available on how to create a problem class for PETSc in DOLFIN? Or is there any advise on my problem class? Any advice is appreciated.

  • output
2-norm of F:  7.313617309588403
snes iter: 0, error: 7.313617309588403
ksp iter: 0, error: 6.1512154267884736
ksp iter: 1, error: 0.11526454588477297
ksp iter: 2, error: 0.0013454238357176631
ksp iter: 3, error: 2.066540868355655e-05
ksp iter: 4, error: 5.599560254996777e-07
ksp iter: 5, error: 1.483005982258138e-08
2-norm of F:  0.0
snes iter: 1, error: 0.0

  • code
from dolfin import *
from petsc4py import PETSc
from mpi4py import MPI as mpi
from mpi4py.MPI import COMM_WORLD, COMM_SELF
import numpy as np
import math

# base for the output
base_path = "/path/to/output/"

 # Define the dimensions of the rectangle
x0, y0 = 0.0, -1.0  # Lower left corner
x1, y1 = 5.0, 1.0  # Upper right corner
nx, ny = 25, 10    # Number of cells in x and y directions

# Create the rectangle mesh
msh = RectangleMesh(Point(x0, y0), Point(x1, y1), nx, ny)


class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
    
class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0)
    
class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],5.)


## Create MeshFunction to mark the boundaries ##
boundaries = MeshFunction("size_t", msh, msh.topology().dim()-1)
boundaries.set_all(3)
Boundary().mark(boundaries, 0)
Inflow().mark(boundaries, 1)
Outflow().mark(boundaries, 2)

## Setup Taylor Hood Function Space ##
CG2 = VectorElement("CG", msh.ufl_cell(), 2)
CG1 = FiniteElement("CG", msh.ufl_cell(), 1)
W1 = MixedElement([CG2, CG1])
W = FunctionSpace(msh, W1)

## Create Trial, Test, and Solutions Functions ##
U  = Function(W)
Te  = TestFunction(W)
Tr = TrialFunction(W)
u, p = split(U)
v, q = split(Te)
du, dp = split(Tr)

## Set up No-slip BC on the walls ##
bcs = []
zero_1d = Constant(0.0)
zero_2d = Constant((0.0,0.0))
bcs.append(DirichletBC(W.sub(0), zero_2d, boundaries, 0))

## Set Lid-driven flow BC ##
inflow_exp = Expression(("1-x[1]*x[1]","0.0"),degree=2)
inflow = Function(W)
bcs.append(DirichletBC(W.sub(0), inflow_exp, boundaries, 1))

## Set do nothing on the outflow ##
bcs.append(DirichletBC(W.sub(1), zero_1d, boundaries, 2))

## Define residule ##
f = Constant((0.0, 0.0))
nu = Constant(0.02)

F_ufl = inner(dot(grad(u), u), v) * dx + nu * inner(grad(u), grad(v)) * dx - p * div(v) * dx - div(u) * q * dx + inner(f, v) * dx
J_ufl = derivative(F_ufl, U, Tr)

# potentially problematic
class NonLinearProblem():
    def __init__(self, U, F_ufl, J_ufl, P_ufl, bcs):
        self.U = U
        self.F_ufl = F_ufl
        self.J_ufl = J_ufl
        self.P_ufl = P_ufl
        self.bcs = bcs

        self.F = None
        self.J = None
        self.P = None
        
    def F_func(self, snes, x, F):
        local = as_backend_type(self.U.vector()).vec().getArray()
        local[:] = x
        as_backend_type(self.U.vector()).vec().ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        J_mat = assemble(self.J_ufl)
        F_vec = assemble(self.F_ufl)
        for bc in self.bcs:
            bc.apply(J_mat)
            bc.apply(F_vec)
        
        self.J = as_backend_type(J_mat).mat()
        self.J.assemblyBegin()
        self.J.assemblyEnd()

        self.b = as_backend_type(F_vec).vec()
        self.b.assemblyBegin()
        self.b.assemblyEnd()

        F.copy(self.b)

        print("2-norm of F: ", F.norm(PETSc.NormType.NORM_2))

    def J_func(self, snes, x, J, P):
        J.copy(self.J)

problem = NonLinearProblem(U, F_ufl, J_ufl, None, bcs)

# set snes
snes = PETSc.SNES().create(COMM_WORLD)
snes.setType("newtonls")
J_mat, F_vec = assemble_system(J_ufl, F_ufl, bcs)
snes.setFunction(problem.F_func, as_backend_type(F_vec).vec())
snes.setJacobian(problem.J_func, as_backend_type(J_mat).mat())
snes.setTolerances(rtol=1e-3, atol=1e-6, max_it=25)
snes.setMonitor(lambda snes, its, res: print(f"snes iter: {its}, error: {res}"))
snes.setFromOptions()

# set ksp
ksp = snes.getKSP()
ksp.setMonitor(lambda ksp, its, res: print(f"ksp iter: {its}, error: {res}"))
ksp.setTolerances(rtol=1e-8, atol=1e-9)
ksp.setType("gmres")
ksp.getPC().setType("hypre")
ksp.getPC().setHYPREType('boomeramg')

# solve and print #
X = Function(W)
x = as_backend_type(X.vector()).vec()
snes.solve(None, x)
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

u_, p_ = X.split()
ufile_pvd = File(f"{base_path}/result_u.pvd")
ufile_pvd << u_
ufile_pvd = File(f"{base_path}/result_p.pvd")
ufile_pvd << p_

del snes

The best publicly available problem class is presented by Prof. Francesco Ballarin in tutorial_navier_stokes

Thank you for your advice. I’ve tried that in DOLFINX, and it works well. However, I’m finding it challenging to translate the same approach into DOLFIN. Specifically, I don’t know exactly what the following lines of code do, nor have I been able to find a way to achieve similar functionality in DOLFIN:

dolfinx.fem.petsc.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)

dolfinx.fem.petsc.set_bc(F_vec, self._bcs, x, -1.0)

The source code of apply_lifting says that it “Modify b such that: b ← b - scale * A_j (g_j - x0_j)”, but there is no g or b elsewhere. Maybe it stands for the value of bc and x?

The source code of set_bc seems to set the value in F_vec to be scalar * bcs on bc, but I’m not sure either. I think there is a better tutorial (than source code and documentation) for them lies somewhere.

I’ve also search for the fenics version of multiphenicsx. Although they are good demos for DOLFIN, it seems doesn’t involves creating a problem for PETSc.

I have experience with dolfinx/multiphenicsx not with dolfin/multiphenics, I’d recommend looking for the BlockPETScSNESSolver used in 02_navier_stokes.

As for:

  • dolfinx.fem.petsc.apply_lifting adjusts the RHS vector by taking into account at block index j of matrix A; Dirichlet BC value g, using the specified scale and an optional initial guess x0.
  • dolfinx.fem.petsc.set_bc can be seen in DirichletBC::set

Thank you for your advise, and I would be grateful if you could explain the detail of dolfinx.fem.petsc.apply_lifting.

I found a tutorial on dolfinx (Application of Dirichlet boundary conditions — FEniCS Tutorial @ Sorbonne), and the equation 12 of it (a.k.a.

\begin{pmatrix} A_{d,d} & 0\\ 0 & I \end{pmatrix}\begin{pmatrix} u_{d}\\ u_{bc} \end{pmatrix} =\begin{pmatrix} b_{d} -A_{d,bc} g\\ g \end{pmatrix}

) seems describes the action of

dolfinx.fem.petsc.apply_lifting(b, [A], [g(or bcs)], x0=[0], scale=+1.0)

. Am I right? Then does

dolfinx.fem.petsc.apply_lifting(b, [A], [g(or bcs)], x0=[x], scale=-1.0)

stands for

\begin{pmatrix} b_{d} +A_{d,bc}( g-x_{0})\\ g \end{pmatrix}

or

\begin{pmatrix} b_{d} +A_{d,bc}( g-x_{0})\\ g-x_{0} \end{pmatrix}

or negative the second term, or neither of them?

I tried

# potentially problematic
class NonLinearProblem():
    def __init__(self, U, F_ufl, J_ufl, P_ufl, bcs):
        self.U = U
        self.F_ufl = - F_ufl
        self.J_ufl = J_ufl
        self.P_ufl = P_ufl
        self.bcs = bcs

        self.g = Function(W)
        for bc in bcs:
            bc.apply(self.g.vector())

        self.x_g = Function(W)
        self.bcs_x_g = [DirichletBC(W.sub(1), self.x_g.sub(1), boundaries, 2), DirichletBC(W.sub(0), self.x_g.sub(0), boundaries, 1), DirichletBC(W.sub(0), self.x_g.sub(0), boundaries, 0)]

        self.F = None
        self.J = None
        self.P = None
        
    def F_func(self, snes, x, F):
        as_backend_type(self.U.vector()).vec()[:] = x
        as_backend_type(self.U.vector()).vec().ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        as_backend_type(self.x_g.vector()).vec().copy(as_backend_type(self.U.vector()).vec())
        as_backend_type(self.x_g.vector()).vec().axpy(-1, as_backend_type(self.g.vector()).vec())

        J_mat, F_vec = assemble_system(self.J_ufl, self.F_ufl, self.bcs_x_g)

        self.J = as_backend_type(J_mat).mat()
        self.J.assemblyBegin()
        self.J.assemblyEnd()

        self.b = as_backend_type(F_vec).vec()
        self.b.assemblyBegin()
        self.b.assemblyEnd()

        F.copy(self.b)

        print("2-norm of F: ", F.norm(PETSc.NormType.NORM_2))

    def J_func(self, snes, x, J, P):
        J.copy(self.J)

, but the output is exactly the same. Does assemble_system in dolfin do the same thing as dolfinX’s assemble_vector_block? Or how can I apply_lifting in DOLFIN?

Something new happened.
I found that there was some mistake on my understanding of PETSc grammar. After I fixed it, the class becomes

class NonLinearProblem():
    def __init__(self, U, F_ufl, J_ufl, P_ufl, bcs):
        self.U = U
        self.F_ufl = - F_ufl
        self.J_ufl = J_ufl
        self.P_ufl = P_ufl
        self.bcs = bcs

        self.g = Function(W)
        for bc in bcs:
            bc.apply(self.g.vector())

        self.x_g = Function(W)
        self.bcs_x_g = [DirichletBC(W.sub(1), self.x_g.sub(1), boundaries, 2), DirichletBC(W.sub(0), self.x_g.sub(0), boundaries, 1), DirichletBC(W.sub(0), self.x_g.sub(0), boundaries, 0)]
        # self.bcs_x_g = [DirichletBC(W.sub(1), self.U.sub(1), boundaries, 2), DirichletBC(W.sub(0), self.U.sub(0), boundaries, 1), DirichletBC(W.sub(0), self.U.sub(0), boundaries, 0)]
        self.F = None
        self.J = None
        self.P = None
        
    def F_func(self, snes, x, F):
        x.copy(as_backend_type(self.U.vector()).vec())
        as_backend_type(self.U.vector()).vec().ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        as_backend_type(self.U.vector()).vec().copy(as_backend_type(self.x_g.vector()).vec())
        as_backend_type(self.x_g.vector()).vec().axpy(-1, as_backend_type(self.g.vector()).vec())
        as_backend_type(self.x_g.vector()).vec().ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        J_mat, F_vec = assemble_system(self.J_ufl, self.F_ufl, self.bcs_x_g)

        self.J = as_backend_type(J_mat).mat()
        self.J.assemblyBegin()
        self.J.assemblyEnd()

        self.b = as_backend_type(F_vec).vec()
        
        self.b.assemblyBegin()
        self.b.assemblyEnd()
        self.b.copy(F)
        
        print("2-norm of F: ", F.norm(PETSc.NormType.NORM_2))
        
    def J_func(self, snes, x, J, P):
        self.J.copy(J)

But strange thing still happens: the pressure space is fixed to zero by something. I tried to set the pressure of outflow by one, i.e.

## Set up No-slip BC on the walls ##
bcs = []
zero_1d = Constant(1.)
...
## Set do nothing on the outflow ##
bcs.append(DirichletBC(W.sub(1), zero_1d, boundaries, 2))
...

, and the phenomenon is shown in a more obvious way (the upper half shows the pressure, and the bottom half shows the velocity):


Any advice is appreciated.