Error while trying to solve coupled PDEs

variational formulation
1

Hi all,
I’m trying to solve the following coupled problem

- \nabla \cdot \nabla \phi = -\psi^2 + \rho_{fx} \\ - \nabla \cdot \nabla \psi +\psi^{4/3} \psi + \phi \psi = 0

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)
2

This was a pain to debug. The fix is the following:

w = df.Function(W)
pv, sv = df.TestFunctions(W)
phi, psi = df.split(w)

In debugging I noticed that the matrix you’re assembling is empty. This implies that the derivative isn’t correctly formed. And indeed this lead to finding that w was defined as a Function and not actually used in your nonlinear formulation. Therefore the Gateaux derviative is zero, and no matrix is assembled.

The fix will allow you to assemble the problem; however, with regards to getting the solver to work I refer you to some tips and tricks. The nonlinear problem you’ve written is not trivial given the fractional powers.

3

Thank you so much! So fenics keep track of all the functions I declared. I thought that “Function” modeled a generic field defined on the domain.

For what concerns convergence, is it really connected to the nonlinear term? Because I tried to remove the nonlinear term or even to comment all the second line of F (so that it should simplify to a Poisson equation) and still get the same error.

4

Keep in mind that the python layer of DOLFIN works with UFL so that many objects have an algebraic representation by inheriting the functionality of UFL. So when you solve the nonlinear problem, the Gateaux derivative is computed using UFL with the UFL representation of the DOLFIN objects (i.e. your Function).

Regarding convergence: The original error you were getting has nothing to do with convergence. I just noticed that running the code with fix above and the residual and Jacobian as defined for your formulation will require more care/work, since you have fractional powers.

5

I see, I’m just trying to run the same example with your correction and now I get

*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached

I thought that it was connected to the nonlinearity you mentioned, but even if I remove it I get the same error.