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
