Skip to content

Commit

Permalink
Updating the He/V ratio to MLB-EOS for loop punching (#80, #116).
Browse files Browse the repository at this point in the history
  • Loading branch information
Sophie Blondel committed Mar 4, 2024
1 parent 6d02126 commit 86b1ca3
Show file tree
Hide file tree
Showing 10 changed files with 209 additions and 66 deletions.
75 changes: 59 additions & 16 deletions test/core/network/PSINetworkTester.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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);
Expand Down Expand Up @@ -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();

Expand All @@ -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);
Expand Down Expand Up @@ -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
Expand All @@ -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());
Expand Down Expand Up @@ -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
Expand All @@ -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);

Expand Down Expand Up @@ -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
Expand All @@ -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);

Expand Down Expand Up @@ -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
Expand All @@ -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());
Expand Down
14 changes: 14 additions & 0 deletions xolotl/core/include/xolotl/core/Constants.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
31 changes: 5 additions & 26 deletions xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double, 30> 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 <typename TSpeciesEnum>
class PSIClusterGenerator :
Expand Down Expand Up @@ -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
Expand Down
14 changes: 14 additions & 0 deletions xolotl/core/include/xolotl/core/network/detail/ClusterData.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,7 @@ struct ClusterDataCommon
DEPTH,
TAU_BURSTING,
F_BURSTING,
TEMPERATURE,
NUM_FLOAT_VALS
};

Expand Down Expand Up @@ -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
{
Expand Down
46 changes: 31 additions & 15 deletions xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,11 @@ PSIClusterGenerator<TSpeciesEnum>::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())
{
}

Expand All @@ -38,7 +42,11 @@ PSIClusterGenerator<TSpeciesEnum>::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())
{
}

Expand All @@ -49,9 +57,9 @@ PSIClusterGenerator<TSpeciesEnum>::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;
Expand Down Expand Up @@ -150,17 +158,23 @@ PSIClusterGenerator<TSpeciesEnum>::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<Species>) {
if (region[Species::D].length() <
util::max((double)(_groupingWidthA + 1),
Expand Down Expand Up @@ -302,9 +316,9 @@ KOKKOS_INLINE_FUNCTION
bool
PSIClusterGenerator<TSpeciesEnum>::select(const Region& region) const
{
using psi::getMaxHePerV;
using psi::hasDeuterium;
using psi::hasTritium;
using util::getMaxHePerVLoop;

// Remove 0
auto isZeroPoint = [](const Region& reg) {
Expand Down Expand Up @@ -401,8 +415,9 @@ PSIClusterGenerator<TSpeciesEnum>::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();
Expand All @@ -411,16 +426,17 @@ PSIClusterGenerator<TSpeciesEnum>::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)) {
return false;
}

// Too many helium
if (lo[Species::He] > (AmountType)getMaxHePerV(hiV)) {
if (lo[Species::He] >
(AmountType)util::getMaxHePerVLoop(hiV, _lattice, _temperature)) {
return false;
}

Expand Down
Loading

0 comments on commit 86b1ca3

Please sign in to comment.