Do not apply polynomial inside KKS solve (#85)

* h(phi) should be computed outside
ORNL · Jan 4, 2025 · 6892d5d · 6892d5d
1 parent 4faaae6
commit 6892d5d
Show file tree

Hide file tree

Showing 4 changed files with 14 additions and 24 deletions.
diff --git a/src/CALPHADFreeEnergyFunctionsTernary.cc b/src/CALPHADFreeEnergyFunctionsTernary.cc
@@ -526,8 +526,8 @@ void CALPHADFreeEnergyFunctionsTernary::computePhasesFreeEnergies(
 //-----------------------------------------------------------------------
 // output: x
 int CALPHADFreeEnergyFunctionsTernary::computePhaseConcentrations(
-    const double temperature, const double* const conc, const double* const phi,
-    double* x)
+    const double temperature, const double* const conc,
+    const double* const hphi, double* x)
 {
     // assert(conc[0] == conc[0]);
     // assert(conc[1] == conc[1]);
@@ -554,11 +554,9 @@ int CALPHADFreeEnergyFunctionsTernary::computePhaseConcentrations(
         L_AB_S, L_AC_S, L_BC_S, L_ABC_S, fA, fB, fC);
     // assert(fC[0] == fC[0]);
 
-    const double hphi = interp_func(conc_interp_func_type_, phi[0]);
-
     CALPHADConcSolverTernary solver;
-    solver.setup(conc[0], conc[1], hphi, RTinv, L_AB_L, L_AC_L, L_BC_L, L_AB_S,
-        L_AC_S, L_BC_S, L_ABC_L, L_ABC_S, fA, fB, fC);
+    solver.setup(conc[0], conc[1], hphi[0], RTinv, L_AB_L, L_AC_L, L_BC_L,
+        L_AB_S, L_AC_S, L_BC_S, L_ABC_L, L_ABC_S, fA, fB, fC);
     int ret = solver.ComputeConcentration(x, newton_tol_, newton_maxits_);
 #ifndef HAVE_OPENMP_OFFLOAD
     if (ret == -1)

diff --git a/src/KKSFreeEnergyFunctionDiluteBinary.cc b/src/KKSFreeEnergyFunctionDiluteBinary.cc
@@ -227,8 +227,8 @@ void KKSFreeEnergyFunctionDiluteBinary::computePhasesFreeEnergies(
 //-----------------------------------------------------------------------
 
 int KKSFreeEnergyFunctionDiluteBinary::computePhaseConcentrations(
-    const double temperature, const double* const conc, const double* const phi,
-    double* x)
+    const double temperature, const double* const conc,
+    const double* const hphi, double* x)
 
 {
 #ifndef HAVE_OPENMP_OFFLOAD
@@ -239,15 +239,13 @@ int KKSFreeEnergyFunctionDiluteBinary::computePhaseConcentrations(
     assert(maxiters_ > 1);
 #endif
 
-    const double hphi = interp_func(conc_interp_func_type_, phi[0]);
-
     const double fB = computeFB(temperature);
 
     // conc could be outside of [0.,1.] in a trial step
     double c0 = conc[0] >= 0. ? conc[0] : 0.;
     c0        = c0 <= 1. ? c0 : 1.;
     KKSdiluteBinaryConcSolver solver;
-    solver.setup(c0, hphi, 1. - hphi, fA_, fB);
+    solver.setup(c0, hphi[0], 1. - hphi[0], fA_, fB);
     int ret = solver.ComputeConcentration(x, tol_, maxiters_);
 
     return ret;

diff --git a/src/ParabolicFreeEnergyFunctionsBinary.cc b/src/ParabolicFreeEnergyFunctionsBinary.cc
@@ -203,16 +203,14 @@ void ParabolicFreeEnergyFunctionsBinary::computePhasesFreeEnergies(
 //-----------------------------------------------------------------------
 
 int ParabolicFreeEnergyFunctionsBinary::computePhaseConcentrations(
-    const double temperature, const double* const conc, const double* const phi,
-    double* x)
+    const double temperature, const double* const conc,
+    const double* const hphi, double* x)
 {
-    const double hphi0 = interp_func(conc_interp_func_type_, phi[0]);
-    const double hphi1 = interp_func(conc_interp_func_type_, phi[1]);
-
     // solve system of equations to find (cl,cs) given conc[0] and hphi
     // x: initial guess and solution
     ParabolicConcSolverBinary solver;
-    solver.setup(conc[0], hphi0, hphi1, temperature - Tref_, coeffL_, coeffA_);
+    solver.setup(
+        conc[0], hphi[0], hphi[1], temperature - Tref_, coeffL_, coeffA_);
     int ret = solver.ComputeConcentration(x);
 #if 0
     if (ret == -1)

diff --git a/src/QuadraticFreeEnergyFunctionsBinaryThreePhase.cc b/src/QuadraticFreeEnergyFunctionsBinaryThreePhase.cc
@@ -170,20 +170,16 @@ void QuadraticFreeEnergyFunctionsBinaryThreePhase::computePhasesFreeEnergies(
 //-----------------------------------------------------------------------
 
 int QuadraticFreeEnergyFunctionsBinaryThreePhase::computePhaseConcentrations(
-    const double temperature, const double* const conc, const double* const phi,
-    double* x)
+    const double temperature, const double* const conc,
+    const double* const hphi, double* x)
 {
     (void)temperature;
 
-    const double hphi0 = interp_func(conc_interp_func_type_, phi[0]);
-    const double hphi1 = interp_func(conc_interp_func_type_, phi[1]);
-    const double hphi2 = interp_func(conc_interp_func_type_, phi[2]);
-
     // solve system of equations to find (cl,cs) given conc[0] and hphi
     // x: initial guess and solution
     QuadraticConcSolverBinaryThreePhase solver;
     solver.setup(
-        conc[0], hphi0, hphi1, hphi2, Al_, ceql_, Aa_, ceqa_, Ab_, ceqb_);
+        conc[0], hphi[0], hphi[1], hphi[2], Al_, ceql_, Aa_, ceqa_, Ab_, ceqb_);
     int ret = solver.ComputeConcentration(x);
 #if 0
     if (ret == -1)