root/trunk/matml/src/ternary/spinodal.c

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

  The spinodal calculation function.
  ***************************************/


#include "freenergy.h"
#include "gibbs.h"
#include <stdlib.h>  /*+ For malloc() and free() +*/


/*++++++++++++++++++++++++++++++++++++++
  This function uses dissection and secant iteration to estimate the root of
  the spinodal function on a line in ternary space between points A and B.

  inline double interpolate Returns the fraction B whose weighted average with
  A minimizes the spinodal function.

  double A2 C2 variable at point A.

  double A3 C3 variable at point A.

  double AS Spinodal function value at point A.

  double A2 C2 variable at point B.

  double A3 C3 variable at point B.

  double AS Spinodal function value at point B.

  double T Temperature.

  double P Pressure (ignored for now).

  energy_params *eparams Energy parameter functions (all of them).

  int efunc Index indicating which energy function to calculate the spinodal
  for.

  int iters Number of secant iterations to perform.
++++++++++++++++++++++++++++++++++++++*/

static inline double interpolate
(double A2, double A3, double AS, double B2, double B3, double BS,
 double T, double P, energy_params *eparams, int efunc, int iters)
{
  int i, binaries=10;
  energy_point X;
  double interp=0., intfrac=1., XS;

  for (i=0; i<binaries; i++)
    {
      X.C2 = 0.5 * (A2+B2);
      X.C3 = 0.5 * (A3+B3);
      X.efunc = -1;
      free_energies (&X, 1, T,P, eparams,efunc, WITH_DERIVATIVES_NO_INFINITY);
      XS = X.G22 * X.G33 - X.G23 * X.G23;

      if (AS * XS > 0)
        {
          interp += intfrac * 0.5;
          intfrac *= 0.5;
          A2 = X.C2;
          A3 = X.C3;
          AS = XS;
        }
      else
        {
          intfrac *= 0.5;
          B2 = X.C2;
          B3 = X.C3;
          BS = XS;
        }
    }

  for (i=0; i<iters; i++)
    {
      X.C2 = (AS*B2 - BS*A2) / (AS - BS);
      X.C3 = (AS*B3 - BS*A3) / (AS - BS);
      X.efunc = -1;
      free_energies (&X, 1, T,P, eparams,efunc, WITH_DERIVATIVES_NO_INFINITY);
      XS = X.G22 * X.G33 - X.G23 * X.G23;

      if (AS * XS > 0)
        {
          interp += intfrac * AS / (AS - BS);
          intfrac *= -BS / (AS-BS);
          A2 = X.C2;
          A3 = X.C3;
          AS = XS;
        }
      else
        {
          intfrac *= AS / (AS-BS);
          B2 = X.C2;
          B3 = X.C3;
          BS = XS;
        }
    }

  return interp + intfrac * AS / (AS - BS);
}


/*++++++++++++++++++++++++++++++++++++++
  This calculates and returns all of the spinodal curves for a set of points
  and triangles.  The spinodal criterion in a ternary system is:
  +latex+\begin{equation}
  +latex+  \label{eq:ternaryspinodal}
  +latex+  G_{22}G_{33} - G_{23}^2 = 0.
  +latex+\end{equation}
  +latex+where $G_{ij}$ is the second derivative of free energy with respect
  +latex+to the $i$ and $j$ concentration variables.
  +html+ <center><i>G</i><sub>22</sub> <i>G</i><sub>33</sub> -
  +html+ <i>G</i><sub>23</sub><sup>2</sup>,</center>
  +html+ where <i>G<sub>ij</sub></i> is the second derivative of free energy
  +html+ with respect to the <i>i</i> and <i>j</i> concentration variables.
  This calculates the criterion at each vertex of every triangle, and if it
  crosses zero on a line segment in the triangle, adds any new end points to
  the points list, and the line to the returned phase boundary structure.
  +latex+\par
  +html+ <p>
  This is second order accurate in the initial resolution, which might not be
  sufficient for some applications.  To improve the resolution, one can call
  this with
  +latex+{\tt spinreturn} set to {\tt NULL},
  +html+ <tt>spinreturn</tt> set to <tt>NULL</tt>,
  so it generates new points close to the spinodal, any number of times, then
  call it with a pointer to a valid
  +latex+{\tt phase\_boundary}
  +html+ <tt>phase_boundary</tt>
  structure to get a "final" spinodal estimate.  This is analogous to secant
  method iteration.

  int calc_spinodal It returns zero or an error code.

  energy_point **points Points in the array, this function adds points to the
  end of the array.

  int *n_points Pointer to the number of points in the points array, which this
  function modifies as it adds new points.

  int *triangle_indices 3*n_triangles integers indicating the points which make
  up each triangle.

  int n_triangles Number of triangles.

  double T Temperature.

  double P Pressure (ignored for now).

  energy_params *eparams Free energy function parameters (all of them).

  int efunc Index indicating which set of energy parameters to use.

  phase_boundary *spinreturn Pointer to phase boundary structure which contains
  the spinodal on return, NULL to simply add points close to the spinodal.
  ++++++++++++++++++++++++++++++++++++++*/

int calc_spinodal
(energy_point **points, int *n_points, int *triangle_indices, int n_triangles,
 double T, double P, energy_params *eparams, int efunc,
 phase_boundary *spinreturn)
{
  int tri, numsegs=0, thesegs_size=20, *thesegs, numnewpoints=0,
    thenewpoints_size=20;
  energy_point *thenewpoints;

  //printf ("Calculating free energies at all points\n");
  //if (tri = free_energies (*points, *n_points, T, P, eparams, efunc,
  //                       WITH_DERIVATIVES_NO_INFINITY))
  //  return tri;

  // Initial allocation of new points and segments
#ifdef DEBUG
  printf ("Allocating new points and segments arrays\n");
#endif
  if (!(thenewpoints = (energy_point *) malloc
        (thenewpoints_size * sizeof (energy_point))))
    return 1;
  if (!(thesegs = (int *) malloc (2 * thesegs_size * sizeof (int))))
    { free (thenewpoints); return 1; }

  // Main loop over triangles to calculate spinodal intersections with each
#ifdef DEBUG
  printf ("Entering main spinodal loop, %d points, first triangle points %d, %d, %d\n", *n_points, triangle_indices[0], triangle_indices[1], triangle_indices[2]);
#endif
  for (tri=0; tri<n_triangles; tri++)
    {
      double spin[3], A2=-1, A3, B2=-1, B3, interp;
      int total;
      energy_point P0=(*points)[triangle_indices[3*tri]],
        P1=(*points)[triangle_indices[3*tri+1]],
        P2=(*points)[triangle_indices[3*tri+2]];

      if (P0.efunc == efunc)
        {
          // The spinodal criterion
          spin[0] = P0.G22 * P0.G33 - P0.G23 * P0.G23;
          spin[1] = P1.G22 * P1.G33 - P1.G23 * P1.G23;
          spin[2] = P2.G22 * P2.G33 - P2.G23 * P2.G23;

          /*+ For each triangle, this sets a "total" score according to the
            signs of the spinodal criterion values at the corners. +*/
          total = ((spin[0] > 0.) ? 2 : ((spin[0] == 0.) ? 1 : 0));
          total += ((spin[1] > 0.) ? 8 : ((spin[1] == 0.) ? 4 : 0));
          total += ((spin[2] > 0.) ? 32 : ((spin[2] == 0.) ? 16 : 0));

#ifdef DEBUG
          printf ("Tri %d (%d,%d,%d): spins %g,%g,%g; total %d\n",
                  tri, triangle_indices[3*tri], triangle_indices[3*tri+1],
                  triangle_indices[3*tri+2], spin[0],spin[1],spin[2], total);
#endif
        }
      else
        total = 0;

      // First decide whether to allocate a new segment
      switch (total)
        {
        case 0:  // neg neg neg
        case 21: //  0   0   0  -- degenerate, probably flat
        case 42: // pos pos pos
        case 1:  //  0  neg neg
        case 4:  // neg  0  neg
        case 16: // neg neg  0
        case 41: //  0  pos pos
        case 38: // pos  0  pos
        case 26: // pos pos  0
        case 5:  // 0 0 neg
        case 37: // 0 0 pos
        case 17: // 0 neg 0
        case 25: // 0 pos 0
        case 20: // neg 0 0
        case 22: // pos 0 0
          {
            A2 = B2 = -1; // No new points for these cases
            break; // No new seg for these cases either
          }
        default:
          {
#ifdef DEBUG
            printf ("Making a new segment...");
#endif
            if (++numsegs >= thesegs_size)
              if (!(thesegs = (int *) realloc
                    (thesegs, 2*(thesegs_size*=2)*sizeof (int))))
                { free (thenewpoints); free (thesegs); return 1; }
          }
        }

      switch (total)
        {
          // New Points between 0 and 1, 0 and 2
        case 2:  // pos neg neg
        case 40: // neg pos pos
          {
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P0.C2 + interp * P1.C2;
            A3 = (1-interp) * P0.C3 + interp * P1.C3;
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            B2 = (1-interp) * P0.C2 + interp * P2.C2;
            B3 = (1-interp) * P0.C3 + interp * P2.C3;
            thesegs[2*numsegs-2] = *n_points + numnewpoints;
            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1;
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
          // New Points between 0 and 1, 1 and 2
        case 8:  // neg pos neg
        case 34: // pos neg pos
          {
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P0.C2 + interp * P1.C2;
            A3 = (1-interp) * P0.C3 + interp * P1.C3;
            interp = interpolate
              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            B2 = (1-interp) * P1.C2 + interp * P2.C2;
            B3 = (1-interp) * P1.C3 + interp * P2.C3;
            thesegs[2*numsegs-2] = *n_points + numnewpoints;
            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1;
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
          // New Points between 0 and 2, 1 and 2
        case 10: // pos pos neg
        case 32: // neg neg pos
          {
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P0.C2 + interp * P2.C2;
            A3 = (1-interp) * P0.C3 + interp * P2.C3;
            interp = interpolate
              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            B2 = (1-interp) * P1.C2 + interp * P2.C2;
            B3 = (1-interp) * P1.C3 + interp * P2.C3;
            thesegs[2*numsegs-2] = *n_points + numnewpoints;
            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1;
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
          /* Degenerate cases: two zeroes mean no new points
        case 5:  // 0 0 neg
        case 37: // 0 0 pos
          {
            A2 = B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri];
            thesegs[2*numsegs-1] = triangle_indices[3*tri+1];
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
        case 17: // 0 neg 0
        case 25: // 0 pos 0
          {
            A2 = B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri];
            thesegs[2*numsegs-1] = triangle_indices[3*tri+2];
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
        case 20: // neg 0 0
        case 22: // pos 0 0
          {
            A2 = B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri+1];
            thesegs[2*numsegs-1] = triangle_indices[3*tri+2];
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
            }*/
          // Degenerate cases: one zero one pos one neg means one new point
        case 9:  // 0 pos neg
        case 33: // 0 neg pos
          {
            interp = interpolate
              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P1.C2 + interp * P2.C2;
            A3 = (1-interp) * P1.C3 + interp * P2.C3;
            B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri];
            thesegs[2*numsegs-1] = *n_points + numnewpoints;
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
        case 6:  // pos 0 neg
        case 36: // neg 0 pos
          {
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P0.C2 + interp * P2.C2;
            A3 = (1-interp) * P0.C3 + interp * P2.C3;
            B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri+1];
            thesegs[2*numsegs-1] = *n_points + numnewpoints;
#ifdef DEBUG
            printf ("%d, %d\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1]);
#endif
            break;
          }
        case 18: // pos neg 0
        case 24: // neg pos 0
          {
            interp = interpolate
              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9);
            A2 = (1-interp) * P0.C2 + interp * P1.C2;
            A3 = (1-interp) * P0.C3 + interp * P1.C3;
            B2 = -1;
            thesegs[2*numsegs-2] = triangle_indices[3*tri+2];
            thesegs[2*numsegs-1] = *n_points + numnewpoints;
#ifdef DEBUG
            printf ("%d, %d (%d)\n", thesegs[2*numsegs-2], thesegs[2*numsegs-1], total);
#endif
            break;
          }
        }

      if (A2 > -1.)
        {
          if (++numnewpoints >= thenewpoints_size)
            if (!(thenewpoints = (energy_point *)
                  realloc (thenewpoints,
                           (thenewpoints_size*=2)*sizeof (energy_point))))
              { free (thenewpoints); free (thesegs); return 1; }
          thenewpoints [numnewpoints-1] .C2 = A2;
          thenewpoints [numnewpoints-1] .C3 = A3;
          thenewpoints [numnewpoints-1] .efunc = efunc;
          thenewpoints [numnewpoints-1] .T = -1;
          thenewpoints [numnewpoints-1] .P = -1;
#ifdef DEBUG
          printf ("New point A%d: %g, %g\n", *n_points+numnewpoints-1, A2, A3);
#endif
        }

      if (B2 > -1.)
        {
          if (++numnewpoints >= thenewpoints_size)
            if (!(thenewpoints = (energy_point *)
                  realloc (thenewpoints,
                           (thenewpoints_size*=2)*sizeof (energy_point))))
              { free (thenewpoints); free (thesegs); return 1; }
          thenewpoints [numnewpoints-1] .C2 = B2;
          thenewpoints [numnewpoints-1] .C3 = B3;
          thenewpoints [numnewpoints-1] .efunc = efunc;
          thenewpoints [numnewpoints-1] .T = -1;
          thenewpoints [numnewpoints-1] .P = -1;
#ifdef DEBUG
          printf ("New point B%d: %g, %g\n", *n_points+numnewpoints-1, B2, B3);
#endif
        }
    }

  // Calculate new points' free energies
  if (tri=free_energies (thenewpoints, numnewpoints, T, P, eparams, efunc,
                         WITH_DERIVATIVES_NO_INFINITY))
    {
      free(thenewpoints);
      free(thesegs);
      return tri;
    }

  // Look for duplicate segs and points among numnewpoints
  // *** THIS IS A TODO... ***


  // Add new points to old points array
  if (!(*points = (energy_point *) realloc
        (*points, (*n_points+numnewpoints)*sizeof (energy_point))))
    { free (thenewpoints); free (thesegs); return 1; }
  for (tri=0; tri<numnewpoints; tri++)
    (*points) [*n_points+tri] = thenewpoints [tri];
  *n_points += numnewpoints;

  // Fill in spinreturn
  if (spinreturn)
    {
      spinreturn->compos = 2;
      if (numsegs)
        {
          if (!(thesegs = (int *) realloc
                (thesegs, 2*numsegs*sizeof (int))))
            { free (thenewpoints); free (thesegs); return 1; }
        }
      else
        {
          free (thesegs);
          thesegs = NULL;
        }
      spinreturn->edges = thesegs;
#ifdef DEBUG
      printf ("First edge: %d %d\n", spinreturn->edges [0], spinreturn->edges [1]);
#endif
      spinreturn->n_edges = numsegs;
      if (!(spinreturn->phases = (int *) malloc (sizeof (int))))
        { free (thenewpoints); free (thesegs); return 1; }
      *(spinreturn->phases) = efunc;
      spinreturn->n_phases = 1;
      spinreturn->type = SPINODAL;
    }

  free(thenewpoints);
  return 0;
}