RCS Header: /cvsroot/petscgraphics/cahnhill.c,v 1.8 2003/05/21 22:52:29 hazelsct Exp
This file contains the heart of the Cahn-Hilliard formulation, in particular
the functions which build the finite difference residuals and Jacobian.
Included Files
- #include "cahnhill.h"
- #include </usr/lib/petsc/include/petscts.h>
- #include </usr/lib/petsc/include/petscda.h>
Preprocessor definitions
#define PSIPRIME_FLOPS 5
#define PSIDOUBLEPRIME_FLOPS 8
#define RESIDUAL_FLOPS_2D
#define RESIDUAL_FLOPS_3D
This computes the Jacobian matrix at each iteration, starting with the alpha
term which is pre-computed at the beginning by
ch_jacobian_alpha_2d().
int ch_jacobian_2d ( Vec X, Mat* A, Mat* B, MatStructure* flag, void* ptr )
- int ch_jacobian_2d
- It returns 0 or an error code.
- Vec X
- The vector of unknowns.
- Mat* A
- The Jacobian matrix returned to PETSc.
- Mat* B
- The matrix preconditioner, in this case the same matrix.
- MatStructure* flag
- Flag saying the nonzeroes are the same each time.
- void* ptr
- Application context structure.
This computes the Jacobian matrix at each iteration, starting with the alpha
term which is pre-computed at the beginning by
ch_jacobian_alpha_3d().
int ch_jacobian_3d ( Vec X, Mat* A, Mat* B, MatStructure* flag, void* ptr )
- int ch_jacobian_3d
- It returns 0 or an error code.
- Vec X
- The vector of unknowns.
- Mat* A
- The Jacobian matrix returned to PETSc.
- Mat* B
- The matrix preconditioner, in this case the same matrix.
- MatStructure* flag
- Flag saying the nonzeroes are the same each time.
- void* ptr
- Application context structure.
This creates the initial alpha and J matrices, where alpha is the alpha term
component of the Jacobian. Since the alpha term is linear, this part of the
Jacobian need only be calculated once.
int ch_jacobian_alpha_2d ( AppCtx* user )
- int ch_jacobian_alpha_2d
- It returns zero or an error message.
- AppCtx* user
- The application context structure pointer.
This creates the initial alpha and J matrices, where alpha is the alpha term
component of the Jacobian. Since the alpha term is linear, this part of the
Jacobian need only be calculated once.
int ch_jacobian_alpha_3d ( AppCtx* user )
- int ch_jacobian_alpha_3d
- It returns zero or an error message.
- AppCtx* user
- The application context structure pointer.
This evaluates the nonlinear finite difference approximation to the residuals
Ri.
Note that it loops on the interior points and the boundary separately, to
avoid conditional statements within the double loop over the local grid
indices.
int ch_residual_vector_2d ( Vec X, Vec F, void* ptr )
- int ch_residual_vector_2d
- It returns zero or an error value.
- Vec X
- The pre-allocated local vector of unknowns.
- Vec F
- The pre-allocated local vector of residuals, filled by this function.
- void* ptr
- Data passed in the application context.
This evaluates the nonlinear finite difference approximation to the residuals
Ri.
Note that it loops on the interior points and the boundary separately, to
avoid conditional statements within the double loop over the local grid
indices.
Thus far, only periodic boundary conditions work.
int ch_residual_vector_3d ( Vec X, Vec F, void* ptr )
- int ch_residual_vector_3d
- It returns zero or an error value.
- Vec X
- The pre-allocated local vector of unknowns.
- Vec F
- The pre-allocated local vector of residuals, filled by this function.
- void* ptr
- Data passed in the application context.
This calculates the second derivative of homogeneous free energy with respect
to concentration, for insertion into the Jacobian matrix. See the
ch_psiprime() function for the definition of Psi.
static inline PetscScalar ch_psidoubleprime ( PetscScalar C, PetscScalar mparam )
- PetscScalar ch_psidoubleprime
- It returns the second derivative
Psi''(C).
- PetscScalar C
- The concentration.
- PetscScalar mparam
- The model parameter
m.
This calculates the derivative of homogeneous free energy with respect to
concentration, which is a component of the chemical potential. It currently
uses
Psi' = C (1-C) (0.5+m+C)
which gives (meta)stable equilibria at
C=0 and 1 and an unstable equilibrium at C = 1/2 +
m; if m>0
then the 0 phase is stable and vice versa.
static inline PetscScalar ch_psiprime ( PetscScalar C, PetscScalar mparam )
- PetscScalar ch_psiprime
- It returns the derivative itself.
- PetscScalar C
- The concentration.
- PetscScalar mparam
- The model parameter
m.
This function computes the residual from indices to points in the
concentration array. ``Up'' refers to the positive
y-direction, ``down'' to negative y, ``left'' to
negative x and ``right'' to positive x.
static inline PetscScalar ch_residual_2d ( PetscScalar* conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam, PetscScalar hx, PetscScalar hy, int upup, int upleft, int up, int upright, int leftleft, int left, int current, int right, int rightright, int downleft, int down, int downright, int downdown )
- PetscScalar ch_residual_2d
- Returns the residual itself
- PetscScalar* conc
- Array of concentrations
- PetscScalar alpha
- Model parameter
kappa alpha
- PetscScalar beta
- Model parameter
kappa beta
- PetscScalar mparam
- Model parameter
m
- PetscScalar hx
- Inverse square
x-spacing hx-2
- PetscScalar hy
- Inverse square
y-spacing hy-2
- int upup
- Index to array position two cells up from current
- int upleft
- Index to array position one cell up and one left from current
- int up
- Index to array position one cell up from current
- int upright
- Index to array position one cell up and one right from current
- int leftleft
- Index to array position two cells left of current
- int left
- Index to array position one cell left of current
- int current
- Index to current cell array position
- int right
- Index to array position one cell right of current
- int rightright
- Index to array position two cells right of current
- int downleft
- Index to array position one cell down and one left from current
- int down
- Index to array position one cell down from current
- int downright
- Index to array position one cell down and one right from current
- int downdown
- Index to array position two cells down from current
This calculates the
beta term, kappa*beta times the laplacian of Psi prime(C),
then subtracts the
alpha term, kappa*alpha times the Laplacian of the Laplacian of C.
This function computes the residual from indices to points in the
concentration array. ``Front refers to the positive
z-direction, ``back'' to negative z, ``up'' to positive
y, ``down'' to negative y, ``left'' to
negative x and ``right'' to positive x.
static inline PetscScalar ch_residual_3d ( PetscScalar* conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam, PetscScalar hx, PetscScalar hy, PetscScalar hz, int frontfront, int frontup, int frontleft, int front, int frontright, int frontdown, int upup, int upleft, int up, int upright, int leftleft, int left, int current, int right, int rightright, int downleft, int down, int downright, int downdown, int backup, int backleft, int back, int backright, int backdown, int backback )
- PetscScalar ch_residual_3d
- Returns the residual itself
- PetscScalar* conc
- Array of concentrations
- PetscScalar alpha
- Model parameter
kappa alpha
- PetscScalar beta
- Model parameter
kappa beta
- PetscScalar mparam
- Model parameter
m
- PetscScalar hx
- Inverse square
x-spacing hx-2
- PetscScalar hy
- Inverse square
y-spacing hy-2
- PetscScalar hz
- Inverse square
z-spacing hz-2
- int frontfront
- Index to array position two cells in front of current
- int frontup
- Index to array position one cell front and one up from current
- int frontleft
- Index to array position one cell front and one left from current
- int front
- Index to array position one cell in front of current
- int frontright
- Index to array position one cell front and one right from current
- int frontdown
- Index to array position one cell front and one down from current
- int upup
- Index to array position two cells up from current
- int upleft
- Index to array position one cell up and one left from current
- int up
- Index to array position one cell up from current
- int upright
- Index to array position one cell up and one right from current
- int leftleft
- Index to array position two cells left of current
- int left
- Index to array position one cell left of current
- int current
- Index to current cell array position
- int right
- Index to array position one cell right of current
- int rightright
- Index to array position two cells right of current
- int downleft
- Index to array position one cell down and one left from current
- int down
- Index to array position one cell down from current
- int downright
- Index to array position one cell down and one right from current
- int downdown
- Index to array position two cells down from current
- int backup
- Index to array position one cell back and one up from current
- int backleft
- Index to array position one cell back and one left from current
- int back
- Index to array position one cell in back of current
- int backright
- Index to array position one cell back and one right from current
- int backdown
- Index to array position one cell back and one down from current
- int backback
- Index to array position two cells in back of current
This calculates the
beta term, kappa*beta times the laplacian of Psi prime(C),
then subtracts the
alpha term, kappa*alpha times the Laplacian of the Laplacian of C.