diff --git a/content/cahn-hilliard/cahn_hilliard.py b/content/cahn-hilliard/cahn_hilliard.py
new file mode 100644
index 0000000..5dcced0
--- /dev/null
+++ b/content/cahn-hilliard/cahn_hilliard.py
@@ -0,0 +1,191 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import sys
+import os
+import math
+import random
+
+random.seed(42) # for reproducibility
+
+sys.path.append("..")
+
+import modsimtools as util
+import ug4py.pyugcore as ug4
+import ug4py.pyconvectiondiffusion as cd
+import ug4py.pylimex as limex
+
+
+# Grid and domain
+# gridName = "grids/laplace_sample_grid_2d_l50.ugx"
+gridName = "grids/fish_grid_with_eye.ugx"
+numRefs = 3
+mandatorySubsets = ["Inner", "Boundary"]
+
+# Create domain (modsimtools handles loading and refinement)
+dom = util.CreateDomain(gridName, numRefs, mandatorySubsets)
+
+
+# Approximation space: two unknowns `mu` and `c` (piecewise linear)
+approxSpace = ug4.ApproximationSpace2d(dom)
+approxSpace.add_fct("mu", "Lagrange", 1)
+approxSpace.add_fct("c", "Lagrange", 1)
+approxSpace.init_levels()
+approxSpace.init_top_surface()
+approxSpace.print_statistic()
+
+
+# Parameters
+class CahnHilliardModel:
+ M = 10000.0
+ gamma = 2.0 # ~ layer thickness
+ cinf = 1.0
+ f0 = lambda c: 0.25*(c**2 - CahnHilliardModel.cinf**2)**2
+ f1 = lambda c: c**3 - CahnHilliardModel.cinf * c
+ f2 = lambda c: 3*c**2 - CahnHilliardModel.cinf
+ ka = 0.03
+ ki = 0.1 # \approx 3*ka, to have a non-trivial steady state at c = 0.25 (since ka*(1-c) - ki*c = 0 => c = ka/(ka+ki) = 0.01/0.04 = 0.25)
+
+
+
+
+# Create element-wise discretizations (FE)
+elemDisc = {}
+elemDisc["mu"] = cd.ConvectionDiffusionFE2d("mu", "Inner")
+elemDisc["c"] = cd.ConvectionDiffusionFE2d("c", "Inner")
+
+
+
+
+# Eqn 1: \partial_t c - M \triangle \mu = q(c)
+elemDisc["mu"].set_mass_scale(0.0)
+elemDisc["mu"].set_mass(elemDisc["c"].value())
+elemDisc["mu"].set_diffusion(CahnHilliardModel.M)
+elemDisc["mu"].set_reaction(0.0)
+elemDisc["mu"].set_source(0.0)
+
+
+# Optional: Add a source term q(c) to the first equation, e.g., to model production/degradation of the species.
+# Turing source term: q(c) = ka * (1 - c) - ki * c
+if False:
+ qc = ug4.PythonUserFunction2d(lambda c: CahnHilliardModel.ka * (1 - c) - CahnHilliardModel.ki * c, 1)
+ qc.set_input_and_deriv(0, elemDisc["c"].value(), lambda c: -CahnHilliardModel.ki - CahnHilliardModel.ka) # Derivative with respect to c is -ki - ka
+ elemDisc["mu"].set_source(qc)
+
+
+# Eqn 2: \mu = f'(c) - \lambda \triangle c
+
+# create Python user function and bind input/derivative to elemDisc["c"].value()
+fPrime = ug4.PythonUserFunction2d(CahnHilliardModel.f1, 1) # f'(c) = c^3 - c
+fPrime.set_input_and_deriv(0, elemDisc["c"].value(), CahnHilliardModel.f2) # f''(c) = 3c^2 - 1
+
+elemDisc["c"].set_stationary()
+elemDisc["c"].set_diffusion(CahnHilliardModel.gamma)
+elemDisc["c"].set_reaction(fPrime)
+elemDisc["c"].set_source(elemDisc["mu"].value())
+
+
+
+# Domain discretization
+domainDisc = ug4.DomainDiscretization2dCPU1(approxSpace)
+domainDisc.add(elemDisc["mu"])
+domainDisc.add(elemDisc["c"])
+
+
+# Debug writer (optional)
+try:
+ dbgWriter = ug4.GridFunctionDebugWriter2dCPU1(approxSpace)
+except Exception:
+ dbgWriter = None
+
+
+# Solution gridfunction and initialization
+u = ug4.GridFunction2dCPU1(approxSpace)
+
+
+# u.set(0.0)
+# Initialize with random perturbation around 0.33 for c, and 0 for mu, as in original Lua
+
+# TODO: Determine the (three) steady states?
+# Which are unstable and which are stable?
+
+# Initial values
+def MyInitialValueMu(x, y, t, si):
+ return 0.0
+
+def MyInitialValueC(x, y, t, si):
+ return 0.01 * (random.random()-0.5)
+
+ug4.Interpolate(ug4.PythonUserNumber2d(MyInitialValueMu), u, "mu")
+ug4.Interpolate(ug4.PythonUserNumber2d(MyInitialValueC), u, "c")
+
+
+
+# TODO: Determine characteristic time scale and adjust endTime accordingly.
+
+startTime = 0.0
+charLength = 100.0 # Assuming a characteristic length scale of 100.0
+charTimeDiff = charLength**4/(CahnHilliardModel.M * CahnHilliardModel.gamma) # Time scale for Cahn-Hilliard
+charTimeFeat = charLength**3 # Time scale based on feature evolution,
+# assuming a characteristic velocity of 1.0 (this is a rough estimate and can be adjusted based on the specific problem setup)
+
+print (f"Characteristic time scale based on diffusion: {charTimeDiff}")
+print (f"Characteristic time scale based on reaction: {charTimeFeat}")
+
+dt = min(charTimeDiff, charTimeFeat)*1e-4
+endTime = max(charTimeDiff, charTimeFeat)*100
+
+
+
+
+
+limexDesc = util.GetLimexDefaultDesc()
+print("LIMEX config:")
+print(limexDesc)
+limexDesc["TOL"]=1e-2
+limex = util.CreateLimexIntegrator(domainDisc, limexDesc, dt)
+
+limex.set_dt_min(dt*1e-3)
+limex.set_dt_max((endTime-startTime)/50)
+
+# Create callback observer for measuring mass.
+ti = []
+mi = []
+ji = []
+def MyMassCallback(usol, step, time, dt) :
+ m=ug4.Integral(usol, "c", "Inner")
+ j=ug4.H1SemiNorm(usol, "mu", 4, "Inner")
+ j=-CahnHilliardModel.M*j*j
+ print (f"Step {step}, time {time:.4f}, mass {m:.6f}, energy change {j:.6f}")
+
+ ti.append(time)
+ mi.append(m)
+ ji.append(j)
+massobserver = ug4.PythonCallbackObserver2dCPU1(MyMassCallback)
+limex.attach_observer(massobserver)
+
+
+
+# Create callback observer for file I/O.
+def MyVTKCallback(usol, step, time, dt) :
+ ug4.WriteGridFunctionToVTK(usol, "vtk/CahnHilliardFish"+str(int(step)).zfill(5)+".vtu")
+vtkobserver = ug4.PythonCallbackObserver2dCPU1(MyVTKCallback)
+limex.attach_observer(vtkobserver)
+
+
+# Run problem!
+limex.apply(u, endTime, u, startTime)
+
+
+if True:
+ import numpy as np
+ import matplotlib.pyplot as plt
+
+ plt.plot(ti, ji, 'o-') # circles connected by lines
+ plt.xlabel('time')
+ plt.ylabel('mass')
+ plt.grid()
+ plt.show()
+
+
+print("done")
diff --git a/content/cahn-hilliard/cahn_hilliard_equation.tex b/content/cahn-hilliard/cahn_hilliard_equation.tex
new file mode 100644
index 0000000..7f2d023
--- /dev/null
+++ b/content/cahn-hilliard/cahn_hilliard_equation.tex
@@ -0,0 +1,222 @@
+
+\documentclass[11pt]{article}
+\usepackage{amsmath,amssymb}
+\usepackage[a4paper,margin=1in]{geometry}
+
+\title{Introduction to the Cahn--Hilliard Equation}
+\date{}
+
+\begin{document}
+\maketitle
+
+\section{Physical Idea}
+
+The Cahn--Hilliard equation describes phase separation in a binary mixture while conserving mass.
+Let
+\[
+c(\mathbf{x},t)
+\]
+denote the concentration of one component.
+
+Mass conservation requires
+\[
+\frac{d}{dt}\int_\Omega c\,dV = 0.
+\]
+
+\section{Free Energy Functional}
+
+The free energy is
+\[
+F[c]
+=
+\int_\Omega
+\left(
+f(c)
++
+\frac{\kappa}{2}|\nabla c|^2
+\right)
+dV.
+\]
+
+The bulk free energy is often chosen as
+\[
+f(c)=\frac{1}{4}(c^2-1)^2,
+\]
+which has minima at \(c=\pm1\).
+
+\section{Chemical Potential}
+
+The chemical potential is the variational derivative of the free energy:
+\[
+\mu=\frac{\delta F}{\delta c}.
+\]
+
+This yields
+\[
+\mu=f'(c)-\kappa\nabla^2 c.
+\]
+
+For the double-well potential,
+\[
+f'(c)=c^3-c,
+\]
+so
+\[
+\mu=c^3-c-\kappa\nabla^2 c.
+\]
+
+\section{Mass Conservation and Flux}
+
+The conservation law is
+\[
+\frac{\partial c}{\partial t}+\nabla\cdot\mathbf{J}=0.
+\]
+
+Assuming a diffusive flux
+\[
+\mathbf{J}=-M\nabla\mu,
+\]
+where \(M\) is the mobility, gives
+\[
+\frac{\partial c}{\partial t}
+=
+\nabla\cdot(M\nabla\mu).
+\]
+
+\section{Classical Cahn--Hilliard Equation}
+
+Combining the above relations gives
+\[
+\frac{\partial c}{\partial t}
+=
+\nabla\cdot
+\left(
+M\nabla
+\left(
+f'(c)-\kappa\nabla^2 c
+\right)
+\right).
+\]
+
+For constant mobility,
+\[
+\frac{\partial c}{\partial t}
+=
+M\nabla^2
+\left(
+f'(c)-\kappa\nabla^2 c
+\right).
+\]
+
+Using the quartic potential,
+\[
+\frac{\partial c}{\partial t}
+=
+M\nabla^2(c^3-c-\kappa\nabla^2 c).
+\]
+
+\section{Fourth-Order Nature}
+
+Since
+\[
+\nabla^2(\nabla^2 c)=\nabla^4 c,
+\]
+the equation is fourth order in space.
+
+A common finite-element formulation introduces
+\[
+\mu=f'(c)-\kappa\nabla^2 c,
+\]
+and solves
+\[
+\frac{\partial c}{\partial t}
+=
+\nabla\cdot(M\nabla\mu).
+\]
+
+\section{Energy Dissipation}
+
+The free energy decreases monotonically:
+\[
+\frac{dF}{dt}
+=
+-\int_\Omega M|\nabla\mu|^2\,dV
+\le 0.
+\]
+
+\section{Spinodal Decomposition}
+
+For a nearly homogeneous state
+\[
+c=c_0+\epsilon,
+\]
+perturbations grow when
+\[
+f''(c_0)<0.
+\]
+
+This leads to spontaneous phase separation.
+
+\section{Linear Stability Analysis}
+
+The growth rate of a Fourier mode is
+\[
+\lambda(k)
+=
+-Mk^2
+\left(
+f''(c_0)+\kappa k^2
+\right).
+\]
+
+Long wavelengths can grow when \(f''(c_0)<0\), while short wavelengths are stabilized by the gradient term.
+
+\section{Boundary Conditions}
+
+Typical no-flux conditions are
+\[
+\nabla\mu\cdot n=0,
+\qquad
+\nabla c\cdot n=0.
+\]
+
+\section{Relation to Allen--Cahn}
+
+Allen--Cahn:
+\[
+\partial_t c=-M\mu.
+\]
+
+Cahn--Hilliard:
+\[
+\partial_t c=\nabla\cdot(M\nabla\mu).
+\]
+
+Allen--Cahn does not conserve mass, while Cahn--Hilliard does.
+
+\section{Coupling to Mechanics}
+
+A chemo-mechanical free energy can be written as
+\[
+F=
+\int
+\left(
+f(c)
++\frac{\kappa}{2}|\nabla c|^2
++\psi_{\mathrm{elastic}}(\varepsilon,c)
+\right)dV.
+\]
+
+The chemical potential becomes
+\[
+\mu
+=
+\frac{\partial f}{\partial c}
+-\kappa\nabla^2 c
++
+\frac{\partial \psi_{\mathrm{elastic}}}{\partial c}.
+\]
+
+Applications include alloys, battery electrodes, hydrides, and phase-field models.
+
+\end{document}
diff --git a/content/cahn-hilliard/grids/fish_grid.ugx b/content/cahn-hilliard/grids/fish_grid.ugx
new file mode 100644
index 0000000..28451af
--- /dev/null
+++ b/content/cahn-hilliard/grids/fish_grid.ugx
@@ -0,0 +1,26 @@
+
+
+ -100 0 0 100 0 0 0 100 0 0 -100 0 70 20 0 70 -20 0 110 -30 0 110 30 0 20 70 0 20 -70 0 50 90 0 50 -90 0 60 120 0 60 -120 0 0 0 0 -40 35 0 -40 -35 0 -50 50 0 -50 -50 0 -50 0 0 10 -35 0 10 35 0
+ 7 1 6 1 5 6 4 7 4 1 1 5 2 8 9 3 12 2 10 8 12 10 11 9 13 11 13 3 10 2 3 11 8 4 5 9 4 5 5 14 14 4 0 15 15 8 15 2 15 14 0 16 16 9 16 3 16 14 0 17 17 2 17 15 0 18 18 3 18 16 14 19 19 0 19 15 19 16 9 20 20 14 20 5 20 16 8 21 21 14 21 4 21 15
+ 12 10 2 2 8 10 9 3 11 3 13 11 1 5 6 4 1 7 4 5 1 14 5 4 2 15 8 3 16 9 15 2 17 15 17 0 16 3 18 16 18 0 15 14 19 15 19 0 16 14 19 16 19 0 5 14 20 5 20 9 16 9 20 16 20 14 4 8 21 4 21 14 15 8 21 15 21 14
+
+
+ 14 15 16 19 20 21
+ 4 5 6 7 14 15 18 19 20 21 22 23 24 25 26 27 28 31 34 35 36 37 38 39 40 41 42 43 44 45 46
+ 0 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
+
+
+ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 17 18
+ 0 1 2 3 8 9 10 11 12 13 16 17 29 30 32 33
+
+
+
+
+
+
+
+
+ 0 0
+
+
+
diff --git a/content/cahn-hilliard/grids/fish_grid_with_eye.ugx b/content/cahn-hilliard/grids/fish_grid_with_eye.ugx
new file mode 100644
index 0000000..a829225
--- /dev/null
+++ b/content/cahn-hilliard/grids/fish_grid_with_eye.ugx
@@ -0,0 +1,30 @@
+
+
+ -100 0 0 100 0 0 0 100 0 0 -100 0 70 20 0 70 -20 0 110 -30 0 110 30 0 20 70 0 20 -70 0 50 90 0 50 -90 0 60 120 0 60 -120 0 0 0 0 -40 -35 0 -50 -50 0 -50 0 0 10 -35 0 10 35 0 -50 27 0 -43.937822064277718 23.4999998108059813 0 -43.9378223919713733 16.4999996216119733 0 -50 13 0 -56.0621782634158876 16.5000007567760925 0 -56.0621782634158876 23.4999992432239075 0 -50 50 0 -75 25 0 -25 75 0
+ 7 1 6 1 5 6 4 7 4 1 1 5 2 8 9 3 12 2 10 8 12 10 11 9 13 11 13 3 10 2 3 11 8 4 5 9 4 5 5 14 14 4 0 15 15 9 15 3 15 14 0 16 16 3 16 15 14 17 17 0 17 15 9 18 18 14 18 5 18 15 8 19 19 14 19 4 21 20 22 21 23 22 24 23 25 24 20 25 24 0 17 24 23 17 0 25 0 27 27 26 26 28 28 2 27 25 27 20 20 26 21 26 22 17 14 22 22 19 19 21 19 26 28 19 28 8
+ 12 10 2 2 8 10 9 3 11 3 13 11 1 5 6 4 1 7 4 5 1 14 5 4 3 15 9 15 3 16 15 16 0 15 14 17 15 17 0 5 14 18 5 18 9 15 9 18 15 18 14 4 8 19 4 19 14 24 0 17 24 23 17 25 24 0 0 27 25 20 25 27 20 26 27 21 26 20 22 17 23 22 17 14 22 19 14 21 22 19 26 21 19 19 26 28 19 28 8 28 8 2
+
+
+ 14 15 17 18 19
+ 4 5 6 7 14 15 18 19 20 21 22 23 24 27 28 29 30 31 32 33 34 35 36 37 54 55 60 61 62 59 58 57 56 46 45 44 47 52 53
+ 0 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
+
+
+ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 16 26 27 28
+ 0 1 2 3 8 9 10 11 12 13 16 17 25 26 48 49 50 51
+
+
+ 20 25 21 22 24 23
+ 43 38 42 41 40 39
+
+
+
+
+
+
+
+
+ 0 0
+
+
+
diff --git a/content/cahn-hilliard/grids/laplace_sample_grid_2d_l50.ugx b/content/cahn-hilliard/grids/laplace_sample_grid_2d_l50.ugx
new file mode 100644
index 0000000..728698f
--- /dev/null
+++ b/content/cahn-hilliard/grids/laplace_sample_grid_2d_l50.ugx
@@ -0,0 +1,30 @@
+
+
+ 100 100 0 100 -100 0 -100 -100 0 -100 100 0 0 100 0 -100 0 0 0 -100 0 100 0 0 0 0 0
+ 0 4 4 3 3 5 5 2 2 6 6 1 1 7 7 0 4 8 8 7 5 8 6 8
+ 0 4 8 7 3 5 8 4 2 6 8 5 1 7 8 6
+
+
+ 8
+ 8 9 10 11
+ 0 1 2 3
+
+
+ 0 3 2 1 4 5 6 7
+ 0 1 2 3 4 5 6 7
+
+
+
+
+
+
+
+ 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1
+ 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10 1 11 1
+ 0 1 1 1 2 1 3 1
+
+
+ 0 0
+
+
+