Mixed-dimensional nonlinear solver minimal working example

General
1

I’ve been trying to work with the mixed-dimensional branch by @cdaversin, and have been struggling get a minimal working example of a solution to a nonlinear coupled equation using this branch.

Starting from the example script in /demo in the latest branch of `mixed-dimensional’, I am trying merely to change the script to use the nonlinear solver instead of the linear one (I realize this is useless, I just want to make a minimal working example, it’s doesn’t matter what the problem is).

The code listed below generates this 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:  68d9ed51bee6723163e7c4fbfffe31e6e3079bca
*** -------------------------------------------------------------------------

Is there something I should do with any preconditioners? Is something wrong with how the boundaries and internal constraint is set up? Why does using solve(a==L) but otherwise using the exact same code work just fine?

All help would be greatly appreciated! :slight_smile:

Here is the code:

from dolfin import *

### Step 1 : Build the meshes ###

# Parent mesh (Omega) : very basic unit square
n = 30
mesh = UnitSquareMesh(n, n)

# Submesh (Gamma, for the Lagrange multiplier) is the line x=0.5
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
    marker[f] = 0.5 - DOLFIN_EPS < f.midpoint().x() < 0.5 + DOLFIN_EPS
gamma_mesh = MeshView.create(marker, 1)

### Step 2 : Define the function spaces ###

V = FunctionSpace(mesh, "CG", 1)
LM = FunctionSpace(gamma_mesh, "DG", 0)
W = MixedFunctionSpace(V,LM)
u = Function(W.sub_space(0))
l = Function(W.sub_space(1))

v,e = TestFunctions(W)

### Step 3 : The variational formulation ###

# Homogeneous Dirichlet on left and right (x=0 and x=1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
bc = DirichletBC(V, Constant(0.0), boundary)

# Define explicitely the integration domains we will use in the bilinear and linear forms
dV = Measure("dx", domain=W.sub_space(0).mesh()) # To integrate on the 2D domain
dL = Measure("dx", domain=W.sub_space(1).mesh()) # To integrate on the 1D submesh (gamma)
# Function used in the RHS
c = 0.1
fu = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
##fu = Expression("10*exp(-50*(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)))", degree=2)
fl = Constant(c) # u = c on Gamma
# Variational formulation
a = inner(grad(u),grad(v))*dV + v*l*dL + u*e*dL 
L = fu*v*dV + fl*e*dL
F = a - L

### Step 4 : Solve the problem ###
sol = Function(W)
solve(F == 0, sol, bc)
2

I determined what was going on. I had redefined the solution as sol = Function(W), but earlier I had defined

u = Function(W.sub_space(0))
l = Function(W.sub_space(1))

and used u and l in the formulation of F.

1 Like