root/trunk/rheoplast/cahnhill.c

/***************************************
  $Header$

  This provides a simple Cahn-Hilliard module for Rheoplast.  The free energy
  goes as
  +latex+\begin{equation}
  +latex+  \label{ch_energy}
  +latex+  {\cal F} = \int \left(\beta\Psi(C) + \frac{\alpha}{2}|\nabla C|^2
  +latex+  \right)dV,
  +latex+\end{equation}
  +latex+where $\Psi(C)$ is the homogeneous free energy, given here as
  +latex+$C^2 (1-C)^2$, or with the {\tt -ch\_polymer} option, a Flory-Huggins
  +latex+polymer solution free energy given by
  +latex+\begin{equation}
  +latex+  \label{ch_polymer_energy}
  +latex+  \beta\Psi(C) = \frac{C}{m}\ln C + (1-C)\ln(1-C) + \chi C(1-C).
  +latex+\end{equation}
  +html+ the integral of a gradient penalty and a homogeneous free energy, the
  +html+ latter of which is <i>C</i><sup>2</sup> (1-<i>C</i>)<sup>2</sup>, or
  +html+ in the polymer option selected by <tt>-ch_polymer</tt>, is given by:
  +html+ <center>beta Psi(<i>C</i>) = (1/<i>m</i>)<i>C</i> ln <i>C</i> +
  +html+ (1-<i>C</i>) ln (1-<i>C</i>) + chi <i>C</i>(1-<i>C</i>).</center>
  The chemical potential is then the variation of this free energy, given by
  +latex+\begin{equation}
  +latex+  \label{ch_potential}
  +latex+  \mu = \frac{\delta\cal F}{\delta C} = \beta\Psi'(C) -
  +latex+  \alpha\nabla^2 C.
  +latex+\end{equation}
  +html+ the homogeneous free energy derivative term and a laplacian term
  +html+ corresponding to the gradient penalty.
  ***************************************/

#include "rheoplast.h"


#undef __FUNCT__
#define __FUNCT__ "psiprime"

/*++++++++++++++++++++++++++++++++++++++
  This abstracts out the function for
  +latex+$\Psi'(C)$,
  +html+ Psi'(<i>C</i>),
  the derivative of homogeneous free energy, so it can be easily modified.
  +latex+Since $\Psi(C) = C^2 (1-C)^2 = C^4 - 2C^3 + C^2$, this returns
  +latex+$4C^3 - 6C^2 + 2C$.
  +html+ Since Psi(<i>C</i>) = <i>C</i><sup>2</sup> (1-<i>C</i>)<sup>2</sup> =
  +html+ <i>C</i><sup>4</sup> - 4<i>C</i><sup>3</sup> + <i>C</i><sup>2</sup>,
  +html+ this returns 4<i>C</i><sup>3</sup> - 6<i>C</i><sup>2</sup>
  +html+ + 2<i>C</i>.
  +latex+\par
  +html+ <p>
  There is also a polymer thermo option based on Flory-Huggins thermodynamics
  +latex+given in equation \ref{ch_polymer_energy}.
  +html+ described above.
  Enable it using option
  +latex+{\tt -ch\_polymer}
  +html+ <tt>-ch_polymer</tt>
  and control it using options
  +latex+{\tt -ch\_polymer\_chi} and {\tt -ch\_polymer\_m}
  +html+ <tt>-ch_polymer_chi</tt> and <tt>-ch_polymer_m</tt>
  (defaults: 0.58, 640).  Note that with
  +latex+$m=1$
  +html+ <i>m</i>=1
  this reduces to the regular solution model.

  PetscScalar psiprime It returns the derivative of homogeneous free
  energy.

  PetscScalar C The
  +latex+$C$
  +html+ <i>C</i>
  parameter it's a function of.

  chparm *thecahnhill Cahn-Hilliard parameter structure.
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar psiprime (PetscScalar C, chparm *thecahnhill)
{
  PetscTruth polymer_solution;
  int ierr;
  ierr = PetscOptionsHasName (PETSC_NULL, "-ch_polymer", &polymer_solution);
  CHKERRQ (ierr);
  
  if (polymer_solution)
    {
      return 1.0/thecahnhill->m*log(C)-log(1.0-C)+thecahnhill->chi*(1.0-2.0*C)+1.0/thecahnhill->m-1.0;
    }

  else
    return 4.*C*C*C - 6.*C*C + 2.*C;
}


#undef __FUNCT__
#define __FUNCT__ "psidoubleprime"

/*++++++++++++++++++++++++++++++++++++++
  This abstracts out the function for
  +latex+$\Psi''(C)$,
  +html+ Psi''(<i>C</i>),
  the second derivative of homogeneous free energy, so it can be easily
  modified.
  +latex+Since $\Psi(C) = C^2 (1-C)^2 = C^4 - 2C^3 + C^2$, this returns
  +latex+$12C^2 - 12C + 2$.
  +html+ Since Psi(<i>C</i>) = <i>C</i><sup>2</sup> (1-<i>C</i>)<sup>2</sup> =
  +html+ <i>C</i><sup>4</sup> - 4<i>C</i><sup>3</sup> + <i>C</i><sup>2</sup>,
  +html+ this returns 12<i>C</i><sup>2</sup> - 12<i>C</i> + 2.
  +latex+\par
  +html+ <p>
  There is also a polymer thermo option based on Flory-Huggins thermodynamics
  +latex+given in equation \ref{ch_polymer_energy}.
  +html+ described above.
  Enable it using option
  +latex+{\tt -ch\_polymer}
  +html+ <tt>-ch_polymer</tt>
  and control it using options
  +latex+{\tt -ch\_polymer\_chi} and {\tt -ch\_polymer\_m}
  +html+ <tt>-ch_polymer_chi</tt> and <tt>-ch_polymer_m</tt>
  (defaults: 0.58, 640).  Note that with
  +latex+$m=1$
  +html+ <i>m</i>=1
  this reduces to the regular solution model.

  PetscScalar psidoubleprime It returns the second derivative of homogeneous
  free energy.

  PetscScalar C The
  +latex+$C$
  +html+ <i>C</i>
  parameter it's a function of.

  chparm *thecahnhill Cahn-Hilliard parameter structure.
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar psidoubleprime (PetscScalar C, chparm *thecahnhill)
{
  PetscTruth polymer_solution;
  int ierr;
  ierr = PetscOptionsHasName (PETSC_NULL, "-ch_polymer", &polymer_solution);
  CHKERRQ (ierr);
  
  if (polymer_solution)
    return 1.0/thecahnhill->m/C + 1./(1.0-C) - 2.*thecahnhill->chi;
  else
    return 12.*C*C - 12.*C + 2;
}


#undef __FUNCT__
#define __FUNCT__ "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
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  struct typedef.

  int cahnhill_first_setup Returns zero (or an error code).

  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
  +latex+{\tt AppCtx}
  +html+ <tt>AppCtx</tt>
  struct typedef, into whose
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  structure this inserts parameters from the command line.
  ++++++++++++++++++++++++++++++++++++++*/

int cahnhill_first_setup
(PetscTruth threedee, int *vars, int *tempvars, int *stencilwid, AppCtx *data)
{
  chparm *thecahnhill = &data->thecahnhill;
  int ierr;

  thecahnhill->Cvar  = (*vars)++;
  thecahnhill->muvar = (*tempvars)++;

  /*+ The standard model temporarily sets
    +latex+$\alpha$ and $\beta$ to $\epsilon/\delta x$ and $\sigma$,
    +html+ <i>alpha</i> and <i>beta</i> to <i>epsilon</i>/<i>dx</i> and
    +html+ <i>sigma</i>,
    and mobility to the dummy value -1, so that if not overridden by the user
    parameter setting, its default value of
    +latex+$\epsilon^2$ can be set in {\tt cahnhill\_labels\_initcond()}.
    +html+ <i>epsilon</i><sup>2</sup> can be set in
    +html+ <tt>cahnhill_labels_initcond()</tt>.
    +*/
  ierr = PetscOptionsHasName (PETSC_NULL, "-ch_polymer",
                              &thecahnhill->polymer_solution); CHKERRQ (ierr);
  if (thecahnhill->polymer_solution)
    {
      thecahnhill->chi=0.58;
      thecahnhill->m=640;
      ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_polymer_chi",
                                    &thecahnhill->chi, PETSC_NULL);
      CHKERRQ (ierr);
      ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_polymer_m",
                                    &thecahnhill->m, PETSC_NULL);
      CHKERRQ (ierr);
    }

  *stencilwid = PetscMax (*stencilwid, 2);

  return 0;
}


#undef __FUNCT__
#define __FUNCT__ "cahnhill_labels_initcond"

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

  int cahnhill_labels_initcond Returns zero (or an error code).

  PetscScalar *globalarray The global field array.

  int nx Overall
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the global array.

  int ny Overall
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the global array.

  int nz Overall
  +latex+$z$-width
  +html+ <i>z</i>-width
  of the global array.

  int xm The
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the local part of the array.

  int ym The
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the local part of the array.

  int zm The
  +latex+$z$-width
  +html+ <i>z</i>-width
  of the local part of the array.

  int xs The (integer)
  +latex+$x$-coordinate
  +html+ <i>x</i>-coordinate
  of the start of the local part of the array.

  int ys The (integer)
  +latex+$y$-coordinate
  +html+ <i>y</i>-coordinate
  of the start of the local part of the array.

  int zs The (integer)
  +latex+$z$-coordinate
  +html+ <i>z</i>-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
  +latex+{\tt AppCtx}
  +html+ <tt>AppCtx</tt>
  struct typedef, whose
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  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.
  ++++++++++++++++++++++++++++++++++++++*/

int 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)
{
  chparm *thecahnhill = &data->thecahnhill;
  int ierr;
  PetscScalar surftens, intwidth, amplitude,cal_Diffusivity,t_scale,
    delta = PetscMin (data->xwid/nx, data->ywid/ny);

  PetscScalar voltage, rhoM;

  if (data->threedee)
    delta = PetscMin (delta, data->zwid/nz);

  /*+ 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
    +latex+{\tt -ch\_intwidth} and {\tt -ch\_surftens}
    +html+ <tt>-ch_intwidth</tt> and <tt>-ch_surftens</tt>
    (intwidth is multiplied by
    +latex+$\Delta x$),
    +html+ <i>dx</i>),
    or by setting the homogeneous and gradient penalty coefficients themselves
    using
    +latex+{\tt -ch\_alpha} and {\tt -ch\_beta}.
    +html+ <tt>-ch_alpha</tt> and <tt>-ch_beta</tt>.
    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).
    +*/
  thecahnhill->intwidth = 1.;
  thecahnhill->surftens = 1.;
  amplitude=1.;
  voltage=1.;
  rhoM=100000.0;

  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_intwidth", &thecahnhill->intwidth,
                                PETSC_NULL); CHKERRQ (ierr);
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_surftens", &thecahnhill->surftens,
                                PETSC_NULL); CHKERRQ (ierr);
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_amplitude", &amplitude,
                                PETSC_NULL); CHKERRQ (ierr);
  amplitude *= delta;
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-voltage", &voltage,
                                PETSC_NULL); CHKERRQ (ierr);
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-rhoM", &rhoM,
                              PETSC_NULL); CHKERRQ (ierr);

  thecahnhill->intwidth *= delta;         /* True thickness */
  thecahnhill->alpha = sqrt(18.)*thecahnhill->surftens * (thecahnhill->intwidth/3.1);
  thecahnhill->beta = sqrt(18.)*thecahnhill->surftens / (thecahnhill->intwidth/3.1);
  DPRINTF ("Interface thickness = %g, tension = %g\n", thecahnhill->intwidth, thecahnhill->surftens);
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_alpha", &thecahnhill->alpha,
                                PETSC_NULL); CHKERRQ (ierr);
  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_beta", &thecahnhill->beta,
                                PETSC_NULL); CHKERRQ (ierr);
  DPRINTF ("Final alpha = %g, beta = %g , amplitude/delta = %g\n",thecahnhill->alpha,
           thecahnhill->beta, amplitude/delta);

  /*+By default,
    +latex+$\kappa$ is set to $\epsilon^2/\beta$,
    +html+ <i>kappa</i> is set to <i>epsilon</i><sup>2</sup>/<i>beta</i>,
    such that the diffusion timescale over the interface thickness is constant;
    this is controlled by command line option
    +latex+{\tt -ch\_mobility}.
    +html+ <tt>-ch_mobility</tt>.
    +*/
  thecahnhill->mobility = thecahnhill->intwidth*thecahnhill->intwidth;
  if (data->electra)
    {
      thecahnhill->mobility =0.5; /*1/<psi''> ~ 1/2*/
      thecahnhill->Pe=96484.*rhoM*voltage/thecahnhill->beta;
      DPRINTF("Peclet Number = %g\n",thecahnhill->Pe);
    }
  cal_Diffusivity=thecahnhill->mobility*thecahnhill->beta*1.;

  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_mobility",
                                &thecahnhill->mobility, PETSC_NULL); 
    
  CHKERRQ (ierr);

  /* Set DA labels, constraints and symmetry types */
  data->label[thecahnhill->Cvar] = "Cahn-Hilliard-C";
  data->timestyle[thecahnhill->Cvar] = TIMESTEP_ONLY_VAR+thecahnhill->Cvar;
  data->plot_types[thecahnhill->Cvar] = FIELD_SCALAR;
  data->symmtypes[thecahnhill->Cvar] =
    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMIN_START |
    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMIN_START |
    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMIN_START |
    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMAX_START |
    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMAX_START |
    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMAX_START;

  /*+Default polymer Flory-Huggins parameters are
    +latex+$\chi=0.58$ and $m=640$,
    +html+ <i>chi</i>=0.58, <i>m</i>=640,
    as discussed in
    +latex+appendix \ref{func_psiprime_cahnhill.c}, page
    +latex+\pageref{func_psiprime_cahnhill.c}.
    +html+ the <a href="#func-psiprime"><tt>psiprime()</tt></a> function in
    +html+ this file. 
    +*/
  ierr = PetscOptionsHasName (PETSC_NULL, "-ch_polymer",
                              &thecahnhill->polymer_solution); CHKERRQ (ierr);
  if (thecahnhill->polymer_solution)
    {
      thecahnhill->chi=0.58;
      thecahnhill->m=640.;
      ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_polymer_chi",
                                    &thecahnhill->chi, PETSC_NULL);
      CHKERRQ (ierr);
      ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_polymer_m",
                                    &thecahnhill->m, PETSC_NULL);
      CHKERRQ (ierr);
    }

  /*+ The explicit finite difference timestep size used here is
    +latex+$(\Delta x)^3/40\kappa$,
    +html+ delta x/40 kappa,
    which works for the fourth-order polynomial free energy in 2-D when nx=ny,
    +latex+$\epsilon=\Delta x$ ({\tt intwidth}=1), and $\sigma=1$.
    +html+ <i>epsilon</i>=<i>dx</i> ({\tt intwidth}=1), and <i>sigma</i>=1.
    +*/
  *max_explicit_deltat = PetscMin
    (*max_explicit_deltat, (data->threedee ? .1/6 : .025)/sqrt(18) *
     delta*delta*delta / thecahnhill->mobility);
  DPRINTF ("Kappa = %g, D= %g, calculated deltat = %g\n", thecahnhill->mobility,cal_Diffusivity, *max_explicit_deltat);

  if (!data->load_data) /* Make sure we're not stepping on loaded data */
    {
      PetscTruth layers, trilayer, particles, sinusoidal, randoms;
      int xmin=0,xmax=0,ymin=0,ymax=0,zmin=0,zmax=1, ix,iy,iz;
      PetscScalar fluct=0.;

      DPRINTF ("Setting Cahn-Hilliard initial condition\n",0);

      /* Uniform empty starting point */
      for (iz=0; iz<(data->threedee ? zm: 1); iz++)
        for (iy=0; iy<ym; iy++)
          for (ix=0; ix<xm; ix++)
            globalarray [((iz*ym + iy)*xm + ix)*vars + thecahnhill->Cvar] =
              ((thecahnhill->polymer_solution) ? 0.002: 0.);
         
      ierr = PetscOptionsHasName (PETSC_NULL, "-ch_layers", &layers);
      CHKERRQ (ierr);
      ierr = PetscOptionsHasName (PETSC_NULL, "-ch_trilayer", &trilayer);
      CHKERRQ (ierr);
      ierr = PetscOptionsHasName (PETSC_NULL, "-ch_particles", &particles);
      CHKERRQ (ierr);
      ierr = PetscOptionsHasName (PETSC_NULL, "-ch_sincathode", &sinusoidal);
      CHKERRQ (ierr);
      ierr = PetscOptionsHasName (PETSC_NULL, "-ch_random_center", &randoms);
      CHKERRQ (ierr);

      if (layers || trilayer || particles || randoms)
        {
          ierr = PetscPrintf (PETSC_COMM_WORLD,
                              "The -ch_layers, -ch_trilayer, -ch_particles and -ch_random* options are gone.\nPlease see the documentation on initcond.c for general initial conditions.\n");
          CHKERRQ (ierr);
        }
      else if (sinusoidal)
        {
          PetscScalar cathode;
          ierr = PetscOptionsGetScalar (PETSC_NULL, "-ch_sincathode",
                                        &cathode, PETSC_NULL);

          DPRINTF ("Sinusoidal cathode with mean %g, amplitude %d\n", cathode, amplitude);

          /*Sinusoidal Perturbation */
          xmin=0;
          xmax=nx;
          zmin=0;
          zmax=nz;            
          ymin=0;

          for (iz = (data->threedee ? PetscMax (zs, zmin) : 0);
               iz < (data->threedee ? PetscMin (zs+zm, zmax): 1); iz++)
            for (ix=PetscMax(xs,xmin); ix<PetscMin(xs+xm,xmax); ix++)
              {
                PetscScalar ymax_scalar = cathode*ny + (amplitude/delta)*
                  (sin(2.* PETSC_PI* ix/nx) + sin(2.* PETSC_PI* iz/nz));
                for (iy=PetscMax(ys,ymin); iy<PetscMin(ys+ym,(int)ymax_scalar);
                     iy++)
                  globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
                               thecahnhill->Cvar] = 1.;
                if (iy+1>ymax_scalar && iy<ymax_scalar && iy<ys+ym)
                  {
                    globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
                                 thecahnhill->Cvar] = ymax_scalar - iy;
                  }
              }
        }
    }

  return 0;
}


#define C(point) (x [vars*(point) + thecahnhill->Cvar])
#define Cfunc(point) (func [vars*(point) + thecahnhill->Cvar])
#define mu(point) (temp [tempvars*(point) + thecahnhill->muvar])

#define vu(point) (x [vars*(point) + thevortex->uvar])
#define vv(point) (x [vars*(point) + thevortex->vvar])
#define pu(point) (x [vars*(point) + thepressure->uvar])
#define pv(point) (x [vars*(point) + thepressure->vvar])
#define pw(point) (x [vars*(point) + thepressure->wvar])

#define V(point) (x [vars*(point) + thepotential->Vvar])
#define sigma(point) (temp [tempvars*(point) + thepotential->sigmavar])
#define sigmaprime(point) (temp [tempvars*(point) + thepotential->sigmaprimevar])
#ifdef VIEW_SIGMA
#define sigma_eff(point) (x [vars*(point) + thepotential->sigma_effvar])
#define sigma_efffunc(point) (func [vars*(point) + thepotential->sigma_effvar])
#else
#define sigma_eff(point) (temp [tempvars*(point) + thepotential->sigma_effvar])
#endif


#undef __FUNCT__
#define __FUNCT__ "cahnhill_temp_parameters_line"

/*++++++++++++++++++++++++++++++++++++++
  This calculates the Cahn-Hilliard temporary parameter
  +latex+$\mu$,
  +html+ <i>mu</i>,
  which is the chemical potential given by
  +latex+equation \ref{ch_potential}.
  -latex-the equation in the intro to this file.

  int cahnhill_temp_parameters_line Returns zero (or an error code).

  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
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$y$-increment.
  +html+ <i>y</i>-increment.

  int gym The
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$z$-increment.
  +html+ <i>z</i>-increment.

  PetscScalar xmin First node
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar xmax Last node plus one
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar ycoord This line
  +latex+$y$-coordinate.
  +html+ <i>y</i>-coordinate.

  PetscScalar zcoord This line
  +latex+$z$-coordinate
  +html+ <i>z</i>-coordinate.

  PetscScalar time Current simulation time.

  AppCtx *data Pointer to the main simulation parameter structure, which
  includes the
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  struct typedef, from which this gets needed parameters.
  ++++++++++++++++++++++++++++++++++++++*/

int 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)
{
  chparm *thecahnhill = &data->thecahnhill;
  PetscScalar deltax_m2 = data->deltax_m2, deltay_m2 = data->deltay_m2,
    deltaz_m2 = data->deltaz_m2;
  int i, vars = data->vars, tempvars = data->tempvars;

  for (i=0; i<points; i++)
    mu(i) = ( thecahnhill->beta * psiprime (C(i), thecahnhill) -
      thecahnhill->alpha * (deltax_m2 * (C(i-1) - 2.*C(i) + C(i+1)) +
                            deltay_m2 * (C(i-gxm) - 2.*C(i) + C(i+gxm))) )/thecahnhill->beta;

  if (data->threedee)
    for (i=0; i<points; i++)
      mu(i) -= ( thecahnhill->alpha * deltaz_m2 *
        (C(i-gxm*gym) - 2.*C(i) + C(i+gxm*gym)) )/thecahnhill->beta;

  return 0;
}


#undef __FUNCT__
#define __FUNCT__ "cahnhill_temp_parameters_boundary_line"

/*++++++++++++++++++++++++++++++++++++++
  This does the same thing as cahnhill_temp_parameters_line but on the
  boundary.

  int cahnhill_temp_parameters_boundary_line Returns zero (or an error code).

  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
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$y$-increment.
  +html+ <i>y</i>-increment.

  int gym The
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$z$-increment.
  +html+ <i>z</i>-increment.

  PetscScalar xmin First node
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar xmax Last node plus one
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar ycoord This line
  +latex+$y$-coordinate.
  +html+ <i>y</i>-coordinate.

  PetscScalar zcoord This line
  +latex+$z$-coordinate
  +html+ <i>z</i>-coordinate.

  PetscScalar time Current simulation time.

  AppCtx *data Pointer to the main simulation parameter structure, which
  includes the
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  struct typedef, from which this gets needed parameters.

  int side Side on which to calculate the temporary fields.
  ++++++++++++++++++++++++++++++++++++++*/

int 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)
{
  chparm *thecahnhill = &data->thecahnhill;
  PetscScalar deltax_m2 = data->deltax_m2, deltay_m2 = data->deltay_m2,
    deltaz_m2 = data->deltaz_m2;
  int i, vars = data->vars, tempvars = data->tempvars;

  if (side==YMIN_BOUNDARY)
    {
      /*Replace C(i-gxm) with C(i+gxm)*/
      for (i=0; i<points; i++)
        
        mu(i) = ( thecahnhill->beta * psiprime (C(i), thecahnhill) -
                  thecahnhill->alpha * (deltax_m2 * (C(i-1) - 2.*C(i) + C(i+1)) +
                                        deltay_m2 * (C(i+gxm) - 2.*C(i) + C(i+gxm))) )/thecahnhill->beta;
      
      if (data->threedee)
        for (i=0; i<points; i++)
          mu(i) -= ( thecahnhill->alpha * deltaz_m2 *
                     (C(i-gxm*gym) - 2.*C(i) + C(i+gxm*gym)) )/thecahnhill->beta;
    }
  if (side==YMAX_BOUNDARY)
    {
      /*Replace C(i+gxm) with C(i-gxm)*/
      for (i=0; i<points; i++)
        {
          mu(i) = ( thecahnhill->beta * psiprime (C(i), thecahnhill) -
                    thecahnhill->alpha * (deltax_m2 * (C(i-1) - 2.*C(i) + C(i+1)) +
                                          deltay_m2 * (C(i-gxm) - 2.*C(i) + C(i-gxm))) )/thecahnhill->beta;
          
          if (data->threedee)
            for (i=0; i<points; i++)
              mu(i) -= ( thecahnhill->alpha * deltaz_m2 *
                         (C(i-gxm*gym) - 2.*C(i) + C(i+gxm*gym)) )/thecahnhill->beta;
        }
    }

  return 0;
}


#undef __FUNCT__
#define __FUNCT__ "cahnhill_interior_line_function"

/*++++++++++++++++++++++++++++++++++++++
  This calculates the time derivative of
  +latex+$C$
  +html+ <i>C</i>
  using the divergence of flux, which goes down the gradient of
  chemical potential
  +latex+given by equation \ref{ch_potential}.  That time derivative is thus
  +latex+given by
  +latex+\begin{equation}
  +latex+  \label{ch_dynamics}
  +latex+  \frac{\partial C}{\partial t} = \nabla\cdot(\kappa\nabla\mu).
  +latex+\end{equation}
  +latex+where $\kappa$ is the mobility.
  +html+ described above, with mobility <i>kappa</i>.

  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
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$y$-increment.
  +html+ <i>y</i>-increment.

  int gym The
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$z$-increment.
  +html+ <i>z</i>-increment.

  PetscScalar xmin First node
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar xmax Last node plus one
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar ycoord This line
  +latex+$y$-coordinate.
  +html+ <i>y</i>-coordinate.

  PetscScalar zcoord This line
  +latex+$z$-coordinate.
  +html+ <i>z</i>-coordinate.

  PetscScalar time Current simulation time.

  AppCtx *data Pointer to the main simulation parameter structure, which
  includes the
  +latex+{\tt chparm}
  +html+ <tt>chparm</tt>
  struct typedef, from which this gets needed parameters.
  ++++++++++++++++++++++++++++++++++++++*/

int 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)
{
  chparm *thecahnhill = &data->thecahnhill;
  PetscScalar deltax_m2 = data->deltax_m2, deltay_m2 = data->deltay_m2,
    deltaz_m2 = data->deltaz_m2, deltax_m1 = sqrt (deltax_m2),
    deltay_m1 = sqrt (deltay_m2), deltaz_m1 = sqrt (deltaz_m2);
  int i, vars = data->vars, tempvars = data->tempvars;

  PetscScalar echem_coeff;int ierr;
 
  for (i=0; i<points; i++)
    {
      /*+For now this assumes uniform 
        +latex+$\kappa$.
        -latex-kappa. +*/
      Cfunc(i) = thecahnhill->mobility *
        (deltax_m2 * (mu(i-1) - 2.*mu(i) + mu(i+1)) +
         deltay_m2 * (mu(i-gxm) - 2.*mu(i) + mu(i+gxm)));
    }

  if (data->threedee)
    for (i=0; i<points; i++)
      Cfunc(i) += thecahnhill->mobility * deltaz_m2 *
        (mu(i-gxm*gym) - 2.*mu(i) + mu(i+gxm*gym));

  if (data->electra)
    {
      echemparm *thepotential = &data->thepotential;
   
      for (i=0; i<points; i++)
        {
                          
          Cfunc(i) += thecahnhill->Pe * ( sigma(i)*( deltax_m2*(V(i-1) -2.*V(i) +V(i+1)) + deltay_m2*(V(i-gxm) -2.*V(i) +V(i+gxm)) )+ 
                                      0.5*deltax_m2*(sigma(i+1)*(V(i+1)-V(i)) -sigma(i-1)*(V(i)-V(i-1)))+
                                      0.5*deltay_m2*(sigma(i+gxm)*(V(i+gxm)-V(i)) -sigma(i-gxm)*(V(i)-V(i-gxm))) ); 
        } 
      if(data->threedee)
          for (i=0; i<points; i++)
            {
              Cfunc(i)+=  thecahnhill->Pe * ( sigma(i)*( deltaz_m2*(V(i-gxm*gym) -2.*V(i) +V(i+gxm*gym)) )+                           
                                         0.5*deltaz_m2*(sigma(i+gxm*gym)*(V(i+gxm*gym)-V(i)) -sigma(i-gxm*gym)*(V(i)-V(i-gxm*gym))) ); 
            } 
    }
   
  
  /*+It includes a convective term with first-order upwinding for
    velocity-vorticity flow.+*/
  if (data->vortflow)
    {
      vortparm *thevortex = &data->thevortex;
      PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
      
      for (i=0; i<points; i++)
        {
          /* convective coefficients for upwinding */
          vx = vu(i); avx = PetscAbsScalar(vx);
          vxp = .5*(vx+avx); vxm = .5*(vx-avx);
          vy = vv(i); avy = PetscAbsScalar(vy);
          vyp = .5*(vy+avy); vym = .5*(vy-avy);
          
          Cfunc(i) -= deltax_m1 * ( vxp*(C(i) - C(i-1)) + vxm*(C(i+1) - C(i)) )
            + deltay_m1 * ( vyp*(C(i) - C(i-gxm)) + vym*(C(i+gxm) - C(i)) );
        }
    }

  /*+And a convective term with first-order upwinding for velocity-pressure
    flow.+*/
  if (data->pressflow)
    {
      pressparm *thepressure = &data->thepressure;
      
      for (i=0; i<points; i++)
        Cfunc(i) -= deltax_m1 *
          ((pu(i)+pu(i+1)>0) ? pu(i) * (C(i)-C(i-1)) : pu(i+1) * (C(i+1)-C(i)))
          + deltay_m1 *
          ((pv(i)+pv(i+gxm)>0) ? pv(i)*(C(i)-C(i-gxm)): pv(i+gxm)*(C(i+gxm)-C(i)));
      if (data->threedee)
        for (i=0; i<points; i++)
          Cfunc(i) -= deltaz_m1 *
              ((pw(i)+pw(i+gxm*gym)>0) ? pw(i)*(C(i)-C(i-gxm*gym)) :
              pw(i+gxm*gym)*(C(i+gxm*gym)-C(i)));
    }

  return 0;
}


#undef __FUNCT__
#define __FUNCT__ "cahnhill_boundary_line_function"

/*++++++++++++++++++++++++++++++++++++++
  This calculates the time derivatives and constraint functions for the
  Cahn-Hilliard equations for a boundary line.
  Note the
  +latex+{\tt shearstrain}
  +html+ <tt>shearstrain</tt>
  module interaction programmed here;

  int cahnhill_interior_line_function Returns zero (or an error code).

  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
  +latex+$x$-width
  +html+ <i>x</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$y$-increment.
  +html+ <i>y</i>-increment.

  int gym The
  +latex+$y$-width
  +html+ <i>y</i>-width
  of the ``local'' vector's array, including shadow nodes, for the
  +latex+$z$-increment.
  +html+ <i>z</i>-increment.

  PetscScalar xmin First node
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar xmax Last node plus one
  +latex+$x$-coordinate.
  +html+ <i>x</i>-coordinate.

  PetscScalar ycoord This line
  +latex+$y$-coordinate.
  +html+ <i>y</i>-coordinate.

  PetscScalar zcoord This line
  +latex+$z$-coordinate.
  +html+ <i>z</i>-coordinate.

  PetscScalar time Current simulation time.

  AppCtx *data Pointer to the main simulation parameter structure, which
  includes the
  +latex+{\tt vortparm}
  +html+ <tt>vortparm</tt>
  struct typedef, from which this gets needed parameters.

  int side Side on which to calculate the function values.
  ++++++++++++++++++++++++++++++++++++++*/

int 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)
{
  chparm *thecahnhill = &data->thecahnhill;
  PetscScalar deltax_m2 = data->deltax_m2, deltay_m2 = data->deltay_m2,
    deltaz_m2 = data->deltaz_m2, deltax_m1 = sqrt (deltax_m2),
    deltay_m1 = sqrt (deltay_m2), deltaz_m1 = sqrt (deltaz_m2);
  int i, vars = data->vars, tempvars = data->tempvars;
  
  if (side==YMIN_BOUNDARY)
  {
    /* At ymin, (i-gxm) is replaced by (i+gxm). */
    for (i=0; i<points; i++)
      {
        Cfunc(i) = thecahnhill->mobility *
                (deltax_m2 * (mu(i-1) - 2.*mu(i) + mu(i+1)) +
                 deltay_m2 * (mu(i+gxm) - 2.*mu(i) + mu(i+gxm)));
      }
    if (data->threedee)
      for (i=0; i<points; i++)
        Cfunc(i) += thecahnhill->mobility * deltaz_m2 *
          (mu(i+gxm*gym) - 2.*mu(i) + mu(i-gxm*gym));

    if (data->electra)//replace (i-gxm) with i for first derivative terms
      {
        echemparm *thepotential = &data->thepotential;

        for