Changeset 481 for trunk/matml

Show
Ignore:
Timestamp:
10/06/2009 04:36:11 PM (3 years ago)
Author:
powell
Message:

Interpolator for more accurate spinodal calculation.

Files:
1 modified

Legend:

Unmodified
Added
Removed

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

    r480 r481  
    99#include "gibbs.h" 
    1010#include <stdlib.h>  /*+ For malloc() and free() +*/ 
     11 
     12 
     13/*++++++++++++++++++++++++++++++++++++++ 
     14  This function uses secant iteration to estimate the root of the spinodal 
     15  function on a line in ternary space between points A and B. 
     16 
     17  inline double interpolate Returns the fraction B whose weighted average with 
     18  A minimizes the spinodal function. 
     19 
     20  double A2 C2 variable at point A. 
     21 
     22  double A3 C3 variable at point A. 
     23 
     24  double AS Spinodal function value at point A. 
     25 
     26  double A2 C2 variable at point B. 
     27 
     28  double A3 C3 variable at point B. 
     29 
     30  double AS Spinodal function value at point B. 
     31 
     32  double T Temperature. 
     33 
     34  double P Pressure (ignored for now). 
     35 
     36  energy_params *eparams Energy parameter functions (all of them). 
     37 
     38  int efunc Index indicating which energy function to calculate the spinodal 
     39  for. 
     40 
     41  int iters Number of secant iterations to perform. 
     42++++++++++++++++++++++++++++++++++++++*/ 
     43 
     44static inline double interpolate 
     45(double A2, double A3, double AS, double B2, double B3, double BS, 
     46 double T, double P, energy_params *eparams, int efunc, int iters) 
     47{ 
     48  int i; 
     49  energy_point X; 
     50  double interp=0., intfrac=1., XS; 
     51 
     52  for (i=0; i<iters; i++) 
     53    { 
     54      X.C2 = (AS*B2 - BS*A2) / (AS - BS); 
     55      X.C3 = (AS*B3 - BS*A3) / (AS - BS); 
     56 
     57      X.efunc = -1; 
     58      free_energies (&X, 1, T,P, eparams,efunc, WITH_DERIVATIVES_NO_INFINITY); 
     59      XS = X.G22 * X.G33 - X.G23 * X.G23; 
     60 
     61      printf ("A=%g,%g->%g B=%g,%g->%g, X=%g,%g->%g\n", A2,A3,AS, B2,B3,BS, 
     62              X.C2,X.C3,XS); 
     63 
     64      if (AS * XS > 0) 
     65        { 
     66          interp += intfrac * AS / (AS - BS); 
     67          intfrac *= -BS / (AS-BS); 
     68          A2 = X.C2; 
     69          A3 = X.C3; 
     70          AS = XS; 
     71        } 
     72      else 
     73        { 
     74          intfrac *= AS / (AS-BS); 
     75          B2 = X.C2; 
     76          B3 = X.C3; 
     77          BS = XS; 
     78        } 
     79    } 
     80 
     81  return interp + intfrac * AS / (AS - BS); 
     82} 
    1183 
    1284 
    … …  
    89161  for (tri=0; tri<n_triangles; tri++) 
    90162    { 
    91       double spin[3], A2=-1, A3, B2=-1, B3; 
     163      double spin[3], A2=-1, A3, B2=-1, B3, interp; 
    92164      int total; 
    93165      energy_point P0=(*points)[triangle_indices[3*tri]], 
    … …  
    157229        case 40: // neg pos pos 
    158230          { 
    159             A2 = (spin[0]*P1.C2 - spin[1]*P0.C2)/(spin[0]-spin[1]); 
    160             A3 = (spin[0]*P1.C3 - spin[1]*P0.C3)/(spin[0]-spin[1]); 
    161             B2 = (spin[0]*P2.C2 - spin[2]*P0.C2)/(spin[0]-spin[2]); 
    162             B3 = (spin[0]*P2.C3 - spin[2]*P0.C3)/(spin[0]-spin[2]); 
     231            interp = interpolate 
     232              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9); 
     233            A2 = (1-interp) * P0.C2 + interp * P1.C2; 
     234            A3 = (1-interp) * P0.C3 + interp * P1.C3; 
     235            interp = interpolate 
     236              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     237            B2 = (1-interp) * P0.C2 + interp * P2.C2; 
     238            B3 = (1-interp) * P0.C3 + interp * P2.C3; 
    163239            thesegs[2*numsegs-2] = *n_points + numnewpoints; 
    164240            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1; 
    … …  
    172248        case 34: // pos neg pos 
    173249          { 
    174             A2 = (spin[0]*P1.C2 - spin[1]*P0.C2)/(spin[0]-spin[1]); 
    175             A3 = (spin[0]*P1.C3 - spin[1]*P0.C3)/(spin[0]-spin[1]); 
    176             B2 = (spin[1]*P2.C2 - spin[2]*P1.C2)/(spin[1]-spin[2]); 
    177             B3 = (spin[1]*P2.C3 - spin[2]*P1.C3)/(spin[1]-spin[2]); 
     250            interp = interpolate 
     251              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9); 
     252            A2 = (1-interp) * P0.C2 + interp * P1.C2; 
     253            A3 = (1-interp) * P0.C3 + interp * P1.C3; 
     254            interp = interpolate 
     255              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     256            B2 = (1-interp) * P1.C2 + interp * P2.C2; 
     257            B3 = (1-interp) * P1.C3 + interp * P2.C3; 
    178258            thesegs[2*numsegs-2] = *n_points + numnewpoints; 
    179259            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1; 
    … …  
    187267        case 32: // neg neg pos 
    188268          { 
    189             A2 = (spin[0]*P2.C2 - spin[2]*P0.C2)/(spin[0]-spin[2]); 
    190             A3 = (spin[0]*P2.C3 - spin[2]*P0.C3)/(spin[0]-spin[2]); 
    191             B2 = (spin[1]*P2.C2 - spin[2]*P1.C2)/(spin[1]-spin[2]); 
    192             B3 = (spin[1]*P2.C3 - spin[2]*P1.C3)/(spin[1]-spin[2]); 
     269            interp = interpolate 
     270              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     271            A2 = (1-interp) * P0.C2 + interp * P2.C2; 
     272            A3 = (1-interp) * P0.C3 + interp * P2.C3; 
     273            interp = interpolate 
     274              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     275            B2 = (1-interp) * P1.C2 + interp * P2.C2; 
     276            B3 = (1-interp) * P1.C3 + interp * P2.C3; 
    193277            thesegs[2*numsegs-2] = *n_points + numnewpoints; 
    194278            thesegs[2*numsegs-1] = *n_points + numnewpoints + 1; 
    … …  
    236320        case 33: // 0 neg pos 
    237321          { 
    238             A2 = (spin[1]*P2.C2 - spin[2]*P1.C2)/(spin[1]-spin[2]); 
    239             A3 = (spin[1]*P2.C3 - spin[2]*P1.C3)/(spin[1]-spin[2]); 
     322            interp = interpolate 
     323              (P1.C2,P1.C3,spin[1], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     324            A2 = (1-interp) * P1.C2 + interp * P2.C2; 
     325            A3 = (1-interp) * P1.C3 + interp * P2.C3; 
    240326            B2 = -1; 
    241327            thesegs[2*numsegs-2] = triangle_indices[3*tri]; 
    … …  
    249335        case 36: // neg 0 pos 
    250336          { 
    251             A2 = (spin[0]*P2.C2 - spin[2]*P0.C2)/(spin[0]-spin[2]); 
    252             A3 = (spin[0]*P2.C3 - spin[2]*P0.C3)/(spin[0]-spin[2]); 
     337            interp = interpolate 
     338              (P0.C2,P0.C3,spin[0], P2.C2,P2.C3,spin[2], T,P,eparams,efunc, 9); 
     339            A2 = (1-interp) * P0.C2 + interp * P2.C2; 
     340            A3 = (1-interp) * P0.C3 + interp * P2.C3; 
    253341            B2 = -1; 
    254342            thesegs[2*numsegs-2] = triangle_indices[3*tri+1]; 
    … …  
    262350        case 24: // neg pos 0 
    263351          { 
    264             A2 = (spin[0]*P1.C2 - spin[1]*P0.C2)/(spin[0]-spin[1]); 
    265             A3 = (spin[0]*P1.C3 - spin[1]*P0.C3)/(spin[0]-spin[1]); 
     352            interp = interpolate 
     353              (P0.C2,P0.C3,spin[0], P1.C2,P1.C3,spin[1], T,P,eparams,efunc, 9); 
     354            A2 = (1-interp) * P0.C2 + interp * P1.C2; 
     355            A3 = (1-interp) * P0.C3 + interp * P1.C3; 
    266356            B2 = -1; 
    267357            thesegs[2*numsegs-2] = triangle_indices[3*tri+2];