Segmentation fault while solving Stokes equation

Errors
1

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)
2

The problem is that you have not defined your preconditioner as a matrix, but as a vector

It should be:

dv1, dv2 = split(dv)  # splitting dv
J_pc = inner(grad(dv1), grad(v1))*dx + dv2*v2*dx

Please note that your system is not converging, even with:

      PETScOptions.set("ksp_type", "preonly")
      PETScOptions.set("pc_type", "lu")
      PETScOptions.set("ksp_view")
3

Thank you @dokken! With modifications, the code works although the solutions are not converging. I think the issue is the sign of the pressure term in the weak form.

F = beta*inner(grad(u1),grad(v1))*dx - div(v1)*u2*dx + v2*div(u1)*dx

If I chose it to be positive as shown in the demo of the Stokes equation on the FEniCS website then the solutions are convergent.

F = beta*inner(grad(u1),grad(v1))*dx + div(v1)*u2*dx + v2*div(u1)*dx

Although I’m not sure if changing the sign of the pressure term as shown in the demo implies right physical picture.