Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins
# as discussed on 28/11/2019
# Example can be found at
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/
# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1
# Runs with fenics 2019.01
# The resulting .pvd file can be opened using default settings in ParaView
import matplotlib.pyplot as plt
import mshr as ms
import numpy as np
from scipy.optimize import curve_fit
import random
from dolfin import *
# Class representing the intial conditions
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
random.seed(2 + MPI.rank(MPI.comm_world))
super().__init__(**kwargs)
def eval(self, values, x):
if (x[0] > -0.1) & (x[0] < 0.1) & (x[1] > -0.1) & (x[1] < 0.1) & (x[2] > -0.1) & (x[2] < 0.1):
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[0] = 0.82
else:
values[0] = 0.3
values[1] = 0.0
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
# Class for interfacing with the Newton solver
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)
# Model parameters
kappa = 0.25*10e-04 # surface parameter
dt = 10.0e-02 # time step
X = 2.2 # Chi Flory Huggins parameter
# time stepping family, e.g.:
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta = 0.5
# Form compiler optionsmai
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and build function space
x_mesh = 10*96
y_mesh = 1
# mesh = UnitSquareMesh.create(x_mesh, y_mesh, CellType.Type.quadrilateral)
# mesh = RectangleMesh(Point(0, 0), Point(1, 0.02), x_mesh, y_mesh)
domain = ms.Sphere(Point(0, 0, 0), 1.0)
mesh = ms.generate_mesh(domain, 20)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Define trial and test functions
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)
# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
# Compute the chemical potential df/dc
c = variable(c)
dfdc = ln(c/(1-c))+(1-2*X*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*c*(1-c)*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(grad(c), grad(v))*dx
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 10*dt
while (t < T):
file << (u.split()[0], t)
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())