I am using the mixed-dimensional branch from @cdaversin via the docker container. I define a MixedNonlinearVariationalSolveraccording to the information in this post Mixed-dimensional branch - how to define NonlinearVariationalProblem?. I get a runtime error when I use the manually defined solver but not when I use the solve() function.
Definition of MixedNonlinearVariationalProblem
Below is the code I used to define the MixedNonlinearVariationalSolver, where F is the variational form und U the function.
F_blocks = df.extract_blocks(F)
Jacobian = []
for Fi in F_blocks:
for ui in U.split():
d = df.derivative(Fi, ui)
Jacobian.append(d)
problem = df.MixedNonlinearVariationalProblem(F_blocks, U.split(), bcs, J=Jacobian)
solver = df.MixedNonlinearVariationalSolver(problem)
MWE
Following MWE solves fine with solve() but gives a runtime error for the manually defined MixedNonlinearVariationalSolver.
import dolfin as df
# %% Meshes
# Create parent mesh
L1 = 1
L2 = 1
L = L1+L2
Nx = int(L*100)
mesh = df.IntervalMesh(Nx, 0, L)
# Create submeshes
marker = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for f in df.facets(mesh):
cond_left = 0 - df.DOLFIN_EPS < f.midpoint().x() < 1 + df.DOLFIN_EPS
# cond_right = 1 - df.DOLFIN_EPS < f.midpoint().x() < 2 + df.DOLFIN_EPS
if cond_left:
marker[f] = 1
else:
marker[f] = 2
submesh_left = df.MeshView.create(marker, 1)
submesh_right = df.MeshView.create(marker, 2)
# %% Markers
boundary_markers_0 = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
boundary_markers_1 = df.MeshFunction('size_t', submesh_left, submesh_left.topology().dim()-1, 0)
boundary_markers_2 = df.MeshFunction('size_t', submesh_right, submesh_right.topology().dim()-1, 0)
class Left(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], 0, df.DOLFIN_EPS)
class Right(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], 2, df.DOLFIN_EPS)
left = Left()
left.mark(boundary_markers_0, 1)
left.mark(boundary_markers_1, 1)
right = Right()
right.mark(boundary_markers_0, 2)
right.mark(boundary_markers_2, 2)
# %% Function spaces and measures
V0 = df.FunctionSpace(mesh, 'P', 2)
V1 = df.FunctionSpace(submesh_left, 'P', 2)
V2 = df.FunctionSpace(submesh_right, 'P', 2)
W = df.MixedFunctionSpace(V0, V1, V2)
dx0 = df.Measure("dx", domain=W.sub_space(0).mesh())
dx1 = df.Measure("dx", domain=W.sub_space(1).mesh())
dx2 = df.Measure("dx", domain=W.sub_space(2).mesh())
# %% Boundary conditions
bc_0 = df.DirichletBC(V0, df.Constant(0), boundary_markers_0, 1)
bc_1 = df.DirichletBC(V1, df.Constant(0), boundary_markers_1, 1)
bc_2 = df.DirichletBC(V2, df.Constant(0), boundary_markers_2, 2)
bcs = [bc_0, bc_1, bc_2]
# %% Variational problem
U = df.Function(W)
u0, u1, u2 = U.sub(0), U.sub(1), U.sub(2)
v0, v1, v2 = df.TestFunctions(W)
f1 = df.Constant(-20)
f2 = df.Constant(20)
# Uncoupled (works for both variants)
# F0 = u0.dx(0)*v0.dx(0)*dx0 - f1*v0*dx0
# Coupled (only works for solve())
F0 = u0.dx(0)*v0.dx(0)*dx0 - u1*v0*dx1 - u2*v0*dx2
F1 = u1.dx(0)*v1.dx(0)*dx1 -f1*v1*dx1
F2 = u2.dx(0)*v2.dx(0)*dx2 -f2*v2*dx2
F = F0+F1+F2
# %% Compute solution
# works for coupled and uncoupled equations
# df.solve(F == 0, U, bcs)
# gives Error for coupled equations
F_blocks = df.extract_blocks(F)
Jacobian = []
for Fi in F_blocks:
for ui in U.split():
d = df.derivative(Fi, ui)
Jacobian.append(d)
problem = df.MixedNonlinearVariationalProblem(F_blocks, U.split(), bcs, J=Jacobian)
solver = df.MixedNonlinearVariationalSolver(problem)
solver.solve()
Error Message
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 set given (local) rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------