### SIGSEGV fault when using PETSc FAS solver

39

views

0

I am solving a nonlinear variational problem using FEniCS and was experimenting with using the PETSc FAS solver. When I do so, I receive a segmentation fault. The problem is unique to FAS, as I am able to solve the problem with other PETSc SNES solvers. Code and output are posted below, and GDB backtracking output is available if that might be helpful. I've tried the code on a Linux operating system and on a Macbook with the same result.

```
D = Circle(Point(0.0, 0.0), 1.0, 100) # outer radius is 1.0
w = Circle(Point(0.0, 0.0), 0.4, 3000) # inner radius is 0.4
D = D - w
D = D - Rectangle(Point(-1.0, -1.0), Point(1.0, 0.0))
D = D - w
class Omega:
def inside(self, x):
return x[0] < DOLFIN_EPS and between(x[0]**2/(x[0]**2 + x[1]**2), (0.0, 0.5))
mu, rho, omega, Omega =0.1, 1.0, 500.0, Omega()
R1 = Expression((('0', '1'), \
('1', '0')), degree=2) # interchange x and y coordinates (reflect about y = x)
R2 = Expression((('0', '-1'), \
('1', '0')), degree=2) # rotate 90 degrees counterclockwise
mesh = generate_mesh(D, 100)
cc = CellFunction("bool", mesh, False)
circle = AutoSubDomain(lambda x: Omega.inside(x) or x[0]*x[0] + x[1]*x[1] < .2)
circle.mark(cc, True)
mesh = refine(mesh, cc)
n = FacetNormal(mesh)
class inflow_condition_0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and not Omega.inside(x) and (x[0]**2 + x[1]**2) < 0.2
class inflow_condition_1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and Omega.inside(x) and (x[0]**2 + x[1]**2) < 0.2
class radial_flow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] < DOLFIN_EPS_LARGE
class outflow_condition(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] ** 2 + x[1] ** 2 > .5 and (abs(x[1]) > DOLFIN_EPS or abs(x[0]) > 1.0 - DOLFIN_EPS)
outflow_condition = outflow_condition()
inflow_condition_0 = inflow_condition_0()
inflow_condition_1 = inflow_condition_1()
radial_flow = radial_flow()
class inflow_condition_2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and not inflow_condition_0.inside(x, on_boundary) and not
inflow_condition_1.inside(x, on_boundary) and (x[0] ** 2 + x[1] ** 2) < 0.3
inflow_condition_2 = inflow_condition_2()
boundary_markers = MeshFunction('size_t', mesh, 1) #create boundary markers for the four different boundaries
ds = Measure('ds', subdomain_data=boundary_markers)
outflow_condition.mark(boundary_markers, 1)
inflow_condition_0.mark(boundary_markers, 2)
inflow_condition_1.mark(boundary_markers, 3)
radial_flow.mark(boundary_markers, 4)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)
a = Expression(('0', '0'), degree=2)
inflow1 = DirichletBC(W.sub(0), a, inflow_condition_0, method="topological")
inflow3 = DirichletBC(W.sub(0), a, inflow_condition_2, method="topological")
a = Expression(('x[0]/.4', 'x[1]/.4'), degree=2)
inflow2 = DirichletBC(W.sub(0), a, inflow_condition_1, method="topological")
a = Expression('0.0', degree=1)
outflow = DirichletBC(W.sub(1), a, outflow_condition, method="topological")
a = Expression(('0', '0'), degree=2)
radialflow = DirichletBC(W.sub(0), a, radial_flow, method="topological")
bcu = [inflow1, inflow2, inflow3, radialflow]
bcp = [outflow]
bcs = bcu + bcp
p1 = Expression('alpha_U', alpha_U=1.0 degree=1, domain=mesh)
p2 = Expression('alpha_L', alpha_L=0.0, degree=1, domain=mesh)
class alpha(Expression):
def eval(self, values, x):
if not Omega.inside(x):
values[0] = p1(x)
else:
values[0] = p2(x)
alpha = alpha(degree = 1)
s = Function(W)
u, p = split(s)
v, q = TestFunctions(W)
r = Expression(('x[0]', 'x[1]'), degree=2, element=W.sub(0).ufl_element())
a = mu*inner(sym(grad(u)), sym(grad(v)))*dx + rho*dot(grad(u)*u, v)*dx \
+ alpha*inner(u, v)*dx - 2*rho*inner(omega*R2*u, v)*dx - p*div(v)*dx + q*div(u)*dx
L = rho*omega*omega*inner(r, v)*dx
F = a - L
tstart = time()
J = derivative(F, s)
problem = NonlinearVariationalProblem(F, s, bcs=bcs, J=J)
PETScOptions().set("snes_type", 'fas')
solver = NonlinearVariationalSolver(problem)
solver.parameters['nonlinear_solver'] = 'snes'
solver.parameters['snes_solver']['linear_solver'] = 'mumps'
PETScOptions().set('snes_fas_monitor')
PETScOptions().set('fas_levels_snes_', 'ngs')
PETScOptions().set('npc_fas_levels_snes_gs_sweeps', 3)
solver.solve()
```

Community: FEniCS Project

Please login to add an answer/comment or follow this question.