Why this simple code is not working?

General
1

Hi Fenics community,
I’m new to Fenics, trying to reproduce the Cahn-Hilliard Equation, but encountered a solver exception. I’m using the docker image of docker.io/nielsbohr/fenics-notebook:latest

Here’s the snippet

import numpy
import random

# from mpi4py import MPI
from petsc4py import PETSc

from fenics import *

mesh = UnitSquareMesh(81, 81)
V = FiniteElement("Lagrange", triangle, 4)
ME = FunctionSpace(mesh, V*V)

du    = TrialFunction(ME)
q,v  = TestFunctions(ME)
# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L # problem itself
        self.a = a # jacobian of u
        # self.bc = bc
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        # Jacobian a
        assemble(self.a, tensor=A)

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())  # c
        values[1] = 0.0  # mu 
    def value_shape(self):
        return (2,)
u_init = InitialConditions() # time consuming
u = interpolate(u_init, ME)
u0 = interpolate(u_init, ME)

# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1

a = derivative(L, u, du)

problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
# solver = KrylovSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())

And I got the following error

Blockquote
RuntimeError Traceback (most recent call last)
Input In [14], in <cell line: 8>()
9 t += dt
10 u0.vector()[:] = u.vector()
—> 11 solver.solve(problem, u.vector())
12 file_c << (u.split()[0], t)
13 file_mu << (u.split()[1], t)

RuntimeError:
*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1653063982664/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: ba376b6aebd7a9bc089be46b50bdd9f5c548fb91

2

Please start by using a standard DOLFIN docker image, like: Package fenics-gmsh · GitHub as I have no clue how the nielsbohr’s docker image looks like.

Secondly, you are redefining u at

i.e. changing it to:

u.interpolate(u_init)
u0.interpolate(u_init)

resolves the issue.