From c6e525d08ad5b18e2f23c3aacd28fa5db93bb175 Mon Sep 17 00:00:00 2001 From: "Fattebert J.-L." Date: Thu, 17 Apr 2025 16:15:21 -0400 Subject: [PATCH 1/6] New solver StochioAB --- ...rgyBinaryMultiOrderThreePhasesStochioAB.cc | 71 +++++ ...ergyBinaryMultiOrderThreePhasesStochioAB.h | 56 ++++ ...ntrationsMultiOrderThreePhasesStochioAB.cc | 100 +++++++ ...entrationsMultiOrderThreePhasesStochioAB.h | 37 +++ source/CMakeLists.txt | 5 + source/FreeEnergyStrategyFactory.h | 38 ++- source/MultiOrderBinaryDrivingForce.h | 4 +- .../MultiOrderBinaryThreePhasesDrivingForce.h | 4 +- ...rBinaryThreePhasesDrivingForceStochioAB.cc | 278 ++++++++++++++++++ ...erBinaryThreePhasesDrivingForceStochioAB.h | 41 +++ ...riumPhaseConcentrationsBinaryMultiOrder.cc | 5 +- ...briumPhaseConcentrationsBinaryMultiOrder.h | 2 - ...iumThreePhasesBinaryMultiOrderStochioAB.cc | 83 ++++++ ...riumThreePhasesBinaryMultiOrderStochioAB.h | 36 +++ ...licFreeEnergyMultiOrderBinaryThreePhase.cc | 11 +- ...olicFreeEnergyMultiOrderBinaryThreePhase.h | 7 +- ...ergyMultiOrderBinaryThreePhaseStochioAB.cc | 218 ++++++++++++++ ...nergyMultiOrderBinaryThreePhaseStochioAB.h | 73 +++++ source/PhaseConcentrationsStrategyFactory.h | 46 ++- source/QuatModelParameters.cc | 2 + source/QuatModelParameters.h | 4 + 21 files changed, 1076 insertions(+), 45 deletions(-) create mode 100644 source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc create mode 100644 source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h create mode 100644 source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc create mode 100644 source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h create mode 100644 source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc create mode 100644 source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h create mode 100644 source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc create mode 100644 source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h create mode 100644 source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc create mode 100644 source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc new file mode 100644 index 00000000..09fa305d --- /dev/null +++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc @@ -0,0 +1,71 @@ +// 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 "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h" + +#include +#include "Database2JSON.h" +namespace pt = boost::property_tree; + +#include "SAMRAI/tbox/InputManager.h" +#include "SAMRAI/pdat/CellData.h" +#include "SAMRAI/hier/Index.h" +#include "SAMRAI/math/HierarchyCellDataOpsReal.h" + +using namespace SAMRAI; + +#include + +//======================================================================= + +CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB:: + CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB( + const int norderp_A, boost::property_tree::ptree calphad_pt, + std::shared_ptr newton_db, + const Thermo4PFM::ConcInterpolationType conc_interp_func_type, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id) + : CALPHADFreeEnergyBinaryMultiOrderThreePhases< + Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase>( + calphad_pt, newton_db, conc_interp_func_type, norderp_A, mvstrategy, + conc_l_id, conc_a_id, conc_b_id) +{ + tbox::plog << "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB..." + << std::endl; + setup(calphad_pt, newton_db); + + d_multiorder_driving_force.reset( + new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, norderp_A)); +} + +//======================================================================= + +void CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::setup( + pt::ptree calphad_pt, std::shared_ptr newton_db) +{ + pt::ptree newton_pt; + if (newton_db) copyDatabase(newton_db, newton_pt); + // set looser tol in solver + newton_pt.put("tol", 1.e-6); + d_ceq_fenergy.reset(new Thermo4PFM::CALPHADFreeEnergyFunctionsBinary( + calphad_pt, newton_pt, d_energy_interp_func_type, + d_conc_interp_func_type)); +} + + +bool CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::computeCeqT( + const double temperature, const Thermo4PFM::PhaseIndex pi0, + const Thermo4PFM::PhaseIndex pi1, double* ceq) +{ + std::cout << "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::" + "computeCeqT..." + << std::endl; + return d_ceq_fenergy->computeCeqT(temperature, &ceq[0], 50, true); +} diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h new file mode 100644 index 00000000..2c363430 --- /dev/null +++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h @@ -0,0 +1,56 @@ +// 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_CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB +#define included_CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB + +#include "CALPHADFreeEnergyBinaryMultiOrderThreePhases.h" +#include "InterpolationType.h" +#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h" + +#include "CALPHADFreeEnergyFunctionsBinaryThreePhase.h" +#include "CALPHADFreeEnergyFunctionsBinary.h" + +#include "SAMRAI/pdat/CellData.h" +#include "SAMRAI/tbox/Database.h" +#include "SAMRAI/hier/Box.h" +class MolarVolumeStrategy; + +#include +#include +#include + +class CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB + : public CALPHADFreeEnergyBinaryMultiOrderThreePhases< + Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase> +{ + public: + CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB( + const int norderp_A, boost::property_tree::ptree calphad_db, + std::shared_ptr newton_db, + const Thermo4PFM::ConcInterpolationType conc_interp_func_type, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id); + + bool computeCeqT(const double temperature, const Thermo4PFM::PhaseIndex pi0, + const Thermo4PFM::PhaseIndex pi1, double* ceq) override; + + private: + void setup(boost::property_tree::ptree calphad_pt, + std::shared_ptr newton_db); + + std::shared_ptr d_ceq_fenergy; + + + std::shared_ptr + d_multiorder_driving_force; +}; + +#endif diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc new file mode 100644 index 00000000..7517fba9 --- /dev/null +++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc @@ -0,0 +1,100 @@ +// 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 "Database2JSON.h" + +#include +namespace pt = boost::property_tree; + +#include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h" + +CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: + CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB( + 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), + d_model_parameters(model_parameters) +{ + tbox::plog << "CALPHADequilPhaseConcMultiOrderThreePhasesStochioAB..." + << std::endl; + std::shared_ptr conc_calphad_db = + conc_db->getDatabase("Calphad"); + std::string calphad_filename = conc_calphad_db->getString("filename"); + + std::shared_ptr calphad_db; + boost::property_tree::ptree calphad_pt; + + if (calphad_filename.compare(calphad_filename.size() - 4, 4, "json") == 0) { + boost::property_tree::read_json(calphad_filename, calphad_pt); + } else { + calphad_db.reset(new tbox::MemoryDatabase("calphad_db")); + tbox::InputManager::getManager()->parseInputFile(calphad_filename, + calphad_db); + copyDatabase(calphad_db, calphad_pt); + } +} + +int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: + computePhaseConcentrations(const double temp, double* c, double* hphi, + double* x) +{ + assert(!std::isnan(x[0])); + + const double epsilon1 = 1.e-4; + const double epsilon2 = 2.e-4; + + // initialize to NaN to trigger error if used when not set + double xkks = tbox::IEEE::getSignalingNaN(); + double xeq = tbox::IEEE::getSignalingNaN(); + + const double cA = d_model_parameters.getStochioA(); + const double cB = d_model_parameters.getStochioB(); + + if (hphi[0] >= epsilon1) { + // solve explicit KKS problem + xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / hphi[0]; + +#if 0 + for (short i = 0; i < 3; i++) + std::cerr << hphi[i] << ", "; + std::cerr << ", c=" << c[0] << std::endl; + std::cerr << "x = " << x[0] << std::endl; + std::cerr << "xkks = " << xkks << std::endl; + std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = " + << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl; +#endif + } + + if (hphi[0] < epsilon2) { + xeq = d_model_parameters.ceq_liquid(temp); + // std::cout << "xeq = " << xeq << std::endl; + assert(!std::isnan(xeq)); + } + + if (hphi[0] >= epsilon2) { + x[0] = xkks; + } else if (hphi[0] <= epsilon1) { + x[0] = xeq; + } else { + // map epsilon1 < hphi < epsilon2 to (0,1) + double h = (hphi[0] - epsilon1) / (epsilon2 - epsilon1); + // mix equilibrium compositions and KKS solution + double f = + Thermo4PFM::interp_func(Thermo4PFM::EnergyInterpolationType::PBG, h); + x[0] = (1. - f) * xeq + f * xkks; + } + x[1] = cA; + x[2] = cB; + + return 0; +} diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h new file mode 100644 index 00000000..83a2649a --- /dev/null +++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h @@ -0,0 +1,37 @@ +// 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_CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB +#define included_CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB + +#include "EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.h" +#include "InterpolationType.h" +#include "QuatModelParameters.h" + +#include "SAMRAI/tbox/InputManager.h" + +class CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB + : public EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases +{ + public: + CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB( + 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); + + ~CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB() {} + + private: + QuatModelParameters d_model_parameters; + + int computePhaseConcentrations(double, double*, double*, double*) override; +}; + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 145ec860..cd468efb 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -30,6 +30,7 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryDrivingForce.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForce.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForceStochioB.cc + ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc ${CMAKE_SOURCE_DIR}/source/CahnHilliardDoubleWell.cc ${CMAKE_SOURCE_DIR}/source/SinteringUWangRHSStrategy.cc ${CMAKE_SOURCE_DIR}/source/WangSinteringCompositionRHSStrategy.cc @@ -79,12 +80,14 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderBinaryThreePhase.cc ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderTernaryThreePhase.cc ${CMAKE_SOURCE_DIR}/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc + ${CMAKE_SOURCE_DIR}/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsStrategyMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsStrategyMultiOrderThreePhases.cc ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsThreePhasesStochioB.cc ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB.cc + ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumThreePhasesTernaryMultiOrder.cc ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumPhaseConcentrationsBinary.cc @@ -93,7 +96,9 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryFolchPlapp.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryFolchPlappStochioB.cc ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc + ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.cc + ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyTernary.cc ${CMAKE_SOURCE_DIR}/source/tools.cc ${CMAKE_SOURCE_DIR}/source/FreeEnergyStrategy.cc diff --git a/source/FreeEnergyStrategyFactory.h b/source/FreeEnergyStrategyFactory.h index 03e3dd3e..0a388daf 100644 --- a/source/FreeEnergyStrategyFactory.h +++ b/source/FreeEnergyStrategyFactory.h @@ -20,6 +20,7 @@ #include "CALPHADFreeEnergyStrategyBinaryFolchPlapp.h" #include "CALPHADFreeEnergyBinaryMultiOrderThreePhases.h" #include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.h" +#include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h" #include "CALPHADFreeEnergyStrategyBinaryFolchPlappStochioB.h" #include "ParabolicFreeEnergyBinary.h" #include "ParabolicFreeEnergyMultiOrderBinary.h" @@ -27,6 +28,7 @@ #include "QuadraticFreeEnergyMultiOrderBinary.h" #include "QuadraticFreeEnergyMultiOrderTernaryThreePhase.h" #include "ParabolicFreeEnergyMultiOrderBinaryThreePhase.h" +#include "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h" #include "KKSdiluteBinary.h" #include "BiasDoubleWellBeckermannFreeEnergyStrategy.h" #include "BiasDoubleWellUTRCFreeEnergyStrategy.h" @@ -51,6 +53,8 @@ class FreeEnergyStrategyFactory std::shared_ptr free_energy_strategy; if (model_parameters.with_concentration()) { + const double cB = model_parameters.getStochioB(); + const double cA = model_parameters.getStochioA(); if (model_parameters.isConcentrationModelCALPHAD()) { std::shared_ptr calphad_db; @@ -82,7 +86,16 @@ class FreeEnergyStrategyFactory if (model_parameters.withMultipleOrderP()) { tbox::plog << "MultiOrder..." << std::endl; if (conc_b_scratch_id >= 0) { - if (model_parameters.getStochioB() >= 0.) { + if (cA >= 0. && cB >= 0.) { + tbox::plog << "StochioAB..." << std::endl; + free_energy_strategy.reset( + new CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB( + model_parameters.norderpA(), calphad_pt, + newton_db, + model_parameters.conc_interp_func_type(), + mvstrategy, conc_l_scratch_id, + conc_a_scratch_id, conc_b_scratch_id)); + } else if (model_parameters.getStochioB() >= 0.) { tbox::plog << "StochioB..." << std::endl; free_energy_strategy.reset( new CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB( @@ -208,13 +221,22 @@ class FreeEnergyStrategyFactory tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase." ".." << std::endl; - free_energy_strategy.reset( - new ParabolicFreeEnergyMultiOrderBinaryThreePhase( - conc_db->getDatabase("Parabolic"), - model_parameters.conc_interp_func_type(), - model_parameters.norderpA(), mvstrategy, - conc_l_scratch_id, conc_a_scratch_id, - conc_b_scratch_id)); + if (cA > 0. && cB > 0.) { + free_energy_strategy.reset( + new ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB( + conc_db->getDatabase("Parabolic"), + model_parameters.norderpA(), mvstrategy, + conc_l_scratch_id, conc_a_scratch_id, + conc_b_scratch_id)); + + } else { + free_energy_strategy.reset( + new ParabolicFreeEnergyMultiOrderBinaryThreePhase( + conc_db->getDatabase("Parabolic"), + model_parameters.norderpA(), mvstrategy, + conc_l_scratch_id, conc_a_scratch_id, + conc_b_scratch_id)); + } } else { tbox::plog << "ParabolicFreeEnergyMultiOrderBinary" << std::endl; diff --git a/source/MultiOrderBinaryDrivingForce.h b/source/MultiOrderBinaryDrivingForce.h index 4f5fd0ab..c1960256 100644 --- a/source/MultiOrderBinaryDrivingForce.h +++ b/source/MultiOrderBinaryDrivingForce.h @@ -15,14 +15,12 @@ #include "SAMRAI/hier/Patch.h" -using namespace SAMRAI; - class MultiOrderBinaryDrivingForce { public: MultiOrderBinaryDrivingForce(FreeEnergyStrategyBinary* fenergy_strategy); - void addDrivingForce(hier::Patch& patch, const int temperature_id, + void addDrivingForce(SAMRAI::hier::Patch& patch, const int temperature_id, const int phase_id, const int conc_l_id, const int conc_a_id, const int f_l_id, const int f_a_id, const int rhs_id); diff --git a/source/MultiOrderBinaryThreePhasesDrivingForce.h b/source/MultiOrderBinaryThreePhasesDrivingForce.h index 294f77b6..aaa02e29 100644 --- a/source/MultiOrderBinaryThreePhasesDrivingForce.h +++ b/source/MultiOrderBinaryThreePhasesDrivingForce.h @@ -15,15 +15,13 @@ #include "SAMRAI/hier/Patch.h" -using namespace SAMRAI; - class MultiOrderBinaryThreePhasesDrivingForce { public: MultiOrderBinaryThreePhasesDrivingForce( FreeEnergyStrategyBinary* fenergy_strategy, const int norderp_A); - void addDrivingForce(hier::Patch& patch, const int temperature_id, + void addDrivingForce(SAMRAI::hier::Patch& patch, const int temperature_id, const int phase_id, const int conc_l_id, const int conc_a_id, const int conc_b_id, const int f_l_id, const int f_a_id, const int f_b_id, diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc new file mode 100644 index 00000000..3f94bbab --- /dev/null +++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc @@ -0,0 +1,278 @@ +#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h" + +#include "functions.h" +#include "SAMRAI/pdat/CellData.h" + + +using namespace SAMRAI; + +MultiOrderBinaryThreePhasesDrivingForceStochioAB:: + MultiOrderBinaryThreePhasesDrivingForceStochioAB( + FreeEnergyStrategyBinary* fenergy_strategy, const int norderp_A) + : d_fenergy_strategy(fenergy_strategy), d_norderp_A(norderp_A) +{ +} + +void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( + hier::Patch& patch, const int temperature_id, const int phase_id, + const int conc_l_id, const int conc_a_id, const int conc_b_id, + const int f_l_id, const int f_a_id, const int f_b_id, const int rhs_id) +{ + assert(phase_id >= 0); + assert(f_l_id >= 0); + assert(f_a_id >= 0); + assert(f_b_id >= 0); + assert(rhs_id >= 0); + assert(conc_l_id >= 0); + assert(conc_a_id >= 0); + assert(conc_b_id >= 0); + assert(temperature_id >= 0); + + std::shared_ptr > phase( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(phase_id))); + assert(phase); + assert(phase->getDepth() > 1); + + std::shared_ptr > temperature( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(temperature_id))); + assert(temperature); + + std::shared_ptr > fl( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(f_l_id))); + assert(fl); + + std::shared_ptr > fa( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(f_a_id))); + assert(fa); + + std::shared_ptr > fb( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(f_b_id))); + assert(fb); + + std::shared_ptr > cl( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(conc_l_id))); + assert(cl); + + std::shared_ptr > ca( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(conc_a_id))); + assert(ca); + + std::shared_ptr > cb( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(conc_b_id))); + assert(cb); + + std::shared_ptr > rhs( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(rhs_id))); + + assert(fl->getGhostCellWidth()[0] == fa->getGhostCellWidth()[0]); + assert(cl->getGhostCellWidth()[0] == ca->getGhostCellWidth()[0]); + assert(phase->getDepth() > 1); + assert(phase->getDepth() == rhs->getDepth()); + + const hier::Box& pbox(patch.getBox()); + const int norderp = phase->getDepth(); + + double* ptr_temp = temperature->getPointer(); + double* ptr_f_l = fl->getPointer(); + double* ptr_f_a = fa->getPointer(); + double* ptr_f_b = fb->getPointer(); + double* ptr_c_l = cl->getPointer(); + double* ptr_c_a = ca->getPointer(); + double* ptr_c_b = cb->getPointer(); + + const hier::Box& rhs_gbox = 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 = 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 = phase->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 = fl->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 = cl->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_local(norderp); + + std::vector ptr_rhs(norderp); + for (short i = 0; i < norderp; i++) + ptr_rhs[i] = rhs->getPointer(i); + + std::vector ptr_phi(norderp); + for (short i = 0; i < norderp; i++) + ptr_phi[i] = phase->getPointer(i); + + // tbox::plog<<"d_norderp_A = "<computeMuL(t, cl); + + // + // see Moelans, Acta Mat 59 (2011) + // + + // driving forces, assuming muL = muA = muB + double dfa = (fa - muL * ca); + double dfb = (fb - muL * cb); + double dfl = (fl - muL * cl); + assert(!std::isnan(dfa)); + + // interpolation polynomials + double hphiA = 0.; + for (short i = 0; i < d_norderp_A; i++) { + const double phiA = std::max(0., ptr_phi[i][idx_pf]); + hphiA += phiA * phiA; + } + assert(!std::isnan(hphiA)); + + double hphiB = 0.; + for (short i = d_norderp_A; i < norderp - 1; i++) { + const double phiB = std::max(0., ptr_phi[i][idx_pf]); + hphiB += phiB * phiB; + } + assert(!std::isnan(hphiB)); + + const double phiL = std::max(0., ptr_phi[norderp - 1][idx_pf]); + double hphil = phiL * phiL; + + // std::cout<<"h = "< 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_local[i] = 2. * std::max(0., 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_local[i] = 2. * std::max(0., ptr_phi[i][idx_pf]) * + ((1. - hphiB) * dfb - hphil * dfl - hphiA * dfa) * + sum2inv; + // liquid phase order parameter + rhs_local[norderp - 1] = + 2. * std::max(0., 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_local[i])); + + // damp driving force when close to 100% phase B + const double epsilon1 = 1.e-4; + const double epsilon2 = 2.e-4; + double factor = 1.; + if (hphil < epsilon1) { + factor = 0.; + } else { + if (hphil < epsilon2) { + // map epsilon1 < hphi < epsilon2 to (0,1) + double h = (hphil - epsilon1) / (epsilon2 - epsilon1); + factor = Thermo4PFM::interp_func( + Thermo4PFM::EnergyInterpolationType::PBG, h); + } + } + + for (short i = 0; i < norderp; i++) + ptr_rhs[i][idx_rhs] -= factor * (rhs_local[i]); + } + } + } +} diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h new file mode 100644 index 00000000..a04753ea --- /dev/null +++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h @@ -0,0 +1,41 @@ +// 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 MultiOrderBinaryThreePhasesDrivingForceStochioAB_H +#define MultiOrderBinaryThreePhasesDrivingForceStochioAB_H + +#include "FreeEnergyStrategyBinary.h" + +#include "SAMRAI/hier/Patch.h" + +using namespace SAMRAI; + +class MultiOrderBinaryThreePhasesDrivingForceStochioAB +{ + public: + MultiOrderBinaryThreePhasesDrivingForceStochioAB( + FreeEnergyStrategyBinary* fenergy_strategy, const int norderp_A); + + void addDrivingForce(hier::Patch& patch, const int temperature_id, + const int phase_id, const int conc_l_id, + const int conc_a_id, const int conc_b_id, + const int f_l_id, const int f_a_id, const int f_b_id, + const int rhs_id); + + private: + /*! + * Class to evaluate chemical potential in liquid phase + */ + FreeEnergyStrategyBinary* d_fenergy_strategy; + + int d_norderp_A; +}; + +#endif diff --git a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc index 85c95811..f85dcb29 100644 --- a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc +++ b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc @@ -16,8 +16,6 @@ ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder:: ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder( const int conc_l_id, const int conc_a_id, - const Thermo4PFM::EnergyInterpolationType energy_interp_func_type, - const Thermo4PFM::ConcInterpolationType conc_interp_func_type, std::shared_ptr conc_db) : EquilibriumPhaseConcentrationsBinaryMultiOrder(conc_l_id, conc_a_id, conc_db) @@ -33,5 +31,6 @@ ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder:: double Tref = input_db->getDouble("Tref"); d_fenergy.reset(new Thermo4PFM::ParabolicFreeEnergyFunctionsBinary( - Tref, coeffL, coeffA, energy_interp_func_type, conc_interp_func_type)); + Tref, coeffL, coeffA, Thermo4PFM::EnergyInterpolationType::LINEAR, + Thermo4PFM::ConcInterpolationType::LINEAR)); } diff --git a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h index a429a5e3..edfd4e30 100644 --- a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h +++ b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h @@ -24,8 +24,6 @@ class ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder public: ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder( const int conc_l_id, const int conc_a_id, - const Thermo4PFM::EnergyInterpolationType energy_interp_func_type, - const Thermo4PFM::ConcInterpolationType conc_interp_func_type, std::shared_ptr conc_db); ~ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder() {} diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc new file mode 100644 index 00000000..84af0581 --- /dev/null +++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc @@ -0,0 +1,83 @@ +// 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 "ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h" +#include "FuncFort.h" +#include "ParabolicTools.h" + +ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB:: + ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB( + 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), + d_model_parameters(model_parameters) +{ +} + +ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB:: + ~ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(){}; + +int ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB:: + computePhaseConcentrations(const double temp, double* c, double* hphi, + double* x) +{ + assert(!std::isnan(x[0])); + + const double epsilon1 = 1.e-4; + const double epsilon2 = 2.e-4; + + // initialize to NaN to trigger error if used when not set + double xkks = tbox::IEEE::getSignalingNaN(); + double xeq = tbox::IEEE::getSignalingNaN(); + + const double cA = d_model_parameters.getStochioA(); + const double cB = d_model_parameters.getStochioB(); + + if (hphi[0] >= epsilon1) { + // solve explicit KKS problem + xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / hphi[0]; + +#if 0 + for (short i = 0; i < 3; i++) + std::cerr << hphi[i] << ", "; + std::cerr << ", c=" << c[0] << std::endl; + std::cerr << "x = " << x[0] << std::endl; + std::cerr << "xkks = " << xkks << std::endl; + std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = " + << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl; +#endif + } + + if (hphi[0] < epsilon2) { + xeq = d_model_parameters.ceq_liquid(temp); + // std::cout << "xeq = " << xeq << std::endl; + assert(!std::isnan(xeq)); + } + + if (hphi[0] >= epsilon2) { + x[0] = xkks; + } else if (hphi[0] <= epsilon1) { + x[0] = xeq; + } else { + // map epsilon1 < hphi < epsilon2 to (0,1) + double h = (hphi[0] - epsilon1) / (epsilon2 - epsilon1); + // mix equilibrium compositions and KKS solution + double f = + Thermo4PFM::interp_func(Thermo4PFM::EnergyInterpolationType::PBG, h); + x[0] = (1. - f) * xeq + f * xkks; + } + x[1] = cA; + x[2] = cB; + + return 0; +} diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h new file mode 100644 index 00000000..bc176310 --- /dev/null +++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h @@ -0,0 +1,36 @@ +// 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_ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB +#define included_ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB + +#include "EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.h" +#include "QuatModelParameters.h" + +class ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB + : public EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases +{ + public: + ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB( + 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); + + ~ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(); + + protected: + virtual int computePhaseConcentrations(const double t, double* c, + double* hphi, double* x); + + private: + const QuatModelParameters d_model_parameters; +}; + +#endif diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc index 18b08472..b6467ce1 100644 --- a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc @@ -24,13 +24,12 @@ using namespace SAMRAI; ParabolicFreeEnergyMultiOrderBinaryThreePhase:: ParabolicFreeEnergyMultiOrderBinaryThreePhase( - std::shared_ptr input_db, - const Thermo4PFM::ConcInterpolationType conc_interp_func_type, - const short norderp_A, MolarVolumeStrategy* mvstrategy, - const int conc_l_id, const int conc_a_id, const int conc_b_id) + std::shared_ptr input_db, const short norderp_A, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id) : FreeEnergyStrategyBinary(Thermo4PFM::EnergyInterpolationType::LINEAR, - conc_interp_func_type, conc_l_id, conc_a_id, - conc_b_id, false), + Thermo4PFM::ConcInterpolationType::LINEAR, + conc_l_id, conc_a_id, conc_b_id, false), d_mv_strategy(mvstrategy) { tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase..." diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h index 79adf021..5f3d4ab6 100644 --- a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h @@ -21,10 +21,9 @@ class ParabolicFreeEnergyMultiOrderBinaryThreePhase { public: ParabolicFreeEnergyMultiOrderBinaryThreePhase( - std::shared_ptr input_db, - const Thermo4PFM::ConcInterpolationType conc_interp_func_type, - const short norderp_A, MolarVolumeStrategy* mvstrategy, - const int conc_l_id, const int conc_a_id, const int conc_b_id); + std::shared_ptr input_db, const short norderp_A, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id); ~ParabolicFreeEnergyMultiOrderBinaryThreePhase(); diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc new file mode 100644 index 00000000..66c6e570 --- /dev/null +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc @@ -0,0 +1,218 @@ +// 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 "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h" +#include "ParabolicTools.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 + +//======================================================================= + +ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB( + std::shared_ptr input_db, const short norderp_A, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id) + : FreeEnergyStrategyBinary(Thermo4PFM::EnergyInterpolationType::LINEAR, + Thermo4PFM::ConcInterpolationType::LINEAR, + conc_l_id, conc_a_id, conc_b_id, false), + d_mv_strategy(mvstrategy) +{ + tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB..." + << std::endl; + + assert(norderp_A > 0); + assert(conc_b_id >= 0); + + double coeffL[3][2]; + readParabolicData(input_db, "Liquid", coeffL); + + double coeffA[3][2]; + readParabolicData(input_db, "PhaseA", coeffA); + + double coeffB[3][2]; + readParabolicData(input_db, "PhaseB", coeffB); + + double Tref = input_db->getDouble("Tref"); + + d_parabolic_fenergy.reset( + new Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase( + Tref, coeffL, coeffA, coeffB, + Thermo4PFM::EnergyInterpolationType::LINEAR, + Thermo4PFM::ConcInterpolationType::LINEAR)); + + d_multiorder_driving_force.reset( + new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, norderp_A)); + + // 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 << "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:" << + // 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; +} + +ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + ~ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(){}; + +bool ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeCeqT( + const double temperature, const Thermo4PFM::PhaseIndex pi0, + const Thermo4PFM::PhaseIndex pi1, double* ceq) +{ + TBOX_ERROR( + "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeCeqT() " + "not " + "implemented"); + return false; +} + + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + computeFreeEnergy(const double temperature, double* c_i, + const Thermo4PFM::PhaseIndex pi, const bool gp) +{ + assert(d_mv_strategy != nullptr); + + double f = d_parabolic_fenergy->computeFreeEnergy(temperature, c_i, pi, gp); + f *= d_mv_strategy->computeInvMolarVolume(temperature, c_i, pi); + return f; +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + computeDerivFreeEnergy(const double temperature, double* c_i, + const Thermo4PFM::PhaseIndex pi) +{ + double deriv; + d_parabolic_fenergy->computeDerivFreeEnergy(temperature, c_i, pi, &deriv); + deriv *= d_mv_strategy->computeInvMolarVolume(temperature, c_i, pi); + return deriv; +} + + +void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::addDrivingForce( + const double time, hier::Patch& patch, const int temperature_id, + const int phase_id, const int conc_id, const int f_l_id, const int f_a_id, + const int f_b_id, const int rhs_id) +{ + (void)time; + (void)f_b_id; + + d_multiorder_driving_force->addDrivingForce(patch, temperature_id, phase_id, + d_conc_l_id, d_conc_a_id, + d_conc_b_id, f_l_id, f_a_id, + f_b_id, rhs_id); +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuL( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseL, + &mu); + return mu * + d_mv_strategy->computeInvMolarVolume(t, &c, + Thermo4PFM::PhaseIndex::phaseL); +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuA( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseA, + &mu); + return mu * + d_mv_strategy->computeInvMolarVolume(t, &c, + Thermo4PFM::PhaseIndex::phaseA); +} + +//======================================================================= + +double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuB( + const double t, const double c0) +{ + double c = c0; + double mu; + d_parabolic_fenergy->computeDerivFreeEnergy(t, &c, + Thermo4PFM::PhaseIndex::phaseB, + &mu); + return mu * + d_mv_strategy->computeInvMolarVolume(t, &c, + Thermo4PFM::PhaseIndex::phaseB); +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + computeSecondDerivativeEnergyPhaseL(const double temp, + 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_mv_strategy->computeInvMolarVolume( + temp, &c_l[0], Thermo4PFM::PhaseIndex::phaseL); +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + computeSecondDerivativeEnergyPhaseA(const double temp, + 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_mv_strategy->computeInvMolarVolume( + temp, &c_a[0], Thermo4PFM::PhaseIndex::phaseA); +} + +//======================================================================= + +void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:: + computeSecondDerivativeEnergyPhaseB(const double temp, + 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_mv_strategy->computeInvMolarVolume( + temp, &c_b[0], Thermo4PFM::PhaseIndex::phaseB); +} diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h new file mode 100644 index 00000000..d9e6b79d --- /dev/null +++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h @@ -0,0 +1,73 @@ +// 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_ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB +#define included_ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB + +#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h" +#include "FreeEnergyStrategyBinary.h" +#include "MolarVolumeStrategy.h" +#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h" + +class ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB + : public FreeEnergyStrategyBinary +{ + public: + ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB( + std::shared_ptr input_db, const short norderp_A, + MolarVolumeStrategy* mvstrategy, const int conc_l_id, + const int conc_a_id, const int conc_b_id); + + ~ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(); + + void addDrivingForce(const double time, hier::Patch& patch, + const int temperature_id, const int phase_id, + const int conc_id, const int f_l_id, const int f_a_id, + const int f_b_id, const int rhs_id) override; + + void preRunDiagnostics(const double temperature){}; + + bool computeCeqT(const double temperature, const Thermo4PFM::PhaseIndex pi0, + const Thermo4PFM::PhaseIndex pi1, double* ceq); + + protected: + double computeFreeEnergy(const double temperature, double* c_i, + const Thermo4PFM::PhaseIndex pi, + const bool gp) override; + + double computeDerivFreeEnergy(const double temperature, double* c_i, + const Thermo4PFM::PhaseIndex pi) override; + + private: + MolarVolumeStrategy* d_mv_strategy; + + std::shared_ptr + d_parabolic_fenergy; + + std::shared_ptr + d_multiorder_driving_force; + + void computeSecondDerivativeEnergyPhaseL( + const double temp, const std::vector& c, + std::vector& d2fdc2, const bool use_internal_units) override; + void computeSecondDerivativeEnergyPhaseA( + const double temp, const std::vector& c, + std::vector& d2fdc2, const bool use_internal_units) override; + void computeSecondDerivativeEnergyPhaseB(const double temp, + const std::vector& c, + std::vector& d2fdc2, + const bool use_internal_units); + + double computeMuL(const double t, const double c); + double computeMuA(const double t, const double c); + double computeMuB(const double t, const double c); +}; + +#endif diff --git a/source/PhaseConcentrationsStrategyFactory.h b/source/PhaseConcentrationsStrategyFactory.h index c8bbdbe4..ff527063 100644 --- a/source/PhaseConcentrationsStrategyFactory.h +++ b/source/PhaseConcentrationsStrategyFactory.h @@ -20,6 +20,7 @@ #include "QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder.h" #include "QuadraticEquilibriumThreePhasesTernaryMultiOrder.h" #include "ParabolicEquilibriumThreePhasesBinaryMultiOrder.h" +#include "ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h" #include "PartitionPhaseConcentrationsStrategy.h" #include "PhaseIndependentConcentrationsStrategy.h" #include "Database2JSON.h" @@ -29,14 +30,12 @@ #include "CALPHADequilibriumPhaseConcentrationsStrategyMultiOrder.h" #include "CALPHADequilibriumPhaseConcentrationsStrategyMultiOrderThreePhases.h" #include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB.h" +#include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h" #include "CALPHADFunctions.h" #include #include -#include -#include - class PhaseConcentrationsStrategyFactory { public: @@ -68,6 +67,11 @@ class PhaseConcentrationsStrategyFactory std::shared_ptr phase_conc_strategy; + const double cA = model_parameters.getStochioA(); + const double cB = model_parameters.getStochioB(); + tbox::plog << "cA = " << cA << std::endl; + tbox::plog << "cB = " << cB << std::endl; + if (model_parameters.kks_phase_concentration()) { tbox::plog << "Phase concentration determined by KKS" << std::endl; if (model_parameters.isConcentrationModelCALPHAD()) { @@ -75,16 +79,21 @@ class PhaseConcentrationsStrategyFactory if (ncompositions == 1) { tbox::plog << "Binary alloy..." << std::endl; const bool subl = Thermo4PFM::checkSublattice(calphad_pt); - const double cB = model_parameters.getStochioB(); - tbox::plog << "cB = " << cB << std::endl; if (conc_b_scratch_id >= 0) { tbox::plog << "Three phases..." << std::endl; // three phases if (model_parameters.withMultipleOrderP()) { // multi-order parameters model tbox::plog << "Multi-order parameters..." << std::endl; - if (cB >= 0.) { - tbox::plog << "Stochiometric case..." << std::endl; + if (cB >= 0. && cA >= 0.) { + tbox::plog << "AB stochiometric case..." << std::endl; + phase_conc_strategy.reset( + new CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB( + model_parameters.norderpA(), conc_l_scratch_id, + conc_a_scratch_id, conc_b_scratch_id, + model_parameters, conc_db)); + } else if (cB >= 0.) { + tbox::plog << "B stochiometric case..." << std::endl; phase_conc_strategy.reset( new CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB( model_parameters.norderpA(), conc_l_scratch_id, @@ -213,18 +222,23 @@ class PhaseConcentrationsStrategyFactory 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)); + if (cB >= 0. && cA >= 0.) { + phase_conc_strategy.reset( + new ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB( + 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 ParabolicEquilibriumThreePhasesBinaryMultiOrder( + 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 ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder( - conc_l_scratch_id, conc_a_scratch_id, - model_parameters.energy_interp_func_type(), - model_parameters.conc_interp_func_type(), - conc_db)); + conc_l_scratch_id, conc_a_scratch_id, conc_db)); } } else { phase_conc_strategy.reset( diff --git a/source/QuatModelParameters.cc b/source/QuatModelParameters.cc index 5006ab49..d93fc97b 100644 --- a/source/QuatModelParameters.cc +++ b/source/QuatModelParameters.cc @@ -90,6 +90,7 @@ QuatModelParameters::QuatModelParameters() : d_moving_frame_velocity(def_val) d_liquidus_slope = def_val; d_average_concentration = def_val; d_stochio_cB = -1.; + d_stochio_cA = -1.; for (short i = 0; i < 3; i++) { d_ceq0_solidA[i] = def_val; @@ -433,6 +434,7 @@ void QuatModelParameters::readConcDB(std::shared_ptr conc_db) conc_db->getBoolWithDefault("init_phase_conc_eq", true); d_stochio_cB = conc_db->getDoubleWithDefault("stochio_cB", -1.); + d_stochio_cA = conc_db->getDoubleWithDefault("stochio_cA", -1.); readEquilibriumCompositions(conc_db); } diff --git a/source/QuatModelParameters.h b/source/QuatModelParameters.h index 31bbc23c..b7798414 100644 --- a/source/QuatModelParameters.h +++ b/source/QuatModelParameters.h @@ -141,7 +141,10 @@ class QuatModelParameters { return d_initc_in_phase[d_ncompositions + index]; } + double getStochioB() const { return d_stochio_cB; } + double getStochioA() const { return d_stochio_cA; } + double meltingT() const { return d_meltingT; } double interfaceMobility() const { return d_interface_mobility; } double rescale_factorT() const { return d_rescale_factorT; } @@ -626,6 +629,7 @@ class QuatModelParameters // stochio composition of phase B (if set to value >= 0) double d_stochio_cB; + double d_stochio_cA; // free energy parameters: // f(phi) = d_phase_well_scale * g(phi) From dd81595626bfbd7704f9d3243c0254b2063d8128 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 8 May 2025 14:22:02 -0400 Subject: [PATCH 2/6] Limit changes in cl --- source/CMakeLists.txt | 1 + source/LimitDifference.cc | 72 +++++++++++++++++++++++++++++++++++++++ source/LimitDifference.h | 37 ++++++++++++++++++++ source/QuatModel.cc | 18 ++++++++-- source/QuatModel.h | 2 ++ 5 files changed, 128 insertions(+), 2 deletions(-) create mode 100644 source/LimitDifference.cc create mode 100644 source/LimitDifference.h diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index cd468efb..e6bf9752 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -26,6 +26,7 @@ set(SOURCES_WPENALTY ) set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc + ${CMAKE_SOURCE_DIR}/source/LimitDifference.cc ${CMAKE_SOURCE_DIR}/source/ParabolicTools.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryDrivingForce.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForce.cc diff --git a/source/LimitDifference.cc b/source/LimitDifference.cc new file mode 100644 index 00000000..9a62df7d --- /dev/null +++ b/source/LimitDifference.cc @@ -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. + +#include "LimitDifference.h" + +#include "FuncFort.h" + +#include "SAMRAI/pdat/CellIterator.h" +#include "SAMRAI/pdat/CellData.h" + +using namespace SAMRAI; + +LimitDifference::LimitDifference( + const std::shared_ptr hierarchy) + : d_hierarchy(hierarchy) +{ +} + +void LimitDifference::apply(const int cell_data_id, const int ref_cell_data_id, + const double diff) + +{ + tbox::pout << "Limit difference to " << diff << std::endl; + const int maxln = d_hierarchy->getFinestLevelNumber(); + for (int ln = 0; ln <= maxln; ln++) { + std::shared_ptr patch_level = + d_hierarchy->getPatchLevel(ln); + + apply(patch_level, cell_data_id, ref_cell_data_id, diff); + } +} + +void LimitDifference::apply(const std::shared_ptr level, + const int cell_data_id, const int ref_cell_data_id, + const double diff) + +{ + for (hier::PatchLevel::Iterator p(level->begin()); p != level->end(); ++p) { + std::shared_ptr patch = *p; + + std::shared_ptr > data( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(cell_data_id))); + assert(data); + const hier::Box& gbox = data->getGhostBox(); + + std::shared_ptr > ref_data( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(ref_cell_data_id))); + assert(ref_data); + + pdat::CellIterator iend(pdat::CellGeometry::end(gbox)); + for (pdat::CellIterator i(pdat::CellGeometry::begin(gbox)); i != iend; + ++i) { + pdat::CellIndex cell = *i; + const double d = (*data)(cell) - (*ref_data)(cell); + if (std::abs(d) > diff) { + if (d > 0.) + (*data)(cell) += diff; + else + (*data)(cell) -= diff; + } + } + } +} diff --git a/source/LimitDifference.h b/source/LimitDifference.h new file mode 100644 index 00000000..6d8c4f04 --- /dev/null +++ b/source/LimitDifference.h @@ -0,0 +1,37 @@ +// 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_LimitDifference +#define included_LimitDifference + +#include "SAMRAI/hier/PatchHierarchy.h" +#include "SAMRAI/hier/PatchLevel.h" + +class LimitDifference +{ + public: + LimitDifference( + const std::shared_ptr hierarchy); + + /*! + * Change values in cell_data so that they differ from ref_cell_data + * by less than "diff" + */ + void apply(const int cell_data_id, const int ref_cell_data_id, + const double diff); + + private: + const std::shared_ptr d_hierarchy; + + void apply(const std::shared_ptr level, + const int cell_data_id, const int ref_cell_data_id, + const double diff); +}; + +#endif diff --git a/source/QuatModel.cc b/source/QuatModel.cc index d4464fde..6ea87da2 100644 --- a/source/QuatModel.cc +++ b/source/QuatModel.cc @@ -50,6 +50,7 @@ #include "computeQDiffs.h" #include "HierarchyStencilOps.h" #include "WangRigidBodyForces.h" +#include "LimitDifference.h" #include "Database2JSON.h" namespace pt = boost::property_tree; @@ -2323,7 +2324,8 @@ double QuatModel::Advance(void) bool update_solution = false; if (d_model_parameters.with_concentration() && - d_model_parameters.isConcentrationModelCALPHAD()) + (d_model_parameters.isConcentrationModelCALPHAD() || + d_model_parameters.getStochioA() >= 0.)) resetRefPhaseConcentrations(); if (d_model_parameters.with_orientation()) { @@ -5021,6 +5023,10 @@ void QuatModel::computePhaseConcentrations( d_phase_scratch_id, d_conc_scratch_id); + if (d_model_parameters.getStochioA() >= 0.) { + limitClDifference(1.e-4); + } + t_phase_conc_timer->stop(); } @@ -5274,7 +5280,7 @@ void QuatModel::resetRefPhaseConcentrations() assert(d_conc_l_ref_id >= 0); assert(d_conc_a_ref_id >= 0); - // tbox::pout << "QuatModel::resetRefPhaseConcentrations()" << std::endl; + tbox::pout << "QuatModel::resetRefPhaseConcentrations()" << std::endl; math::HierarchyCellDataOpsReal cellops(d_patch_hierarchy); @@ -5284,6 +5290,14 @@ void QuatModel::resetRefPhaseConcentrations() cellops.copyData(d_conc_b_ref_id, d_conc_b_id, false); } +//======================================================================= + +void QuatModel::limitClDifference(const double delta) +{ + LimitDifference limit(d_patch_hierarchy); + limit.apply(d_conc_l_id, d_conc_l_ref_id, delta); +} + //======================================================================= /* void QuatModel::setPhaseConcentrationsToEquilibrium(const double* const ceq) diff --git a/source/QuatModel.h b/source/QuatModel.h index 806619f9..cbb06278 100644 --- a/source/QuatModel.h +++ b/source/QuatModel.h @@ -298,6 +298,8 @@ class QuatModel : public PFModel void setDiffusionVisit( const std::shared_ptr hierarchy); + void limitClDifference(const double delta); + private: void setAuxilliaryCompositions(); void registerPhaseConcentrationDiffusionVariables(); From afaf46c17dfd4c40681acc5e0ef56afbda5e194f Mon Sep 17 00:00:00 2001 From: "Fattebert J.-L" Date: Thu, 8 May 2025 15:40:45 -0400 Subject: [PATCH 3/6] use dt to set limit on delta cl --- source/LimitDifference.cc | 3 +++ source/QuatIntegrator.cc | 2 +- source/QuatModel.cc | 10 ++++++++-- source/QuatModel.h | 2 +- 4 files changed, 13 insertions(+), 4 deletions(-) diff --git a/source/LimitDifference.cc b/source/LimitDifference.cc index 9a62df7d..ae9c79c3 100644 --- a/source/LimitDifference.cc +++ b/source/LimitDifference.cc @@ -27,6 +27,9 @@ void LimitDifference::apply(const int cell_data_id, const int ref_cell_data_id, const double diff) { + assert(cell_data_id >= 0); + assert(ref_cell_data_id >= 0); + tbox::pout << "Limit difference to " << diff << std::endl; const int maxln = d_hierarchy->getFinestLevelNumber(); for (int ln = 0; ln <= maxln; ln++) { diff --git a/source/QuatIntegrator.cc b/source/QuatIntegrator.cc index 43b8f43e..e8e05a35 100644 --- a/source/QuatIntegrator.cc +++ b/source/QuatIntegrator.cc @@ -3206,7 +3206,7 @@ int QuatIntegrator::evaluateRHSFunction(double time, SundialsAbstractVector* y, // compute phase concentrations again if they depend on velocity // tbox::pout<<"Evaluate c_L, c_S..."<computePhaseConcentrations(hierarchy); + d_quat_model->computePhaseConcentrations(time, hierarchy); } } while (need_iterate); diff --git a/source/QuatModel.cc b/source/QuatModel.cc index 6ea87da2..fc40403d 100644 --- a/source/QuatModel.cc +++ b/source/QuatModel.cc @@ -1214,7 +1214,8 @@ void QuatModel::registerPhaseConcentrationVariables() assert(d_conc_b_var); } - if (d_model_parameters.isConcentrationModelCALPHAD()) { + if (d_model_parameters.isConcentrationModelCALPHAD() || + d_model_parameters.getStochioA() >= 0.) { d_conc_l_ref_var.reset( new pdat::CellVariable(tbox::Dimension(NDIM), "conc_l_ref", d_ncompositions)); @@ -4986,6 +4987,7 @@ void QuatModel::evaluateEnergy( //======================================================================= void QuatModel::computePhaseConcentrations( +const double time, const std::shared_ptr hierarchy) { assert(d_phase_conc_strategy != nullptr); @@ -5024,7 +5026,9 @@ void QuatModel::computePhaseConcentrations( d_conc_scratch_id); if (d_model_parameters.getStochioA() >= 0.) { - limitClDifference(1.e-4); + const double dt = time - d_current_time; + + d_quat_model->limitClDifference(dt * 1.e5); } t_phase_conc_timer->stop(); @@ -5294,6 +5298,8 @@ void QuatModel::resetRefPhaseConcentrations() void QuatModel::limitClDifference(const double delta) { + assert(d_conc_l_ref_id >= 0); + LimitDifference limit(d_patch_hierarchy); limit.apply(d_conc_l_id, d_conc_l_ref_id, delta); } diff --git a/source/QuatModel.h b/source/QuatModel.h index cbb06278..3376ea75 100644 --- a/source/QuatModel.h +++ b/source/QuatModel.h @@ -224,7 +224,7 @@ class QuatModel : public PFModel void computeVectorWeightsDiagnostics( std::shared_ptr hierarchy); - void computePhaseConcentrations( + void computePhaseConcentrations(const double time, std::shared_ptr hierarchy); void evaluateEnergy(const std::shared_ptr hierarchy, From 5a494d81cd3725f0a2475fdb8c8735e8f8600151 Mon Sep 17 00:00:00 2001 From: "Fattebert J.-L" Date: Mon, 12 May 2025 10:18:08 -0400 Subject: [PATCH 4/6] Use liquid phase parameter in setting limit --- source/LimitDifference.cc | 29 +++++++++++++++++++++-------- source/LimitDifference.h | 4 ++-- source/QuatIntegrator.cc | 6 ++++-- source/QuatModel.cc | 17 ++++++++--------- source/QuatModel.h | 4 ++-- 5 files changed, 37 insertions(+), 23 deletions(-) diff --git a/source/LimitDifference.cc b/source/LimitDifference.cc index ae9c79c3..d3444f84 100644 --- a/source/LimitDifference.cc +++ b/source/LimitDifference.cc @@ -24,6 +24,7 @@ LimitDifference::LimitDifference( } void LimitDifference::apply(const int cell_data_id, const int ref_cell_data_id, + const int weight_id, const int depth, const double diff) { @@ -36,12 +37,14 @@ void LimitDifference::apply(const int cell_data_id, const int ref_cell_data_id, std::shared_ptr patch_level = d_hierarchy->getPatchLevel(ln); - apply(patch_level, cell_data_id, ref_cell_data_id, diff); + apply(patch_level, cell_data_id, ref_cell_data_id, weight_id, depth, + diff); } } void LimitDifference::apply(const std::shared_ptr level, const int cell_data_id, const int ref_cell_data_id, + const int weight_id, const int depth, const double diff) { @@ -59,17 +62,27 @@ void LimitDifference::apply(const std::shared_ptr level, patch->getPatchData(ref_cell_data_id))); assert(ref_data); + std::shared_ptr > weight( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(weight_id))); + assert(weight); + + const double inv1pd = 1. / (1. + diff); pdat::CellIterator iend(pdat::CellGeometry::end(gbox)); for (pdat::CellIterator i(pdat::CellGeometry::begin(gbox)); i != iend; ++i) { pdat::CellIndex cell = *i; - const double d = (*data)(cell) - (*ref_data)(cell); - if (std::abs(d) > diff) { - if (d > 0.) - (*data)(cell) += diff; - else - (*data)(cell) -= diff; - } + assert((*ref_data)(cell) >= 0.); + assert((*ref_data)(cell) <= 1.); + const double delta = (*data)(cell) - (*ref_data)(cell); + const double w = std::max(0., (*weight)(cell, depth)); + const double x = 2. * w; + const double t = x > 1. ? 1. : x * x * x * (10. - 15. * x + 6. * x); + const double tt = (t + diff) * inv1pd; + // const double tt = std::abs(delta)>diff ? diff/delta : 1.; + (*data)(cell) = (*ref_data)(cell) + tt * delta; + if ((*data)(cell) > 1.) (*data)(cell) = 1.; + if ((*data)(cell) < 0.) (*data)(cell) = 0.; } } } diff --git a/source/LimitDifference.h b/source/LimitDifference.h index 6d8c4f04..5c0cee2d 100644 --- a/source/LimitDifference.h +++ b/source/LimitDifference.h @@ -24,14 +24,14 @@ class LimitDifference * by less than "diff" */ void apply(const int cell_data_id, const int ref_cell_data_id, - const double diff); + const int wright_id, const int depth, const double diff); private: const std::shared_ptr d_hierarchy; void apply(const std::shared_ptr level, const int cell_data_id, const int ref_cell_data_id, - const double diff); + const int wright_id, const int depth, const double diff); }; #endif diff --git a/source/QuatIntegrator.cc b/source/QuatIntegrator.cc index e8e05a35..fb964d4c 100644 --- a/source/QuatIntegrator.cc +++ b/source/QuatIntegrator.cc @@ -3077,7 +3077,8 @@ void QuatIntegrator::setCoefficients( if (d_with_concentration && d_model_parameters.concentrationModelNeedsPhaseConcentrations()) { - d_quat_model->computePhaseConcentrations(hierarchy); + d_quat_model->computePhaseConcentrations(time - d_current_time, + hierarchy); } // mobilities may depend on cl and cs, thus they should be computed after @@ -3206,7 +3207,8 @@ int QuatIntegrator::evaluateRHSFunction(double time, SundialsAbstractVector* y, // compute phase concentrations again if they depend on velocity // tbox::pout<<"Evaluate c_L, c_S..."<computePhaseConcentrations(time, hierarchy); + d_quat_model->computePhaseConcentrations(time - d_current_time, + hierarchy); } } while (need_iterate); diff --git a/source/QuatModel.cc b/source/QuatModel.cc index fc40403d..e6b996f4 100644 --- a/source/QuatModel.cc +++ b/source/QuatModel.cc @@ -735,6 +735,7 @@ void QuatModel::Initialize(std::shared_ptr& input_db, const double cLref = d_model_parameters.concL_ref(); if (cLref >= 0.) { tbox::plog << "Set reference cL to " << cLref << std::endl; + assert(cLref <= 1.); if (d_conc_l_ref_id >= 0) cellops.setToScalar(d_conc_l_ref_id, cLref, false); cellops.setToScalar(d_conc_l_id, cLref, false); @@ -4942,9 +4943,7 @@ void QuatModel::evaluateEnergy( } if (d_model_parameters.concentrationModelNeedsPhaseConcentrations()) { assert(d_phase_conc_strategy != nullptr); - d_phase_conc_strategy->computePhaseConcentrations( - hierarchy, d_temperature_scratch_id, d_phase_scratch_id, - d_conc_scratch_id); + computePhaseConcentrations(0., hierarchy); } } @@ -4987,8 +4986,7 @@ void QuatModel::evaluateEnergy( //======================================================================= void QuatModel::computePhaseConcentrations( -const double time, - const std::shared_ptr hierarchy) + const double dt, const std::shared_ptr hierarchy) { assert(d_phase_conc_strategy != nullptr); @@ -5026,9 +5024,7 @@ const double time, d_conc_scratch_id); if (d_model_parameters.getStochioA() >= 0.) { - const double dt = time - d_current_time; - - d_quat_model->limitClDifference(dt * 1.e5); + limitClDifference(dt * 1.e5); } t_phase_conc_timer->stop(); @@ -5301,7 +5297,8 @@ void QuatModel::limitClDifference(const double delta) assert(d_conc_l_ref_id >= 0); LimitDifference limit(d_patch_hierarchy); - limit.apply(d_conc_l_id, d_conc_l_ref_id, delta); + limit.apply(d_conc_l_id, d_conc_l_ref_id, d_phase_id, + d_model_parameters.norderp() - 1, delta); } //======================================================================= @@ -5356,6 +5353,8 @@ void QuatModel::setRefPhaseConcentrationsToEquilibrium(const double* const ceq) { assert(d_conc_l_ref_id >= 0); assert(d_conc_a_ref_id >= 0); + assert(ceq[0] >= 0.); + assert(ceq[0] <= 1.); tbox::pout << "QuatModel::setRefPhaseConcentrationsToEquilibrium(ceq)" << std::endl; diff --git a/source/QuatModel.h b/source/QuatModel.h index 3376ea75..8d107cc1 100644 --- a/source/QuatModel.h +++ b/source/QuatModel.h @@ -224,8 +224,8 @@ class QuatModel : public PFModel void computeVectorWeightsDiagnostics( std::shared_ptr hierarchy); - void computePhaseConcentrations(const double time, - std::shared_ptr hierarchy); + void computePhaseConcentrations( + const double time, std::shared_ptr hierarchy); void evaluateEnergy(const std::shared_ptr hierarchy, const double time, double& total_energy, From 5dd8c36ed7a130bf85fb515a0b2bb1ca06f05d57 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Sat, 17 May 2025 20:47:03 -0400 Subject: [PATCH 5/6] Stabilize computation of cl --- ...ntrationsMultiOrderThreePhasesStochioAB.cc | 39 ++++---- source/CMakeLists.txt | 1 - source/LimitDifference.cc | 88 ------------------- source/LimitDifference.h | 37 -------- ...rBinaryThreePhasesDrivingForceStochioAB.cc | 17 +--- source/QuatModel.cc | 18 +--- source/QuatModel.h | 2 - source/fortran/ConcFort.h | 11 +++ 8 files changed, 29 insertions(+), 184 deletions(-) delete mode 100644 source/LimitDifference.cc delete mode 100644 source/LimitDifference.h diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc index 7517fba9..a22e23fb 100644 --- a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc +++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc @@ -50,36 +50,29 @@ int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: { assert(!std::isnan(x[0])); - const double epsilon1 = 1.e-4; - const double epsilon2 = 2.e-4; - - // initialize to NaN to trigger error if used when not set - double xkks = tbox::IEEE::getSignalingNaN(); - double xeq = tbox::IEEE::getSignalingNaN(); + const double epsilon1 = 1.e-3; + const double epsilon2 = 2.e-3; const double cA = d_model_parameters.getStochioA(); const double cB = d_model_parameters.getStochioB(); - if (hphi[0] >= epsilon1) { - // solve explicit KKS problem - xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / hphi[0]; + const double xeq = d_model_parameters.ceq_liquid(temp); + // std::cout << "xeq = " << xeq << std::endl; + assert(!std::isnan(xeq)); + // solve explicit KKS problem + const double a = std::max(epsilon1, hphi[0]); + double xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / a; + xkks = xeq + (1. - xeq) * std::tanh(xkks - xeq); #if 0 - for (short i = 0; i < 3; i++) - std::cerr << hphi[i] << ", "; - std::cerr << ", c=" << c[0] << std::endl; - std::cerr << "x = " << x[0] << std::endl; - std::cerr << "xkks = " << xkks << std::endl; - std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = " - << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl; + for (short i = 0; i < 3; i++) + std::cerr << hphi[i] << ", "; + std::cerr << ", c=" << c[0] << std::endl; + std::cerr << "x = " << x[0] << std::endl; + std::cerr << "xkks = " << xkks << std::endl; + std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = " + << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl; #endif - } - - if (hphi[0] < epsilon2) { - xeq = d_model_parameters.ceq_liquid(temp); - // std::cout << "xeq = " << xeq << std::endl; - assert(!std::isnan(xeq)); - } if (hphi[0] >= epsilon2) { x[0] = xkks; diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index e6bf9752..cd468efb 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -26,7 +26,6 @@ set(SOURCES_WPENALTY ) set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc - ${CMAKE_SOURCE_DIR}/source/LimitDifference.cc ${CMAKE_SOURCE_DIR}/source/ParabolicTools.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryDrivingForce.cc ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForce.cc diff --git a/source/LimitDifference.cc b/source/LimitDifference.cc deleted file mode 100644 index d3444f84..00000000 --- a/source/LimitDifference.cc +++ /dev/null @@ -1,88 +0,0 @@ -// 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 "LimitDifference.h" - -#include "FuncFort.h" - -#include "SAMRAI/pdat/CellIterator.h" -#include "SAMRAI/pdat/CellData.h" - -using namespace SAMRAI; - -LimitDifference::LimitDifference( - const std::shared_ptr hierarchy) - : d_hierarchy(hierarchy) -{ -} - -void LimitDifference::apply(const int cell_data_id, const int ref_cell_data_id, - const int weight_id, const int depth, - const double diff) - -{ - assert(cell_data_id >= 0); - assert(ref_cell_data_id >= 0); - - tbox::pout << "Limit difference to " << diff << std::endl; - const int maxln = d_hierarchy->getFinestLevelNumber(); - for (int ln = 0; ln <= maxln; ln++) { - std::shared_ptr patch_level = - d_hierarchy->getPatchLevel(ln); - - apply(patch_level, cell_data_id, ref_cell_data_id, weight_id, depth, - diff); - } -} - -void LimitDifference::apply(const std::shared_ptr level, - const int cell_data_id, const int ref_cell_data_id, - const int weight_id, const int depth, - const double diff) - -{ - for (hier::PatchLevel::Iterator p(level->begin()); p != level->end(); ++p) { - std::shared_ptr patch = *p; - - std::shared_ptr > data( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(cell_data_id))); - assert(data); - const hier::Box& gbox = data->getGhostBox(); - - std::shared_ptr > ref_data( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(ref_cell_data_id))); - assert(ref_data); - - std::shared_ptr > weight( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(weight_id))); - assert(weight); - - const double inv1pd = 1. / (1. + diff); - pdat::CellIterator iend(pdat::CellGeometry::end(gbox)); - for (pdat::CellIterator i(pdat::CellGeometry::begin(gbox)); i != iend; - ++i) { - pdat::CellIndex cell = *i; - assert((*ref_data)(cell) >= 0.); - assert((*ref_data)(cell) <= 1.); - const double delta = (*data)(cell) - (*ref_data)(cell); - const double w = std::max(0., (*weight)(cell, depth)); - const double x = 2. * w; - const double t = x > 1. ? 1. : x * x * x * (10. - 15. * x + 6. * x); - const double tt = (t + diff) * inv1pd; - // const double tt = std::abs(delta)>diff ? diff/delta : 1.; - (*data)(cell) = (*ref_data)(cell) + tt * delta; - if ((*data)(cell) > 1.) (*data)(cell) = 1.; - if ((*data)(cell) < 0.) (*data)(cell) = 0.; - } - } -} diff --git a/source/LimitDifference.h b/source/LimitDifference.h deleted file mode 100644 index 5c0cee2d..00000000 --- a/source/LimitDifference.h +++ /dev/null @@ -1,37 +0,0 @@ -// 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_LimitDifference -#define included_LimitDifference - -#include "SAMRAI/hier/PatchHierarchy.h" -#include "SAMRAI/hier/PatchLevel.h" - -class LimitDifference -{ - public: - LimitDifference( - const std::shared_ptr hierarchy); - - /*! - * Change values in cell_data so that they differ from ref_cell_data - * by less than "diff" - */ - void apply(const int cell_data_id, const int ref_cell_data_id, - const int wright_id, const int depth, const double diff); - - private: - const std::shared_ptr d_hierarchy; - - void apply(const std::shared_ptr level, - const int cell_data_id, const int ref_cell_data_id, - const int wright_id, const int depth, const double diff); -}; - -#endif diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc index 3f94bbab..65282e7f 100644 --- a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc +++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc @@ -255,23 +255,8 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( for (short i = 0; i < norderp; i++) assert(!std::isnan(rhs_local[i])); - // damp driving force when close to 100% phase B - const double epsilon1 = 1.e-4; - const double epsilon2 = 2.e-4; - double factor = 1.; - if (hphil < epsilon1) { - factor = 0.; - } else { - if (hphil < epsilon2) { - // map epsilon1 < hphi < epsilon2 to (0,1) - double h = (hphil - epsilon1) / (epsilon2 - epsilon1); - factor = Thermo4PFM::interp_func( - Thermo4PFM::EnergyInterpolationType::PBG, h); - } - } - for (short i = 0; i < norderp; i++) - ptr_rhs[i][idx_rhs] -= factor * (rhs_local[i]); + ptr_rhs[i][idx_rhs] -= rhs_local[i]; } } } diff --git a/source/QuatModel.cc b/source/QuatModel.cc index e6b996f4..4be259bb 100644 --- a/source/QuatModel.cc +++ b/source/QuatModel.cc @@ -50,7 +50,6 @@ #include "computeQDiffs.h" #include "HierarchyStencilOps.h" #include "WangRigidBodyForces.h" -#include "LimitDifference.h" #include "Database2JSON.h" namespace pt = boost::property_tree; @@ -5023,10 +5022,6 @@ void QuatModel::computePhaseConcentrations( d_phase_scratch_id, d_conc_scratch_id); - if (d_model_parameters.getStochioA() >= 0.) { - limitClDifference(dt * 1.e5); - } - t_phase_conc_timer->stop(); } @@ -5280,7 +5275,7 @@ void QuatModel::resetRefPhaseConcentrations() assert(d_conc_l_ref_id >= 0); assert(d_conc_a_ref_id >= 0); - tbox::pout << "QuatModel::resetRefPhaseConcentrations()" << std::endl; + // tbox::pout << "QuatModel::resetRefPhaseConcentrations()" << std::endl; math::HierarchyCellDataOpsReal cellops(d_patch_hierarchy); @@ -5290,17 +5285,6 @@ void QuatModel::resetRefPhaseConcentrations() cellops.copyData(d_conc_b_ref_id, d_conc_b_id, false); } -//======================================================================= - -void QuatModel::limitClDifference(const double delta) -{ - assert(d_conc_l_ref_id >= 0); - - LimitDifference limit(d_patch_hierarchy); - limit.apply(d_conc_l_id, d_conc_l_ref_id, d_phase_id, - d_model_parameters.norderp() - 1, delta); -} - //======================================================================= /* void QuatModel::setPhaseConcentrationsToEquilibrium(const double* const ceq) diff --git a/source/QuatModel.h b/source/QuatModel.h index 8d107cc1..888b51cd 100644 --- a/source/QuatModel.h +++ b/source/QuatModel.h @@ -298,8 +298,6 @@ class QuatModel : public PFModel void setDiffusionVisit( const std::shared_ptr hierarchy); - void limitClDifference(const double delta); - private: void setAuxilliaryCompositions(); void registerPhaseConcentrationDiffusionVariables(); diff --git a/source/fortran/ConcFort.h b/source/fortran/ConcFort.h index 427022e9..51d23e8e 100644 --- a/source/fortran/ConcFort.h +++ b/source/fortran/ConcFort.h @@ -183,6 +183,17 @@ void ADD_FLUX_4TH(const int& ifirst0, const int& ilast0, const int& ifirst1, #endif const int& ngflux, const int* physb); +void ZERO_FLUX_PHYSB(const int& ifirst0, const int& ilast0, const int& ifirst1, + const int& ilast1, +#if (NDIM == 3) + const int& ifirst2, const int& ilast2, +#endif + const double* flux0, const double* flux1, +#if (NDIM == 3) + const double* flux2, +#endif + const int& ngflux, const int* physb); + void CONCENTRATION_FLUX_SPINODAL( const int& ifirst0, const int& ilast0, const int& ifirst1, const int& ilast1, From 40bcb1f25a183974925a75a025ed612e6ee1f21a Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Sun, 25 May 2025 16:04:07 -0400 Subject: [PATCH 6/6] More changes to stabilize computation of cl --- ...rgyBinaryMultiOrderThreePhasesStochioAB.cc | 18 ++++++++- ...ergyBinaryMultiOrderThreePhasesStochioAB.h | 7 +++- ...ntrationsMultiOrderThreePhasesStochioAB.cc | 28 ++++++------- source/FreeEnergyStrategyFactory.h | 6 ++- ...rBinaryThreePhasesDrivingForceStochioAB.cc | 39 ++++++++++++++++--- ...erBinaryThreePhasesDrivingForceStochioAB.h | 11 +++--- ...ergyMultiOrderBinaryThreePhaseStochioAB.cc | 3 +- source/QuatModel.cc | 11 +++--- 8 files changed, 83 insertions(+), 40 deletions(-) diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc index 09fa305d..d1e791c1 100644 --- a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc +++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc @@ -31,7 +31,7 @@ CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB:: std::shared_ptr newton_db, const Thermo4PFM::ConcInterpolationType conc_interp_func_type, MolarVolumeStrategy* mvstrategy, const int conc_l_id, - const int conc_a_id, const int conc_b_id) + const int conc_a_id, const int conc_b_id, const int conc_id) : CALPHADFreeEnergyBinaryMultiOrderThreePhases< Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase>( calphad_pt, newton_db, conc_interp_func_type, norderp_A, mvstrategy, @@ -42,7 +42,8 @@ CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB:: setup(calphad_pt, newton_db); d_multiorder_driving_force.reset( - new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, norderp_A)); + new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, conc_id, + norderp_A)); } //======================================================================= @@ -69,3 +70,16 @@ bool CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::computeCeqT( << std::endl; return d_ceq_fenergy->computeCeqT(temperature, &ceq[0], 50, true); } + +void CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::addDrivingForce( + const double time, hier::Patch& patch, const int temperature_id, + const int phase_id, const int conc_id, const int f_l_id, const int f_a_id, + const int f_b_id, const int rhs_id) +{ + (void)time; + + d_multiorder_driving_force->addDrivingForce(patch, temperature_id, phase_id, + d_conc_l_id, d_conc_a_id, + d_conc_b_id, f_l_id, f_a_id, + f_b_id, rhs_id); +} diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h index 2c363430..8b4d7256 100644 --- a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h +++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h @@ -37,11 +37,16 @@ class CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB std::shared_ptr newton_db, const Thermo4PFM::ConcInterpolationType conc_interp_func_type, MolarVolumeStrategy* mvstrategy, const int conc_l_id, - const int conc_a_id, const int conc_b_id); + const int conc_a_id, const int conc_b_id, const int conc_id); bool computeCeqT(const double temperature, const Thermo4PFM::PhaseIndex pi0, const Thermo4PFM::PhaseIndex pi1, double* ceq) override; + void addDrivingForce(const double time, hier::Patch& patch, + const int temperature_id, const int phase_id, + const int conc_id, const int f_l_id, const int f_a_id, + const int f_b_id, const int rhs_id) override; + private: void setup(boost::property_tree::ptree calphad_pt, std::shared_ptr newton_db); diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc index a22e23fb..6d21f4e5 100644 --- a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc +++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc @@ -50,8 +50,7 @@ int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: { assert(!std::isnan(x[0])); - const double epsilon1 = 1.e-3; - const double epsilon2 = 2.e-3; + const double epsilon = 1e-1; const double cA = d_model_parameters.getStochioA(); const double cB = d_model_parameters.getStochioB(); @@ -61,9 +60,15 @@ int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: assert(!std::isnan(xeq)); // solve explicit KKS problem - const double a = std::max(epsilon1, hphi[0]); - double xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / a; - xkks = xeq + (1. - xeq) * std::tanh(xkks - xeq); + const double a = hphi[0] + epsilon * std::exp(-hphi[0] / epsilon); + const double dkks = (c[0] - hphi[0] * xeq - hphi[1] * cA - hphi[2] * cB) / a; + + const double xmin = 0.7857; + const double xmax = 0.999; + const double d = dkks > 0. ? (xmax - xeq) : xeq - xmin; + const double t = std::tanh(dkks / d); + const double xkks = xeq + d * t; + #if 0 for (short i = 0; i < 3; i++) std::cerr << hphi[i] << ", "; @@ -74,18 +79,7 @@ int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB :: << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl; #endif - if (hphi[0] >= epsilon2) { - x[0] = xkks; - } else if (hphi[0] <= epsilon1) { - x[0] = xeq; - } else { - // map epsilon1 < hphi < epsilon2 to (0,1) - double h = (hphi[0] - epsilon1) / (epsilon2 - epsilon1); - // mix equilibrium compositions and KKS solution - double f = - Thermo4PFM::interp_func(Thermo4PFM::EnergyInterpolationType::PBG, h); - x[0] = (1. - f) * xeq + f * xkks; - } + x[0] = xkks; x[1] = cA; x[2] = cB; diff --git a/source/FreeEnergyStrategyFactory.h b/source/FreeEnergyStrategyFactory.h index 0a388daf..6fb29675 100644 --- a/source/FreeEnergyStrategyFactory.h +++ b/source/FreeEnergyStrategyFactory.h @@ -43,7 +43,8 @@ class FreeEnergyStrategyFactory static std::shared_ptr create( QuatModelParameters& model_parameters, const int ncompositions, const int conc_l_scratch_id, const int conc_a_scratch_id, - const int conc_b_scratch_id, MolarVolumeStrategy* mvstrategy, + const int conc_b_scratch_id, const int conc_scratch_id, + MolarVolumeStrategy* mvstrategy, MeltingTemperatureStrategy* meltingT_strategy, const double Tref, std::shared_ptr conc_db) { @@ -94,7 +95,8 @@ class FreeEnergyStrategyFactory newton_db, model_parameters.conc_interp_func_type(), mvstrategy, conc_l_scratch_id, - conc_a_scratch_id, conc_b_scratch_id)); + conc_a_scratch_id, conc_b_scratch_id, + conc_scratch_id)); } else if (model_parameters.getStochioB() >= 0.) { tbox::plog << "StochioB..." << std::endl; free_energy_strategy.reset( diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc index 65282e7f..334e0250 100644 --- a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc +++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc @@ -8,8 +8,11 @@ using namespace SAMRAI; MultiOrderBinaryThreePhasesDrivingForceStochioAB:: MultiOrderBinaryThreePhasesDrivingForceStochioAB( - FreeEnergyStrategyBinary* fenergy_strategy, const int norderp_A) - : d_fenergy_strategy(fenergy_strategy), d_norderp_A(norderp_A) + FreeEnergyStrategyBinary* fenergy_strategy, const int conc_id, + const int norderp_A) + : d_fenergy_strategy(fenergy_strategy), + d_conc_id(conc_id), + d_norderp_A(norderp_A) { } @@ -27,6 +30,7 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( assert(conc_a_id >= 0); assert(conc_b_id >= 0); assert(temperature_id >= 0); + assert(d_conc_id >= 0); std::shared_ptr > phase( SAMRAI_SHARED_PTR_CAST, hier::PatchData>( @@ -34,6 +38,11 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( assert(phase); assert(phase->getDepth() > 1); + std::shared_ptr > conc( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch.getPatchData(d_conc_id))); + assert(conc); + std::shared_ptr > temperature( SAMRAI_SHARED_PTR_CAST, hier::PatchData>( patch.getPatchData(temperature_id))); @@ -75,6 +84,8 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( assert(fl->getGhostCellWidth()[0] == fa->getGhostCellWidth()[0]); assert(cl->getGhostCellWidth()[0] == ca->getGhostCellWidth()[0]); + if (conc->getGhostCellWidth()[0] != phase->getGhostCellWidth()[0]) + std::cerr << "Wrong number of ghosts!!!" << std::endl; assert(phase->getDepth() > 1); assert(phase->getDepth() == rhs->getDepth()); @@ -88,6 +99,7 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( double* ptr_c_l = cl->getPointer(); double* ptr_c_a = ca->getPointer(); double* ptr_c_b = cb->getPointer(); + double* ptr_c = conc->getPointer(); const hier::Box& rhs_gbox = rhs->getGhostBox(); int imin_rhs = rhs_gbox.lower(0); @@ -211,15 +223,15 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( // interpolation polynomials double hphiA = 0.; for (short i = 0; i < d_norderp_A; i++) { - const double phiA = std::max(0., ptr_phi[i][idx_pf]); - hphiA += phiA * phiA; + const double phi = std::max(0., ptr_phi[i][idx_pf]); + hphiA += phi * phi; } assert(!std::isnan(hphiA)); double hphiB = 0.; for (short i = d_norderp_A; i < norderp - 1; i++) { - const double phiB = std::max(0., ptr_phi[i][idx_pf]); - hphiB += phiB * phiB; + const double phi = std::max(0., ptr_phi[i][idx_pf]); + hphiB += phi * phi; } assert(!std::isnan(hphiB)); @@ -238,6 +250,21 @@ void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce( assert(!std::isnan(hphiA)); assert(!std::isnan(hphiB)); + // muL ill-defined for hphil very small leads to instabilities + // replace driving force by stabilization term in this case + // that minimizes ||hphil*cl-hphiA*cA-hphiB*cB-c|| + // const double epsilon = 1.e-4; + // if (hphil < epsilon) { + //const double factor = 1.e4; + //double c = ptr_c[idx_pf]; + // dfl = 0.; + //double r = hphil * cl + hphiA * ca + hphiB * cb - c; + // dfl += factor * r * cl; + // dfa += factor * r * ca; + // dfb += factor * r * cb; + // std::cout<<"dfb = "<