From 86b1ca30ef0bb87b3431ecf137225961fc9b9b5c Mon Sep 17 00:00:00 2001 From: Sophie Blondel Date: Mon, 4 Mar 2024 10:53:17 -0500 Subject: [PATCH] Updating the He/V ratio to MLB-EOS for loop punching (#80, #116). --- test/core/network/PSINetworkTester.cpp | 75 +++++++++++++++---- xolotl/core/include/xolotl/core/Constants.h | 14 ++++ .../xolotl/core/network/PSIClusterGenerator.h | 31 ++------ .../xolotl/core/network/detail/ClusterData.h | 14 ++++ .../core/network/impl/PSIClusterGenerator.tpp | 46 ++++++++---- .../xolotl/core/network/impl/PSIReaction.tpp | 18 +++-- .../core/network/impl/PSIReactionNetwork.tpp | 12 ++- .../core/network/impl/ReactionNetwork.tpp | 1 + xolotl/core/src/network/PSINetworkHandler.cpp | 8 +- xolotl/util/include/xolotl/util/MathUtils.h | 56 ++++++++++++++ 10 files changed, 209 insertions(+), 66 deletions(-) diff --git a/test/core/network/PSINetworkTester.cpp b/test/core/network/PSINetworkTester.cpp index 4cf7274ab..4e57cb3de 100644 --- a/test/core/network/PSINetworkTester.cpp +++ b/test/core/network/PSINetworkTester.cpp @@ -34,7 +34,8 @@ BOOST_AUTO_TEST_CASE(fullyRefined) std::string parameterFile = "param.txt"; std::ofstream paramFile(parameterFile); paramFile << "netParam=8 1 1 1 1" << std::endl - << "process=reaction" << std::endl; + << "process=reaction" << std::endl + << "tempParam=1000" << std::endl; paramFile.close(); // Create a fake command line to read the options @@ -47,11 +48,17 @@ BOOST_AUTO_TEST_CASE(fullyRefined) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; - // Get the boundaries from the options + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType::AmountType maxD = 2.0 / 3.0 * (double)maxHe; NetworkType::AmountType maxT = 2.0 / 3.0 * (double)maxHe; NetworkType network({maxHe, maxD, maxT, maxV, maxI}, 1, opts); @@ -581,6 +588,7 @@ BOOST_AUTO_TEST_CASE(reducedMatrixMethod) std::ofstream paramFile(parameterFile); paramFile << "netParam=8 1 1 1 1" << std::endl << "process=reaction" << std::endl + << "tempParam=1000" << std::endl << "petscArgs=-snes_mf_operator" << std::endl; paramFile.close(); @@ -594,11 +602,17 @@ BOOST_AUTO_TEST_CASE(reducedMatrixMethod) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; - // Get the boundaries from the options + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType::AmountType maxD = 2.0 / 3.0 * (double)maxHe; NetworkType::AmountType maxT = 2.0 / 3.0 * (double)maxHe; NetworkType network({maxHe, maxD, maxT, maxV, maxI}, 1, opts); @@ -853,7 +867,8 @@ BOOST_AUTO_TEST_CASE(HeliumSpeciesList) std::string parameterFile = "param.txt"; std::ofstream paramFile(parameterFile); paramFile << "netParam=8 0 0 2 2" << std::endl - << "process=reaction" << std::endl; + << "process=reaction" << std::endl + << "tempParam=1000" << std::endl; paramFile.close(); // Create a fake command line to read the options @@ -866,11 +881,17 @@ BOOST_AUTO_TEST_CASE(HeliumSpeciesList) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; - // Get the boundaries from the options + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType network({maxHe, maxV, maxI}, 1, opts); BOOST_REQUIRE(!network.hasDeuterium()); @@ -1145,7 +1166,8 @@ BOOST_AUTO_TEST_CASE(DeuteriumSpeciesList) std::string parameterFile = "param.txt"; std::ofstream paramFile(parameterFile); paramFile << "netParam=8 1 0 1 1" << std::endl - << "process=reaction" << std::endl; + << "process=reaction" << std::endl + << "tempParam=1000" << std::endl; paramFile.close(); // Create a fake command line to read the options @@ -1158,11 +1180,17 @@ BOOST_AUTO_TEST_CASE(DeuteriumSpeciesList) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; - // Get the boundaries from the options + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType::AmountType maxD = 2.0 / 3.0 * (double)maxHe; NetworkType network({maxHe, maxD, maxV, maxI}, 1, opts); @@ -1490,7 +1518,8 @@ BOOST_AUTO_TEST_CASE(TritiumSpeciesList) std::string parameterFile = "param.txt"; std::ofstream paramFile(parameterFile); paramFile << "netParam=8 0 1 1 1" << std::endl - << "process=reaction" << std::endl; + << "process=reaction" << std::endl + << "tempParam=1000" << std::endl; paramFile.close(); // Create a fake command line to read the options @@ -1503,11 +1532,17 @@ BOOST_AUTO_TEST_CASE(TritiumSpeciesList) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; - // Get the boundaries from the options + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType::AmountType maxT = 2.0 / 3.0 * (double)maxHe; NetworkType network({maxHe, maxT, maxV, maxI}, 1, opts); @@ -2040,7 +2075,8 @@ BOOST_AUTO_TEST_CASE(LargeBubble) std::string parameterFile = "param.txt"; std::ofstream paramFile(parameterFile); paramFile << "netParam=8 0 0 5 2" << std::endl - << "process=reaction largeBubble" << std::endl; + << "process=reaction largeBubble" << std::endl + << "tempParam=1000" << std::endl; paramFile.close(); // Create a fake command line to read the options @@ -2053,10 +2089,17 @@ BOOST_AUTO_TEST_CASE(LargeBubble) using Spec = NetworkType::Species; using Composition = NetworkType::Composition; + // Get the temperature and lattice constant from the options + double latticeConst = opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter(); + double temperature = opts.getTempParam(); + // Get the boundaries from the options NetworkType::AmountType maxV = opts.getMaxV(); NetworkType::AmountType maxI = opts.getMaxI(); - NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV); + NetworkType::AmountType maxHe = + util::getMaxHePerVLoop(maxV, latticeConst, temperature); NetworkType network({maxHe, maxV, maxI}, 1, opts); BOOST_REQUIRE(!network.hasDeuterium()); diff --git a/xolotl/core/include/xolotl/core/Constants.h b/xolotl/core/include/xolotl/core/Constants.h index e9f574b46..42cac5c4d 100644 --- a/xolotl/core/include/xolotl/core/Constants.h +++ b/xolotl/core/include/xolotl/core/Constants.h @@ -105,5 +105,19 @@ constexpr double tungstenDensity = 62.8; // V value for which the grouping behavior changes for PSI constexpr IdType psiVThreshold = 150; +// Parameters for W Benedict EOS +constexpr double gammaEOS = 2.67; +constexpr double gEOS = 156.1e9; +constexpr double kEOS = 1.380649e-23; +constexpr double A_11 = -5.5991e-27; +constexpr double A01 = 1.74e-26; +constexpr double A21 = 4.9833e-30; +constexpr double A02 = -4.4658e-24; +constexpr double A22 = -8.7825e-27; +constexpr double A03 = 1.7595e-22; +constexpr double A23 = 1.7608e-23; +constexpr double A_13 = -3.2615e-21; +constexpr double A_23 = 3.1524e-20; + } /* end namespace core */ } /* end namespace xolotl */ diff --git a/xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h b/xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h index aa00d232f..db4423cf1 100644 --- a/xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h +++ b/xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h @@ -8,32 +8,6 @@ namespace core { namespace network { -namespace psi -{ -KOKKOS_INLINE_FUNCTION -double -getMaxHePerV(double amtV) noexcept -{ - using AmountType = IReactionNetwork::AmountType; - - /** - * The maximum number of helium atoms that can be combined with a - * vacancy cluster with size equal to the index i. - * It could support a mixture of up to nine - * helium atoms with one vacancy. - */ - constexpr Kokkos::Array maxHePerV = {0, 9, 14, 18, 20, 27, 30, - 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 98, 100, 101, 103, - 105, 107, 109, 110, 112, 116}; - - if ((int)amtV < maxHePerV.size()) { - return maxHePerV[(int)amtV]; - } - return util::max((4.0 * amtV), - (double)maxHePerV[maxHePerV.size() - 1] + amtV - - (double)maxHePerV.size() + 1.0); -} -} // namespace psi template class PSIClusterGenerator : @@ -105,6 +79,11 @@ class PSIClusterGenerator : AmountType _groupingMin; AmountType _groupingWidthA; AmountType _groupingWidthB; + + // The temperature + double _temperature{933.0}; + // The lattice parameter + double _lattice{xolotl::core::tungstenLatticeConstant}; }; } // namespace network } // namespace core diff --git a/xolotl/core/include/xolotl/core/network/detail/ClusterData.h b/xolotl/core/include/xolotl/core/network/detail/ClusterData.h index 075e00463..6693ca987 100644 --- a/xolotl/core/include/xolotl/core/network/detail/ClusterData.h +++ b/xolotl/core/include/xolotl/core/network/detail/ClusterData.h @@ -151,6 +151,7 @@ struct ClusterDataCommon DEPTH, TAU_BURSTING, F_BURSTING, + TEMPERATURE, NUM_FLOAT_VALS }; @@ -356,6 +357,19 @@ struct ClusterDataCommon setVal(_floatVals, F_BURSTING, val); } + KOKKOS_INLINE_FUNCTION + double + getTemperature() const + { + return _floatVals[TEMPERATURE]; + } + + void + setTemperature(double val) + { + setVal(_floatVals, TEMPERATURE, val); + } + int transitionSize() const { diff --git a/xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp b/xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp index 6fc1f4264..36207a51e 100644 --- a/xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp @@ -21,7 +21,11 @@ PSIClusterGenerator::PSIClusterGenerator( _maxI(opts.getMaxI()), _groupingMin(opts.getGroupingMin()), _groupingWidthA(opts.getGroupingWidthA()), - _groupingWidthB(opts.getGroupingWidthB()) + _groupingWidthB(opts.getGroupingWidthB()), + _temperature(opts.getTempParam()), + _lattice(opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter()) { } @@ -38,7 +42,11 @@ PSIClusterGenerator::PSIClusterGenerator( _maxI(opts.getMaxI()), _groupingMin(opts.getGroupingMin()), _groupingWidthA(opts.getGroupingWidthA()), - _groupingWidthB(opts.getGroupingWidthB()) + _groupingWidthB(opts.getGroupingWidthB()), + _temperature(opts.getTempParam()), + _lattice(opts.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + opts.getLatticeParameter()) { } @@ -49,9 +57,9 @@ PSIClusterGenerator::refine( const Region& region, BoolArray& result) const { using detail::toIndex; - using psi::getMaxHePerV; using psi::hasDeuterium; using psi::hasTritium; + using util::getMaxHePerVLoop; for (auto& r : result) { r = true; @@ -150,17 +158,23 @@ PSIClusterGenerator::refine( } // Else refine around the edge - auto maxDPerV = [maxD = _maxD](AmountType amtV) { - return (2.0 / 3.0) * getMaxHePerV(amtV) * (maxD > 0); + auto maxDPerV = [lattice = _lattice, temperature = _temperature, + maxD = _maxD](AmountType amtV) { + return (2.0 / 3.0) * + util::getMaxHePerVLoop(amtV, lattice, temperature) * (maxD > 0); }; - auto maxTPerV = [maxT = _maxT](AmountType amtV) { - return (2.0 / 3.0) * getMaxHePerV(amtV) * (maxT > 0); + auto maxTPerV = [lattice = _lattice, temperature = _temperature, + maxT = _maxT](AmountType amtV) { + return (2.0 / 3.0) * + util::getMaxHePerVLoop(amtV, lattice, temperature) * (maxT > 0); }; if (hi[Species::V] > 1 and lo[Species::V] <= _maxV) { double factor = 1.0e-1; - if (lo[Species::He] <= getMaxHePerV(hi[Species::V] - 1) && - hi[Species::He] - 1 >= getMaxHePerV(lo[Species::V] - 1)) { + if (lo[Species::He] <= util::getMaxHePerVLoop(hi[Species::V] - 1, + _lattice, _temperature) && + hi[Species::He] - 1 >= util::getMaxHePerVLoop(lo[Species::V] - 1, + _lattice, _temperature)) { if constexpr (hasDeuterium) { if (region[Species::D].length() < util::max((double)(_groupingWidthA + 1), @@ -302,9 +316,9 @@ KOKKOS_INLINE_FUNCTION bool PSIClusterGenerator::select(const Region& region) const { - using psi::getMaxHePerV; using psi::hasDeuterium; using psi::hasTritium; + using util::getMaxHePerVLoop; // Remove 0 auto isZeroPoint = [](const Region& reg) { @@ -401,8 +415,9 @@ PSIClusterGenerator::select(const Region& region) const } } - auto maxHPerV = [](AmountType amtV) { - return (2.0 / 3.0) * getMaxHePerV(amtV); + auto maxHPerV = [lattice = _lattice, temperature = _temperature]( + AmountType amtV) { + return (2.0 / 3.0) * util::getMaxHePerVLoop(amtV, lattice, temperature); }; Composition lo = region.getOrigin(); @@ -411,8 +426,8 @@ PSIClusterGenerator::select(const Region& region) const // The edge if (region[Species::V].end() > 1) { auto hiV = util::min(hi[Species::V] - 1, _maxV); - auto hiHe = - util::min(hi[Species::He] - 1, (AmountType)getMaxHePerV(_maxV)); + auto hiHe = util::min(hi[Species::He] - 1, + (AmountType)util::getMaxHePerVLoop(_maxV, _lattice, _temperature)); // No impurities above maxV if (lo[Species::V] > _maxV and !lo.isOnAxis(Species::V)) { @@ -420,7 +435,8 @@ PSIClusterGenerator::select(const Region& region) const } // Too many helium - if (lo[Species::He] > (AmountType)getMaxHePerV(hiV)) { + if (lo[Species::He] > + (AmountType)util::getMaxHePerVLoop(hiV, _lattice, _temperature)) { return false; } diff --git a/xolotl/core/include/xolotl/core/network/impl/PSIReaction.tpp b/xolotl/core/include/xolotl/core/network/impl/PSIReaction.tpp index d5e7b992e..c4c47e270 100644 --- a/xolotl/core/include/xolotl/core/network/impl/PSIReaction.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/PSIReaction.tpp @@ -95,8 +95,11 @@ PSIProductionReaction::computeFlux( auto avV = this->_clusterData->bubbleAvV(); // Sigmoid - double x = - ((comp[Species::He] + avHe) / psi::getMaxHePerV(avV)) - 1.0; + double x = ((comp[Species::He] + avHe) / + util::getMaxHePerVLoop(avV, + this->_clusterData->latticeParameter(), + cl.getTemperature(gridIndex))) - + 1.0; double sigmo = 1.0 / (1.0 + std::exp(-50.0 * x)); double nI = 1.0; @@ -289,8 +292,11 @@ PSIProductionReaction::computePartialDerivatives( auto avV = this->_clusterData->bubbleAvV(); // Sigmoid - double x = - ((comp[Species::He] + avHe) / psi::getMaxHePerV(avV)) - 1.0; + double x = ((comp[Species::He] + avHe) / + util::getMaxHePerVLoop(avV, + this->_clusterData->latticeParameter(), + cl.getTemperature(gridIndex))) - + 1.0; double sigmo = 1.0 / (1.0 + std::exp(-50.0 * x)); double nI = 1.0; @@ -731,8 +737,10 @@ PSIDissociationReaction::computeBindingEnergy(double time) Composition comp(clReg.getOrigin()); AmountType lowerV = 16, higherV = 31; AmountType minV = 1; + for (auto i = 1; i < higherV; i++) { - auto maxHe = psi::getMaxHePerV(i); + auto maxHe = util::getMaxHePerVLoop(i, + this->_clusterData->latticeParameter(), cl.getTemperature(0)); if (comp[Species::He] > maxHe) minV = i; } diff --git a/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp b/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp index 2ee24bff5..40f346de4 100644 --- a/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp @@ -752,7 +752,9 @@ PSIReactionGenerator::addLargeBubbleReactions( this->getCluster(this->largestClusterId).getRegion(); Composition hiLargest = largestReg.getUpperLimitPoint(); auto largestV = hiLargest[Species::V] - 1; - auto largestHe = psi::getMaxHePerV(largestV); + auto largestHe = + util::getMaxHePerVLoop(largestV, this->_clusterData.latticeParameter(), + this->_clusterData.getTemperature()); // General case constexpr auto species = NetworkType::getSpeciesRange(); @@ -780,10 +782,14 @@ PSIReactionGenerator::addLargeBubbleReactions( if (lo1.isOnAxis(Species::He) or lo2.isOnAxis(Species::He)) { // Is it trap mutation? if (bounds[Species::He].first > - psi::getMaxHePerV(bounds[Species::V].first)) { + util::getMaxHePerVLoop(bounds[Species::V].first, + this->_clusterData.latticeParameter(), + this->getCluster(i).getTemperature(0))) { AmountType iSize = 1; while (bounds[Species::He].first > - psi::getMaxHePerV(bounds[Species::V].first + iSize)) { + util::getMaxHePerVLoop(bounds[Species::V].first + iSize, + this->_clusterData.latticeParameter(), + this->getCluster(i).getTemperature(0))) { iSize++; } // Get the corresponding I cluster diff --git a/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp b/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp index 9b5f81519..bcb092108 100644 --- a/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp @@ -50,6 +50,7 @@ ReactionNetwork::ReactionNetwork(const Subpaving& subpaving, this->setTauBursting(opts.getBurstingDepth()); this->setFBursting(opts.getBurstingFactor()); _clusterData.h_view().setTransitionSize(opts.getTransitionSize()); + _clusterData.h_view().setTemperature(opts.getTempParam(0)); auto map = opts.getProcesses(); this->setEnableStdReaction(map["reaction"]); this->setEnableReSolution(map["resolution"]); diff --git a/xolotl/core/src/network/PSINetworkHandler.cpp b/xolotl/core/src/network/PSINetworkHandler.cpp index c0f5ed277..1e03c17f2 100644 --- a/xolotl/core/src/network/PSINetworkHandler.cpp +++ b/xolotl/core/src/network/PSINetworkHandler.cpp @@ -28,11 +28,17 @@ generatePSIReactionNetwork(const options::IOptions& options) { using AmountType = IReactionNetwork::AmountType; + // Get the temperature and lattice constant from the options + double latticeConst = options.getLatticeParameter() <= 0.0 ? + xolotl::core::tungstenLatticeConstant : + options.getLatticeParameter(); + double temperature = options.getTempParam(); + // Get the boundaries from the options AmountType maxV = options.getMaxPureV(); AmountType maxHeV = options.getMaxV(); AmountType maxI = options.getMaxI(); - AmountType maxHe = psi::getMaxHePerV(options.getMaxV()); + AmountType maxHe = util::getMaxHePerVLoop(maxV, latticeConst, temperature); AmountType maxD = 2.0 / 3.0 * (double)maxHe; AmountType maxT = 2.0 / 3.0 * (double)maxHe; AmountType groupingWidthHe = options.getGroupingWidthA(); diff --git a/xolotl/util/include/xolotl/util/MathUtils.h b/xolotl/util/include/xolotl/util/MathUtils.h index 05c2c3349..3530e6b1c 100644 --- a/xolotl/util/include/xolotl/util/MathUtils.h +++ b/xolotl/util/include/xolotl/util/MathUtils.h @@ -16,6 +16,8 @@ #include +#include + namespace xolotl { namespace util @@ -71,6 +73,60 @@ equal(double a, double b) return std::fabs(b - a) < epsilon; } +KOKKOS_INLINE_FUNCTION +double +computeBenedictF1(const double temp) noexcept +{ + return xolotl::core::A_11 / sqrt(temp) + xolotl::core::A01 + + xolotl::core::A21 * temp; +} + +KOKKOS_INLINE_FUNCTION +double +computeBenedictF2(const double temp) noexcept +{ + return xolotl::core::A02 + xolotl::core::A22 * temp; +} + +KOKKOS_INLINE_FUNCTION +double +computeBenedictF3(const double temp) noexcept +{ + return xolotl::core::A_23 / temp + xolotl::core::A_13 / sqrt(temp) + + xolotl::core::A03 + xolotl::core::A23 * temp; +} + +KOKKOS_INLINE_FUNCTION +double +getMaxHePerVLoop(double amtV, double latticeParameter, double temp) noexcept +{ + constexpr Kokkos::Array maxHePerV = { + 0, 9, 14, 18, 20, 27, 30, 35, 40, 45, 50, 55, 60, 65, 70}; + + if (amtV < maxHePerV.size()) { + return maxHePerV[(int)amtV]; + } + + double omega = + 0.5 * latticeParameter * latticeParameter * latticeParameter * 1.0e-27; + double bEOS = latticeParameter * 0.5 * sqrt(3.0) * 1.0e-9; + double nume = amtV * omega; + double Fs = pow(3.0 * omega / (4.0 * xolotl::core::pi), 1.0 / 3.0); + double term1 = computeBenedictF1(temp) * pow(Fs, 1.0 / 3.0) * + pow((double)amtV, 1.0 / 9.0) * + pow(2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS, + -1.0 / 3.0); + double term2 = computeBenedictF2(temp) * pow(Fs, 2.0 / 3.0) * + pow((double)amtV, 2.0 / 9.0) * + pow(2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS, + -2.0 / 3.0); + double term3 = computeBenedictF3(temp) * Fs * pow((double)amtV, 1.0 / 3.0) / + (2.0 * xolotl::core::gammaEOS + xolotl::core::gEOS * bEOS); + + return util::max((nume / (term1 + term2 + term3)), + maxHePerV[maxHePerV.size() - 1] + amtV - maxHePerV.size() + 1.0); +} + /** * This operation computes the Legendre polynomial, P_n (x), of degree n. *