New solver StochioAB

source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc

@@ -0,0 +1,85 @@
+// 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 <boost/property_tree/json_parser.hpp>
+#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 <cassert>
+
+//=======================================================================
+
+CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::
+    CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB(
+        const int norderp_A, boost::property_tree::ptree calphad_pt,
+        std::shared_ptr<tbox::Database> 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_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, conc_id,
+                                                            norderp_A));
+}
+
+//=======================================================================
+
+void CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::setup(
+    pt::ptree calphad_pt, std::shared_ptr<tbox::Database> 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);
+}
+
+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);
+}

source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h

@@ -0,0 +1,61 @@
+// 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 <string>
+#include <vector>
+#include <boost/property_tree/ptree.hpp>
+
+class CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB
+    : public CALPHADFreeEnergyBinaryMultiOrderThreePhases<
+          Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase>
+{
+ public:
+   CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB(
+       const int norderp_A, boost::property_tree::ptree calphad_db,
+       std::shared_ptr<tbox::Database> 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_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<tbox::Database> newton_db);
+
+   std::shared_ptr<Thermo4PFM::CALPHADFreeEnergyFunctionsBinary> d_ceq_fenergy;
+
+
+   std::shared_ptr<MultiOrderBinaryThreePhasesDrivingForceStochioAB>
+       d_multiorder_driving_force;
+};
+
+#endif

source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc

@@ -0,0 +1,87 @@
+// 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 <boost/property_tree/json_parser.hpp>
+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<tbox::Database> 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<tbox::Database> conc_calphad_db =
+       conc_db->getDatabase("Calphad");
+   std::string calphad_filename = conc_calphad_db->getString("filename");
+
+   std::shared_ptr<tbox::MemoryDatabase> 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 epsilon = 1e-1;
+
+   const double cA = d_model_parameters.getStochioA();
+   const double cB = d_model_parameters.getStochioB();
+
+   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 = 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] << ", ";
+   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
+
+   x[0] = xkks;
+   x[1] = cA;
+   x[2] = cB;
+
+   return 0;
+}

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<tbox::Database> conc_db);
+
+   ~CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB() {}
+
+ private:
+   QuatModelParameters d_model_parameters;
+
+   int computePhaseConcentrations(double, double*, double*, double*) override;
+};
+
+#endif

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

source/FreeEnergyStrategyFactory.h

@@ -20,13 +20,15 @@
 #include "CALPHADFreeEnergyStrategyBinaryFolchPlapp.h"
 #include "CALPHADFreeEnergyBinaryMultiOrderThreePhases.h"
 #include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.h"
+#include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h"
 #include "CALPHADFreeEnergyStrategyBinaryFolchPlappStochioB.h"
 #include "ParabolicFreeEnergyBinary.h"
 #include "ParabolicFreeEnergyMultiOrderBinary.h"
 #include "QuadraticFreeEnergyBinary.h"
 #include "QuadraticFreeEnergyMultiOrderBinary.h"
 #include "QuadraticFreeEnergyMultiOrderTernaryThreePhase.h"
 #include "ParabolicFreeEnergyMultiOrderBinaryThreePhase.h"
+#include "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h"
 #include "KKSdiluteBinary.h"
 #include "BiasDoubleWellBeckermannFreeEnergyStrategy.h"
 #include "BiasDoubleWellUTRCFreeEnergyStrategy.h"
@@ -41,7 +43,8 @@ class FreeEnergyStrategyFactory
    static std::shared_ptr<FreeEnergyStrategy> 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<tbox::Database> conc_db)
    {
@@ -51,6 +54,8 @@ class FreeEnergyStrategyFactory
       std::shared_ptr<FreeEnergyStrategy> 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<tbox::MemoryDatabase> calphad_db;
@@ -82,7 +87,17 @@ 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,
+                                conc_scratch_id));
+                     } else if (model_parameters.getStochioB() >= 0.) {
                         tbox::plog << "StochioB..." << std::endl;
                         free_energy_strategy.reset(
                             new CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB(
@@ -208,13 +223,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;