I’m solving NS using dolfin. I try to use PETSc directly, but the dofs in pressure space doesn’t change for no reason, making the solver not converge. I turn the grid size small, and set pressure to 1 at the outlet to make the phenomenon obvious.
Any advice would be appreciated.
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
# Define the dimensions of the rectangle
x0, y0 = 0.0, 0.0 # Lower left corner
x1, y1 = 1.0, 1.0 # Upper right corner
nx, ny = 2, 2 # 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.)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0],1.)
## 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(), 1)
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(1.)
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(("4*(1-x[1])*x[1]","0.0"),degree=2)
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(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)
dof_map_p = np.array(W.sub(1).dofmap().dofs(), dtype=np.int32)
P_sub_perm = PETSc.IS().createGeneral(dof_map_p, comm=PETSc.COMM_SELF)
# question class
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
def F_func(self, snes, x, F):
x_p = x.getSubVector(P_sub_perm)
x_p.assemble()
x_p.view()
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)
J_mat, F_vec = assemble_system(self.J_ufl, self.F_ufl, self.bcs, self.U.vector())
self.J = as_backend_type(J_mat).mat()
self.J.assemble()
self.b = as_backend_type(F_vec).vec()
self.b.assemble()
self.b.copy(F)
print("2-norm of F: ", F.norm(PETSc.NormType.NORM_2))
def J_func(self, snes, x, J, P):
J.zeroEntries()
self.J.copy(J)
J.assemble()
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, max_it = 100)
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
output
Vec Object: 1 MPI process
type: seq
0.
0.
0.
0.
0.
0.
0.
0.
0.
2-norm of F: 4.815064324019405
snes iter: 0, error: 4.815064324019405
ksp iter: 0, error: 5.473145554945983
ksp iter: 1, error: 0.4251872374915494
ksp iter: 2, error: 5.156831346466621e-16
Vec Object: 1 MPI process
type: seq
0.
0.
0.
0.
0.
1.
0.
1.
1.
2-norm of F: 3.1176301986533823
snes iter: 1, error: 3.1176301986533823
ksp iter: 0, error: 7.435868578230604
ksp iter: 1, error: 1.057990131117207
ksp iter: 2, error: 0.05106642125695564
ksp iter: 3, error: 1.2664526980287813e-14
Vec Object: 1 MPI process
type: seq
0.
0.
0.
0.
0.
1.
0.
1.
1.
2-norm of F: 0.9382567294507129
snes iter: 2, error: 0.9382567294507129
ksp iter: 0, error: 1.9018504441501722
ksp iter: 1, error: 0.04137687969950201
ksp iter: 2, error: 0.03872479111499394
ksp iter: 3, error: 2.1429443878265707e-15
Vec Object: 1 MPI process
type: seq
0.
0.
0.
0.
0.
1.
0.
1.
1.
2-norm of F: 0.6499460229653232
snes iter: 3, error: 0.6499460229653232
ksp iter: 0, error: 0.700620163254973
ksp iter: 1, error: 0.06945619829053755
ksp iter: 2, error: 0.0002784835850249214
ksp iter: 3, error: 6.01204069947197e-18
Vec Object: 1 MPI process
type: seq
0.
0.
0.
0.
0.
1.
0.
1.
1.
2-norm of F: 0.7919481372965301
...(continue to linesearch until it get to the max iteration)