Generalize interface to binary solvers (#63)

Use phi0 and phi1 instead of phi and 1-phi

drivers/loopCALPHADConcSolverBinary.cc

@@ -123,7 +123,7 @@ int main(int argc, char* argv[])
             double hphi      = 0.5 + (i % 100) * deviation;
             double c0        = 0.3;
             Thermo4PFM::CALPHADConcSolverBinary solver;
-            solver.setup(c0, hphi, RTinv, Lmix_L, Lmix_A, fA, fB);
+            solver.setup(c0, hphi, 1. - hphi, RTinv, Lmix_L, Lmix_A, fA, fB);
             nits[i] = solver.ComputeConcentration(&xhost[2 * i], 1.e-8, 50);
         }
         auto t2 = Clock::now();

src/CALPHADConcSolverBinary.cc

@@ -37,14 +37,14 @@ void CALPHADConcSolverBinary::computeXi(
 
 //=======================================================================
 
-// solve for c=(c_L, c_A, c_B)
+// solve for c=(c_L, c_A)
 void CALPHADConcSolverBinary::RHS(const double* const c, double* const fvec)
 {
     double xi[2] = { 0., 0. };
 
     computeXi(c, xi);
 
-    fvec[0] = -c0_ + (1.0 - hphi_) * c[0] + hphi_ * c[1];
+    fvec[0] = -c0_ + hphi1_ * c[0] + hphi0_ * c[1];
     fvec[1] = xlogx_deriv(c[0]) - xlogx_deriv(1. - c[0]) - xlogx_deriv(c[1])
               + xlogx_deriv(1. - c[1]) + (xi[0] - xi[1]);
 }
@@ -74,8 +74,8 @@ void CALPHADConcSolverBinary::Jacobian(
     double dxidc[2];
     computeDxiDc(c, dxidc);
 
-    fjac[0][0] = (1.0 - hphi_);
-    fjac[0][1] = hphi_;
+    fjac[0][0] = hphi1_;
+    fjac[0][1] = hphi0_;
 
     fjac[1][0] = dxidc[0] + xlogx_deriv2(c[0]) + xlogx_deriv2(1. - c[0]);
 
@@ -84,13 +84,14 @@ void CALPHADConcSolverBinary::Jacobian(
 
 // set values of internal variables used to evaluate
 // terms in Newton iterations
-void CALPHADConcSolverBinary::setup(const double c0, const double hphi,
-    const double RTinv, const CalphadDataType* const Lmix_L,
+void CALPHADConcSolverBinary::setup(const double c0, const double hphi0,
+    const double hphi1, const double RTinv, const CalphadDataType* const Lmix_L,
     const CalphadDataType* const Lmix_A, const CalphadDataType* const fA,
     const CalphadDataType* const fB)
 {
     c0_    = c0;
-    hphi_  = hphi;
+    hphi0_ = hphi0;
+    hphi1_ = hphi1;
     RTinv_ = RTinv;
 
     for (int ii = 0; ii < 4; ii++)

src/CALPHADConcSolverBinary.h

@@ -25,8 +25,8 @@ class CALPHADConcSolverBinary
     /// setup model paramater values to be used by solver,
     /// including composition "c0" and phase fraction "hphi"
     /// to solve for
-    void setup(const double c0, const double hphi, const double RTinv,
-        const CalphadDataType* const Lmix_L_,
+    void setup(const double c0, const double hphi0, const double hphi1,
+        const double RTinv, const CalphadDataType* const Lmix_L_,
         const CalphadDataType* const Lmix_A_, const CalphadDataType* const fA,
         const CalphadDataType* const fB);
 
@@ -50,7 +50,8 @@ class CALPHADConcSolverBinary
     ///
     /// phase fraction to solve for
     ///
-    double hphi_;
+    double hphi0_;
+    double hphi1_;
 
     double RTinv_;
 

src/CALPHADFreeEnergyFunctionsBinary.cc

@@ -269,7 +269,7 @@ void CALPHADFreeEnergyFunctionsBinary::computePhasesFreeEnergies(
     double RTinv = 1.0 / (GASCONSTANT_R_JPKPMOL * temperature);
 
     CALPHADConcSolverBinary solver;
-    solver.setup(conc, hphi[0], RTinv, Lmix_L, Lmix_A, fA, fB);
+    solver.setup(conc, hphi[0], 1. - hphi[0], RTinv, Lmix_L, Lmix_A, fA, fB);
     int ret = solver.ComputeConcentration(
         c, newton_tol_, newton_maxits_, newton_alpha_);
     if (ret < 0)
@@ -311,7 +311,7 @@ int CALPHADFreeEnergyFunctionsBinary::computePhaseConcentrations(
     // solve system of equations to find (cl,cs) given conc[0] and hphi
     // x: initial guess and solution
     CALPHADConcSolverBinary solver;
-    solver.setup(conc[0], hphi, RTinv, Lmix_L, Lmix_A, fA, fB);
+    solver.setup(conc[0], hphi, 1. - hphi, RTinv, Lmix_L, Lmix_A, fA, fB);
     int ret = solver.ComputeConcentration(
         x, newton_tol_, newton_maxits_, newton_alpha_);
 #if 0

src/KKSFreeEnergyFunctionDiluteBinary.cc

@@ -203,7 +203,7 @@ void KKSFreeEnergyFunctionDiluteBinary::computePhasesFreeEnergies(
     double fB = computeFB(temperature);
 
     KKSdiluteBinaryConcSolver solver;
-    solver.setup(conc, *hphi, fA_, fB);
+    solver.setup(conc, *hphi, 1. - *hphi, fA_, fB);
     int ret = solver.ComputeConcentration(c, tol_, maxiters_, alpha_);
 
 #ifndef HAVE_OPENMP_OFFLOAD
@@ -247,7 +247,7 @@ int KKSFreeEnergyFunctionDiluteBinary::computePhaseConcentrations(
     double c0 = conc[0] >= 0. ? conc[0] : 0.;
     c0        = c0 <= 1. ? c0 : 1.;
     KKSdiluteBinaryConcSolver solver;
-    solver.setup(c0, hphi, fA_, fB);
+    solver.setup(c0, hphi, 1. - hphi, fA_, fB);
     int ret = solver.ComputeConcentration(x, tol_, maxiters_);
 
     return ret;

src/KKSdiluteBinaryConcSolver.cc

@@ -9,7 +9,7 @@ namespace Thermo4PFM
 #endif
 void KKSdiluteBinaryConcSolver::RHS(const double* const c, double* const fvec)
 {
-    fvec[0] = -c0_ + (1.0 - hphi_) * c[0] + hphi_ * c[1];
+    fvec[0] = -c0_ + hphi1_ * c[0] + hphi0_ * c[1];
     fvec[1] = xlogx_deriv(c[0]) - xlogx_deriv(1. - c[0]) - xlogx_deriv(c[1])
               + xlogx_deriv(1. - c[1]) - (fA_ - fB_);
 }
@@ -19,20 +19,21 @@ void KKSdiluteBinaryConcSolver::RHS(const double* const c, double* const fvec)
 void KKSdiluteBinaryConcSolver::Jacobian(
     const double* const c, JacobianDataType** const fjac)
 {
-    fjac[0][0] = (1.0 - hphi_);
-    fjac[0][1] = hphi_;
+    fjac[0][0] = hphi1_;
+    fjac[0][1] = hphi0_;
 
     fjac[1][0] = xlogx_deriv2(c[0]) + xlogx_deriv2(1. - c[0]);
     fjac[1][1] = -xlogx_deriv2(c[1]) - xlogx_deriv2(1. - c[1]);
 }
 
-void KKSdiluteBinaryConcSolver::setup(
-    const double c0, const double hphi, const double fA, const double fB)
+void KKSdiluteBinaryConcSolver::setup(const double c0, const double hphi0,
+    const double hphi1, const double fA, const double fB)
 {
-    c0_   = c0;
-    hphi_ = hphi;
-    fA_   = fA;
-    fB_   = fB;
+    c0_    = c0;
+    hphi0_ = hphi0;
+    hphi1_ = hphi1;
+    fA_    = fA;
+    fB_    = fB;
 }
 #ifdef HAVE_OPENMP_OFFLOAD
 #pragma omp end declare target

src/KKSdiluteBinaryConcSolver.h

@@ -25,8 +25,8 @@ class KKSdiluteBinaryConcSolver
     /// setup model paramater values to be used by solver,
     /// including composition "c0" and phase fraction "hphi"
     /// to solve for
-    void setup(
-        const double c0, const double hphi, const double fA, const double fB);
+    void setup(const double c0, const double hphi0, const double hphi1,
+        const double fA, const double fB);
 
     /// evaluate RHS of the system of eqautions to solve for
     /// specific to this solver
@@ -46,9 +46,10 @@ class KKSdiluteBinaryConcSolver
     double c0_;
 
     ///
-    /// phase fraction to solve for
+    /// phase fractiona to solve for
     ///
-    double hphi_;
+    double hphi0_;
+    double hphi1_;
 
     ///
     /// model parameters

tests/testCALPHADJacobianBinary.cc

@@ -40,7 +40,7 @@ TEST_CASE("Jacobian binary", "[Jacobian binary]")
         std::cout << "Test CALPHADConcSolverBinary...\n";
         {
             Thermo4PFM::CALPHADConcSolverBinary solver;
-            solver.setup(c0, hphi, RTinv, Lmix_L, Lmix_S, fA, fB);
+            solver.setup(c0, hphi, 1. - hphi, RTinv, Lmix_L, Lmix_S, fA, fB);
 
             solver.RHS(x, fvec1);