The PSI phasespace is defined by Hammond-MLB equations, effective rad…

…ius is used for He+bubble reactions (#131).

test/core/network/PSINetworkTester.cpp

@@ -51,7 +51,8 @@ BOOST_AUTO_TEST_CASE(fullyRefined)
 	// Get the boundaries from the options
 	NetworkType::AmountType maxV = opts.getMaxV();
 	NetworkType::AmountType maxI = opts.getMaxI();
-	NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV, opts.getHeVRatio());
+	NetworkType::AmountType maxHe =
+		util::getMaxHePerV(maxV, opts.getHeVRatio());
 	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);
@@ -602,7 +603,8 @@ BOOST_AUTO_TEST_CASE(reducedMatrixMethod)
 	// Get the boundaries from the options
 	NetworkType::AmountType maxV = opts.getMaxV();
 	NetworkType::AmountType maxI = opts.getMaxI();
-	NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV, opts.getHeVRatio());
+	NetworkType::AmountType maxHe =
+		util::getMaxHePerV(maxV, opts.getHeVRatio());
 	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);
@@ -874,7 +876,8 @@ BOOST_AUTO_TEST_CASE(HeliumSpeciesList)
 	// Get the boundaries from the options
 	NetworkType::AmountType maxV = opts.getMaxV();
 	NetworkType::AmountType maxI = opts.getMaxI();
-	NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV, opts.getHeVRatio());
+	NetworkType::AmountType maxHe =
+		util::getMaxHePerV(maxV, opts.getHeVRatio());
 	NetworkType network({maxHe, maxV, maxI}, 1, opts);
 
 	BOOST_REQUIRE(!network.hasDeuterium());
@@ -1170,7 +1173,8 @@ BOOST_AUTO_TEST_CASE(DeuteriumSpeciesList)
 	// Get the boundaries from the options
 	NetworkType::AmountType maxV = opts.getMaxV();
 	NetworkType::AmountType maxI = opts.getMaxI();
-	NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV, opts.getHeVRatio());
+	NetworkType::AmountType maxHe =
+		util::getMaxHePerV(maxV, opts.getHeVRatio());
 	NetworkType::AmountType maxD = 2.0 / 3.0 * (double)maxHe;
 	NetworkType network({maxHe, maxD, maxV, maxI}, 1, opts);
 
@@ -1519,7 +1523,8 @@ BOOST_AUTO_TEST_CASE(TritiumSpeciesList)
 	// Get the boundaries from the options
 	NetworkType::AmountType maxV = opts.getMaxV();
 	NetworkType::AmountType maxI = opts.getMaxI();
-	NetworkType::AmountType maxHe = psi::getMaxHePerV(maxV, opts.getHeVRatio());
+	NetworkType::AmountType maxHe =
+		util::getMaxHePerV(maxV, opts.getHeVRatio());
 	NetworkType::AmountType maxT = 2.0 / 3.0 * (double)maxHe;
 	NetworkType network({maxHe, maxT, maxV, maxI}, 1, opts);
 

xolotl/core/include/xolotl/core/Constants.h

@@ -101,5 +101,19 @@ constexpr double basalTransitionSize = 325;
 // Tungsten density in nm^-3
 constexpr double tungstenDensity = 62.8;
 
+// 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 */

xolotl/core/include/xolotl/core/network/PSIClusterGenerator.h

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <xolotl/core/Constants.h>
 #include <xolotl/util/MathUtils.h>
 
 namespace xolotl
@@ -8,33 +9,6 @@ namespace core
 {
 namespace network
 {
-namespace psi
-{
-KOKKOS_INLINE_FUNCTION
-IReactionNetwork::AmountType
-getMaxHePerV(IReactionNetwork::AmountType amtV, double ratio) 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<AmountType, 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 (amtV < maxHePerV.size()) {
-		return maxHePerV[amtV];
-	}
-	return util::max((AmountType)(ratio * amtV),
-		maxHePerV[maxHePerV.size() - 1] + amtV - (AmountType)maxHePerV.size() +
-			1);
-}
-} // namespace psi
-
 template <typename TSpeciesEnum>
 class PSIClusterGenerator :
 	public plsm::refine::Detector<PSIClusterGenerator<TSpeciesEnum>>
@@ -108,7 +82,11 @@ class PSIClusterGenerator :
 	AmountType _groupingMin;
 	AmountType _groupingWidthA;
 	AmountType _groupingWidthB;
-	double _hevRatio{4.0};
+
+	// The temperature
+	double _temperature{933.0};
+	// The lattice parameter
+	double _lattice{xolotl::core::tungstenLatticeConstant};
 };
 
 // template <typename TSpeciesEnum>

xolotl/core/include/xolotl/core/network/impl/PSIClusterGenerator.tpp

@@ -21,7 +21,10 @@ PSIClusterGenerator<TSpeciesEnum>::PSIClusterGenerator(
 	_groupingMin(opts.getGroupingMin()),
 	_groupingWidthA(opts.getGroupingWidthA()),
 	_groupingWidthB(opts.getGroupingWidthB()),
-	_hevRatio(opts.getHeVRatio())
+	_temperature(opts.getTempParam()),
+	_lattice(opts.getLatticeParameter() <= 0.0 ?
+			xolotl::core::tungstenLatticeConstant :
+			opts.getLatticeParameter())
 {
 }
 
@@ -38,7 +41,10 @@ PSIClusterGenerator<TSpeciesEnum>::PSIClusterGenerator(
 	_groupingMin(opts.getGroupingMin()),
 	_groupingWidthA(opts.getGroupingWidthA()),
 	_groupingWidthB(opts.getGroupingWidthB()),
-	_hevRatio(opts.getHeVRatio())
+	_temperature(opts.getTempParam()),
+	_lattice(opts.getLatticeParameter() <= 0.0 ?
+			xolotl::core::tungstenLatticeConstant :
+			opts.getLatticeParameter())
 {
 }
 
@@ -49,9 +55,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;
@@ -137,18 +143,23 @@ PSIClusterGenerator<TSpeciesEnum>::refine(
 	}
 
 	// Else refine around the edge
-	auto maxDPerV = [hevRatio = _hevRatio, maxD = _maxD](AmountType amtV) {
-		return (2.0 / 3.0) * getMaxHePerV(amtV, hevRatio) * (maxD > 0);
+	auto maxDPerV = [lattice = _lattice, temperature = _temperature,
+						maxD = _maxD](AmountType amtV) {
+		return (2.0 / 3.0) * getMaxHePerVLoop(amtV, lattice, temperature) *
+			(maxD > 0);
 	};
-	auto maxTPerV = [hevRatio = _hevRatio, maxT = _maxT](AmountType amtV) {
-		return (2.0 / 3.0) * getMaxHePerV(amtV, hevRatio) * (maxT > 0);
+	auto maxTPerV = [lattice = _lattice, temperature = _temperature,
+						maxT = _maxT](AmountType amtV) {
+		return (2.0 / 3.0) * getMaxHePerVLoop(amtV, lattice, temperature) *
+			(maxT > 0);
 	};
 	if (hi[Species::V] > 1) {
 		double factor = 1.0e-1;
 
-		if (lo[Species::He] <= getMaxHePerV(hi[Species::V] - 1, _hevRatio) &&
+		if (lo[Species::He] <=
+				getMaxHePerVLoop(hi[Species::V] - 1, _lattice, _temperature) &&
 			hi[Species::He] - 1 >=
-				getMaxHePerV(lo[Species::V] - 1, _hevRatio)) {
+				getMaxHePerVLoop(lo[Species::V] - 1, _lattice, _temperature)) {
 			if constexpr (hasDeuterium<Species>) {
 				if (region[Species::D].length() <
 					util::max((double)(_groupingWidthA + 1),
@@ -287,9 +298,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) {
@@ -382,20 +393,21 @@ PSIClusterGenerator<TSpeciesEnum>::select(const Region& region) const
 		}
 	}
 
-	auto maxHPerV = [hevRatio = _hevRatio](AmountType amtV) {
-		return (2.0 / 3.0) * getMaxHePerV(amtV, hevRatio);
+	auto maxHPerV = [lattice = _lattice, temperature = _temperature](
+						AmountType amtV) {
+		return (2.0 / 3.0) * getMaxHePerVLoop(amtV, lattice, temperature);
 	};
 
 	// The edge
 	if (region[Species::V].end() > 1) {
 		Composition lo = region.getOrigin();
 		Composition hi = region.getUpperLimitPoint();
 		auto hiV = util::min(hi[Species::V] - 1, _maxV);
-		auto hiHe =
-			util::min(hi[Species::He] - 1, getMaxHePerV(_maxV, _hevRatio));
+		auto hiHe = util::min(hi[Species::He] - 1,
+			(AmountType)getMaxHePerVLoop(_maxV, _lattice, _temperature));
 
 		// Too many helium
-		if (lo[Species::He] > getMaxHePerV(hiV, _hevRatio)) {
+		if (lo[Species::He] > getMaxHePerVLoop(hiV, _lattice, _temperature)) {
 			return false;
 		}
 
@@ -673,6 +685,7 @@ PSIClusterGenerator<TSpeciesEnum>::getReactionRadius(
 {
 	const auto& reg = cluster.getRegion();
 	double radius = 0.0;
+
 	if (reg.isSimplex()) {
 		Composition comp(reg.getOrigin());
 		if (comp.isOnAxis(Species::I)) {