Add parabolic free energies

* to be used for quadratic polynomial fits
* linear solvers use matrices in couble precision

src/CMakeLists.txt

@@ -29,6 +29,8 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/src/InterpolationType.cc
             ${CMAKE_SOURCE_DIR}/src/QuadraticFreeEnergyFunctionsBinaryThreePhase.cc
             ${CMAKE_SOURCE_DIR}/src/QuadraticConcSolverBinaryThreePhase.cc
             ${CMAKE_SOURCE_DIR}/src/QuadraticFreeEnergyFunctionsTernaryThreePhase.cc
+            ${CMAKE_SOURCE_DIR}/src/ParabolicFreeEnergyFunctionsBinaryThreePhase.cc
+            ${CMAKE_SOURCE_DIR}/src/ParabolicConcSolverBinaryThreePhase.cc
             ${CMAKE_SOURCE_DIR}/src/xlogx.cc)
 
 if(${SHARED_LIB})

src/Determinant.cc

@@ -128,4 +128,5 @@ ScalarType evalDeterminant(ScalarType** mat)
 
 template double evalDeterminant<5, double>(double**);
 template float evalDeterminant<5, float>(float**);
+template float evalDeterminant<3, float>(float**);
 }

src/LinearSolver.cc

@@ -1,3 +1,4 @@
+#include "ParabolicConcSolverBinaryThreePhase.h"
 #include "QuadraticConcSolverBinaryThreePhase.h"
 
 #include "Determinant.h"
@@ -109,8 +110,9 @@ int LinearSolver<Dimension, SolverType,
     return 0;
 }
 
-template class LinearSolver<3, QuadraticConcSolverBinaryThreePhase,
-    JacobianDataType>;
+template class LinearSolver<3, QuadraticConcSolverBinaryThreePhase, double>;
+template class LinearSolver<3, ParabolicConcSolverBinaryThreePhase, double>;
+
 #ifdef HAVE_OPENMP_OFFLOAD
 #pragma omp end declare target
 #endif

src/LinearSolver.h

@@ -3,7 +3,7 @@
 
 namespace Thermo4PFM
 {
-template <unsigned int Dimension, class SolverType, typename JacobianDataType>
+template <unsigned int Dimension, class SolverType, typename MatrixDataType>
 class LinearSolver
 {
 public:
@@ -22,8 +22,8 @@ class LinearSolver
 #ifdef HAVE_OPENMP_OFFLOAD
 #pragma omp declare target
 #endif
-    void UpdateSolution(double* const x, const double* const fvec,
-        JacobianDataType** const fjac);
+    void UpdateSolution(
+        double* const x, const double* const fvec, MatrixDataType** const fjac);
 
     int ComputeSolutionInternal(double* const conc);
 
@@ -35,7 +35,7 @@ class LinearSolver
     template <typename ScalarType>
     void CopyMatrix(ScalarType** const dst, ScalarType** const src);
 
-    void internalJacobian(const double* const x, JacobianDataType** const fjac)
+    void internalJacobian(const double* const x, MatrixDataType** const fjac)
     {
         static_cast<SolverType*>(this)->Jacobian(x, fjac);
     }

src/ParabolicConcSolverBinaryThreePhase.cc

@@ -0,0 +1,95 @@
+#include "ParabolicConcSolverBinaryThreePhase.h"
+#include "Determinant.h"
+
+#include <iostream>
+
+#define EQSCALING 1.e-4
+
+namespace Thermo4PFM
+{
+//=======================================================================
+
+#ifdef HAVE_OPENMP_OFFLOAD
+#pragma omp declare target
+#endif
+
+//=======================================================================
+
+// solve for c=(c_L, c_A, c_B)
+void ParabolicConcSolverBinaryThreePhase::RHS(
+    const double* const c, double* const fvec)
+{
+    fvec[0] = -c0_ + hphi0_ * c[0] + hphi1_ * c[1] + hphi2_ * c[2];
+
+    const double t   = temperature_;
+    const double muL = (aL_[1] * t + aL_[0]) * c[0] + (bL_[1] * t + bL_[0]);
+    const double muA = (aA_[1] * t + aA_[0]) * c[1] + (bA_[1] * t + bA_[0]);
+    const double muB = (aB_[1] * t + aB_[0]) * c[2] + (bB_[1] * t + bB_[0]);
+    fvec[1]          = EQSCALING * (muA - muL);
+    fvec[2]          = EQSCALING * (muB - muL);
+}
+
+//=======================================================================
+
+void ParabolicConcSolverBinaryThreePhase::Jacobian(
+    const double* const c, double** const fjac)
+{
+    (void)c;
+
+    const double t = temperature_;
+    fjac[0][0]     = hphi0_;
+    fjac[0][1]     = hphi1_;
+    fjac[0][2]     = hphi2_;
+
+    fjac[1][0] = EQSCALING * (-1. * (aL_[1] * t + aL_[0]));
+    fjac[1][1] = EQSCALING * ((aA_[1] * t + aA_[0]));
+    fjac[1][2] = 0.;
+
+    fjac[2][0] = EQSCALING * (-1. * (aL_[1] * t + aL_[0]));
+    fjac[2][1] = 0.;
+    fjac[2][2] = EQSCALING * (aB_[1] * t + aB_[0]);
+#if 0
+    std::cout << "Matrix:" << std::endl;
+    for (int i = 0; i < 3; i++)
+    {
+        for (int j = 0; j < 3; j++)
+            std::cout << fjac[i][j] << " ";
+        std::cout << std::endl;
+    }
+    double d = evalDeterminant<3, double>(fjac);
+    std::cout << "Determinant: " << d << std::endl;
+#endif
+}
+
+// set values of internal variables
+void ParabolicConcSolverBinaryThreePhase::setup(const double c0,
+    const double hphi0, const double hphi1, const double hphi2,
+    const double temperature, const double coeffL[][2],
+    const double coeffA[][2], const double coeffB[][2])
+{
+    c0_          = c0;
+    hphi0_       = hphi0;
+    hphi1_       = hphi1;
+    hphi2_       = hphi2;
+    temperature_ = temperature;
+
+    aL_[0] = coeffL[0][0];
+    aL_[1] = coeffL[0][1];
+    bL_[0] = coeffL[1][0];
+    bL_[1] = coeffL[1][1];
+
+    aA_[0] = coeffA[0][0];
+    aA_[1] = coeffA[0][1];
+    bA_[0] = coeffA[1][0];
+    bA_[1] = coeffA[1][1];
+
+    aB_[0] = coeffB[0][0];
+    aB_[1] = coeffB[0][1];
+    bB_[0] = coeffB[1][0];
+    bB_[1] = coeffB[1][1];
+}
+
+#ifdef HAVE_OPENMP_OFFLOAD
+#pragma omp end declare target
+#endif
+}

src/ParabolicConcSolverBinaryThreePhase.h

@@ -0,0 +1,79 @@
+#ifndef Thermo4PFM_included_ParabolicConcSolverBinaryThreePhase
+#define Thermo4PFM_included_ParabolicConcSolverBinaryThreePhase
+
+#include "LinearSolver.h"
+#include "datatypes.h"
+
+namespace Thermo4PFM
+{
+class ParabolicConcSolverBinaryThreePhase
+    : public LinearSolver<3, ParabolicConcSolverBinaryThreePhase, double>
+{
+public:
+#ifdef HAVE_OPENMP_OFFLOAD
+#pragma omp declare target
+#endif
+    /// compute "internal" concentrations cL, cS1, cS2 by solving KKS
+    /// equations
+    /// conc: solution (concentration in each phase)
+    int ComputeConcentration(double* const conc)
+    {
+        return LinearSolver::ComputeSolution(conc);
+    }
+
+    /// 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 hphi0, const double hphi1,
+        const double hphi2, const double temperature, const double coeffL[][2],
+        const double coeffA[][2], const double coeffB[][2]);
+
+    /// evaluate RHS of the system of eqautions to solve for
+    /// specific to this solver
+    void RHS(const double* const x, double* const fvec);
+
+    /// evaluate Jacobian of system of equations
+    /// specific to this solver
+    void Jacobian(const double* const x, double** const fjac);
+#ifdef HAVE_OPENMP_OFFLOAD
+#pragma omp end declare target
+#endif
+
+private:
+    ///
+    /// Nominal composition to solve for
+    ///
+    double c0_;
+
+    ///
+    /// phase fractions to solve for
+    ///
+    double hphi0_;
+    double hphi1_;
+    double hphi2_;
+
+    double temperature_;
+
+    /*
+     * Phase L coefficients
+     * g(c,T) = (aL1_*T+aL0_)*c*c + (bL1_*T+bL0_)*c + cL1_*T+cL0_
+     */
+    double aL_[2];
+    double bL_[2];
+
+    /*
+     * Phase A coefficients
+     * g(c,T) = (aL1_*T+aL0_)*c*c + (bL1_*T+bL0_)*c + cL1_*T+cL0_
+     */
+    double aA_[2];
+    double bA_[2];
+
+    /*
+     * Phase B coefficients
+     * g(c,T) = (aL1_*T+aL0_)*c*c + (bL1_*T+bL0_)*c + cL1_*T+cL0_
+     */
+    double aB_[2];
+    double bB_[2];
+};
+}
+#endif