Hi all,

I’m trying to solve the following coupled problem

Reading the documentation and some older questions I come up with the following code. However, when I try to run it I get the error Unable to successfully call PETSc function ‘MatSetValuesLocal’ .

Where I am wrong?

```
import numpy as np
import dolfin as df
import matplotlib.pyplot as plt
Lx = 50
Ly = 10
Nx = 200
Ny = 10
# Create mesh and define function space
mesh = df.RectangleMesh(df.Point(0, 0), df.Point(Lx, Ly), Nx, Ny)
P1 = df.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = df.FunctionSpace(mesh, df.MixedElement([P1, P1]))
class DirichletBoundaryPhi(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], 0, df.DOLFIN_EPS) and on_boundary
class DirichletBoundaryPsi(df.SubDomain):
def inside(self, x, on_boundary):
return (
df.near(x[0], 0, df.DOLFIN_EPS)
or df.near(x[0], Lx, df.DOLFIN_EPS)
or df.near(x[1], 0, df.DOLFIN_EPS)
or df.near(x[1], Ly, df.DOLFIN_EPS)
) and on_boundary
bc = [
df.DirichletBC(W.sub(0), df.Constant(-1.0), DirichletBoundaryPhi()),
df.DirichletBC(W.sub(1), df.Constant(+0.0), DirichletBoundaryPsi())
]
(phi, psi) = df.Function(W)
(pv, sv) = df.TestFunction(W)
# Define a variable for the solution
w = df.Function(W)
rhof = df.Expression("exp(-pow(x[0] - 10, 2))", degree=2) # Fixed charge
# Assemble the differential form
F = (
df.inner(df.grad(phi), df.grad(pv)) + (psi**2 + rhof) * pv
+ df.inner(df.grad(psi), df.grad(sv)) + abs(psi)**(3/4)*psi*sv + phi * psi * sv
) * df.dx
# Solve
df.solve(F == 0, w, bc)
# Extract the solution
phi, psi = w.split()
fig, ax = plt.subplots()
im = df.plot(phi)
fig.colorbar(im)
```