diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4a2887b..e9d1c88 100644 --- a/src/CMakeLists.txt +++ b/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}) diff --git a/src/Determinant.cc b/src/Determinant.cc index 44c911b..034b4af 100644 --- a/src/Determinant.cc +++ b/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**); } diff --git a/src/LinearSolver.cc b/src/LinearSolver.cc index fcb4b3e..2d67522 100644 --- a/src/LinearSolver.cc +++ b/src/LinearSolver.cc @@ -1,3 +1,4 @@ +#include "ParabolicConcSolverBinaryThreePhase.h" #include "QuadraticConcSolverBinaryThreePhase.h" #include "Determinant.h" @@ -109,8 +110,9 @@ int LinearSolver; +template class LinearSolver<3, QuadraticConcSolverBinaryThreePhase, double>; +template class LinearSolver<3, ParabolicConcSolverBinaryThreePhase, double>; + #ifdef HAVE_OPENMP_OFFLOAD #pragma omp end declare target #endif diff --git a/src/LinearSolver.h b/src/LinearSolver.h index 10a5804..0403210 100644 --- a/src/LinearSolver.h +++ b/src/LinearSolver.h @@ -3,7 +3,7 @@ namespace Thermo4PFM { -template +template 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 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(this)->Jacobian(x, fjac); } diff --git a/src/ParabolicConcSolverBinaryThreePhase.cc b/src/ParabolicConcSolverBinaryThreePhase.cc new file mode 100644 index 0000000..02d5e09 --- /dev/null +++ b/src/ParabolicConcSolverBinaryThreePhase.cc @@ -0,0 +1,95 @@ +#include "ParabolicConcSolverBinaryThreePhase.h" +#include "Determinant.h" + +#include + +#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 +} diff --git a/src/ParabolicConcSolverBinaryThreePhase.h b/src/ParabolicConcSolverBinaryThreePhase.h new file mode 100644 index 0000000..1222cc6 --- /dev/null +++ b/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 diff --git a/src/ParabolicFreeEnergyFunctionsBinaryThreePhase.cc b/src/ParabolicFreeEnergyFunctionsBinaryThreePhase.cc new file mode 100644 index 0000000..e7a759e --- /dev/null +++ b/src/ParabolicFreeEnergyFunctionsBinaryThreePhase.cc @@ -0,0 +1,345 @@ +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" +#include "ParabolicConcSolverBinaryThreePhase.h" +#include "functions.h" + +#include +#include + +namespace Thermo4PFM +{ + +ParabolicFreeEnergyFunctionsBinaryThreePhase:: + ParabolicFreeEnergyFunctionsBinaryThreePhase(const double Tref, + const double coeffL[][2], const double coeffA[][2], + const double coeffB[][2], + const EnergyInterpolationType energy_interp_func_type, + const ConcInterpolationType conc_interp_func_type) + : Tref_(Tref), + energy_interp_func_type_(energy_interp_func_type), + conc_interp_func_type_(conc_interp_func_type) +{ + aL_[0] = coeffL[0][0]; + aL_[1] = coeffL[0][1]; + bL_[0] = coeffL[1][0]; + bL_[1] = coeffL[1][1]; + cL_[0] = coeffL[2][0]; + cL_[1] = coeffL[2][1]; + + aA_[0] = coeffA[0][0]; + aA_[1] = coeffA[0][1]; + bA_[0] = coeffA[1][0]; + bA_[1] = coeffA[1][1]; + cA_[0] = coeffA[2][0]; + cA_[1] = coeffA[2][1]; + + aB_[0] = coeffB[0][0]; + aB_[1] = coeffB[0][1]; + bB_[0] = coeffB[1][0]; + bB_[1] = coeffB[1][1]; + cB_[0] = coeffB[2][0]; + cB_[1] = coeffB[2][1]; + + std::string fenergy_diag_filename("energy.vtk"); + fenergy_diag_filename_ = new char[fenergy_diag_filename.length() + 1]; + strcpy(fenergy_diag_filename_, fenergy_diag_filename.c_str()); +} + +//----------------------------------------------------------------------- +#ifdef HAVE_OPENMP_OFFLOAD +#pragma omp declare target +#endif +double ParabolicFreeEnergyFunctionsBinaryThreePhase::computeFreeEnergy( + const double temperature, const double* const conc, const PhaseIndex pi, + const bool gp) +{ + double a0, a1, b0, b1, c0, c1; + switch (pi) + { + case PhaseIndex::phaseL: + a0 = aL_[0]; + a1 = aL_[1]; + b0 = bL_[0]; + b1 = bL_[1]; + c0 = cL_[0]; + c1 = cL_[1]; + break; + case PhaseIndex::phaseA: + a0 = aA_[0]; + a1 = aA_[1]; + b0 = bA_[0]; + b1 = bA_[1]; + c0 = cA_[0]; + c1 = cA_[1]; + break; + case PhaseIndex::phaseB: + a0 = aB_[0]; + a1 = aB_[1]; + b0 = bB_[0]; + b1 = bB_[1]; + c0 = cB_[0]; + c1 = cB_[1]; + break; + default: + return 0.; + } + + const double c = conc[0]; + const double t = (temperature - Tref_); + double fe = 0.5 * (a1 * t + a0) * c * c + (b1 * t + b0) * c + (c1 * t + c0); + + // subtract -mu*c to get grand potential + if (gp) + { + double deriv; + computeDerivFreeEnergy(temperature, conc, pi, &deriv); + fe -= deriv * conc[0]; + } + return fe; +} + +//======================================================================= + +void ParabolicFreeEnergyFunctionsBinaryThreePhase::computeDerivFreeEnergy( + const double temperature, const double* const conc, const PhaseIndex pi, + double* deriv) +{ + double a0, a1, b0, b1; + switch (pi) + { + case PhaseIndex::phaseL: + a0 = aL_[0]; + a1 = aL_[1]; + b0 = bL_[0]; + b1 = bL_[1]; + break; + case PhaseIndex::phaseA: + a0 = aA_[0]; + a1 = aA_[1]; + b0 = bA_[0]; + b1 = bA_[1]; + break; + case PhaseIndex::phaseB: + a0 = aB_[0]; + a1 = aB_[1]; + b0 = bB_[0]; + b1 = bB_[1]; + break; + default: +#ifndef HAVE_OPENMP_OFFLOAD + abort(); +#endif + } + + const double c = conc[0]; + const double t = (temperature - Tref_); + double mu = (a1 * t + a0) * c + (b1 * t + b0); + + deriv[0] = mu; +} + +//======================================================================= + +void ParabolicFreeEnergyFunctionsBinaryThreePhase:: + computeSecondDerivativeFreeEnergy(const double temp, + const double* const conc, const PhaseIndex pi, double* d2fdc2) +{ + (void)conc; + + double a0, a1; + switch (pi) + { + case PhaseIndex::phaseL: + a0 = aL_[0]; + a1 = aL_[1]; + break; + case PhaseIndex::phaseA: + a0 = aA_[0]; + a1 = aA_[1]; + break; + case PhaseIndex::phaseB: + a0 = aB_[0]; + a1 = aB_[1]; + break; + default: +#ifndef HAVE_OPENMP_OFFLOAD + abort(); +#endif + } + + d2fdc2[0] = (a1 * (temp - Tref_) + a0); +} + +//======================================================================= + +void ParabolicFreeEnergyFunctionsBinaryThreePhase::computePhasesFreeEnergies( + const double temperature, const double* const hphi, const double conc, + double& fl, double& fa, double& fb) +{ + // std::cout<<"ParabolicFreeEnergyFunctionsBinary::computePhasesFreeEnergies()"< + +#include +#include +#include +#include + +namespace Thermo4PFM +{ + +class ParabolicFreeEnergyFunctionsBinaryThreePhase +{ +public: + ParabolicFreeEnergyFunctionsBinaryThreePhase(const double Tref, + const double coeffL[][2], const double coeffA[][2], + const double coeffB[][2], + const EnergyInterpolationType energy_interp_func_type, + const ConcInterpolationType conc_interp_func_type); + + ~ParabolicFreeEnergyFunctionsBinaryThreePhase() + { + delete[] fenergy_diag_filename_; + }; + + double computeFreeEnergy(const double temperature, const double* const conc, + const PhaseIndex pi, const bool gp = false); + void computeDerivFreeEnergy(const double temperature, + const double* const conc, const PhaseIndex pi, double*); + void computeSecondDerivativeFreeEnergy(const double temp, + const double* const conc, const PhaseIndex pi, double* d2fdc2); + + int computePhaseConcentrations(const double temperature, const double* conc, + const double* const phi, double* x); + void printEnergyVsComposition(const double temperature, std::ostream& os, + const double cmin, const double cmax, const int npts = 100); + double fchem(const double* const phi, const double* const conc, + const double temperature); + +private: + double Tref_; + /* + * Phase L coefficients + * g(c,T) = (aL1_*T+aL0_)*c*c + (bL1_*T+bL0_)*c + cL1_*T+cL0_ + */ + double aL_[2]; + double bL_[2]; + double cL_[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]; + double cA_[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]; + double cB_[2]; + + EnergyInterpolationType energy_interp_func_type_; + ConcInterpolationType conc_interp_func_type_; + + char* fenergy_diag_filename_; + + void computePhasesFreeEnergies(const double temperature, + const double* const hphi, const double conc, double& fl, double& fa, + double& fb); +}; +} +#endif diff --git a/src/QuadraticConcSolverBinaryThreePhase.cc b/src/QuadraticConcSolverBinaryThreePhase.cc index 97dd6f8..2e3dc38 100644 --- a/src/QuadraticConcSolverBinaryThreePhase.cc +++ b/src/QuadraticConcSolverBinaryThreePhase.cc @@ -24,7 +24,7 @@ void QuadraticConcSolverBinaryThreePhase::RHS( //======================================================================= void QuadraticConcSolverBinaryThreePhase::Jacobian( - const double* const c, JacobianDataType** const fjac) + const double* const c, double** const fjac) { (void)c; diff --git a/src/QuadraticConcSolverBinaryThreePhase.h b/src/QuadraticConcSolverBinaryThreePhase.h index d541483..46fe4b4 100644 --- a/src/QuadraticConcSolverBinaryThreePhase.h +++ b/src/QuadraticConcSolverBinaryThreePhase.h @@ -7,8 +7,7 @@ namespace Thermo4PFM { class QuadraticConcSolverBinaryThreePhase - : public LinearSolver<3, QuadraticConcSolverBinaryThreePhase, - JacobianDataType> + : public LinearSolver<3, QuadraticConcSolverBinaryThreePhase, double> { public: #ifdef HAVE_OPENMP_OFFLOAD @@ -35,7 +34,7 @@ class QuadraticConcSolverBinaryThreePhase /// evaluate Jacobian of system of equations /// specific to this solver - void Jacobian(const double* const x, JacobianDataType** const fjac); + void Jacobian(const double* const x, double** const fjac); #ifdef HAVE_OPENMP_OFFLOAD #pragma omp end declare target #endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ab8cd37..18140ba 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -62,6 +62,10 @@ add_executable(testQuadraticFreeEnergyBinaryThreePhase ${TEST_DIR}/testQuadraticFreeEnergyBinaryThreePhase.cc ) add_executable(testQuadraticFreeEnergyTernaryThreePhase ${TEST_DIR}/testQuadraticFreeEnergyTernaryThreePhase.cc ) +add_executable(testParabolicFreeEnergyBinaryThreePhase + ${TEST_DIR}/testParabolicFreeEnergyBinaryThreePhase.cc ) +add_executable(testParabolicConcSolverBinaryThreePhase + ${TEST_DIR}/testParabolicConcSolverBinaryThreePhase.cc ) if(${WITH_OPENMP_OFFLOAD}) add_executable(testOpenMPoffload @@ -103,7 +107,8 @@ target_link_libraries(testQuadraticBinary thermo4pfm) target_link_libraries(testQuadraticConcSolverBinaryThreePhase thermo4pfm) target_link_libraries(testQuadraticFreeEnergyBinaryThreePhase thermo4pfm) target_link_libraries(testQuadraticFreeEnergyTernaryThreePhase thermo4pfm) - +target_link_libraries(testParabolicFreeEnergyBinaryThreePhase thermo4pfm) +target_link_libraries(testParabolicConcSolverBinaryThreePhase thermo4pfm) target_link_libraries(testLoopCALPHADSpeciesPhaseGibbsEnergy thermo4pfm) target_link_libraries(testLoopCALPHADbinaryKKS thermo4pfm) target_link_libraries(testLoopCALPHADbinaryEquilibrium thermo4pfm) @@ -274,6 +279,16 @@ add_test(NAME testQuadraticFreeEnergyTernaryThreePhase ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/testQuadraticFreeEnergyTernaryThreePhase WORKING_DIRECTORY ${TEST_DIR}) +add_test(NAME testParabolicFreeEnergyBinaryThreePhase + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} + ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testParabolicFreeEnergyBinaryThreePhase + WORKING_DIRECTORY ${TEST_DIR}) +add_test(NAME testParabolicConcSolverBinaryThreePhase + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} + ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testParabolicConcSolverBinaryThreePhase + WORKING_DIRECTORY ${TEST_DIR}) if(${WITH_OPENMP_OFFLOAD}) add_test(NAME testOpenMPoffload diff --git a/tests/testParabolicConcSolverBinaryThreePhase.cc b/tests/testParabolicConcSolverBinaryThreePhase.cc new file mode 100644 index 0000000..52898fa --- /dev/null +++ b/tests/testParabolicConcSolverBinaryThreePhase.cc @@ -0,0 +1,116 @@ +#define CATCH_CONFIG_MAIN + +#include "ParabolicConcSolverBinaryThreePhase.h" + +#include "catch.hpp" + +#include +#include +#include + +TEST_CASE("Parabolic conc solver binary three phase KKS, two-phase consistancy", + "[conc solver binary three phase kks, two-phase consistancy]") +{ + // Parameters from Yang, Wang, Xing, Huang, Mat. Today Comm. 28 (2021) + // 102712 Ti-Al + double coeffA[3][2] = { { 22130.9, -4.4354 }, { -10192.4, 2.6104 }, + { -10276.25, -8.7833 } }; + double coeffB[3][2] + = { { 22090.9, -4.4330 }, { -9841.2, 2.4252 }, { -10439.2, -8.8977 } }; + double coeffL[3][2] + = { { 17183.8, -3.9748 }, { -8659.95, 2.4447 }, { -10439.2, -9.6344 } }; + const double Tref = 1500.; + + const double conc = 0.46; + std::cout << "conc = " << conc << std::endl; + + const double temperature = 1520.; + std::cout << "temperature = " << temperature << std::endl; + + const double phi = 0.5; + + // Inputs to the solver + double hphi0 = 1.0 - phi; + double hphi1 = phi; + double hphi2 = 0.0; + + // Create and set up the solver + Thermo4PFM::ParabolicConcSolverBinaryThreePhase solver; + solver.setup( + conc, hphi0, hphi1, hphi2, temperature - Tref, coeffL, coeffA, coeffB); + + // Run the solver + double sol_test[3] = { 0., 0., 0. }; + std::cout << "First test..." << std::endl; + int ret = solver.ComputeConcentration(sol_test); + std::cout << "Solution: " << sol_test[0] << " " << sol_test[1] << " " + << sol_test[2] << std::endl; + REQUIRE(ret >= 0); + + // Plug the solution back into the RHS + double residual[3]; + solver.RHS(sol_test, residual); + + const double tol = 1.e-8; + REQUIRE(std::abs(residual[0]) < 1.1 * tol); + REQUIRE(std::abs(residual[1]) < 1.1 * tol); + REQUIRE(std::abs(residual[2]) < 1.1 * tol); + + // ---------- + // Now do the same but for the second solid phase + hphi0 = 1.0 - phi; + hphi1 = 0.0; + hphi2 = phi; + + Thermo4PFM::ParabolicConcSolverBinaryThreePhase solver2; + solver2.setup( + conc, hphi0, hphi1, hphi2, temperature - Tref, coeffL, coeffA, coeffB); + + // Run the solver + sol_test[0] = 0.; + sol_test[1] = 0.; + sol_test[2] = 0.; + + residual[0] = 0.; + residual[1] = 0.; + residual[2] = 0.; + + std::cout << "Second test..." << std::endl; + ret = solver2.ComputeConcentration(sol_test); + std::cout << "Solution: " << sol_test[0] << " " << sol_test[1] << " " + << sol_test[2] << std::endl; + REQUIRE(ret >= 0); + + // Plug the solution back into the RHS + solver2.RHS(sol_test, residual); + + REQUIRE(std::abs(residual[0]) < 1.1 * tol); + REQUIRE(std::abs(residual[1]) < 1.1 * tol); + REQUIRE(std::abs(residual[2]) < 1.1 * tol); + + // Inputs to the solver + hphi0 = 0.5; + hphi1 = 0.4; + hphi2 = 0.1; + + // setup the solver + solver.setup( + conc, hphi0, hphi1, hphi2, temperature - Tref, coeffL, coeffA, coeffB); + + // Run the solver + sol_test[0] = 0.; + sol_test[1] = 0.; + sol_test[2] = 0.; + std::cout << "Third test..." << std::endl; + ret = solver.ComputeConcentration(sol_test); + std::cout << "Solution: " << sol_test[0] << " " << sol_test[1] << " " + << sol_test[2] << std::endl; + REQUIRE(ret >= 0); + + // Plug the solution back into the RHS + solver.RHS(sol_test, residual); + + REQUIRE(std::abs(residual[0]) < 1.1 * tol); + REQUIRE(std::abs(residual[1]) < 1.1 * tol); + REQUIRE(std::abs(residual[2]) < 1.1 * tol); +} diff --git a/tests/testParabolicFreeEnergyBinaryThreePhase.cc b/tests/testParabolicFreeEnergyBinaryThreePhase.cc new file mode 100644 index 0000000..ef2438d --- /dev/null +++ b/tests/testParabolicFreeEnergyBinaryThreePhase.cc @@ -0,0 +1,79 @@ +#define CATCH_CONFIG_MAIN + +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" +#include "Phases.h" + +#include "catch.hpp" + +#include +#include +#include + +TEST_CASE("Parabolic conc solver binary three phase KKS, two-phase consistancy", + "[conc solver binary three phase kks, two-phase consistancy]") +{ + // Parameters from Yang, Wang, Xing, Huang, Mat. Today Comm. 28 (2021) + // 102712 Ti-Al + double coeffA[3][2] = { { 22130.9, -4.4354 }, { -10192.4, 2.6104 }, + { -10276.25, -8.7833 } }; + double coeffB[3][2] + = { { 22090.9, -4.4330 }, { -9841.2, 2.4252 }, { -10439.2, -8.8977 } }; + double coeffL[3][2] + = { { 17183.8, -3.9748 }, { -8659.95, 2.4447 }, { -10439.2, -9.6344 } }; + const double Tref = 1500.; + + const double conc = 0.46; + std::cout << "conc = " << conc << std::endl; + + const double temperature = 1520.; + std::cout << "temperature = " << temperature << std::endl; + + // Inputs to the solver + double hphi[3] = { 0.32, 0.3, 0.38 }; + + Thermo4PFM::EnergyInterpolationType energy_interp_func_type + = Thermo4PFM::EnergyInterpolationType::PBG; + Thermo4PFM::ConcInterpolationType conc_interp_func_type + = Thermo4PFM::ConcInterpolationType::LINEAR; + + // Create and set up the solver + Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase qfe(Tref, coeffL, + coeffA, coeffB, energy_interp_func_type, conc_interp_func_type); + + const int npts = 100; + std::ofstream os("FvsC_parabolic.dat", std::ios::out); + std::cout << " Compute energies..." << std::endl; + const double cmin = 0.37; + const double cmax = 0.55; + qfe.printEnergyVsComposition(temperature, os, cmin, cmax, npts); + + // Run the solver + double sol_test[3] = { 0., 0., 0. }; + std::cout << "First test..." << std::endl; + int nit + = qfe.computePhaseConcentrations(temperature, &conc, hphi, sol_test); + std::cout << "nit = " << nit << std::endl; + std::cout << sol_test[0] << " " << sol_test[1] << " " << sol_test[2] + << std::endl; + REQUIRE(nit >= 0); + + // Plug the solution back into the derivative computation + Thermo4PFM::PhaseIndex pl = Thermo4PFM::PhaseIndex::phaseL; + double derivl; + qfe.computeDerivFreeEnergy(temperature, sol_test, pl, &derivl); + std::cout << "dfl/dcl = " << derivl << std::endl; + + Thermo4PFM::PhaseIndex pa = Thermo4PFM::PhaseIndex::phaseA; + double deriva; + qfe.computeDerivFreeEnergy(temperature, sol_test + 1, pa, &deriva); + std::cout << "dfa/dca = " << deriva << std::endl; + + Thermo4PFM::PhaseIndex pb = Thermo4PFM::PhaseIndex::phaseB; + double derivb; + qfe.computeDerivFreeEnergy(temperature, sol_test + 2, pb, &derivb); + std::cout << "dfb/dcb = " << derivb << std::endl; + + const double tol = 1.e-8; + REQUIRE(std::abs(deriva - derivl) < tol); + REQUIRE(std::abs(derivb - derivl) < tol); +}