SIGSEGV fault when using PETSc FAS solver


39
views
0
4 weeks ago by
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.

Similar posts:
Search »