make cahnhilliard example dimensional

This patch changes the cahnhilliard example to use dimensional values.

.github/workflows/test.yaml

@@ -134,7 +134,7 @@ jobs:
         id: install
         run: |
           python -um pip install --upgrade wheel
-          python -um pip install --upgrade treelog stringly matplotlib scipy pillow numpy
+          python -um pip install --upgrade treelog stringly matplotlib scipy pillow numpy git+https://github.com/evalf/nutils-SI.git
           # Install Nutils from `dist` dir created in job `build-python-package`.
           python -um pip install --no-index --find-links ./dist nutils
       - name: Test

examples/cahnhilliard.py

@@ -1,101 +1,203 @@
 #! /usr/bin/env python3
 #
-# In this script we solve the Cahn-Hiilliard equation, which models the
-# unmixing of two phases under the effect of surface tension.
+# 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::
+#
+#     ∂/∂η ∫_Ω [ η (φ - φ0) + .5 dt J(η)·∇η ] = 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.
 
-from nutils import mesh, function, solver, export, cli, testing
+from nutils import mesh, function, solver, numeric, export, cli, testing
 from nutils.expression_v2 import Namespace
-import numpy, treelog, itertools, enum, typing
+from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse
+import numpy, itertools, enum
+import treelog as log
 
 class stab(enum.Enum):
-  none = '0' # for educational purposes only
-  linear = '.5 dc^2 (6 - 6 c^2 + 8 c dc - 3 dc^2) / epsilon^2'
-  optimal = '.5 dc^2 (1 - dc^2 / 12) / epsilon^2'
-
-def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[float],
-         contactangle:float, timestep:float, mtol:float, seed:int, circle:bool, stab:stab):
+  nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)'
+  linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)'
+  none = '0' # not unconditionally stable
+
+def main(size:Length, epsilon:Length, mobility:Time/Density, stens:Tension,
+         wtensn:Tension, wtensp:Tension, nelems:int, etype:str, degree:int,
+         timestep:Time, tol:Energy/Length, endtime:Time, seed:int, circle:bool,
+         stab:stab, showflux:bool):
   '''
   Cahn-Hilliard equation on a unit square/circle.
 
   .. arguments::
 
-     nelems [20]
-       Number of elements along domain edge.
+     size [10cm]
+       Domain size.
+     epsilon [2mm]
+       Interface thickness; defaults to an automatic value based on the
+       configured mesh density if left unspecified.
+     mobility [1mL*s/kg]
+       Mobility.
+     stens [50mN/m]
+       Surface tension.
+     wtensn [30mN/m]
+       Wall surface tension for phase -1.
+     wtensp [20mN/m]
+       Wall surface tension for phase +1.
+     nelems [0]
+       Number of elements along domain edge. When set to zero a value is set
+       automatically based on the configured domain size and epsilon.
      etype [square]
        Type of elements (square/triangle/mixed).
-     btype [std]
-       Type of basis function (std/spline), with availability depending on the
-       configured element type.
      degree [2]
        Polynomial degree.
-     epsilon []
-       Interface thickness; defaults to an automatic value based on the
-       configured mesh density if left unspecified.
-     contactangle [90]
-       Wall contact angle in degrees.
-     timestep [.01]
+     timestep [.5s]
        Time step.
-     mtol [.01]
-       Threshold value for chemical potential peak to peak difference, used as
-       a stop criterion.
+     tol [1nJ/m]
+       Newton tolerance.
+     endtime [1min]
+       End of the simulation.
      seed [0]
        Random seed for the initial condition.
      circle [no]
        Select circular domain as opposed to a unit square.
      stab [linear]
-       Stabilization method (linear/optimal/none).
+       Stabilization method (linear/nonlinear/none).
+     showflux [no]
+       Overlay flux vectors on phase plot
   '''
 
-  mineps = 1./nelems
-  if epsilon is None:
-    treelog.info('setting epsilon={}'.format(mineps))
-    epsilon = mineps
-  elif epsilon < mineps:
-    treelog.warning('epsilon under crititical threshold: {} < {}'.format(epsilon, mineps))
+  nmin = round(size / epsilon)
+  if nelems <= 0:
+    nelems = nmin
+    log.info('setting nelems to {}'.format(nelems))
+  elif nelems < nmin:
+    log.warning('mesh is too coarse, consider increasing nelems to {:.0f}'.format(nmin))
+
+  log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi))
 
   domain, geom = mesh.unitsquare(nelems, etype)
-  bezier = domain.sample('bezier', 5) # sample for plotting
+  if circle:
+    angle = (geom-.5) * (numpy.pi/2)
+    geom = .5 + function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2)
+
+  bezier = domain.sample('bezier', 5) # sample for surface plots
+  grid = domain.locate(geom, numeric.simplex_grid([1,1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5) # sample for quivers
+
+  φbasis = ηbasis = domain.basis('std', degree=degree)
+  ηbasis *= stens / epsilon # basis scaling to give η the required unit
 
   ns = Namespace()
-  if not circle:
-    ns.x = geom
-  else:
-    angle = (geom-.5) * (numpy.pi/2)
-    ns.x = function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2)
+  ns.x = size * geom
   ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-  ns.epsilon = epsilon
-  ns.ewall = .5 * numpy.cos(contactangle * numpy.pi / 180)
-  ns.cbasis = ns.mbasis = domain.basis('std', degree=degree)
-  ns.c = function.dotarg('c', ns.cbasis)
-  ns.dc = ns.c - function.dotarg('c0', ns.cbasis)
-  ns.m = function.dotarg('m', ns.mbasis)
-  ns.F = '.5 (c^2 - 1)^2 / epsilon^2'
-  ns.dF = stab.value
+  ns.ε = epsilon
+  ns.σ = stens
+  ns.φ = function.dotarg('φ', φbasis)
+  ns.σmean = (wtensp + wtensn) / 2
+  ns.σdiff = (wtensp - wtensn) / 2
+  ns.σwall = 'σmean + φ σdiff'
+  ns.φ0 = function.dotarg('φ0', φbasis)
+  ns.dφ = 'φ - φ0'
+  ns.η = function.dotarg('η', ηbasis)
+  ns.ψ = '.25 (φ^2 - 1)^2'
+  ns.δψ = stab.value
+  ns.M = mobility
+  ns.J_i = '-M ∇_i(η)'
   ns.dt = timestep
 
-  nrg_mix = domain.integral('F dV' @ ns, degree=7)
-  nrg_iface = domain.integral('.5 ∇_k(c) ∇_k(c) dV' @ ns, degree=7)
-  nrg_wall = domain.boundary.integral('(abs(ewall) + c ewall) dS' @ ns, degree=7)
-  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 ∇_k(m) ∇_k(m)) dV' @ ns, degree=7)
+  nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=7)
+  nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=7)
+  nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=7)
+  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7)
 
   numpy.random.seed(seed)
-  state = dict(c=numpy.random.normal(0,.5,ns.cbasis.shape), m=numpy.random.normal(0,.5,ns.mbasis.shape)) # initial condition
+  state = dict(φ=numpy.random.normal(0, .5, φbasis.shape)) # initial condition
 
-  with treelog.iter.plain('timestep', itertools.count()) as steps:
+  with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
    for istep in steps:
 
-    E = function.eval([nrg_mix, nrg_iface, nrg_wall], **state)
-    treelog.user('energy: {0:.3f} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(sum(E), 100*numpy.array(E)/sum(E)))
-
-    x, c, m = bezier.eval(['x_i', 'c', 'm'] @ ns, **state)
-    export.triplot('phase.png', x, c, tri=bezier.tri, clim=(-1,1))
-    export.triplot('chempot.png', x, m, tri=bezier.tri)
-
-    if numpy.ptp(m) < mtol:
-      break
-
-    state['c0'] = state['c']
-    state = solver.optimize(['c', 'm'], nrg, arguments=state, tol=1e-10)
+    E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state))
+    log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E)))
+
+    state['φ0'] = state['φ']
+    state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1)
+
+    with export.mplfigure('phase.png') as fig:
+      ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')
+      x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state)
+      im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr')
+      im.set_clim(-1, 1)
+      fig.colorbar(im)
+      if showflux:
+        x, J = grid.eval(['x_i', 'J_i'] @ ns, **state)
+        log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:,0], J[:,1]))))
+        ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r')
+        ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b')
+      ax.autoscale(enable=True, axis='both', tight=True)
 
   return state
 
@@ -114,37 +216,34 @@ def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[f
 class test(testing.TestCase):
 
   def test_initial(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=float('inf'), seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ0'], '''
       eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
       L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoBYgCd/w7NQModNFjLtckB0VA0rDCiM+zKBMzPygjMcMr4yJcy0MsUyXY0OcoxMFo19zE5Np/JMDTG
-      yk7KGstQzFgwvMnDNXk0Msm+Njc3SjZizebJEjbPy1w25zLsNfQzSjURLRk3Qt4uBQ==''')
 
   def test_square(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/5o2nDZANsQ1dTQO1QXQmTaaNjw2vDVVNOrPqM5mNmY2xDWwNPDRNMsZyys2KjYdNQoyh8v9
-      yffJ1zXTNSA0Ks2JyorJicluNWQ1IDJXy+/JNck3yWQ1WjX0MVHL78k3yTnJbYgt7Q==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNpVyU0KgCAABeGrK9SiQCiQCsrwB8GjzLF8uJNvNYzBsuP5yGKmWlhxXKMSmxzchPGceF4iVU55+Ck0
-      masD28JDNQ==''')
-
-  def test_contactangle(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=45, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/zs2bjYbNqs19zNqzRXMazaZNkw25zVxNLXOncxDNnU28zUqNeIwTsvKyhc2TDaANdozSMwp
-      yuXJojXiNYo05c/Oyp3Jb8n7NE81iDK2ywfKO8kayXI03jQjLyrLzskZyfrIkb8vTg==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoNzCEOggAYgNFLvbPguILeQRpWgyQPoAYCuBGhWHVzGshOLegMSvffS1/5EqmF2sUnTKJylbO3l6mZ
-      pb2TZ5jLrDWObmGlUDroDWFjq43H3dcvaqdz9TCGP1tYLVU=''')
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], '''
+      eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk
+      ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''')
+    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], '''
+      eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
+      xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
 
   def test_mixedcircle(self):
-    state = main(nelems=3, etype='mixed', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=True, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/y415zQQNUw1dzVTNB01fTS4My81dDQpNWU03DJENDo0rjEVMCfMj8vjysUzozMmz3E0VCyg
-      1I00MDMWM3bQU8t4ym7Kd8uKzoTPgc7PzmfOYs3ZzJvM6jIozunL7sqezgXMYsUuJg==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoljEkKgEAQxP7/B/EkiCCO27gOuJxy8FEGpSEUoaoTkYWTgyAjHRuTZpGzHGgpKchMMz2JRnN/l6j1
-      uctge2W08W8fdlPF5ccXGqBI7Q==''')
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], '''
+      eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK
+      z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''')
+    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], '''
+      eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv
+      x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''')