next up previous contents
Next: Appendix Up: Illuminator Distributed Visualization Library Previous: Developer's Manual   Contents

Example Transient Cahn-Hilliard chts.c

The example program provided with Illuminator solves the Cahn-Hilliard equation in 3-D. The implementation is split into files cahnhill.c (appendix H) containing all of the free energy functions, derivatives, etc.; and chts.c (appendix G) with main(), initialization, callbacks--essentially, the overhead of the program; and a small header with the data structure in cahnhill.h (currently not included in the documentation because of a bug in the cxref documentation generation system). The calls to Illuminator library functions are contained in chts.c.

The idea behind Cahn-Hilliard is that we have a concentration field $C$, with associated thermodynamic free energy density $f$ given by

\begin{displaymath}
f = \frac{\alpha}{2}\left\vert C_{,i}\right\vert^2 + \beta\Psi(C),
\end{displaymath} (4.1)

where the two terms are the gradient penalty and homogeneous free energy $\Psi$, and $\alpha$ and $\beta$ are model constants. This functional form was first used to explain the finite domain size during the initial phase separation in spinodal decomposition. A typical homogeneous free energy function (and the one used in cahnhill.c) looks like
\begin{displaymath}
\Psi = C^2 (1-C)^2.
\end{displaymath} (4.2)

This functional form gives two stable equilibria at $C=0$ and $C=1$, and an unstable equilibrium at $C=\frac{1}{2}$.

The total free energy of the system occupying the body $\Omega$ is

\begin{displaymath}
{\cal F} = \int_\Omega f dV.
\end{displaymath} (4.3)

The chemical potential $\mu$ is given by
\begin{displaymath}
\mu = \frac{\delta{\cal F}}{\delta C} = -\alpha C_{,ii} + \beta\Psi'(C),
\end{displaymath} (4.4)

which leads to the transport equation
\begin{displaymath}
\dot{C} = \left(\kappa(C)\mu_{,i}\right)_{,i}
\end{displaymath} (4.5)

where $\kappa$ is the mobility.

The interface thickness $\epsilon$ is on the order of $\sqrt{\alpha/\beta}$, and the surface energy $\gamma$ on the order of $\sqrt{\alpha\beta}$, so we can set $\alpha=\gamma\epsilon$ and $\beta=\gamma/\epsilon$, with constants which we'll worry about later.

Returning to equation 4.5, it expands to

\begin{displaymath}
\dot{C} = \kappa_{,i}\mu_{,i} + \kappa\mu_{,ii}
\end{displaymath} (4.6)

and further in terms of $C$ to
\begin{displaymath}
\dot{C} = \kappa_{,i}(-\alpha C_{,jj} + \beta\Psi')_{,i}
- \kappa\alpha C_{,iijj} +
\kappa\beta(\Psi'(C))_{,ii}
\end{displaymath} (4.7)

Note this assumes constant $\alpha$ and $\beta$, which we'll have to relax if surface tension depends on another field.

Now we turn to the PETSc implementation, using distributed arrays and TS timestepping solvers. The nonlinear timestepping code uses runtime-selectable timestepping solvers (explicit, implicit, or default Crank-Nicholson) to solve the equation

\begin{displaymath}
\dot{u_i} = F_i(u_j),
\end{displaymath} (4.8)

where $u_i$ here is the vector of unknowns. This is converted by PETSc into a system of nonlinear equations according to the solver type and then solved using its SNES solvers, but we must provide subroutines to calculate $F_i$ and the derivatives $\frac{\partial F_i}{\partial u_j}$ which are used to calculate the Jacobian of the nonlinear system (unless running matrix-free using the -snes_mf option). Starting for now with constant $\kappa$, $\alpha$ and $\beta$, the $F_i$ are given by
\begin{displaymath}
F_i = \kappa\left(-\alpha C_{,iijj} + \beta(\Psi'(C))_{,ii}\right)
\end{displaymath} (4.9)

Let's suppose $C$ is on a 2-D finite difference mesh with uniform but possibly different spacings $h_x$ and $h_y$ in the $x$- and $y$-directions, so we'll let $C_{x,y}$ be the value of $C$ at coordinates $(xh_x, yh_y)$.

Starting with the $\beta$ term, the Laplacian $(\Psi'(C))_{,ii}$ at $(x,y)$ can be approximated using the standard 5-node finite difference stencil:

\begin{displaymath}
(\Psi'(C))_{,ii} \simeq
\frac{\Psi'(C_{x-1,y})-2\Psi'(C_{x...
...frac{\Psi'(C_{x,y-1})-2\Psi'(C_{x,y})+\Psi'(C_{x,y+1})}{h_y^2}
\end{displaymath} (4.10)

or expressed slightly differently, as the sum of the terms:
\begin{displaymath}
\begin{array}[h]{ccc}
& \Psi'(C_{x,y+1})h_y^{-2} & \\
\P...
...{x+1,y})h_x^{-2} \\
& \Psi'(C_{x,y-1})h_y^{-2} &
\end{array}\end{displaymath} (4.11)

So the product of $\kappa\beta$ and this is one part of the function $F_i$, and that part of the derivative $\frac{\partial F_i}{\partial C_j}$ is simply the appropriate coefficient from equation 4.11 times the second derivative $\Psi''(C)$ evaluated at the appropriate point.

For the $\alpha$-term, the Laplacian of the Laplacian is a bit more messy. Using the notation in equation 4.11 but only in the first quadrant (the coefficients are symmetric), the term will be $-\kappa\alpha$ times:

\begin{displaymath}
\begin{array}[h]{ccccc}
& & C_{x,y+2}h_y^{-4} & & \\
& ....
...h_x^{-4} \\
& ... & ... & ... & \\
& & ... & &
\end{array}\end{displaymath} (4.12)

(ellipsis indicates symmetry). This is quite a bit more complicated, but at least it's linear so the derivatives $\frac{\partial F_i}{\partial C_j}$ are constant.

These $F_i$ functions are calculated by the function ch_residual_2d() (appendix H.1.9) and assembled into a PETSc vector in ch_residual_vector_2d() (appendix H.1.5). The derivatives are calculated in two parts: the $\alpha$ term's derivative matrix is built at the start of the run in ch_jacobian_alpha_2d() (appendix H.1.3), and with each nonlinear iteration, the $\beta$ term's derivative matrix is added to that in ch_jacobian_2d() (appendix H.1.1).

Note that this is all in 2-D, the 3-D version is left as an exercise to the reader; though it's already coded in the corresponding _3d functions in cahnhill.c.


next up previous contents
Next: Appendix Up: Illuminator Distributed Visualization Library Previous: Developer's Manual   Contents
root 2002-06-18