-
Notifications
You must be signed in to change notification settings - Fork 2
/
cahn-hilliard.py
338 lines (306 loc) · 12.9 KB
/
cahn-hilliard.py
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
# Cahn-Hilliard equation
# ======================
#
# This demo is implemented in a single Python file,
# :download:`demo_cahn-hilliard.py`, which contains both the variational
# forms and the solver.
#
# This example demonstrates the solution of a particular nonlinear
# time-dependent fourth-order equation, known as the Cahn-Hilliard
# equation. In particular it demonstrates the use of
#
# * The built-in Newton solver
# * Advanced use of the base class ``NonlinearProblem``
# * Automatic linearisation
# * A mixed finite element method
# * The :math:`\theta`-method for time-dependent equations
# * User-defined Expressions as Python classes
# * Form compiler options
# * Interpolation of functions
#
#
# Equation and problem definition
# -------------------------------
#
# The Cahn-Hilliard equation is a parabolic equation and is typically
# used to model phase separation in binary mixtures. It involves
# first-order time derivatives, and second- and fourth-order spatial
# derivatives. The equation reads:
#
# .. math::
# \frac{\partial c}{\partial t} - \nabla \cdot M \left(\nabla\left(\frac{d f}{d c}
# - \lambda \nabla^{2}c\right)\right) &= 0 \quad {\rm in} \ \Omega, \\
# M\left(\nabla\left(\frac{d f}{d c} - \lambda \nabla^{2}c\right)\right) \cdot n &= 0 \quad {\rm on} \ \partial\Omega, \\
# M \lambda \nabla c \cdot n &= 0 \quad {\rm on} \ \partial\Omega.
#
# where :math:`c` is the unknown field, the function :math:`f` is
# usually non-convex in :math:`c` (a fourth-order polynomial is commonly
# used), :math:`n` is the outward directed boundary normal, and
# :math:`M` is a scalar parameter.
#
#
# Mixed form
# ^^^^^^^^^^
#
# The Cahn-Hilliard equation is a fourth-order equation, so casting it
# in a weak form would result in the presence of second-order spatial
# derivatives, and the problem could not be solved using a standard
# Lagrange finite element basis. A solution is to rephrase the problem
# as two coupled second-order equations:
#
# .. math::
# \frac{\partial c}{\partial t} - \nabla \cdot M \nabla\mu &= 0 \quad {\rm in} \ \Omega, \\
# \mu - \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in} \ \Omega.
#
# The unknown fields are now :math:`c` and :math:`\mu`. The weak
# (variational) form of the problem reads: find :math:`(c, \mu) \in V
# \times V` such that
#
# .. math::
# \int_{\Omega} \frac{\partial c}{\partial t} q \, {\rm d} x + \int_{\Omega} M \nabla\mu \cdot \nabla q \, {\rm d} x
# &= 0 \quad \forall \ q \in V, \\
# \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \, {\rm d} x - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
# &= 0 \quad \forall \ v \in V.
#
#
# Time discretisation
# ^^^^^^^^^^^^^^^^^^^
#
# Before being able to solve this problem, the time derivative must be
# dealt with. Apply the :math:`\theta`-method to the mixed weak form of
# the equation:
#
# .. math::
#
# \int_{\Omega} \frac{c_{n+1} - c_{n}}{dt} q \, {\rm d} x + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
# &= 0 \quad \forall \ q \in V \\
# \int_{\Omega} \mu_{n+1} v \, {\rm d} x - \int_{\Omega} \frac{d f_{n+1}}{d c} v \, {\rm d} x - \int_{\Omega} \lambda \nabla c_{n+1} \cdot \nabla v \, {\rm d} x
# &= 0 \quad \forall \ v \in V
#
# where :math:`dt = t_{n+1} - t_{n}` and :math:`\mu_{n+\theta} =
# (1-\theta) \mu_{n} + \theta \mu_{n+1}`. The task is: given
# :math:`c_{n}` and :math:`\mu_{n}`, solve the above equation to find
# :math:`c_{n+1}` and :math:`\mu_{n+1}`.
#
#
# Demo parameters
# ^^^^^^^^^^^^^^^
#
# The following domains, functions and time stepping parameters are used
# in this demo:
#
# * :math:`\Omega = (0, 1) \times (0, 1)` (unit square)
# * :math:`f = 100 c^{2} (1-c)^{2}`
# * :math:`\lambda = 1 \times 10^{-2}`
# * :math:`M = 1`
# * :math:`dt = 5 \times 10^{-6}`
# * :math:`\theta = 0.5`
#
#
# Implementation
# --------------
#
# This demo is implemented in the :download:`demo_cahn-hilliard.py`
# file.
#
# First, the modules :py:mod:`random` :py:mod:`matplotlib`
# :py:mod:`dolfin` module are imported::
import random
from dolfin import *
# .. index:: Expression
#
# A class which will be used to represent the initial conditions is then
# created::
# 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):
values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
# It is a subclass of :py:class:`Expression
# <dolfin.functions.expression.Expression>`. In the constructor
# (``__init__``), the random number generator is seeded. If the program
# is run in parallel, the random number generator is seeded using the
# rank (process number) to ensure a different sequence of numbers on
# each process. The function ``eval`` returns values for a function of
# dimension two. For the first component of the function, a randomized
# value is returned. The method ``value_shape`` declares that the
# :py:class:`Expression <dolfin.functions.expression.Expression>` is
# vector valued with dimension two.
#
# .. index::
# single: NonlinearProblem; (in Cahn-Hilliard demo)
#
# A class which will represent the Cahn-Hilliard in an abstract from for
# use in the Newton solver is now defined. It is a subclass of
# :py:class:`NonlinearProblem <dolfin.cpp.NonlinearProblem>`. ::
# 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)
# The constructor (``__init__``) stores references to the bilinear
# (``a``) and linear (``L``) forms. These will used to compute the
# Jacobian matrix and the residual vector, respectively, for use in a
# Newton solver. The function ``F`` and ``J`` are virtual member
# functions of :py:class:`NonlinearProblem
# <dolfin.cpp.NonlinearProblem>`. The function ``F`` computes the
# residual vector ``b``, and the function ``J`` computes the Jacobian
# matrix ``A``.
#
# Next, various model parameters are defined::
# 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
# .. index::
# singe: form compiler options; (in Cahn-Hilliard demo)
#
# It is possible to pass arguments that control aspects of the generated
# code to the form compiler. The lines ::
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# tell the form to apply optimization strategies in the code generation
# phase and the use compiler optimization flags when compiling the
# generated C++ code. Using the option ``["optimize"] = True`` will
# generally result in faster code (sometimes orders of magnitude faster
# for certain operations, depending on the equation), but it may take
# considerably longer to generate the code and the generation phase may
# use considerably more memory).
#
# A unit square mesh with 97 (= 96 + 1) vertices in each direction is
# created, and on this mesh a :py:class:`FunctionSpace
# <dolfin.functions.functionspace.FunctionSpace>` ``ME`` is built using
# a pair of linear Lagrangian elements. ::
# Create mesh and build function space
mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Trial and test functions of the space ``ME`` are now defined::
# Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
# .. index:: split functions
#
# For the test functions, :py:func:`TestFunctions
# <dolfin.functions.function.TestFunctions>` (note the 's' at the end)
# is used to define the scalar test functions ``q`` and ``v``. The
# :py:class:`TrialFunction <dolfin.functions.function.TrialFunction>`
# ``du`` has dimension two. Some mixed objects of the
# :py:class:`Function <dolfin.functions.function.Function>` class on
# ``ME`` are defined to represent :math:`u = (c_{n+1}, \mu_{n+1})` and
# :math:`u0 = (c_{n}, \mu_{n})`, and these are then split into
# sub-functions::
# 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)
# The line ``c, mu = split(u)`` permits direct access to the components
# of a mixed function. Note that ``c`` and ``mu`` are references for
# components of ``u``, and not copies.
#
# .. index::
# single: interpolating functions; (in Cahn-Hilliard demo)
#
# Initial conditions are created by using the class defined at the
# beginning of the demo and then interpolating the initial conditions
# into a finite element space::
# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
# The first line creates an object of type ``InitialConditions``. The
# following two lines make ``u`` and ``u0`` interpolants of ``u_init``
# (since ``u`` and ``u0`` are finite element functions, they may not be
# able to represent a given function exactly, but the function can be
# approximated by interpolating it in a finite element space).
#
# .. index:: automatic differentiation
#
# The chemical potential :math:`df/dc` is computed using automated
# differentiation::
# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
# The first line declares that ``c`` is a variable that some function
# can be differentiated with respect to. The next line is the function
# :math:`f` defined in the problem statement, and the third line
# performs the differentiation of ``f`` with respect to the variable
# ``c``.
#
# It is convenient to introduce an expression for :math:`\mu_{n+\theta}`::
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# which is then used in the definition of the variational forms::
# 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
# This is a statement of the time-discrete equations presented as part
# of the problem statement, using UFL syntax. The linear forms for the
# two equations can be summed into one form ``L``, and then the
# directional derivative of ``L`` can be computed to form the bilinear
# form which represents the Jacobian matrix::
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# .. index::
# single: Newton solver; (in Cahn-Hilliard demo)
#
# The DOLFIN Newton solver requires a :py:class:`NonlinearProblem
# <dolfin.cpp.NonlinearProblem>` object to solve a system of nonlinear
# equations. Here, we are using the class ``CahnHilliardEquation``,
# which was declared at the beginning of the file, and which is a
# sub-class of :py:class:`NonlinearProblem
# <dolfin.cpp.NonlinearProblem>`. We need to instantiate objects of both
# ``CahnHilliardEquation`` and :py:class:`NewtonSolver
# <dolfin.cpp.NewtonSolver>`::
# 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
# The string ``"lu"`` passed to the Newton solver indicated that an LU
# solver should be used. The setting of
# ``parameters["convergence_criterion"] = "incremental"`` specifies that
# the Newton solver should compute a norm of the solution increment to
# check for convergence (the other possibility is to use ``"residual"``,
# or to provide a user-defined check). The tolerance for convergence is
# specified by ``parameters["relative_tolerance"] = 1e-6``.
#
# To run the solver and save the output to a VTK file for later visualization,
# the solver is advanced in time from :math:`t_{n}` to :math:`t_{n+1}` until
# a terminal time :math:`T` is reached::
# Output file
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)
# The string ``"compressed"`` indicates that the output data should be
# compressed to reduce the file size. Within the time stepping loop, the
# solution vector associated with ``u`` is copied to ``u0`` at the
# beginning of each time step, and the nonlinear problem is solved by
# calling :py:func:`solver.solve(problem, u.vector())
# <dolfin.cpp.NewtonSolver.solve>`, with the new solution vector
# returned in :py:func:`u.vector() <dolfin.cpp.Function.vector>`. The
# ``c`` component of the solution (the first component of ``u``) is then
# written to file at every time step.