AdaptiveNonlinearVariationalSolver does not pass solver parameters


242
views
1
5 months ago by
Hi,
this seems similar to this resolved issue, however it happens in AdaptiveNonlinearVariationalSolver instead of NonlinearVariationalSolver. Essentially, setting

solver2 = AdaptiveNonlinearVariationalSolver(problem2, M)
solver2.parameters['nonlinear_variational_solver']["newton_solver"]["linear_solver"] = "gmres"

does not change the solver in the simulation and I end up with the UMFPACK out of memory error:

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "mwe.py", line 98, in <module>
    solver2.solve(ATOL)
  File "/usr/lib/python3/dist-packages/dolfin/fem/adaptivesolving.py", line 124, in solve
    cpp.AdaptiveNonlinearVariationalSolver.solve(self, tol)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2017.2.0
*** Git changeset:  4c59bbdb45b95db2f07f4e3fd8985c098615527f
*** -------------------------------------------------------------------------

Aborted (core dumped)​

The best MWE I could come up with based on the problem I'm trying to solve is this:
from fenics import *

parameters["allow_extrapolation"] = True
parameters["refinement_algorithm"] = "plaza_with_parent_facets"

# Parameters
TOL = 1e-6
ATOL = 1e-10
el = tetrahedron
Kk1 = Constant(0.1) # permeability compartment 1
Kk2 = Constant(0.2) # permeability compartment 2
Kk3 = Constant(0.15) # permeability compartment 3
b12 = b13 = b23 = Constant(0.1) # compartment exchange coefficients
s1 = s2 = s3 = Constant(0.0) # porous media source terms
g = Constant(0.0) # Neumann boundary condition 
r1 = Constant(0.0) # decay rate
D = Constant(0.01) # diffusion coefficient
Pec = Constant(2.0)
dt = 1e-3

# Create mesh
nx = 10
mesh = UnitCubeMesh(nx, nx, nx)

# Define function space
P1 = FiniteElement('P', el, 1) # pressure space
Elem = MixedElement([P1, P1, P1])
W1 = FunctionSpace(mesh, Elem)
W2 = FunctionSpace(mesh, 'P', 1) # concentration space

# Define boundary conditions
def inject_p(x, on_boundary):
    return near(x[0], 0) and near(x[1], 0) and near(x[2], 0)

def outlet_p(x, on_boundary):
    return near(x[0], 1) and near(x[1], 1) and near(x[2], 1)

t = 1
pD = Expression('120*t', degree=1, t=t)
bcp = [DirichletBC(W1.sub(0), pD, inject_p, method='pointwise'),
        DirichletBC(W1.sub(0), Constant(0.0), outlet_p, method='pointwise')]
bcc = [DirichletBC(W2, Constant(1.0), inject_p, method='pointwise')]

# Define variational problem
q1, q2, q3 = TestFunctions(W1)
pp = Function(W1)
p1, p2, p3 = split(pp)
v1 = TestFunction(W2)
c1 = Function(W2)
c1_n = interpolate(Constant(0.0), W2)

# Permeability tensor
d = p1.geometric_dimension()
I = Identity(d)
K1 = Kk1*I
K2 = Kk2*I
K3 = Kk3*I

# Variational problem
# Darcy's law compartment 1
F1 = dot(grad(q1), K1*grad(p1))*dx + q1*(b12*(p1-p2) + b13*(p1-p3))*dx - q1*s1*dx + q1*g*ds
# Darcy's law compartment 2
F1 += dot(grad(q2), K2*grad(p2))*dx + q2*(b12*(p2-p1) + b23*(p2-p3))*dx - q2*s2*dx + q2*g*ds
# Darcy's law compartment 3
F1 += dot(grad(q3), K3*grad(p3))*dx + q3*(b13*(p3-p1) + b23*(p3-p2))*dx - q3*s3*dx + q3*g*ds
# Reaction/diffusion/advection kinetics species 1
F2 = v1*(1/dt)*(c1-c1_n)*dx + v1*dot((-K1*grad(p1)), grad(c1))*dx +\
        (1/Pec)*dot(grad(v1), grad(c1))*dx + v1*r1*c1*dx


# Solver for porous problem
J1 = derivative(F1, pp)
problem1 = NonlinearVariationalProblem (F1, pp, J=J1, bcs=bcp)
solver1 = NonlinearVariationalSolver(problem1)
solver1.parameters["newton_solver"]["linear_solver"] = "gmres"
solver1.parameters["newton_solver"]["absolute_tolerance"] = TOL
solver1.parameters["newton_solver"]["relative_tolerance"] = TOL

# Solver for solute tracking
J2 = derivative(F2, c1)
M = c1*dx
problem2 = NonlinearVariationalProblem (F2, c1, J=J2, bcs=bcc)
solver2 = AdaptiveNonlinearVariationalSolver(problem2, M)
solver2.parameters['nonlinear_variational_solver']["newton_solver"]["linear_solver"] = "gmres"
solver2.parameters['nonlinear_variational_solver']["newton_solver"]["absolute_tolerance"] = TOL
solver2.parameters['nonlinear_variational_solver']["newton_solver"]["relative_tolerance"] = TOL

# Solve porous problem
solver1.solve()
p1_n, p2_n, p3_n = pp.split()

solver2.solve(ATOL)​

Is there a workaround that I could possibly use?

Thanks!

Community: FEniCS Project

1 Answer


4
5 months ago by
I have been able to sort this issue out. The adaptive solver has a separate solver for error control, which needs to be set in addition to the problem solver:

solver2 = AdaptiveNonlinearVariationalSolver(problem2, M)
solver2.parameters["nonlinear_variational_solver"]["newton_solver"]["linear_solver"] = "gmres"
solver2.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "gmres"​
Please login to add an answer/comment or follow this question.

Similar posts:
Search »