Hello,
I am running into a Runtime error when attempting to solve a phase-separation problem using the Cahn-Hilliard formalism. The approach is based on the demo: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html.
When running the code below, I obtain the following error:
RuntimeErrorTraceback (most recent call last)
<ipython-input-2-5315e8a12b8d> in <module>
164
165 u_old.vector()[:] = u_new.vector()
--> 166 solve(Res == 0, u_new)
167
168 # # Plot initial condition
/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
218 # tolerance)
219 elif isinstance(args[0], ufl.classes.Equation):
--> 220 _solve_varproblem(*args, **kwargs)
221
222 # Default case, just call the wrapped C++ solve function
/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
264 solver = NonlinearVariationalSolver(problem)
265 solver.parameters.update(solver_parameters)
--> 266 solver.solve()
267
268
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 solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
The code producing this error is:
# In order to use solve this problem we need to import all necessary modules.
from __future__ import print_function
import numpy as np
from fenics import *
from mshr import *
import ufl as ufl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# First let's illustrate the free energy function chosen.
a = 4.
b = 1.
cst = 1/3.
# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT
# on one of the two corners (0, 1) and (1, 0)
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], 0) and near(x[1], 1)) or
(near(x[0], 1) and near(x[1], 0)))) and on_boundary)
def map(self, x, y):
if near(x[0], 1) and near(x[1], 1):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
elif near(x[0], 1):
y[0] = x[0] - 1.
y[1] = x[1]
else: # near(x[1], 1)
y[0] = x[0]
y[1] = x[1] - 1.
# Create mesh and define function space with periodic boundary conditions
# mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
N = 64
mesh = UnitSquareMesh.create(N, N, CellType.Type.triangle)
dim = mesh.topology().dim() # dimensionality of the problem
print("dim = ", dim)
plot(mesh, linewidth=0.2)
plt.show()
degreeElements = 1
V = FiniteElement('Lagrange', mesh.ufl_cell(), degreeElements)
VFS = FunctionSpace(mesh, MixedElement([V,V,V,V]),constrained_domain=PeriodicBoundary())
# Define functions
u_new = Function(VFS) # solution for the next step
u_old = Function(VFS) # solution from previous step
u_new.rename("fields","")
# Split mixed functions
phi1_new, phi2_new, mu1_new, mu2_new = split(u_new)
phi1_old, phi2_old, mu1_old, mu2_old = split(u_old)
# Define test functions
tf = TestFunction(VFS)
q1, q2, v1, v2 = split(tf)
# Class representing the intial conditions
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
np.random.seed(123)
super().__init__(**kwargs)
def eval(self, values, x):
values[0] = cst*np.random.rand(1) # concentration
values[1] = cst - values[0]
values[2:] = 0.0 # chemical potential
def value_shape(self):
return (4,)
# Create intial conditions and interpolate
u_init = InitialConditions(degree=degreeElements)
u_new.interpolate(u_init)
u_old.interpolate(u_init)
phi2_new = cst - phi1_new
phi2_old = cst - phi1_old
# define energy
phi1_new = variable(phi1_new)
phi2_new = variable(phi2_new)
f = - 0.5*a*(phi1_new*phi1_new + phi2_new*phi2_new + (1 - phi1_new - phi2_new)*(1 - phi1_new - phi2_new)) \
+ 0.25*b*(phi1_new*phi1_new + phi2_new*phi2_new + (1 - phi1_new - phi2_new)*(1 - phi1_new - phi2_new))**2
dfdphi1 = diff(f, phi1_new)
dfdphi2 = diff(f, phi2_new)
# parameters
lmbda = 1.0e-04 # surface parameter
dt = 1.0e-6 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# mu_(n+theta)
mu1_mid = (1.0-theta)*mu1_old + theta*mu1_new
mu2_mid = (1.0-theta)*mu2_old + theta*mu2_new
# Weak statement of the equations
Res_0 = (phi1_new - phi1_old)*q1/dt*dx \
+ (phi2_new - phi2_old)*q2/dt*dx \
+ dot(grad(mu1_mid), grad(q1))*dx \
+ dot(grad(mu2_mid), grad(q2))*dx
Res_1 = (mu1_new - dfdphi1)*v1*dx \
+ (mu2_new - dfdphi2)*v2*dx \
- 0.5*lmbda*dot(grad(phi2_new + 2.*phi1_new), grad(v1))*dx \
- 0.5*lmbda*dot(grad(phi1_new + 2.*phi2_new), grad(v2))*dx
Res = Res_0 + Res_1
u_old.vector()[:] = u_new.vector()
solve(Res == 0, u_new)
When reducing the problem to 1 variable, the code runs without error. I can provide this code if needed.
I am running this code in a docker container based on the image quay.io/fenicsproject/stable:current (ID f99661bca197).
I would be grateful for any help, and I’d be happy to edit this post if it doesn’t conform with the posting standards.
Thanks,