RCS Header: /cvsroot/rheoplast/cahnhill.c,v 1.37 2004/08/18 21:19:52 hazelsct Exp
This provides a simple Cahn-Hilliard module for Rheoplast. The free energy
goes as
polymer solution free energy given by
the integral of a gradient penalty and a homogeneous free energy, the
latter of which in the polymer option 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_interior_line_function"
#define __FUNCT__ "cahnhill_interior_line_jacobian"
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().
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.
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.
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:
-
Small random fluctuations are selected using the
-ch_random_center and -ch_random_fluct
options. These produce a uniform distribution centered at the
``center'' value and with width twice the fluctuation value.
-
A two-layer initial condition in the
y-direction is chosen using the -ch_layers
option, whose argument is the fraction of the domain which will be
in the ``bottom''
C=1 region. -
The
-ch_trilayer
option is used for the ``metal-electrolyte-metal'' simulation. This
creates three layers in the
y-direction
with width half that of the domain, which is to say, one-quarter of
that of the whole symmetric domain.
-
The
-ch_particles
option is used for the ``colliding particles'' simulation. This
creates a symmetric simulation with a square particle centered on
the middle of the
x-axis
with width half that of the domain, which is to say, one-quarter of
that of the whole symmetric domain.
-
The default initial condition is a square (cube) with side length
equal to half of the domain width, either in the center, or if all
symmetries are on, then at the origin.
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.
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.
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.
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.