Implementation of a Cahn-Hilliard model

Objective of this 1st tutorial

The objective of this tutorial is to implement a kernel to simulate the Cahn-Hilliard equation (Eq. (508) with chemical potential defined by Eq. (507)). The starting point is the kernel ADE_for_CH-Tutorial which simulates the Advection-Diffusion Equation (ADE) with the lattice Boltzmann method. This tutorial provides all code lines to write in each file for simulating three test cases with the Cahn-Hilliard model: serpentine, spinodal decomposition and nucleation. The verification tests will be carried out with the input files (.ini) contained in the folder Tutorial_Cahn-Hilliard/ in run_training_lbm.

It is assumed that you have already downloaded LBM_Saclay_Rech-Dev and you are working on your own branch (see Before programming: create a new branch from master).

1. Create a new kernel

Copy the ADE kernel template and rename all files with new extension _CH

Copy the folder ADE_for_CH-Tutorial which is contained in kernels/Templates_for_Tutorials, and rename it Cahn-Hilliard_Tutorial in the folder kernels

Commands
$ cd LBM_Saclay_Rech-Dev/src/kernels
$ cp -r Templates_for_Tutorials/ADE_for_CH-Tutorial ./Cahn-Hilliard_Tutorial
Remarks

Your new developments will be performed in the folder Cahn-Hilliard_Tutorial

The folder name Cahn-Hilliard_Tutorial will be used for compilation

All files have the extension _ADE for Advection Diffusion Equation. Rename them with _CH for Cahn-Hilliard.

Edit them and change strings _ADE with _CH

Change strings _ADE with _CH
- Open all files in your favorite editor (e.g. geany or codium).
- Search in all files the strings _ADE and replace them with _CH (Ctrl H with geany).
Change the problem name: in file Setup.h modify the problem name ADE with CH
```
register_problem<Problem, ModelParams::Tag2Quadra>("ADE", "base");
```
Modification with CH
register_problem<Problem, ModelParams::Tag2Quadra>("CH", "base");
Remark

Your new problem will be named CH inside your .ini input file i.e. in section [lbm] your problem must be:
[lbm] problem=CH

Compile your new kernel

Generate your makefile (only once)

Compile

Run with input files of ADE

Execute your new kernel CH of LBM_Saclay either with the simple diffusion of gaussian or with the advected gaussian in the folder run_training_lbm/Tutorial_Cahn-Hilliard

Check the solutions with paraview.

2. Modifications inside each file for Cahn-Hilliard

Advice

During your implementation, it is strongly recommended to compile regularly (make or make -j 22) to detect any typos.

For compilation errors, it is useful to write the terminal outputs inside a file e.g. compil.log:

$ make -j 22 2>&1 | tee compil.log

or alternatively compil2.log

$ make VERBOSE=1 2>&1 | tee compil2.log

File Index_CH.h

Declaration of new fields

In enum ComponentIndex, only the indices for macroscopic fields such as concentration $c$ (IC) and velocity components $u_x$ (IU), $u_y$ (IV) & $u_z$ (IW) are declared:

enum ComponentIndex {
  IU     , /*!< X velocity / momentum index */
  IV     , /*!< Y velocity / momentum index */
  IW     , /*!< Z velocity / momentum index */
  IC     , /*!< Phase field index */
  COMPONENT_SIZE /*!< invalid index, just counting number of fields */
};

Add the new indices for chemical potential IMU and laplacian ILAPLAC.

Add new outputs

In struct index2names only the strings for concentration and velocity components are written:

map[IC ] = "comp";
map[IU ] = "vx";
map[IV ] = "vy";
map[IW ] = "vz";

Add a new output for chemical potential.

Solution
map[IC ] = "comp";
map[IU ] = "vx";
map[IV ] = "vy";
map[IW ] = "vz";
map[IMU] = "mu";
Remark

The keyword mu will be used in the .ini input file to write the chemical potential.

File Models_CH.h

Read and declare Cahn-Hilliard parameters in ModelParams

The Cahn-Hillard equation involves three parameters: the surface tension $\sigma$, the interface thickness $W$ and the mobility $\mathcal{M}_\phi$.

Read and declare Cahn-Hilliard parameters in ModelParams

Add specific functions for Cahn-Hilliard

Add a new function for derivative of double-well g_prime defined by the first term inside the bracket of Eq. (507).

Add a new function M2 to compute the chemical potential $\mu$ (MU) defined by Eq. (507)

Solution function M2
// =======================================================
// second order moment of feq
KOKKOS_INLINE_FUNCTION
       real_t M2(EquationTag1 tag, const LBMState &lbmState) const {
       const real_t mu = sigma * 1.5 / W * (g_prime(lbmState[IC]) - SQR(W) * lbmState[ILAPLAC]);
  return mu;
    }
Remarks

IC & ILAPLAC have been declared in Index_CH.h.

The parameters sigma (surface tension $\sigma$) and W (interface thickness $W$) have been read and declared inside ModelParams (previous box).

The new function g_prime has been defined just above

Add a new function tauM for collision rate using the mobility $\mathcal{M}_\phi$ of Cahn-Hilliard equation

Solution tauM
// =======================================================
// relaxation coef for LBM scheme of phase field equation
KOKKOS_INLINE_FUNCTION
       real_t tauM (EquationTag1 tag, const LBMState &lbmState) const {
            real_t tau = 0.5 + (e2 * mobility * dt / SQR(dx));
  return (tau);
     }
Remark

The parameter mobility has been read and declared in ModelParams.

You could have used the parameter D as the mobility coefficient. In that case, that new function is useless. But take care, especially if a future coupling with a transport equation is planned.

Add the standard hyperbolic tangent function phi0 defined by Eq. (385) (for initial condition)

Add conditions and keywords for several initial conditions

Currently only one condition exits for intial condition:

if (initTypeStr == "gaussian") initType = ADE_GAUSSIAN;

Add the necessary conditions for three Cahn-Hilliard initializations

File LBMScheme_CH.h

Function setup_collider with BGK

The function setup_collider is written for solving Advection-Diffusion Equation. It must be modified for solving Cahn-Hilliard equation.

First, the 2nd order moment, is currently $c$

const real_t M0 = Model.MomentM0(tag, lbmState);
    const RVect<dim> uc = Model.MomentM1<dim>(tag, lbmState);
    const real_t M2 = Model.MomentM2(tag, lbmState);

Modify the call of function to compute $\mu$ instead of $c$

Next, the appropriate function must be called for collision rate to take into account the mobility $\mathcal{M}_\phi$

 // compute collision rate
collider.tau = Model.tau(tag, lbmState);

Finally, the equilibrium distribution function is currently written for ADE (lines 4 & 10)

         // ipop = 0
         collider.f[0]   = Base::get_f_val(tag, IJK, 0);
         collider.S0[0]  = 0.0 ;
         collider.feq[0] = w[0]*M0;

         // ipop > 0
         for (int ipop = 1; ipop < npop; ++ipop) {
         collider.f[ipop]   = this->get_f_val(tag, IJK, ipop);
         collider.S0[ipop]  = 0.0 ;
         collider.feq[ipop] = w[ipop] * (M2 + c_cs2 * Base::compute_scal(ipop, uc));

Modify it to take into account $c$ (moment 0) and $\mu$ (moment 2) (see SMEMaG course)

Function init_macro

Only one initial condition exists for ADE with ADE_GAUSSIAN:

 if (Model.initType == ADE_GAUSSIAN) {
  c  = Model.ampl * exp( - ( SQR(x-Model.x0) / (2 * SQR(Model.sigmax)) + SQR(y-Model.y0) / (2 * SQR(Model.sigmay)) )) ;
  vx = Model.initVX;
  vy = Model.initVY;
}

Write three initial conditions corresponding to keywords SERPENTINE, RANDOM_SPINODAL and RANDOM_NUCLEATION

Function update_macro

Currently, only the concentration and the velocity are updated:

const real_t c  = moment_c;
const real_t vx = Model.initVX;
const real_t vy = Model.initVY;

this->set_lbm_val(IJK, IC , c );
this->set_lbm_val(IJK, IU , vx );
this->set_lbm_val(IJK, IV , vy );

Add for chemical potential and add a condition for velocity of SERPENTINE option

Function update_macro_grad

Currently the function is empty because for ADE, there is no need to compute additional gradients or laplacian with the LB method.

KOKKOS_INLINE_FUNCTION
void update_macro_grad(IVect<dim> IJK) const {

}

This is not the case for Cahn-Hilliard equation because a laplacian is involved in the definition of chemical potential. Add the computation of laplacian.

File InitConditionsTypes.h

Add the keywords of initial conditions

In file InitConditionsTypes.h, two keywords are currently declared:

enum InitCondition{
        PHASE_FIELD_INIT_UNDEFINED,
        ADE_GAUSSIAN
};

Add the three new keywords for Cahn-Hilliard: SERPENTINE, RANDOM_SPINODAL and RANDOM_NUCLEATION.

3. Verifications of your implementation

Three .ini input files are available in the directory Tutorial_Cahn-Hilliard to check your implementation for each initial condition:

Tuto-CH_Test01-Serpentine.ini

Tuto-CH_Test02-Spinodal.ini

Tuto-CH_Test03-Nucleation.ini

It is recommended to start with Tuto-CH_Test01-Serpentine.ini because comparisons can be done with csv files contained in the folder Contours_CAC.

Run Tuto-CH_Test01-Serpentine.ini

Run Tuto-CH_Test01-Serpentine.ini

Commands
Go to the appropriate folder
$ cd run_training_lbm/Tutorial_Cahn-Hilliard
and run LBM_Saclay
$ ../../build_CH/src/LBM_saclay Tuto-CH_Test01-Serpentine.ini

Once the simulations is over, post-process your results with paraview.

Post-process with paraview

Open paraview (e.g. version 5.11) and compare the contours from Cahn-Hilliard .vti files and csv files of subfolder Contours_CAC.

Commands for contours

Open all .vti files and select comp

Ctrl space and Cell Data to Point Data and Apply

Clic on contour and select field comp with value 0.5 and Apply

Commands for csv files

Open csv files (choose option CSV Reader) + OK + Apply

Crtl space + Table To Points + select Points:0 for X Column and Points:2 for Y Column

Play with Point Size in section Styling

Exercise

Make a video with paraview for the two other test cases:

Spinodal decomposition: Tuto-CH_Test02-Spinodal.ini

Nucleation: Tuto-CH_Test03-Nucleation.ini

4. Push your developments on Codev-Tuleap

Save your developments on codev-tuleap.cea.fr

Once your verifications are carried out, save only your implementation (not your output .vti files). See commands on After programming: push your developments

Section author: Alain Cartalade