I am trying to solve the 2D Diffusion equation with a 1D differential equation governing the flux through the bottom wall. To solve this I am using the cecile/mixed-dimensional branch installed through docker (tag 27-03-20) .
Problem:
In 2D medium:
\frac{\partial c}{\partial t} = \frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2}Bottom Wall:
\frac{\partial b}{\partial t} = c_{y=0} (1-b) - b\Gamma_{bottom} : \frac{\partial b}{\partial t} = -\frac{\partial c}{\partial y}Other Walls:
\Gamma_{walls} : \mathbf{n}\cdot \nabla c = 0
Variational Formulation:
F = \int_{\Omega} \left[ \frac{1}{\Delta t}\left(c-c^n \right) v_c +\frac{\partial c}{\partial x}\frac{\partial v_c}{\partial x} + \frac{\partial c}{\partial y}\frac{\partial v_c}{\partial y} \right] d\Omega - \int_{\Gamma_s} \frac{\partial c}{\partial y} v_c ~d\Gamma = 0F \text{ +=} \int_{\Gamma_s} \frac{1}{\Delta t}(b-b^n) w - g(c,b) w ~d\Gamma = 0where,
\frac{\partial c}{\partial y} = -g(c,b) \qquad , \qquad g(c,b) = c(1-b)-b
MWE:
from dolfin import *
mesh = UnitSquareMesh(10, 10)
# Domains and Boundaries
domains = MeshFunction('size_t', mesh, mesh.topology().dim()) ; domains.set_all(0)
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1) ; boundaries.set_all(0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0, DOLFIN_EPS) and on_boundary
Bottom = Bottom() ; Bottom.mark(domains,1) ; Bottom.mark(boundaries,1)
# SubMesh
marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
Bottom.mark(marker,1)
SubMesh = MeshView.create(marker,1)
domains_submesh = MeshFunction('size_t', SubMesh, SubMesh.topology().dim(), 0)
# Define the mixed function space
Q1 = FunctionSpace(mesh, "Lagrange", 1)
Q2 = FunctionSpace(SubMesh, "Lagrange", 1)
W = MixedFunctionSpace(Q1,Q2)
# Measures
dV = Measure("dx", domain=W.sub_space(0).mesh(), subdomain_data = domains)
dL = Measure("dx", domain=W.sub_space(1).mesh(), subdomain_data = domains_submesh)
# Initial Conditions
c_n = interpolate(Constant(1), W.sub_space(0))
b_n = interpolate(Constant(0), W.sub_space(1))
# Define test functions
(v_c , v_b) = TestFunctions(W)
# Functions at Current time
c = Function(W.sub_space(0))
b = Function(W.sub_space(1))
# Constants
t = 0 ; dt = 0.001
Δt = Constant(dt)
def g(c,b): return c*(1-b)-b
# Variational Form
F = (c-c_n)/Δt*v_c*dV + dot(grad(c),grad(v_c))*dV + g(c,b)*v_c*dL
F += (b-b_n)/Δt*v_b*dL - g(c,b)*v_b*dL
sol = Function(W)
# solve(F==0, sol, bcs=[], solver_parameters= {"nonlinear_solver": "newton",
# "newton_solver": {"linear_solver": "mumps"}})
# Manually Define Problem
## Calculate Jacobian Matrix
# from ufl.algorithms import expand_derivatives
F_blocks = extract_blocks(F)
Jacobian = []
for Fi in F_blocks:
for ui in sol.split():
d = derivative(Fi, ui)
# d = expand_derivatives(d) # results in Jacobian =[form[],form[],form[],form[]]
Jacobian.append(d)
## Define Problem and Solver
problem = MixedNonlinearVariationalProblem(F_blocks, sol.split(), bcs=[], J=Jacobian)
solver = MixedNonlinearVariationalSolver(problem)
## Solver Parameters
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['linear_solver'] = 'mumps'
prm['newton_solver']['convergence_criterion'] = 'residual'
solver.solve()
Error:
fenics@73585bea6869:~/shared$ python3 MWE.py
[problem] create list of residual forms OK
Calling FFC just-in-time (JIT) compiler, this may take some time.
[problem] create list of jacobian forms OK, J_list size = 4
[MixedNonlinearVariationalProblem]::check_forms - To be implemented
Solving mixed nonlinear variational problem.
Newton iteration 0: r (abs) = 1.559e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File "MWE.py", line 82, in <module>
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. ...
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------
Additional:
The reason I am using mumps is because of the following error, which arises for the default linear_solver:
UMFPACK V5.7.1 (Oct 10, 2014): ERROR: required argument(s) missing
Traceback (most recent call last):
File "MWE.py", line 81, in <module>
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an 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.
*** -------------------------------------------------------------------------
I’d be very grateful for any help. Thanks in advance!