File cahnhill.c

RCS Header: /cvsroot/rheoplast/cahnhill.c,v 1.55 2006/02/28 21:00:07 wanida Exp


This provides a simple Cahn-Hilliard module for Rheoplast. The free energy goes as the integral of a gradient penalty and a homogeneous free energy, the latter of which is C2 (1-C)2, or in the polymer option selected by -ch_polymer, is given by:

beta Psi(C) = (1/m)C ln C + (1-C) ln (1-C) + chi C(1-C).
The chemical potential is then the variation of this free energy, given by the homogeneous free energy derivative term and a laplacian term corresponding to the gradient penalty.

Included Files

Preprocessor definitions

#define __FUNCT__ "psiprime"

#define __FUNCT__ "psidoubleprime"

#define __FUNCT__ "cahnhill_first_setup"

#define __FUNCT__ "cahnhill_labels_initcond"

#define SMALL_PRIME 1571

#define MEDIUM_PRIME 524287

#define LARGE_PRIME 2147483647

#define MY_RANDOM( pseudo_random )

#define C( point )

#define Cfunc( point )

#define mu( point )

#define vu( point )

#define vv( point )

#define pu( point )

#define pv( point )

#define pw( point )

#define V( point )

#define sigma( point )

#define sigmaprime( point )

#define sigma_eff( point )

#define __FUNCT__ "cahnhill_temp_parameters_line"

#define __FUNCT__ "cahnhill_temp_parameters_boundary_line"

#define __FUNCT__ "cahnhill_interior_line_function"

#define __FUNCT__ "cahnhill_boundary_line_function"

#define __FUNCT__ "cahnhill_interior_line_jacobian"

Global Function cahnhill_boundary_line_function()

This calculates the time derivatives and constraint functions for the Cahn-Hilliard equations for a boundary line. Note the shearstrain module interaction programmed here;

void cahnhill_boundary_line_function ( PetscScalar* x, PetscScalar* func, PetscScalar* temp, PetscTruth** mixed_constraints, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data, int side )

PetscScalar* x
The field variables from which to evaluate the function.
PetscScalar* func
Where to put the evaluated function.
PetscScalar* temp
Array of temporary field variables.
PetscTruth** mixed_constraints
Arrays of boolean variables indicating constraint equations in mixed timestep-constraint fields.
int points
Number of points to evaluate at.
int gxm
The x-width of the ``local'' vector's array, including shadow nodes, for the y-increment.
int gym
The y-width of the ``local'' vector's array, including shadow nodes, for the z-increment.
PetscScalar xmin
First node x-coordinate.
PetscScalar xmax
Last node plus one x-coordinate.
PetscScalar ycoord
This line y-coordinate.
PetscScalar zcoord
This line z-coordinate.
PetscScalar time
Current simulation time.
AppCtx* data
Pointer to the main simulation parameter structure, which includes the vortparm struct typedef, from which this gets needed parameters.
int side
 
void cahnhill_interior_line_function It returns nothing.
It includes a convective term with first-order upwinding for velocity-vorticity flow.
And a convective term with first-order upwinding for velocity-pressure flow.
It includes a convective term with first-order upwinding for velocity-vorticity flow.
And a convective term with first-order upwinding for velocity-pressure flow.

Global Function cahnhill_first_setup()

The basic setup, assigning the number of solved and temporary field variables, the stencil width, and using options to set the parameters in the chparm struct typedef.

void cahnhill_first_setup ( PetscTruth threedee, int* vars, int* tempvars, int* stencilwid, AppCtx* data )

PetscTruth threedee
Request support for 3-D.
int* vars
Pointer to the number of solved field variables.
int* tempvars
Pointer to the number of temporary field variables.
int* stencilwid
Pointer to the stencil width.
AppCtx* data
Pointer to the AppCtx struct typedef, into whose chparm structure this inserts parameters from the command line.
The standard model temporarily sets alpha and beta to epsilon/dx and sigma, and mobility to the dummy value -1, so that if not overridden by the user parameter setting, its default value of epsilon2 can be set in cahnhill_labels_initcond().

Global Function cahnhill_interior_line_function()

This calculates the time derivative of C using the divergence of flux, which goes down the gradient of chemical potential described above, with mobility kappa.

void cahnhill_interior_line_function ( PetscScalar* x, PetscScalar* func, PetscScalar* temp, PetscTruth** mixed_constraints, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )

PetscScalar* x
The field variables from which to evaluate the function.
PetscScalar* func
Where to put the evaluated function.
PetscScalar* temp
Array of temporary field variables.
PetscTruth** mixed_constraints
Arrays of boolean variables indicating constraint equations in mixed timestep-constraint fields.
int points
Number of points to evaluate at.
int gxm
The x-width of the ``local'' vector's array, including shadow nodes, for the y-increment.
int gym
The y-width of the ``local'' vector's array, including shadow nodes, for the z-increment.
PetscScalar xmin
First node x-coordinate.
PetscScalar xmax
Last node plus one x-coordinate.
PetscScalar ycoord
This line y-coordinate.
PetscScalar zcoord
This line z-coordinate.
PetscScalar time
Current simulation time.
AppCtx* data
Pointer to the main simulation parameter structure, which includes the chparm struct typedef, from which this gets needed parameters.
For now this assumes uniform kappa.
It includes a convective term with first-order upwinding for velocity-vorticity flow.
And a convective term with first-order upwinding for velocity-pressure flow.

Global Function cahnhill_interior_line_jacobian()

This calculates the Jacobian of the equations corresponding to the Cahn-Hilliard variables.

void cahnhill_interior_line_jacobian ( PetscScalar* x, PetscScalar* temp, Mat J, int points, int gxm, int gym, int firstrow, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )

PetscScalar* x
The field variables from which to evaluate the Jacobian.
PetscScalar* temp
Array of temporary field variables.
Mat J
Where to put the evaluated Jacobian.
int points
Number of points to evaluate at.
int gxm
The x-width of the ``local'' vector's array, including shadow nodes, for the y-increment.
int gym
The y-width of the ``local'' vector's array, including shadow nodes, for the z-increment.
int firstrow
The matrix row number corresponding to the first point in the line.
PetscScalar xmin
First node x-coordinate.
PetscScalar xmax
Last node plus one x-coordinate.
PetscScalar ycoord
This line y-coordinate.
PetscScalar zcoord
This line z-coordinate.
PetscScalar time
Current simulation time.
AppCtx* data
Pointer to the main simulation parameter structure, which includes the chparm struct typedef, from which this gets needed parameters.
First this calculates the fixed coefficients in the alpha gradient penalty term of the transport equation. The variable jvalue will hold the Jacobian values for insertion into the matrix. These are fixed in eight of the thirteen non-zeroes (18 of the 25 in 3-D); the other five (seven in 3-D) depend on the second derivative of the Psi function with respect to C, as discussed in the psidoubleprime() function in this file. The fixed Jacobian values are stared right away in the jvalue array, and fixed parts of the C-dependent Jacobian values are temporarily stored in the fvalue array for subsequent insertion into jvalue and addition to the variable parts.

Global Function cahnhill_labels_initcond()

This sets up the field variable labels, maximum stable explicit deltat, and initial condition for the Cahn-Hilliard variables.

void cahnhill_labels_initcond ( PetscScalar* globalarray, int nx, int ny, int nz, int xm, int ym, int zm, int xs, int ys, int zs, int vars, AppCtx* data, PetscScalar* max_explicit_deltat )

PetscScalar* globalarray
The global field array.
int nx
Overall x-width of the global array.
int ny
Overall y-width of the global array.
int nz
Overall z-width of the global array.
int xm
The x-width of the local part of the array.
int ym
The y-width of the local part of the array.
int zm
The z-width of the local part of the array.
int xs
The (integer) x-coordinate of the start of the local part of the array.
int ys
The (integer) y-coordinate of the start of the local part of the array.
int zs
The (integer) z-coordinate of the start of the local part of the array.
int vars
Total number of field variables to be solved.
AppCtx* data
Pointer to the AppCtx struct typedef, whose chparm structure this uses for various purposes.
PetscScalar* max_explicit_deltat
Pointer to the largest allowable explicit timestep size for this equation, which this function can set/modify.
The standard model sets the interface thickness to thrice the minimum grid spacing and surface tension to 1 by default. These are controlled by command line options, either -ch_intwidth and -ch_surftens (intwidth is multiplied by dx), or by setting the homogeneous and gradient penalty coefficients themselves using -ch_alpha and -ch_beta. The constants in the code (3.1 and square root of 18) get the surface tension and interface thickness correct for the standard model; different values are needed to get the polymer model correct (though arguably, in the polymer model one should set beta to 1 and use alpha to adjust the interface thickness).
By default, kappa is set to epsilon2/beta, such that the diffusion timescale over the interface thickness is constant; this is controlled by command line option -ch_mobility.
Default polymer Flory-Huggins parameters are chi=0.58, m=640, as discussed in the psiprime() function in this file.
The explicit finite difference timestep size used here is delta x/40 kappa, which works for the fourth-order polynomial free energy in 2-D when nx=ny, epsilon=dx ({\tt intwidth}=1), and sigma=1.
For random fluctuations in this module, it's necessary to generate different random sequences on each process, even though rand() will return the same sequences. So we create a random counter which takes on values between zero and SMALL_PRIME-1; this is initialized to rank times LARGE_PRIME (which should make it sort of random), and incremented by MEDIUM_PRIME then re-moduloed to SMALL_PRIME for each random number generated. This counter, in turn, is divided by SMALL_PRIME and the result added to rand()/RAND_MAX then modulo 1, in order to give a "different random number" (at least at the resolution of SMALL_PRIME) in each process. This is not a great algorithm, but should do for the purpose needed here.
Eventually it might be good to make this resource available to the whole code.
If we're not loading in data as the initial condition, the C field is initiated here. There are several types of initial conditions here which are chosen by command-line parameters:
If the option -ch_random_fluct is selected without -ch_random_center, then random fluctuations with uniform distribution of width twice the fluctuation value are added to any other initial condition present.

Global Function cahnhill_temp_parameters_boundary_line()

void cahnhill_temp_parameters_boundary_line ( PetscScalar* x, PetscScalar* temp, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data, int side )

Global Function cahnhill_temp_parameters_line()

This calculates the Cahn-Hilliard temporary parameter mu, which is the chemical potential given by the equation in the intro to this file.

void cahnhill_temp_parameters_line ( PetscScalar* x, PetscScalar* temp, int points, int gxm, int gym, PetscScalar xmin, PetscScalar xmax, PetscScalar ycoord, PetscScalar zcoord, PetscScalar time, AppCtx* data )

PetscScalar* x
Array with the "real" field variables.
PetscScalar* temp
Array with the temporary field variables.
int points
Number of points at which to calculate the temporary variables.
int gxm
The x-width of the ``local'' vector's array, including shadow nodes, for the y-increment.
int gym
The y-width of the ``local'' vector's array, including shadow nodes, for the z-increment.
PetscScalar xmin
First node x-coordinate.
PetscScalar xmax
Last node plus one x-coordinate.
PetscScalar ycoord
This line y-coordinate.
PetscScalar zcoord
This line z-coordinate.
PetscScalar time
Current simulation time.
AppCtx* data
Pointer to the main simulation parameter structure, which includes the chparm struct typedef, from which this gets needed parameters.

Local Function psidoubleprime()

This abstracts out the function for Psi''(C), the second derivative of homogeneous free energy, so it can be easily modified. Since Psi(C) = C2 (1-C)2 = C4 - 4C3 + C2, this returns 12C2 - 12C + 2.

There is also a polymer thermo option based on Flory-Huggins thermodynamics described above. Enable it using option -ch_polymer and control it using options -ch_polymer_chi and -ch_polymer_m (defaults: 0.58, 640). Note that with m=1 this reduces to the regular solution model.

static inline PetscScalar psidoubleprime ( PetscScalar C, chparm* thecahnhill )

PetscScalar psidoubleprime
It returns the second derivative of homogeneous free energy.
PetscScalar C
The C parameter it's a function of.
chparm* thecahnhill
Cahn-Hilliard parameter structure.

Local Function psiprime()

This abstracts out the function for Psi'(C), the derivative of homogeneous free energy, so it can be easily modified. Since Psi(C) = C2 (1-C)2 = C4 - 4C3 + C2, this returns 4C3 - 6C2 + 2C.

There is also a polymer thermo option based on Flory-Huggins thermodynamics described above. Enable it using option -ch_polymer and control it using options -ch_polymer_chi and -ch_polymer_m (defaults: 0.58, 640). Note that with m=1 this reduces to the regular solution model.

static inline PetscScalar psiprime ( PetscScalar C, chparm* thecahnhill )

PetscScalar psiprime
It returns the derivative of homogeneous free energy.
PetscScalar C
The C parameter it's a function of.
chparm* thecahnhill
Cahn-Hilliard parameter structure.