PDDO-MultiPhysics Examples Phase field and microstructure

Phase-field crystal

Conserved Swift-Hohenberg dynamics, sixth order in space, written as a mixed system: a homogeneous density crystallizes into a periodic lattice.

phase_field_crystal.yaml 32,768 nodes sbdf2, 19 min

1The Problem

The phase-field crystal model describes a crystal by a density field that oscillates on the scale of the atomic spacing. A homogeneous state with small noise is unstable at these parameters and orders into a density wave with wavenumber about one. The mean density is conserved.

$$\frac{\partial\psi}{\partial t} = \nabla^2\mu, \qquad \mu = \big(r + (1 + \nabla^2)^2\big)\psi + \psi^3 = (1+r)\psi + 2\nabla^2\psi + \nabla^4\psi + \psi^3$$

2The Problem File

The sixth-order equation becomes two equations of at most fourth order, in ψ and the chemical potential μ. The SBDF2 scheme treats the linear part implicitly and the cubic term explicitly, with a linear stabilization Cs*psi added implicitly and subtracted explicitly, so the implicit matrix is constant. Jacobi-GMRES struggles with this sixth-order operator, so the solver factorizes it instead (2 factorizations in the whole run). The factorization dominates the first steps; after it each step takes about a second.

The complete problem file, examples/phase_field_crystal.yaml 
# Phase-field crystal (Swift-Hohenberg conserved dynamics, 6th order in space):
#   dt(psi) = lap(mu),   mu = (r + (1 + lap)^2) psi + psi^3
#           = (1 + r) psi + 2 lap(psi) + bih(psi) + psi^3
# written as a mixed system (highest derivative 4: PDDO order 4). A homogeneous state
# psi0 = -0.25 with noise crystallizes into a periodic density wave with wavenumber ~1
# (lattice spacing 2 pi); the mean of psi is conserved. IMEX (SBDF2): the cubic term is
# explicit with a linear stabilization Cs*psi (implicit - explicit), so the implicit matrix
# is constant.
"k">name: phase_field_crystal
"k">domain:
  "k">size: [25.1327, 25.1327, 25.1327]      # 4 wavelengths (8 pi)
  "k">n: [32, 32, 32]
  "k">periodic: true
"k">parameters:
  "k">r: -0.35
  "k">psi0: -0.25
  "k">Cs: 2.0                          # linear stabilization of the explicit cubic term
"k">fields:
  "k">psi: scalar
  "k">mu: scalar
"k">equations:
  "k">psi: "dt(psi) = imp(lap(mu))"
  "k">mu: "mu = imp((1 + r + Cs)*psi) + 2*lap(psi) + bih(psi) + expl(psi**3 - Cs*psi)"
"k">initial_conditions:
  "k">psi: "psi0 + 0.05*randn()"
"k">time:
  "k">end: 150.0
  "k">dt: 0.5
  "k">scheme: sbdf2                     # linear part implicit, constant matrix
"k">output:
  "k">every: 20
"k">monitors:
  "k">psi_mean: "mean(psi)"
  "k">psi_max: "max(psi)"
  "k">psi_min: "min(psi)"
  "k">energy: "integral(1/2*psi*((1 + r)*psi + 2*lap(psi) + bih(psi)) + psi**4/4)"
The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
Figure 1. The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.
The mixed system with the stabilized explicit cubic term.
Figure 2. The mixed system with the stabilized explicit cubic term.

3Results

Nodes
32×32×32 = 32,768
Periodic
xyz
PDDO order
4, horizon 3.015 Δx
Time scheme
sbdf2
Steps
300 × Δt = 0.5
Newton iterations
301
Factorizations
2
Run time
19 min

The mean density stays at -0.2498 (change below 1×10-10). The free energy falls from 788 to 36.3 as the density wave forms, and the density oscillates between -0.908 and 1.08. By t = 150 the box is ordered locally into rods and stripes with the expected spacing, but not into one perfect crystal: the energy is still falling slowly as the defects between the differently oriented regions anneal.

An isosurface of the density at t = 150: ordered rods and stripes with defects between them.
Figure 3. An isosurface of the density at t = 150: ordered rods and stripes with defects between them.
The density on the surface of the periodic box.
Figure 4. The density on the surface of the periodic box.
Free energy against time.
Figure 5. Free energy against time.
Smallest and largest density: the amplitude of the lattice growing.
Figure 6. Smallest and largest density: the amplitude of the lattice growing.

4Run It Yourself

The problem file ships with PDDO-MultiPhysics. Open it in the graphical interface and press Compute, or run it from the command line:

python -m pddo_mp gui examples/phase_field_crystal.yaml
python -m pddo_mp run examples/phase_field_crystal.yaml

Results go to results/phase_field_crystal/ next to the problem file: a ParaView time series (.pvd), the monitors as CSV, and the final state. Any entry can be changed on the command line, for example --set time.end=75 for a run to half the end time. The user manual describes every part of the problem file and of the interface.

Have a coupled problem of your own?

Send us the physics and the geometry — we will write the problem file with you, or build the model you need.

Talk to an engineer About PDDO-MultiPhysics