Fixed the FeCr flux and reaction rates to reproduce SPICES (#100).

xolotl/core/include/xolotl/core/flux/FeCrFitFluxHandler.h

@@ -232,32 +232,35 @@ class FeCrFitFluxHandler : public FluxHandler
 	computeIncidentFlux(
 		double currentTime, double* updatedConcOffset, int xi, int surfacePos)
 	{
-		updatedConcOffset[fluxIndices[0]] += fluxAmplitude * 11.7769; // V1
-		updatedConcOffset[fluxIndices[1]] += fluxAmplitude * 3.6300379; // V2
-		updatedConcOffset[fluxIndices[2]] += fluxAmplitude * 0.94389166; // V3
-		updatedConcOffset[fluxIndices[3]] += fluxAmplitude * 0.76565337; // V4
-		updatedConcOffset[fluxIndices[4]] += fluxAmplitude * 0.38057040; // V5
-		updatedConcOffset[fluxIndices[5]] += fluxAmplitude * 0.34432560; // V9
-		updatedConcOffset[fluxIndices[6]] += fluxAmplitude * 5.9403415; // I1
-		updatedConcOffset[fluxIndices[7]] += fluxAmplitude * 4.3737840; // I2
-		updatedConcOffset[fluxIndices[8]] += fluxAmplitude * 2.1009049; // I3
-		updatedConcOffset[fluxIndices[9]] += fluxAmplitude * 1.7910009; // Free4
+		updatedConcOffset[fluxIndices[0]] += fluxAmplitude * 11.776939094; // V1
+		updatedConcOffset[fluxIndices[1]] += fluxAmplitude * 1.815018938; // V2
+		updatedConcOffset[fluxIndices[2]] +=
+			fluxAmplitude * 0.31463055266666667; // V3
+		updatedConcOffset[fluxIndices[3]] += fluxAmplitude * 0.191413343; // V4
+		updatedConcOffset[fluxIndices[4]] += fluxAmplitude * 0.07611408; // V5
+		updatedConcOffset[fluxIndices[5]] += fluxAmplitude * 0.0382584; // V9
+		updatedConcOffset[fluxIndices[6]] += fluxAmplitude * 5.940341464; // I1
+		updatedConcOffset[fluxIndices[7]] += fluxAmplitude * 2.186892002; // I2
+		updatedConcOffset[fluxIndices[8]] +=
+			fluxAmplitude * 0.70030164666666661; // I3
+		updatedConcOffset[fluxIndices[9]] +=
+			fluxAmplitude * 0.447750235; // Free4
 		updatedConcOffset[fluxIndices[10]] +=
-			fluxAmplitude * 1.3093007; // Free5
+			fluxAmplitude * 0.261860134; // Free5
 		updatedConcOffset[fluxIndices[11]] +=
-			fluxAmplitude * 0.96875166; // Free6
+			fluxAmplitude * 0.16145861; // Free6
 		updatedConcOffset[fluxIndices[12]] +=
-			fluxAmplitude * 0.66476811; // Free7
+			fluxAmplitude * 0.094966872285714293; // Free7
 		updatedConcOffset[fluxIndices[13]] +=
-			fluxAmplitude * 0.18484848; // Free8
+			fluxAmplitude * 0.02310606; // Free8
 		updatedConcOffset[fluxIndices[14]] +=
-			fluxAmplitude * 0.23226550; // Free9
+			fluxAmplitude * 0.025807277555555556; // Free9
 		updatedConcOffset[fluxIndices[15]] +=
-			fluxAmplitude * 0.11597512; // Free12
+			fluxAmplitude * 0.0096645931666666674; // Free12
 		updatedConcOffset[fluxIndices[16]] +=
-			fluxAmplitude * 0.011597512; // Free16
+			fluxAmplitude * 0.00611631; // Free16
 		updatedConcOffset[fluxIndices[17]] +=
-			fluxAmplitude * 0.061616160; // Free20
+			fluxAmplitude * 0.0030808080; // Free20
 
 		return;
 	}

xolotl/core/include/xolotl/core/network/ConstantReaction.h

@@ -311,6 +311,14 @@ class ConstantReaction : public Reaction<TNetwork, TDerived>
 		return 0.0;
 	}
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity)

xolotl/core/include/xolotl/core/network/FeCrReaction.h

@@ -24,6 +24,11 @@ class FeCrProductionReaction :
 	KOKKOS_INLINE_FUNCTION
 	double
 	getRateForProduction(IndexType gridIndex);
+
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex);
 };
 
 class FeCrDissociationReaction :
@@ -86,6 +91,11 @@ class FeCrSinkReaction :
 	KOKKOS_INLINE_FUNCTION
 	double
 	getSinkStrength();
+
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex);
 };
 
 class FeCrTransformReaction :

xolotl/core/include/xolotl/core/network/FeCrReactionNetwork.h

@@ -1,8 +1,8 @@
 #pragma once
 
+#include <xolotl/core/network/ReactionNetwork.h>
 #include <xolotl/core/network/FeCrReaction.h>
 #include <xolotl/core/network/FeCrTraits.h>
-#include <xolotl/core/network/ReactionNetwork.h>
 
 namespace xolotl
 {
@@ -36,6 +36,19 @@ class FeCrReactionNetwork : public ReactionNetwork<FeCrReactionNetwork>
 	IndexType
 	checkLargestClusterId();
 
+	void
+	initializeExtraClusterData(const options::IOptions& options);
+
+	void
+	computeFluxesPreProcess(ConcentrationsView concentrations,
+		FluxesView fluxes, IndexType gridIndex, double surfaceDepth,
+		double spacing);
+
+	void
+	computePartialsPreProcess(ConcentrationsView concentrations,
+		Kokkos::View<double*> values, IndexType gridIndex, double surfaceDepth,
+		double spacing);
+
 private:
 	double
 	checkLatticeParameter(double latticeParameter);

xolotl/core/include/xolotl/core/network/FeCrTraits.h

@@ -90,6 +90,61 @@ struct ReactionNetworkTraits<FeCrReactionNetwork>
 
 	using ClusterGenerator = FeCrClusterGenerator;
 };
+
+namespace detail
+{
+template <typename PlsmContext>
+struct ClusterDataExtra<FeCrReactionNetwork, PlsmContext>
+{
+	using NetworkType = FeCrReactionNetwork;
+
+	template <typename TData>
+	using View = ViewType<TData, PlsmContext>;
+
+	using IndexType = detail::ReactionNetworkIndexType;
+
+	ClusterDataExtra() = default;
+
+	template <typename PC>
+	KOKKOS_INLINE_FUNCTION
+	ClusterDataExtra(const ClusterDataExtra<NetworkType, PC>& data) :
+		netSigma(data.netSigma)
+	{
+	}
+
+	template <typename PC>
+	void
+	deepCopy(const ClusterDataExtra<NetworkType, PC>& data)
+	{
+		if (!data.netSigma.is_allocated()) {
+			return;
+		}
+
+		if (!netSigma.is_allocated()) {
+			netSigma = create_mirror_view(data.netSigma);
+		}
+		deep_copy(netSigma, data.netSigma);
+	}
+
+	std::uint64_t
+	getDeviceMemorySize() const noexcept
+	{
+		std::uint64_t ret = 0;
+
+		ret += netSigma.required_allocation_size(netSigma.extent(0));
+
+		return ret;
+	}
+
+	void
+	initialize(IndexType numClusters)
+	{
+		netSigma = View<double*>("Net Sigma", numClusters);
+	}
+
+	View<double*> netSigma;
+};
+} // namespace detail
 } // namespace network
 } // namespace core
 } // namespace xolotl

xolotl/core/include/xolotl/core/network/IReactionNetwork.h

@@ -418,6 +418,13 @@ class IReactionNetwork
 	getLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex) = 0;
 
+	/**
+	 * @brief Returns the sum of absorption cross sections.
+	 */
+	virtual double
+	getNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex) = 0;
+
 	/**
 	 * Get the diagonal fill for the Jacobian, corresponding to the reactions.
 	 * Also populates the inverse map.

xolotl/core/include/xolotl/core/network/NucleationReaction.h

@@ -113,6 +113,14 @@ class NucleationReaction : public Reaction<TNetwork, TDerived>
 		return 0.0;
 	}
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity);

xolotl/core/include/xolotl/core/network/ReSolutionReaction.h

@@ -108,6 +108,14 @@ class ReSolutionReaction : public Reaction<TNetwork, TDerived>
 	computeLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex);
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity);

xolotl/core/include/xolotl/core/network/Reaction.h

@@ -162,6 +162,18 @@ class Reaction
 			concentrations, clusterId, gridIndex);
 	}
 
+	/**
+	 * \see IReactionNetwork.h
+	 */
+	KOKKOS_INLINE_FUNCTION
+	double
+	contributeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return asDerived()->computeNetSigma(
+			concentrations, clusterId, gridIndex);
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	defineJacobianEntries(Connectivity connectivity)
@@ -337,6 +349,11 @@ class ProductionReaction : public Reaction<TNetwork, TDerived>
 	computeLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex);
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex);
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity);
@@ -452,6 +469,14 @@ class DissociationReaction : public Reaction<TNetwork, TDerived>
 	computeLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex);
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity);

xolotl/core/include/xolotl/core/network/ReactionNetwork.h

@@ -472,6 +472,10 @@ class ReactionNetwork : public ReactionNetworkInterface<TImpl>::Type
 	getLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex) override;
 
+	double
+	getNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex) override;
+
 	IndexType
 	getDiagonalFill(SparseFillMap& fillMap) override;
 

xolotl/core/include/xolotl/core/network/SinkReaction.h

@@ -139,6 +139,11 @@ class SinkReaction : public Reaction<TNetwork, TDerived>
 		return 0.0;
 	}
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex);
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity)

xolotl/core/include/xolotl/core/network/TransformReaction.h

@@ -153,6 +153,14 @@ class TransformReaction : public Reaction<TNetwork, TDerived>
 		return 0.0;
 	}
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity)

xolotl/core/include/xolotl/core/network/TrapMutationReaction.h

@@ -139,6 +139,14 @@ class TrapMutationReaction : public Reaction<TNetwork, TDerived>
 	computeLeftSideRate(ConcentrationsView concentrations, IndexType clusterId,
 		IndexType gridIndex);
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeNetSigma(ConcentrationsView concentrations, IndexType clusterId,
+		IndexType gridIndex)
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	void
 	mapJacobianEntries(Connectivity connectivity)

xolotl/core/include/xolotl/core/network/ZrReaction.h

@@ -51,6 +51,20 @@ class ZrSinkReaction : public SinkReaction<ZrReactionNetwork, ZrSinkReaction>
 
 	using Superclass::Superclass;
 
+	KOKKOS_INLINE_FUNCTION
+	double
+	getSinkBias()
+	{
+		return 1.0;
+	}
+
+	KOKKOS_INLINE_FUNCTION
+	double
+	getSinkStrength()
+	{
+		return 0.0;
+	}
+
 	KOKKOS_INLINE_FUNCTION
 	double
 	computeRate(IndexType gridIndex, double time = 0.0);