This is my modified code:
import random
from dolfin import *
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())
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
mesh = UnitCubeMesh.create(96, 96, 96, CellType.Type.tetrahedron)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1*P1)
du = TrialFunction(ME)
q, v = TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME)
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
mu_mid = (1.0-theta)*mu0 + theta*mu
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.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 50*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
And this is the error message:
Can't add expressions with different shapes.
Traceback (most recent call last):
File "/Users/lorenzomarta/Desktop/Simulazioni_demo_fenics/Cahn-Hilliard/demo_cahn-hilliard.py", line 258, in <module>
f = 100*c**2*(1-c)**2
File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/exproperators.py", line 243, in _rsub
return Sum(o, -self)
File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/algebra.py", line 53, in __new__
error("Can't add expressions with different shapes.")
File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different shapes.