# cahnhilliard.py¶

In this script we solve the Cahn-Hilliard equation, which models the unmixing of two phases (φ=+1 and φ=-1) under the influence of surface tension. It is a mixed equation of two scalar equations for phase φ and chemical potential η:

```dφ/dt = -div(J(η))

η = ψ'(φ) σ / ε - σ ε Δ(φ)
```

along with constitutive relations for the flux vector J = -M ∇η and the double well potential ψ = .25 (φ² - 1)², and subject to boundary conditions ∂ₙφ = -σd / σ ε and ∂ₙη = 0. Parameters are the interface thickness ε, fluid surface tension σ, differential wall surface tension σd, and mobility M.

Cahn-Hilliard is a diffuse interface model, which means that phases do not separate sharply, but instead form a transition zone between the phases. The transition zone has a thickness proportional to ε, as is readily confirmed in one dimension, where a steady state solution on an infinite domain is formed by η(x) = 0, φ(x) = tanh(x / √2 ε).

The main driver of the unmixing process is the double well potential ψ that is proportional to the mixing energy, taking its minima at the pure phases φ=+1 and φ=-1. The interface between the phases adds an energy contribution proportional to its length. At the wall we have a phase-dependent fluid-solid energy. Over time, the system minimizes the total energy:

```E(φ) = ∫_Ω ψ(φ) σ / ε + ∫_Ω .5 σ ε ‖∇φ‖² + ∫_Γ (σm + φ σd)

\                \                  \

mixing energy    interface energy   wall energy
```

Proof: the time derivative of E followed by substitution of the strong form and boundary conditions yields dE/dt = ∫_Ω η dφ/dt = -∫_Ω M ‖∇η‖² ≤ 0. □

Switching to discrete time we set dφ/dt = (φ - φ0) / dt and add a stabilizing perturbation term δψ(φ, φ0) to the doube well potential for reasons outlined below. This yields the following time discrete system:

```φ = φ0 - dt div(J(η))

η = (ψ'(φ) + δψ'(φ, ψ0)) σ / ε - σ ε Δ(φ)
```

with the equivalent weak formulation:

```∂/∂η ∫_Ω [ .5 dt J(η)·∇η - η (φ - φ0) ] = 0

∂/∂φ ∫_Ω [ E(φ) + δψ(φ, φ0) σ / ε - η φ ] = 0
```

For stability we wish for the perturbation δψ to be such that the time discrete system preserves the energy dissipation property E(φ) ≤ E(φ0) for any timestep dt. To derive suitable perturbation terms to this effect, we define without loss of generality δψ’(φ, φ0) = .5 (φ - φ0) f(φ, φ0) and derive the following condition for unconditional stability:

```E(φ) - E(φ0) = ∫_Ω .5 (1 - φ² - .5 (φ + φ0)² - f(φ, φ0)) (φ - φ0)² σ / ε

- ∫_Ω (.5 σ ε ‖∇φ - ∇φ0‖² + dt M ‖∇η‖²) ≤ 0
```

The inequality holds true if the perturbation f is bounded from below such that f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)². To keep the energy minima at the pure phases we additionally impose that f(±1, φ0) = 0, and select 1 - φ² as a suitable upper bound which we will call “nonlinear”.

We next observe that η is linear in φ if f(φ, φ0) = g(φ0) - φ² - (φ + φ0)² for any function g, which dominates if g(φ0) ≥ 1 + .5 (φ + φ0)². While this cannot be made to hold true for all φ, we make it hold for -√2 ≤ φ, φ0 ≤ √2 by defining g(φ0) = 2 + 2 ‖φ0‖ + φ0², which we will call “linear”. This scheme further satisfies a weak minima preservation f(±1, ±‖φ0‖) = 0.

We have thus arrived at the three stabilization schemes implemented here:

• nonlinear: f(φ, φ0) = 1 - φ²

• linear: f(φ, φ0) = 2 + 2 ‖φ0‖ - 2 φ (φ + φ0)

• none: f(φ, φ0) = 0 (not unconditionally stable)

The stab enum in this script defines the schemes in terms of δψ to allow Nutils to construct the residuals through automatic differentiation.

NOTE: This script uses dimensional quantities and requires the nutils.SI module, which is installed separate from the the core nutils.

``` 78from nutils import mesh, function, solver, numeric, export, cli, testing
79from nutils.expression_v2 import Namespace
80from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse
81import numpy
82import itertools
83import enum
84import treelog as log
85
86
87class stab(enum.Enum):
88    nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)'
89    linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)'
90    none = '0'  # not unconditionally stable
91
92
93def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
94         wtensn: Tension, wtensp: Tension, nelems: int, etype: str, degree: int,
95         timestep: Time, tol: Energy/Length, endtime: Time, seed: int, circle: bool,
96         stab: stab, showflux: bool):
97    '''
98    Cahn-Hilliard equation on a unit square/circle.
99
100    .. arguments::
101
102       size [10cm]
103         Domain size.
104       epsilon [2mm]
105         Interface thickness; defaults to an automatic value based on the
106         configured mesh density if left unspecified.
107       mobility [1mL*s/kg]
108         Mobility.
109       stens [50mN/m]
110         Surface tension.
111       wtensn [30mN/m]
112         Wall surface tension for phase -1.
113       wtensp [20mN/m]
114         Wall surface tension for phase +1.
115       nelems 
116         Number of elements along domain edge. When set to zero a value is set
117         automatically based on the configured domain size and epsilon.
118       etype [square]
119         Type of elements (square/triangle/mixed).
120       degree 
121         Polynomial degree.
122       timestep [.5s]
123         Time step.
124       tol [1nJ/m]
125         Newton tolerance.
126       endtime [1min]
127         End of the simulation.
128       seed 
129         Random seed for the initial condition.
130       circle [no]
131         Select circular domain as opposed to a unit square.
132       stab [linear]
133         Stabilization method (linear/nonlinear/none).
134       showflux [no]
135         Overlay flux vectors on phase plot
136    '''
137
138    nmin = round(size / epsilon)
139    if nelems <= 0:
140        nelems = nmin
141        log.info('setting nelems to {}'.format(nelems))
142    elif nelems < nmin:
143        log.warning('mesh is too coarse, consider increasing nelems to {:.0f}'.format(nmin))
144
145    log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi))
146
147    domain, geom = mesh.unitsquare(nelems, etype)
148    if circle:
149        angle = (geom-.5) * (numpy.pi/2)
150        geom = .5 + function.sin(angle) * function.cos(angle)[[1, 0]] / numpy.sqrt(2)
151
152    bezier = domain.sample('bezier', 5)  # sample for surface plots
153    grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5)  # sample for quivers
154
155    φbasis = ηbasis = domain.basis('std', degree=degree)
156    ηbasis *= stens / epsilon  # basis scaling to give η the required unit
157
158    ns = Namespace()
159    ns.x = size * geom
160    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
161    ns.ε = epsilon
162    ns.σ = stens
163    ns.φ = function.dotarg('φ', φbasis)
164    ns.σmean = (wtensp + wtensn) / 2
165    ns.σdiff = (wtensp - wtensn) / 2
166    ns.σwall = 'σmean + φ σdiff'
167    ns.φ0 = function.dotarg('φ0', φbasis)
168    ns.dφ = 'φ - φ0'
169    ns.η = function.dotarg('η', ηbasis)
170    ns.ψ = '.25 (φ^2 - 1)^2'
171    ns.δψ = stab.value
172    ns.M = mobility
173    ns.J_i = '-M ∇_i(η)'
174    ns.dt = timestep
175
176    nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=7)
177    nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=7)
178    nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=7)
179    nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7)
180
181    numpy.random.seed(seed)
182    state = dict(φ=numpy.random.normal(0, .5, φbasis.shape))  # initial condition
183
184    with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
185        for istep in steps:
186
187            E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state))
188            log.user('energy: {0:,.0μJ/m} ({1:.0f}% mixture, {1:.0f}% interface, {1:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E)))
189
190            state['φ0'] = state['φ']
191            state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1)
192
193            with export.mplfigure('phase.png') as fig:
194                ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')
195                x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state)
196                im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr')
197                im.set_clim(-1, 1)
198                fig.colorbar(im)
199                if showflux:
200                    x, J = grid.eval(['x_i', 'J_i'] @ ns, **state)
201                    log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:, 0], J[:, 1]))))
202                    ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r')
203                    ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b')
204                ax.autoscale(enable=True, axis='both', tight=True)
205
206    return state
```

If the script is executed (as opposed to imported), `nutils.cli.run()` calls the main function with arguments provided from the command line.

```212if __name__ == '__main__':
213    cli.run(main)
```

Once a simulation is developed and tested, it is good practice to save a few strategic return values for regression testing. The `nutils.testing` module, which builds on the standard `unittest` framework, facilitates this by providing `nutils.testing.TestCase.assertAlmostEqual64()` for the embedding of desired results as compressed base64 data.

```222class test(testing.TestCase):
223
224    def test_initial(self):
225        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
226                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
227                     etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
228                     endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True)
229        with self.subTest('concentration'):
230            self.assertAlmostEqual64(state['φ0'], '''
231                eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
232                L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''')
233
234    def test_square(self):
235        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
236                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
237                     etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
238                     endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True)
239        with self.subTest('concentration'):
240            self.assertAlmostEqual64(state['φ'], '''
241                eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk
242                ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''')
243        with self.subTest('chemical-potential'):
244            self.assertAlmostEqual64(state['η'], '''
245                eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
246                xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
247
248    def test_mixedcircle(self):
249        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
250                     stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
251                     etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
252                     endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
253        with self.subTest('concentration'):
254            self.assertAlmostEqual64(state['φ'], '''
255                eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK