File vectorphase.c

For vector-valued phase field, there are two field variables p and q, and three temporary fields i>phi, epsilon and tau. Throughout this function, 2-D is assumed, since this formulation is really only good for 2-D at this point (though that may change).

void vectorphase_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 vectorphasers structure this inserts parameters from the command line.

Global Function vectorphase_interior_line_function()

This calculates the time derivatives dp/dt and dq/dt for an interior line. This will eventually use the free energy of Kobayashi, Warren and Carter and resulting dynamics in terms of magnitude phi and angle theta as presented in Physica D 140 (2000), 141.

Unfortunately, there is a singularity in theta diffusivity where its gradient is zero.

Giga and Giga and Kobayashi and Giga have outlined a way to solve this, but we haven't yet implemented it. So for now, I'm going with KWC's original 1998 formulation even though it's not quite right. That formulation begins with a slightly different free energy based on angle theta which is orientation times the rotational symmetry order N₀. The free energy is given by three terms: the homogeneous free energy F(phi;m) (with derivative f), the gradient penalty term in phi with coefficient epsilon²/2, and the gradient penalty term in theta with coefficient phi²mu(phi)/2.

We can then write the variation in energy as two terms A and B which are coefficients on the variations of phi and theta respectively. As is typical for phase field, these variations lead directly to equations for dynamics of phi and theta such that their time derivatives are A and B/phi respectively then both divided by a kinetic parameter tau.

The really neat thing is that we can turn these dynamics into those for the vector field using the simple transformation where p=phi cos(theta) and q=phi sin(theta), which leads to straightforward expressions for dp/dt and dq/dt in terms of the phi and theta time derivatives.

void vectorphase_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 variables (phi, epsilon, tau).
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 vectorphasers struct typedef, from which this gets needed parameters.

The A variation has three terms: the gradient penalty, the homogeneous free energy, and a derivative related to the orientation gradient, which are labeled A₁, A₂ and A₃ respectively. The first is the divergence of the "flux" of phi with epsilon² as its diffusivity.
The second term A₂ is the derivative of the homogeneous free energy with respect to phi. Note that it has local minima at 0 and 1, and if m>0 then the disordered phi=0 phase is stable and the ordered phase metastable, and if m<0 then the ordered phase is stable and the disordered phase metastable.
The third term A₃ is given by the same derivative as above. Note that the vector component spatial derivatives are calculated by central differencing here. The variables gradx and grady represent beta_x phi deltax and beta_y phi delaty respectively.
Because I'm omitting the epsilon' term, the B variation has just one term, given by the divergence of phi mu(phi) beta divided by phi.
These terms come together to produce the time derivatives of the vector components dp/dt and dq/dt.

Global Function vectorphase_labels_initcond()

This sets up the parameters for this equation in the vectorphasers struct typedef, the field variable labels, maximum stable explicit timestep size, and initial condition for the vector phase field variables. Model parameterss include the energy parameter m (-mparam), the ratio epsilon₀/dx (-epsilon), parameters for the phi-theta coupling function mu (-mubar and -muexp), general symmetry order N₀ (-n_general), surface energy anisotropy parameters and kinetic anisotropy parameters which are discussed in vectorphase_temp_parameters_line().

void vectorphase_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 vectorphasers 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 explicit finite difference criterion which seems to work is one-quarter of the square of the mesh spacing.
The initial condition is generally circular nuclei, with some scattered completely randomly (controlled by parameter -nukes), some semi-circles on the edges, for simulations of competitive growth from edges (controlled by parameter -edge_nukes), and some circles limited to the interior for similar simulations (controlled by parameter -interior_nukes).

Global Function vectorphase_temp_parameters_line()

This calculates the temporary field variables in one line of the global array for vector phase field stuff. In this case, we want to calculate the order parameter phi, which is the magnitude of the vector phase field, the orientation-dependent gradient penalty coefficient epsilon which determines interfacial energy and thickness, and orientation-dependent kinetic parameter tau which determines the growth rate.

void vectorphase_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 vectorphasers struct typedef, from which this gets needed parameters.

The order parameter phi is the most straightforward, it's just the magnitude of the vector phase variable.
The gradient penalty coefficient epsilon is quite a bit more complicated, as it's a function of the angle psi which represents the misorientation between the gradient of phi and the crystal orientation theta-bar = theta/N₀. As an intermediate, we'll call psi-bar the angle of the order parameter gradient, as long as there's 2-fold symmetry this can be represented by psi-bar = tan^-1[(dphi/dy)/(dphi/dx)]. (Actually, since we use the C math library atan2 function, there's no 2-fold symmetry restriction in the code.)

The functional form of epsilon given in KWC 1998 is available in the LaTeX-based documentation.

Note that to calculate psi we need nearest-neighbor values of phi, but the function needs nearest neighbor values of epsilon. So this forces us to look to next-nearest-neighbors, expanding the required stencil width for these calculations to 2. Make sure these data exist at the edges, or you'll get some wierd behavior at the interfaces!

Note the abstraction violation here: the code calculating epsilon assumes that phi is already calculated in the negative y direction nearest neighbor(s), but not in the positive direction. This could therefore fail if for some reason the loops go the other way! A longer-term fix would allow two levels of temporary parameters, but it's not clear it's worth this just to save the square roots which are pretty simple. Also note that this code assumes equal spacings in both grid directions.
The kinetic timescale parameter tau is also a function of psi which you can see in the LaTeX-based documentation.