diff --git a/drivers/CMakeLists.txt b/drivers/CMakeLists.txt index bd66545c..34a1b3c4 100644 --- a/drivers/CMakeLists.txt +++ b/drivers/CMakeLists.txt @@ -9,6 +9,11 @@ add_executable(plotCALPHADbinary target_link_libraries(plotCALPHADbinary ${PROJECT_LINK_LIBS_SAMRAI} ${THERMOLIB}) +add_executable(plotEnergyVsCompositionParabolicThreePhase + ${CMAKE_SOURCE_DIR}/drivers/plotEnergyVsCompositionParabolicThreePhase.cc ) +target_link_libraries(plotEnergyVsCompositionParabolicThreePhase ${PROJECT_LINK_LIBS_SAMRAI} + ${THERMOLIB}) + add_executable(computeEquilibriumCompositions computeEquilibriumCompositions.cc ${CMAKE_SOURCE_DIR}/source/Database2JSON.cc @@ -41,6 +46,22 @@ target_link_libraries(computeCALPHADbinaryEquilibrium ${PROJECT_LINK_LIBS_SAMRAI} ${THERMOLIB}) +add_executable(computeParabolic3PhasesEquilibrium + computeParabolic3PhasesEquilibrium.cc + ${CMAKE_SOURCE_DIR}/source/Database2JSON.cc + ${CMAKE_SOURCE_DIR}/source/fortran/functions.f) +target_link_libraries(computeParabolic3PhasesEquilibrium + ${PROJECT_LINK_LIBS_SAMRAI} + ${THERMOLIB}) + +add_executable(computeParabolic3PhasesKKS + computeParabolic3PhasesKKS.cc + ${CMAKE_SOURCE_DIR}/source/Database2JSON.cc + ${CMAKE_SOURCE_DIR}/source/fortran/functions.f) +target_link_libraries(computeParabolic3PhasesKKS + ${PROJECT_LINK_LIBS_SAMRAI} + ${THERMOLIB}) + add_executable(convertDatabase2JSON ${CMAKE_SOURCE_DIR}/source/Database2JSON.cc convertDatabase2JSON.cc ) diff --git a/drivers/computeParabolic3PhasesEquilibrium.cc b/drivers/computeParabolic3PhasesEquilibrium.cc new file mode 100644 index 00000000..20a06e90 --- /dev/null +++ b/drivers/computeParabolic3PhasesEquilibrium.cc @@ -0,0 +1,110 @@ +#include "ParabolicEqConcSolverBinary.h" + +#include "SAMRAI/tbox/SAMRAIManager.h" +#include "SAMRAI/tbox/InputManager.h" +#include "SAMRAI/tbox/Database.h" +#include "SAMRAI/SAMRAI_config.h" +#include "SAMRAI/tbox/SAMRAI_MPI.h" + +#include +#include +#include + +namespace pt = boost::property_tree; +using namespace SAMRAI; + +int main(int argc, char* argv[]) +{ + tbox::SAMRAI_MPI::init(&argc, &argv); + tbox::SAMRAIManager::initialize(); + tbox::SAMRAIManager::startup(); + + int ret = 0; + { + std::string databasename(argv[1]); + double temperature = atof(argv[2]); + std::cout << "Temperature = " << temperature << std::endl; + double c0 = 0.5; + double c1 = 1.5; + if (argc > 3) { + c0 = atof(argv[3]); + c1 = atof(argv[4]); + } + + std::shared_ptr input_db( + new tbox::MemoryDatabase("db")); + std::cout << "Filename = " << databasename << std::endl; + tbox::InputManager::getManager()->parseInputFile(databasename, input_db); + + double coeffL[3][2]; + std::shared_ptr liquid_db = + input_db->getDatabase("Liquid"); + coeffL[0][0] = liquid_db->getDouble("a0"); + coeffL[0][1] = liquid_db->getDouble("a1"); + coeffL[1][0] = liquid_db->getDouble("b0"); + coeffL[1][1] = liquid_db->getDouble("b1"); + coeffL[2][0] = liquid_db->getDouble("c0"); + coeffL[2][1] = liquid_db->getDouble("c1"); + + double coeffA[3][2]; + std::shared_ptr phasea_db = + input_db->getDatabase("PhaseA"); + coeffA[0][0] = phasea_db->getDouble("a0"); + coeffA[0][1] = phasea_db->getDouble("a1"); + coeffA[1][0] = phasea_db->getDouble("b0"); + coeffA[1][1] = phasea_db->getDouble("b1"); + coeffA[2][0] = phasea_db->getDouble("c0"); + coeffA[2][1] = phasea_db->getDouble("c1"); + + double coeffB[3][2]; + std::shared_ptr phaseb_db = + input_db->getDatabase("PhaseB"); + coeffB[0][0] = phaseb_db->getDouble("a0"); + coeffB[0][1] = phaseb_db->getDouble("a1"); + coeffB[1][0] = phaseb_db->getDouble("b0"); + coeffB[1][1] = phaseb_db->getDouble("b1"); + coeffB[2][0] = phaseb_db->getDouble("c0"); + coeffB[2][1] = phaseb_db->getDouble("c1"); + + double Tref = input_db->getDouble("Tref"); + + input_db->printClassData(std::cout); + + Thermo4PFM::ParabolicEqConcSolverBinary solver; + solver.setup(temperature - Tref, coeffL, coeffA); + + const double toln = 1.e-12; + const int max_iters = 100; + const double alpha = 1.; + + double sol[2] = {c0, c1}; + std::cout << "Phases L,A..." << std::endl; + int ret = solver.ComputeConcentration(sol, toln, max_iters, alpha); + std::cout << ret << " iterations" << std::endl; + std::cout << "Solution: " << sol[0] << " " << sol[1] << std::endl; + + sol[0] = c0; + sol[1] = c1; + + std::cout << "Phases L,B..." << std::endl; + solver.setup(temperature - Tref, coeffL, coeffB); + ret = solver.ComputeConcentration(sol, toln, max_iters, alpha); + std::cout << ret << " iterations" << std::endl; + std::cout << "Solution: " << sol[0] << " " << sol[1] << std::endl; + + sol[0] = c0; + sol[1] = c1; + + std::cout << "Phases A,B..." << std::endl; + solver.setup(temperature - Tref, coeffA, coeffB); + ret = solver.ComputeConcentration(sol, toln, max_iters, alpha); + std::cout << ret << " iterations" << std::endl; + std::cout << "Solution: " << sol[0] << " " << sol[1] << std::endl; + } + + tbox::SAMRAIManager::shutdown(); + tbox::SAMRAIManager::finalize(); + tbox::SAMRAI_MPI::finalize(); + + return ret; +} diff --git a/drivers/computeParabolic3PhasesKKS.cc b/drivers/computeParabolic3PhasesKKS.cc new file mode 100644 index 00000000..27d1b436 --- /dev/null +++ b/drivers/computeParabolic3PhasesKKS.cc @@ -0,0 +1,92 @@ +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" + +#include "SAMRAI/tbox/SAMRAIManager.h" +#include "SAMRAI/tbox/InputManager.h" +#include "SAMRAI/tbox/Database.h" +#include "SAMRAI/SAMRAI_config.h" +#include "SAMRAI/tbox/SAMRAI_MPI.h" + +#include +#include +#include + +namespace pt = boost::property_tree; +using namespace SAMRAI; + +int main(int argc, char* argv[]) +{ + tbox::SAMRAI_MPI::init(&argc, &argv); + tbox::SAMRAIManager::initialize(); + tbox::SAMRAIManager::startup(); + + int ret = 0; + { + std::string databasename(argv[1]); + double temperature = atof(argv[2]); + double phiL = atof(argv[3]); + double phiA = atof(argv[4]); + double phiB = atof(argv[5]); + double c = atof(argv[6]); + + std::cout << "Temperature = " << temperature << std::endl; + + std::shared_ptr input_db( + new tbox::MemoryDatabase("db")); + std::cout << "Filename = " << databasename << std::endl; + tbox::InputManager::getManager()->parseInputFile(databasename, input_db); + + double coeffL[3][2]; + std::shared_ptr liquid_db = + input_db->getDatabase("Liquid"); + coeffL[0][0] = liquid_db->getDouble("a0"); + coeffL[0][1] = liquid_db->getDouble("a1"); + coeffL[1][0] = liquid_db->getDouble("b0"); + coeffL[1][1] = liquid_db->getDouble("b1"); + coeffL[2][0] = liquid_db->getDouble("c0"); + coeffL[2][1] = liquid_db->getDouble("c1"); + + double coeffA[3][2]; + std::shared_ptr phasea_db = + input_db->getDatabase("PhaseA"); + coeffA[0][0] = phasea_db->getDouble("a0"); + coeffA[0][1] = phasea_db->getDouble("a1"); + coeffA[1][0] = phasea_db->getDouble("b0"); + coeffA[1][1] = phasea_db->getDouble("b1"); + coeffA[2][0] = phasea_db->getDouble("c0"); + coeffA[2][1] = phasea_db->getDouble("c1"); + + double coeffB[3][2]; + std::shared_ptr phaseb_db = + input_db->getDatabase("PhaseB"); + coeffB[0][0] = phaseb_db->getDouble("a0"); + coeffB[0][1] = phaseb_db->getDouble("a1"); + coeffB[1][0] = phaseb_db->getDouble("b0"); + coeffB[1][1] = phaseb_db->getDouble("b1"); + coeffB[2][0] = phaseb_db->getDouble("c0"); + coeffB[2][1] = phaseb_db->getDouble("c1"); + + double Tref = input_db->getDouble("Tref"); + + input_db->printClassData(std::cout); + + Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase + parabolic_fenergy(Tref, coeffL, coeffA, coeffB, + Thermo4PFM::EnergyInterpolationType::PBG, + Thermo4PFM::ConcInterpolationType::LINEAR); + + double phi[3] = {phiL, phiA, phiB}; + double conc[3] = {c, c, c}; + parabolic_fenergy.computePhaseConcentrations(temperature, &c, phi, conc); + + std::cout << "phi = (" << phiL << "," << phiA << "," << phiB << ")" + << ", c=" << c << std::endl; + std::cout << "KKS solution: cl = " << conc[0] << ", ca=" << conc[1] + << " and cb = " << conc[2] << std::endl; + } + + tbox::SAMRAIManager::shutdown(); + tbox::SAMRAIManager::finalize(); + tbox::SAMRAI_MPI::finalize(); + + return ret; +} diff --git a/drivers/plotEnergyVsCompositionParabolicThreePhase.cc b/drivers/plotEnergyVsCompositionParabolicThreePhase.cc new file mode 100644 index 00000000..9a5a12ef --- /dev/null +++ b/drivers/plotEnergyVsCompositionParabolicThreePhase.cc @@ -0,0 +1,99 @@ +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" + +#include "SAMRAI/tbox/SAMRAIManager.h" +#include "SAMRAI/tbox/InputManager.h" +#include "SAMRAI/tbox/Database.h" +#include "SAMRAI/SAMRAI_config.h" +#include "SAMRAI/tbox/SAMRAI_MPI.h" + +#include +#include +#include + +namespace pt = boost::property_tree; +using namespace SAMRAI; + +int main(int argc, char* argv[]) +{ + tbox::SAMRAI_MPI::init(&argc, &argv); + tbox::SAMRAIManager::initialize(); + tbox::SAMRAIManager::startup(); + + int ret = 0; + { + std::string databasename(argv[1]); + double temperature = atof(argv[2]); + std::cout << "Temperature = " << temperature << std::endl; + double cmin = 0.; + double cmax = 1.; + if (argc > 3) { + cmin = atof(argv[3]); + cmax = atof(argv[4]); + } + + Thermo4PFM::EnergyInterpolationType energy_interp_func_type = + Thermo4PFM::EnergyInterpolationType::PBG; + Thermo4PFM::ConcInterpolationType conc_interp_func_type = + Thermo4PFM::ConcInterpolationType::LINEAR; + + std::shared_ptr input_db( + new tbox::MemoryDatabase("db")); + std::cout << "Filename = " << databasename << std::endl; + tbox::InputManager::getManager()->parseInputFile(databasename, input_db); + + double coeffL[3][2]; + std::shared_ptr liquid_db = + input_db->getDatabase("Liquid"); + coeffL[0][0] = liquid_db->getDouble("a0"); + coeffL[0][1] = liquid_db->getDouble("a1"); + coeffL[1][0] = liquid_db->getDouble("b0"); + coeffL[1][1] = liquid_db->getDouble("b1"); + coeffL[2][0] = liquid_db->getDouble("c0"); + coeffL[2][1] = liquid_db->getDouble("c1"); + + double coeffA[3][2]; + std::shared_ptr phasea_db = + input_db->getDatabase("PhaseA"); + coeffA[0][0] = phasea_db->getDouble("a0"); + coeffA[0][1] = phasea_db->getDouble("a1"); + coeffA[1][0] = phasea_db->getDouble("b0"); + coeffA[1][1] = phasea_db->getDouble("b1"); + coeffA[2][0] = phasea_db->getDouble("c0"); + coeffA[2][1] = phasea_db->getDouble("c1"); + + double coeffB[3][2]; + std::shared_ptr phaseb_db = + input_db->getDatabase("PhaseB"); + coeffB[0][0] = phaseb_db->getDouble("a0"); + coeffB[0][1] = phaseb_db->getDouble("a1"); + coeffB[1][0] = phaseb_db->getDouble("b0"); + coeffB[1][1] = phaseb_db->getDouble("b1"); + coeffB[2][0] = phaseb_db->getDouble("c0"); + coeffB[2][1] = phaseb_db->getDouble("c1"); + + double Tref = input_db->getDouble("Tref"); + + input_db->printClassData(std::cout); + + const int npts = 100; + std::cout << " Compute energies..." << std::endl; + + Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase qfe( + Tref, coeffL, coeffA, coeffB, energy_interp_func_type, + conc_interp_func_type); + + std::stringstream ss; + ss << std::fixed; + ss << std::setprecision(2); + ss << temperature; + std::ofstream os("ParabolicFvsC" + ss.str() + ".csv", std::ios::out); + + qfe.printEnergyVsComposition(temperature, os, cmin, cmax, npts); + } + + tbox::SAMRAIManager::shutdown(); + tbox::SAMRAIManager::finalize(); + tbox::SAMRAI_MPI::finalize(); + + return ret; +} diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 7fb567cf..1124566c 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -72,6 +72,7 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderBinary.cc ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderBinaryThreePhase.cc ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderTernaryThreePhase.cc + ${CMAKE_SOURCE_DIR}/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsStrategyMultiOrder.cc @@ -82,6 +83,7 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumThreePhasesTernaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumPhaseConcentrationsBinary.cc ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumPhaseConcentrationsBinary.cc + ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryThreePhase.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryThreePhaseStochioB.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.cc diff --git a/source/EBSCompositionRHSStrategy.cc b/source/EBSCompositionRHSStrategy.cc index 2232f8d8..c9103e9b 100644 --- a/source/EBSCompositionRHSStrategy.cc +++ b/source/EBSCompositionRHSStrategy.cc @@ -137,7 +137,6 @@ EBSCompositionRHSStrategy::EBSCompositionRHSStrategy( { assert(diffusion_l_id >= 0); assert(temperature_scratch_id >= 0); - assert(d_free_energy_strategy); tbox::plog << "EBSCompositionRHSStrategy" << std::endl; diff --git a/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc b/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc index 2a6613ac..ae2c959f 100644 --- a/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc +++ b/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc @@ -149,6 +149,7 @@ int EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases:: double t = ptr_temp[idx_temp]; // solve KKS equations + assert(!std::isnan(ptr_c_l[idx_c_i])); double x[3] = {ptr_c_l[idx_c_i], ptr_c_a[idx_c_i], ptr_c_b[idx_c_i]}; int ret = computePhaseConcentrations(t, &c, hphi, x); @@ -156,12 +157,17 @@ int EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases:: std::cerr << "EquilibriumPhaseConcentrationsBinaryMultiOrderThre" "ePhases" << std::endl; - std::cerr << "computePhaseConcentrations failed for T=" << t - << ", " - << "hphi = " << hphi[0] << "," << hphi[1] << "," - << hphi[2] << std::endl; + std::cerr << "computePhaseConcentrations failed for T = " << t + << ", c = " << c << ", hphi (L,A,B) = " << hphi[0] + << "," << hphi[1] << "," << hphi[2] << std::endl; + std::cerr << "x = " << x[0] << ", " << x[1] << ", " << x[2] + << std::endl; + std::cerr << "ptr_c_l = " << ptr_c_l[idx_c_i] << std::endl; + std::cerr << "ptr_c_a = " << ptr_c_a[idx_c_i] << std::endl; + std::cerr << "ptr_c_b = " << ptr_c_b[idx_c_i] << std::endl; MPI_Abort(mpi.getCommunicator(), -1); } + assert(!std::isnan(x[0])); ptr_c_l[idx_c_i] = x[0]; ptr_c_a[idx_c_i] = x[1]; ptr_c_b[idx_c_i] = x[2]; diff --git a/source/FreeEnergyStrategyFactory.h b/source/FreeEnergyStrategyFactory.h index 92bf5d0f..f42f7f18 100644 --- a/source/FreeEnergyStrategyFactory.h +++ b/source/FreeEnergyStrategyFactory.h @@ -25,6 +25,7 @@ #include "QuadraticFreeEnergyBinary.h" #include "QuadraticFreeEnergyMultiOrderBinary.h" #include "QuadraticFreeEnergyMultiOrderTernaryThreePhase.h" +#include "ParabolicFreeEnergyMultiOrderBinaryThreePhase.h" #include "KKSdiluteBinary.h" #include "BiasDoubleWellBeckermannFreeEnergyStrategy.h" #include "BiasDoubleWellUTRCFreeEnergyStrategy.h" @@ -205,10 +206,28 @@ class FreeEnergyStrategyFactory tbox::plog << "Using Parabolic model for concentration" << std::endl; tbox::plog << "ParabolicFreeEnergyBinary" << std::endl; - free_energy_strategy.reset(new ParabolicFreeEnergyBinary( - model_parameters.energy_interp_func_type(), - model_parameters.conc_interp_func_type(), mvstrategy, - conc_l_scratch_id, conc_a_scratch_id, conc_db)); + if (model_parameters.norderp() > 1) { + if (conc_b_scratch_id > -1) { + tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase." + ".." + << std::endl; + free_energy_strategy.reset( + new ParabolicFreeEnergyMultiOrderBinaryThreePhase( + conc_db->getDatabase("Parabolic"), + model_parameters.energy_interp_func_type(), + model_parameters.norderpA(), + model_parameters.molar_volume_liquid(), + model_parameters.molar_volume_solid_A(), + model_parameters.molar_volume_solid_B(), + conc_l_scratch_id, conc_a_scratch_id, + conc_b_scratch_id)); + } + } else { + free_energy_strategy.reset(new ParabolicFreeEnergyBinary( + model_parameters.energy_interp_func_type(), + model_parameters.conc_interp_func_type(), mvstrategy, + conc_l_scratch_id, conc_a_scratch_id, conc_db)); + } } else if (model_parameters.isConcentrationModelQuadratic()) { tbox::plog << "Using Quadratic model for concentration" << std::endl; diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc new file mode 100644 index 00000000..e7c2b1f1 --- /dev/null +++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc @@ -0,0 +1,60 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +// +#include "ParabolicEquilibriumThreePhasesBinaryMultiOrder.h" +#include "FuncFort.h" + + +ParabolicEquilibriumThreePhasesBinaryMultiOrder:: + ParabolicEquilibriumThreePhasesBinaryMultiOrder( + const short norderp_A, const int conc_l_id, const int conc_a_id, + const int conc_b_id, const QuatModelParameters& model_parameters, + std::shared_ptr conc_db) + : EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases( + norderp_A, conc_l_id, conc_a_id, conc_b_id, model_parameters, conc_db) +{ + std::shared_ptr input_db = conc_db->getDatabase("Parabolic"); + + double coeffL[3][2]; + std::shared_ptr liquid_db = input_db->getDatabase("Liquid"); + coeffL[0][0] = liquid_db->getDouble("a0"); + coeffL[0][1] = liquid_db->getDouble("a1"); + coeffL[1][0] = liquid_db->getDouble("b0"); + coeffL[1][1] = liquid_db->getDouble("b1"); + coeffL[2][0] = liquid_db->getDouble("c0"); + coeffL[2][1] = liquid_db->getDouble("c1"); + + double coeffA[3][2]; + std::shared_ptr phasea_db = input_db->getDatabase("PhaseA"); + coeffA[0][0] = phasea_db->getDouble("a0"); + coeffA[0][1] = phasea_db->getDouble("a1"); + coeffA[1][0] = phasea_db->getDouble("b0"); + coeffA[1][1] = phasea_db->getDouble("b1"); + coeffA[2][0] = phasea_db->getDouble("c0"); + coeffA[2][1] = phasea_db->getDouble("c1"); + + double coeffB[3][2]; + std::shared_ptr phaseb_db = input_db->getDatabase("PhaseB"); + coeffB[0][0] = phaseb_db->getDouble("a0"); + coeffB[0][1] = phaseb_db->getDouble("a1"); + coeffB[1][0] = phaseb_db->getDouble("b0"); + coeffB[1][1] = phaseb_db->getDouble("b1"); + coeffB[2][0] = phaseb_db->getDouble("c0"); + coeffB[2][1] = phaseb_db->getDouble("c1"); + + double Tref = input_db->getDouble("Tref"); + + d_fenergy.reset(new Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase( + Tref, coeffL, coeffA, coeffB, model_parameters.energy_interp_func_type(), + Thermo4PFM::ConcInterpolationType::LINEAR)); +} + +ParabolicEquilibriumThreePhasesBinaryMultiOrder:: + ~ParabolicEquilibriumThreePhasesBinaryMultiOrder(){}; diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.h b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.h new file mode 100644 index 00000000..fbced498 --- /dev/null +++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.h @@ -0,0 +1,40 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +// +#ifndef included_ParabolicEquilibriumThreePhasesBinaryMultiOrder +#define included_ParabolicEquilibriumThreePhasesBinaryMultiOrder + +#include "EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.h" +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" + +class ParabolicEquilibriumThreePhasesBinaryMultiOrder + : public EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases +{ + public: + ParabolicEquilibriumThreePhasesBinaryMultiOrder( + const short norderp_A, const int conc_l_id, const int conc_a_id, + const int conc_b_id, const QuatModelParameters& model_parameters, + std::shared_ptr conc_db); + + ~ParabolicEquilibriumThreePhasesBinaryMultiOrder(); + + protected: + virtual int computePhaseConcentrations(const double t, double* c, + double* hphi, double* x) + { + return d_fenergy->computePhaseConcentrations(t, c, hphi, x); + } + + private: + std::shared_ptr + d_fenergy; +}; + +#endif diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc new file mode 100644 index 00000000..efeffcf6 --- /dev/null +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc @@ -0,0 +1,455 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +// +#include "ParabolicFreeEnergyMultiOrderBinaryThreePhase.h" + +#include "SAMRAI/tbox/InputManager.h" +#include "SAMRAI/pdat/CellData.h" +#include "SAMRAI/hier/Index.h" +#include "SAMRAI/math/HierarchyCellDataOpsReal.h" + +using namespace SAMRAI; + +#include + +//======================================================================= + +ParabolicFreeEnergyMultiOrderBinaryThreePhase:: + ParabolicFreeEnergyMultiOrderBinaryThreePhase( + std::shared_ptr input_db, + const Thermo4PFM::EnergyInterpolationType energy_interp_func_type, + const short norderp_A, const double vml, const double vma, + const double vmb, const int conc_l_id, const int conc_a_id, + const int conc_b_id) + : FreeEnergyStrategyThreePhase(input_db, vml, vma, vmb, conc_l_id, + conc_a_id, conc_b_id), + d_norderp_A(norderp_A) +{ + tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase..." + << std::endl; + + double coeffL[3][2]; + std::shared_ptr liquid_db = input_db->getDatabase("Liquid"); + coeffL[0][0] = liquid_db->getDouble("a0"); + coeffL[0][1] = liquid_db->getDouble("a1"); + coeffL[1][0] = liquid_db->getDouble("b0"); + coeffL[1][1] = liquid_db->getDouble("b1"); + coeffL[2][0] = liquid_db->getDouble("c0"); + coeffL[2][1] = liquid_db->getDouble("c1"); + + double coeffA[3][2]; + std::shared_ptr phasea_db = input_db->getDatabase("PhaseA"); + coeffA[0][0] = phasea_db->getDouble("a0"); + coeffA[0][1] = phasea_db->getDouble("a1"); + coeffA[1][0] = phasea_db->getDouble("b0"); + coeffA[1][1] = phasea_db->getDouble("b1"); + coeffA[2][0] = phasea_db->getDouble("c0"); + coeffA[2][1] = phasea_db->getDouble("c1"); + + double coeffB[3][2]; + std::shared_ptr phaseb_db = input_db->getDatabase("PhaseB"); + coeffB[0][0] = phaseb_db->getDouble("a0"); + coeffB[0][1] = phaseb_db->getDouble("a1"); + coeffB[1][0] = phaseb_db->getDouble("b0"); + coeffB[1][1] = phaseb_db->getDouble("b1"); + coeffB[2][0] = phaseb_db->getDouble("c0"); + coeffB[2][1] = phaseb_db->getDouble("c1"); + + double Tref = input_db->getDouble("Tref"); + + d_parabolic_fenergy.reset( + new Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase( + Tref, coeffL, coeffA, coeffB, energy_interp_func_type, + Thermo4PFM::ConcInterpolationType::LINEAR)); + + // conversion factor from [J/mol] to [pJ/(mu m)^3] + // vm^-1 [mol/m^3] * 10e-18 [m^3/(mu m^3)] * 10e12 [pJ/J] + // d_jpmol2pjpmumcube = 1.e-6 / d_vm; + + // R = 8.314472 J · K-1 · mol-1 + // tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase:" << + // std::endl; tbox::plog << "Molar volume L =" << vml << std::endl; + // tbox::plog << "Molar volume A =" << vma << std::endl; + // tbox::plog << "jpmol2pjpmumcube=" << d_jpmol2pjpmumcube << std::endl; +} + +ParabolicFreeEnergyMultiOrderBinaryThreePhase:: + ~ParabolicFreeEnergyMultiOrderBinaryThreePhase(){}; + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhase::computeFreeEnergy( + const hier::Box& pbox, std::shared_ptr > cd_temp, + std::shared_ptr > cd_free_energy, + std::shared_ptr > cd_conc_i, + Thermo4PFM::PhaseIndex pi, const double energy_factor) +{ + assert(cd_conc_i->getDepth() == 1); + assert(energy_factor > 0.); + assert(d_parabolic_fenergy); + + double* ptr_temp = cd_temp->getPointer(); + double* ptr_f = cd_free_energy->getPointer(); + double* ptr_c_i = cd_conc_i->getPointer(); + + const hier::Box& temp_gbox = cd_temp->getGhostBox(); + int imin_temp = temp_gbox.lower(0); + int jmin_temp = temp_gbox.lower(1); + int jp_temp = temp_gbox.numberCells(0); + int kmin_temp = 0; + int kp_temp = 0; +#if (NDIM == 3) + kmin_temp = temp_gbox.lower(2); + kp_temp = jp_temp * temp_gbox.numberCells(1); +#endif + + const hier::Box& f_gbox = cd_free_energy->getGhostBox(); + int imin_f = f_gbox.lower(0); + int jmin_f = f_gbox.lower(1); + int jp_f = f_gbox.numberCells(0); + int kmin_f = 0; + int kp_f = 0; +#if (NDIM == 3) + kmin_f = f_gbox.lower(2); + kp_f = jp_f * f_gbox.numberCells(1); +#endif + + const hier::Box& c_i_gbox = cd_conc_i->getGhostBox(); + int imin_c_i = c_i_gbox.lower(0); + int jmin_c_i = c_i_gbox.lower(1); + int jp_c_i = c_i_gbox.numberCells(0); + int kmin_c_i = 0; + int kp_c_i = 0; +#if (NDIM == 3) + kmin_c_i = c_i_gbox.lower(2); + kp_c_i = jp_c_i * c_i_gbox.numberCells(1); +#endif + + int imin = pbox.lower(0); + int imax = pbox.upper(0); + int jmin = pbox.lower(1); + int jmax = pbox.upper(1); + int kmin = 0; + int kmax = 0; +#if (NDIM == 3) + kmin = pbox.lower(2); + kmax = pbox.upper(2); +#endif + + for (int kk = kmin; kk <= kmax; kk++) { + for (int jj = jmin; jj <= jmax; jj++) { + for (int ii = imin; ii <= imax; ii++) { + + const int idx_temp = (ii - imin_temp) + (jj - jmin_temp) * jp_temp + + (kk - kmin_temp) * kp_temp; + + const int idx_f = + (ii - imin_f) + (jj - jmin_f) * jp_f + (kk - kmin_f) * kp_f; + + const int idx_c_i = (ii - imin_c_i) + (jj - jmin_c_i) * jp_c_i + + (kk - kmin_c_i) * kp_c_i; + + double t = ptr_temp[idx_temp]; + double c_i = ptr_c_i[idx_c_i]; + + ptr_f[idx_f] = d_parabolic_fenergy->computeFreeEnergy(t, &c_i, pi); + ptr_f[idx_f] *= energy_factor; + assert(!std::isnan(ptr_f[idx_f])); + } + } + } +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhase::addDrivingForceOnPatch( + std::shared_ptr > cd_rhs, + std::shared_ptr > cd_temperature, + std::shared_ptr > cd_phi, + std::shared_ptr > cd_f_l, + std::shared_ptr > cd_f_a, + std::shared_ptr > cd_f_b, + std::shared_ptr > cd_c_l, + std::shared_ptr > cd_c_a, + std::shared_ptr > cd_c_b, const hier::Box& pbox) +{ + assert(cd_temperature); + assert(cd_f_b); + assert(cd_c_b); + assert(cd_f_l->getGhostCellWidth()[0] == cd_f_a->getGhostCellWidth()[0]); + assert(cd_f_l->getGhostCellWidth()[0] == cd_f_b->getGhostCellWidth()[0]); + assert(cd_c_l->getGhostCellWidth()[0] == cd_c_a->getGhostCellWidth()[0]); + assert(cd_c_l->getGhostCellWidth()[0] == cd_c_b->getGhostCellWidth()[0]); + assert(cd_phi->getDepth() > 1); + assert(cd_phi->getDepth() == cd_rhs->getDepth()); + + const int norderp = cd_phi->getDepth(); + + double* ptr_temp = cd_temperature->getPointer(); + double* ptr_f_l = cd_f_l->getPointer(); + double* ptr_f_a = cd_f_a->getPointer(); + double* ptr_f_b = cd_f_b->getPointer(); + double* ptr_c_l = cd_c_l->getPointer(); + double* ptr_c_a = cd_c_a->getPointer(); + double* ptr_c_b = cd_c_b->getPointer(); + + const hier::Box& rhs_gbox = cd_rhs->getGhostBox(); + int imin_rhs = rhs_gbox.lower(0); + int jmin_rhs = rhs_gbox.lower(1); + int jp_rhs = rhs_gbox.numberCells(0); + int kmin_rhs = 0; + int kp_rhs = 0; +#if (NDIM == 3) + kmin_rhs = rhs_gbox.lower(2); + kp_rhs = jp_rhs * rhs_gbox.numberCells(1); +#endif + + const hier::Box& temp_gbox = cd_temperature->getGhostBox(); + int imin_temp = temp_gbox.lower(0); + int jmin_temp = temp_gbox.lower(1); + int jp_temp = temp_gbox.numberCells(0); + int kmin_temp = 0; + int kp_temp = 0; +#if (NDIM == 3) + kmin_temp = temp_gbox.lower(2); + kp_temp = jp_temp * temp_gbox.numberCells(1); +#endif + + const hier::Box& pf_gbox = cd_phi->getGhostBox(); + int imin_pf = pf_gbox.lower(0); + int jmin_pf = pf_gbox.lower(1); + int jp_pf = pf_gbox.numberCells(0); + int kmin_pf = 0; + int kp_pf = 0; +#if (NDIM == 3) + kmin_pf = pf_gbox.lower(2); + kp_pf = jp_pf * pf_gbox.numberCells(1); +#endif + + // Assuming f_l, f_a, all have same ghost box + const hier::Box& f_i_gbox = cd_f_l->getGhostBox(); + int imin_f_i = f_i_gbox.lower(0); + int jmin_f_i = f_i_gbox.lower(1); + int jp_f_i = f_i_gbox.numberCells(0); + int kmin_f_i = 0; + int kp_f_i = 0; +#if (NDIM == 3) + kmin_f_i = f_i_gbox.lower(2); + kp_f_i = jp_f_i * f_i_gbox.numberCells(1); +#endif + + // Assuming c_l, c_a, all have same ghost box + const hier::Box& c_i_gbox = cd_c_l->getGhostBox(); + int imin_c_i = c_i_gbox.lower(0); + int jmin_c_i = c_i_gbox.lower(1); + int jp_c_i = c_i_gbox.numberCells(0); + int kmin_c_i = 0; + int kp_c_i = 0; +#if (NDIM == 3) + kmin_c_i = c_i_gbox.lower(2); + kp_c_i = jp_c_i * c_i_gbox.numberCells(1); +#endif + + int imin = pbox.lower(0); + int imax = pbox.upper(0); + int jmin = pbox.lower(1); + int jmax = pbox.upper(1); + int kmin = 0; + int kmax = 0; +#if (NDIM == 3) + kmin = pbox.lower(2); + kmax = pbox.upper(2); +#endif + + std::vector rhs(norderp); + + std::vector ptr_rhs(norderp); + for (short i = 0; i < norderp; i++) + ptr_rhs[i] = cd_rhs->getPointer(i); + + std::vector ptr_phi(norderp); + for (short i = 0; i < norderp; i++) + ptr_phi[i] = cd_phi->getPointer(i); + + // tbox::plog<<"d_norderp_A = "< 0.); + const double sum2inv = 1. / sum2; + + hphil *= sum2inv; + hphiA *= sum2inv; + hphiB *= sum2inv; + + assert(!std::isnan(hphiA)); + assert(!std::isnan(hphiB)); + + // solid phase A order parameters + for (short i = 0; i < d_norderp_A; i++) + rhs[i] = 2. * ptr_phi[i][idx_pf] * + ((1. - hphiA) * dfa - hphil * dfl - hphiB * dfb) * + sum2inv; + // solid phase B order parameters + for (short i = d_norderp_A; i < norderp - 1; i++) + rhs[i] = 2. * ptr_phi[i][idx_pf] * + ((1. - hphiB) * dfb - hphil * dfl - hphiA * dfa) * + sum2inv; + // liquid phase order parameter + rhs[norderp - 1] = + 2. * ptr_phi[norderp - 1][idx_pf] * + ((1. - hphil) * dfl - hphiA * dfa - hphiB * dfb) * sum2inv; + for (short i = 0; i < norderp; i++) + assert(!std::isnan(rhs[i])); + + for (short i = 0; i < norderp; i++) + ptr_rhs[i][idx_rhs] -= (rhs[i]); + } + } + } +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhase::computeMuL( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseL, + &mu); + mu *= d_energy_conv_factor_L; + return mu; +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhase::computeMuA( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseA, + &mu); + mu *= d_energy_conv_factor_A; + return mu; +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhase::computeMuB( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseB, + &mu); + mu *= d_energy_conv_factor_B; + return mu; +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhase:: + computeSecondDerivativeEnergyPhaseL(const std::vector& c_l, + std::vector& d2fdc2, + const bool use_internal_units) +{ + d_parabolic_fenergy->computeSecondDerivativeFreeEnergy( + 0., &c_l[0], Thermo4PFM::PhaseIndex::phaseL, &d2fdc2[0]); + if (use_internal_units) + for (short i = 0; i < 3; i++) + d2fdc2[i] *= d_energy_conv_factor_L; +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhase:: + computeSecondDerivativeEnergyPhaseA(const std::vector& c_a, + std::vector& d2fdc2, + const bool use_internal_units) +{ + d_parabolic_fenergy->computeSecondDerivativeFreeEnergy( + 0., &c_a[0], Thermo4PFM::PhaseIndex::phaseA, &d2fdc2[0]); + if (use_internal_units) + for (short i = 0; i < 3; i++) + d2fdc2[i] *= d_energy_conv_factor_A; +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhase:: + computeSecondDerivativeEnergyPhaseB(const std::vector& c_b, + std::vector& d2fdc2, + const bool use_internal_units) +{ + d_parabolic_fenergy->computeSecondDerivativeFreeEnergy( + 0., &c_b[0], Thermo4PFM::PhaseIndex::phaseB, &d2fdc2[0]); + if (use_internal_units) + for (short i = 0; i < 3; i++) + d2fdc2[i] *= d_energy_conv_factor_B; +} diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h new file mode 100644 index 00000000..c385f90e --- /dev/null +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h @@ -0,0 +1,72 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +// +#ifndef included_ParabolicFreeEnergyMultiOrderBinaryThreePhase +#define included_ParabolicFreeEnergyMultiOrderBinaryThreePhase + +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" +#include "FreeEnergyStrategyThreePhase.h" + +class ParabolicFreeEnergyMultiOrderBinaryThreePhase + : public FreeEnergyStrategyThreePhase +{ + public: + ParabolicFreeEnergyMultiOrderBinaryThreePhase( + std::shared_ptr input_db, + const Thermo4PFM::EnergyInterpolationType energy_interp_func_type, + const short norderp_A, const double vml, const double vma, + const double vmb, const int conc_l_id, const int conc_a_id, + const int conc_b_id); + + ~ParabolicFreeEnergyMultiOrderBinaryThreePhase(); + + private: + // + // number of order parameters associated with phase A + // + const short d_norderp_A; + + std::shared_ptr + d_parabolic_fenergy; + + void computeSecondDerivativeEnergyPhaseL( + const std::vector& c, std::vector& d2fdc2, + const bool use_internal_units) override; + void computeSecondDerivativeEnergyPhaseA( + const std::vector& c, std::vector& d2fdc2, + const bool use_internal_units) override; + void computeSecondDerivativeEnergyPhaseB( + const std::vector& c, std::vector& d2fdc2, + const bool use_internal_units) override; + + double computeMuL(const double t, const double c); + double computeMuA(const double t, const double c); + double computeMuB(const double t, const double c); + + void addDrivingForceOnPatch( + std::shared_ptr > cd_rhs, + std::shared_ptr > cd_temperature, + std::shared_ptr > cd_phi, + std::shared_ptr > cd_f_l, + std::shared_ptr > cd_f_a, + std::shared_ptr > cd_f_b, + std::shared_ptr > cd_c_l, + std::shared_ptr > cd_c_a, + std::shared_ptr > cd_c_b, + const hier::Box& pbox) override; + + void computeFreeEnergy( + const hier::Box& pbox, std::shared_ptr > cd_temp, + std::shared_ptr > cd_free_energy, + std::shared_ptr > cd_conc_i, + Thermo4PFM::PhaseIndex pi, const double energy_factor) override; +}; + +#endif diff --git a/source/PhaseConcentrationsStrategyFactory.h b/source/PhaseConcentrationsStrategyFactory.h index 50fce2dc..56639b78 100644 --- a/source/PhaseConcentrationsStrategyFactory.h +++ b/source/PhaseConcentrationsStrategyFactory.h @@ -18,6 +18,7 @@ #include "QuadraticEquilibriumPhaseConcentrationsBinary.h" #include "QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder.h" #include "QuadraticEquilibriumThreePhasesTernaryMultiOrder.h" +#include "ParabolicEquilibriumThreePhasesBinaryMultiOrder.h" #include "PartitionPhaseConcentrationsStrategy.h" #include "PhaseIndependentConcentrationsStrategy.h" #include "Database2JSON.h" @@ -197,44 +198,53 @@ class PhaseConcentrationsStrategyFactory conc_l_scratch_id, conc_a_scratch_id, model_parameters.energy_interp_func_type(), model_parameters.conc_interp_func_type(), conc_db)); - } else { - if (model_parameters.isConcentrationModelParabolic()) { - tbox::plog << "Parabolic..." << std::endl; - assert(conc_b_scratch_id == -1); - phase_conc_strategy.reset( - new ParabolicEquilibriumPhaseConcentrationsBinary( - conc_l_scratch_id, conc_a_scratch_id, - model_parameters.energy_interp_func_type(), - model_parameters.conc_interp_func_type(), conc_db)); - - } else { - if (model_parameters.isConcentrationModelQuadratic()) { - if (model_parameters.norderp() > 1) { - tbox::plog << "Quadratic, MultiOrder..." << std::endl; - if (conc_b_scratch_id >= 0) { - tbox::plog << "Quadratic, MultiOrder, Three phases..." - << std::endl; - phase_conc_strategy.reset( - new QuadraticEquilibriumThreePhasesTernaryMultiOrder( - model_parameters.norderpA(), conc_l_scratch_id, - conc_a_scratch_id, conc_b_scratch_id, - model_parameters, conc_db)); - } else { - phase_conc_strategy.reset( - new QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder( - conc_l_scratch_id, conc_a_scratch_id, - model_parameters, conc_db)); - } + } else { // neither CALPHAD or KKS + if (model_parameters.isConcentrationModelQuadratic()) { + if (model_parameters.norderp() > 1) { + tbox::plog << "Quadratic, MultiOrder..." << std::endl; + if (conc_b_scratch_id >= 0) { + tbox::plog << "Quadratic, MultiOrder, Three phases..." + << std::endl; + phase_conc_strategy.reset( + new QuadraticEquilibriumThreePhasesTernaryMultiOrder( + model_parameters.norderpA(), conc_l_scratch_id, + conc_a_scratch_id, conc_b_scratch_id, + model_parameters, conc_db)); } else { - tbox::plog << "Quadratic..." << std::endl; - assert(conc_b_scratch_id == -1); phase_conc_strategy.reset( - new QuadraticEquilibriumPhaseConcentrationsBinary( + new QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder( conc_l_scratch_id, conc_a_scratch_id, - model_parameters.energy_interp_func_type(), - model_parameters.conc_interp_func_type(), - conc_db)); + model_parameters, conc_db)); + } + } else { + tbox::plog << "Quadratic..." << std::endl; + assert(conc_b_scratch_id == -1); + phase_conc_strategy.reset( + new QuadraticEquilibriumPhaseConcentrationsBinary( + conc_l_scratch_id, conc_a_scratch_id, + model_parameters.energy_interp_func_type(), + model_parameters.conc_interp_func_type(), conc_db)); + } + } else if (model_parameters.isConcentrationModelParabolic()) { + if (model_parameters.norderp() > 1) { + tbox::plog << "Parabolic, MultiOrder..." << std::endl; + if (conc_b_scratch_id >= 0) { + tbox::plog << "Parabolic, MultiOrder, Three phases..." + << std::endl; + phase_conc_strategy.reset( + new ParabolicEquilibriumThreePhasesBinaryMultiOrder( + model_parameters.norderpA(), conc_l_scratch_id, + conc_a_scratch_id, conc_b_scratch_id, + model_parameters, conc_db)); } + } else { + tbox::plog << "Parabolic..." << std::endl; + assert(conc_b_scratch_id == -1); + phase_conc_strategy.reset( + new ParabolicEquilibriumPhaseConcentrationsBinary( + conc_l_scratch_id, conc_a_scratch_id, + model_parameters.energy_interp_func_type(), + model_parameters.conc_interp_func_type(), conc_db)); } } } diff --git a/source/TbasedCompositionDiffusionStrategy.cc b/source/TbasedCompositionDiffusionStrategy.cc index 659234e9..7da216ed 100644 --- a/source/TbasedCompositionDiffusionStrategy.cc +++ b/source/TbasedCompositionDiffusionStrategy.cc @@ -219,7 +219,6 @@ void TbasedCompositionDiffusionStrategy::setDiffusion( { const char interp_func_type = interpChar(); - const hier::Box& pbox = patch->getBox(); const hier::Index& ifirst = pbox.lower(); const hier::Index& ilast = pbox.upper(); @@ -237,6 +236,7 @@ void TbasedCompositionDiffusionStrategy::setDiffusion( assert(phi->getDepth() == d_norderp); { + //tbox::plog<<"d_with_phaseB, d_with3phases = "<getPointer(0) : phi->getPointer(d_norderp - 1); @@ -248,7 +248,6 @@ void TbasedCompositionDiffusionStrategy::setDiffusion( const int nphiA = d_norderpA; const int nphiB = d_norderpB; - // this call assumes the order phiA, phiB, phiL CONCENTRATION_PFMDIFFUSION_OF_TEMPERATURE_THREEPHASES( ifirst(0), ilast(0), ifirst(1), ilast(1), #if (NDIM == 3)