Hi all,

I’m trying to solve Stokes equation using an iterative solver. I want to achieve this in the framework of Newton solver as, later, I want to couple stokes equation with equations of poroelasticity and Darcy flow. I found a demo by @nate which shows how to build a custom Newton solver and add preconditioner to it. When I tried to apply it to the stokes equation, I get segmentation fault. The mesh size I’m using shouldn’t cause any memory issues since it’s not really large mesh. What could be an issue then?

Thanks for all the help! I’m pretty new to the FEniCS and I feel like I’m missing something important here in building a custom solver.

code-

```
#importing required libraries
from fenics import *
class Problem(NonlinearProblem):
def __init__(self, J, J_pc, F, bcs):
self.bilinear_form = J
self.preconditioner_form = J_pc
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)
def F(self, b, x):
pass
def J(self, A, x):
pass
def form(self, A, P, b, x):
assemble(self.linear_form, tensor=b)
assemble(self.bilinear_form, tensor=A)
assemble(self.preconditioner_form, tensor=P)
for bc in self.bcs:
bc.apply(b,x)
bc.apply(A)
bc.apply(P)
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(), PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operators(A,P)
PETScOptions.set("ksp_type","minres")
PETScOptions.set("pc_type","hypre")
PETScOptions.set("ksp_view")
self.linear_solver().set_from_options()
#constants
beta = Constant(1e-8)
#defining the mesh
mesh = RectangleMesh(Point(0,0),Point(10,1),40,20)
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
P2 = VectorElement('P', mesh.ufl_cell(), 2)
element = MixedElement([P2,P1])
V = FunctionSpace(mesh, element)
v1, v2 = TestFunctions(V)
u = Function(V)
dv = TrialFunction(V)
u1, u2 = split(u); #splitting u
#defining boundaries
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
class BoundaryL(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, DOLFIN_EPS)
bxL = BoundaryL()
bxL.mark(boundary_markers, 0)
class BoundaryR(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 10, DOLFIN_EPS)
bxR = BoundaryR()
bxR.mark(boundary_markers, 2)
class BoundaryD(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, DOLFIN_EPS)
bxD = BoundaryD()
bxD.mark(boundary_markers, 1)
class BoundaryU(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1, DOLFIN_EPS)
bxU = BoundaryU()
bxU.mark(boundary_markers, 3)
#boundary conditions
bc2 = DirichletBC(V.sub(1), Constant(0.), bxD)
bc3 = DirichletBC(V.sub(1), Constant(0.), bxL)
bc4 = DirichletBC(V.sub(1), Constant(0.), bxR)
inflow = Expression(("0.","sin(x[0]*pi)"), degree=2)
bc1 = DirichletBC(V.sub(0), inflow, bxU)
bcs = [bc1,bc2,bc3,bc4]
#weak statement of the equations
F = beta*inner(grad(u1),grad(v1))*dx - div(v1)*u2*dx + v2*div(u1)*dx
J = derivative(F,u,dv)
J_pc = inner(grad(u1),grad(v1))*dx + u2*v2*dx
# Define solver
problem = Problem(J, J_pc, F, bcs)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
```

error message-

```
Segmentation fault (core dumped)
```