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)