Saving the left rate information for clusters with sink terms (#112).

xolotl/core/include/xolotl/core/network/PSIReaction.h

@@ -62,6 +62,7 @@ class PSISinkReaction :
 		PSISinkReaction<TSpeciesEnum>>;
 
 	using Superclass::Superclass;
+	using IndexType = typename Superclass::IndexType;
 
 	KOKKOS_INLINE_FUNCTION
 	double
@@ -70,6 +71,10 @@ class PSISinkReaction :
 	KOKKOS_INLINE_FUNCTION
 	double
 	getSinkStrength();
+
+	KOKKOS_INLINE_FUNCTION
+	double
+	computeRate(IndexType gridIndex, double time = 0.0);
 };
 template <typename TSpeciesEnum>
 class PSITrapMutationReaction :

xolotl/core/include/xolotl/core/network/PSITraits.h

@@ -296,12 +296,19 @@ struct ClusterDataExtra<PSIReactionNetwork<TSpeciesEnum>, MemSpace>
 {
 	using NetworkType = PSIReactionNetwork<TSpeciesEnum>;
 
+	template <typename TData>
+	using View = ViewType<TData, MemSpace>;
+
+	using IndexType = detail::ReactionNetworkIndexType;
+
 	ClusterDataExtra() = default;
 
 	template <typename MS>
 	KOKKOS_INLINE_FUNCTION
 	ClusterDataExtra(const ClusterDataExtra<NetworkType, MS>& data) :
-		trapMutationData(data.trapMutationData)
+		trapMutationData(data.trapMutationData),
+		leftSideRates(data.leftSideRates),
+		sinkMap(data.sinkMap)
 	{
 	}
 
@@ -310,17 +317,45 @@ struct ClusterDataExtra<PSIReactionNetwork<TSpeciesEnum>, MemSpace>
 	deepCopy(const ClusterDataExtra<NetworkType, MS>& data)
 	{
 		trapMutationData.deepCopy(data.trapMutationData);
+
+		if (!data.leftSideRates.is_allocated()) {
+			return;
+		}
+
+		if (!leftSideRates.is_allocated()) {
+			leftSideRates = create_mirror_view(data.leftSideRates);
+		}
+		deep_copy(leftSideRates, data.leftSideRates);
+
+		if (!sinkMap.is_allocated()) {
+			sinkMap = create_mirror_view(data.sinkMap);
+		}
+		deep_copy(sinkMap, data.sinkMap);
 	}
 
 	std::uint64_t
 	getDeviceMemorySize() const noexcept
 	{
-		return trapMutationData.getDeviceMemorySize();
+		std::uint64_t ret = 0;
+		ret += trapMutationData.getDeviceMemorySize();
+		ret += leftSideRates.required_allocation_size(leftSideRates.extent(0));
+		ret += sinkMap.required_allocation_size(sinkMap.extent(0));
+		return ret;
+	}
+
+	void
+	initialize(IndexType numClusters)
+	{
+		leftSideRates = View<double*>("Left Side Rates", numClusters);
+		sinkMap = View<IndexType*>("Sink Map", numClusters);
 	}
 
 	using TrapMutationData =
 		TrapMutationClusterData<ClusterDataCommon<MemSpace>>;
 	TrapMutationData trapMutationData;
+
+	View<double*> leftSideRates;
+	View<IndexType*> sinkMap;
 };
 } // namespace detail
 } // namespace network

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

@@ -64,6 +64,13 @@ class SinkReaction : public Reaction<TNetwork, TDerived>
 	double
 	computeRate(IndexType gridIndex, double time = 0.0);
 
+	KOKKOS_INLINE_FUNCTION
+	IndexType
+	getReactantId()
+	{
+		return _reactant;
+	}
+
 private:
 	KOKKOS_INLINE_FUNCTION
 	void

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

@@ -237,9 +237,34 @@ PSISinkReaction<TSpeciesEnum>::getSinkStrength()
 {
 	constexpr double pi = ::xolotl::core::pi;
 	double grainSize = 50000.0; // 50 um
-
 	return 1.0 / (pi * grainSize * grainSize);
 }
+
+template <typename TSpeciesEnum>
+KOKKOS_INLINE_FUNCTION
+double
+PSISinkReaction<TSpeciesEnum>::computeRate(IndexType gridIndex, double time)
+{
+	auto cl = this->_clusterData->getCluster(this->_reactant);
+	double dc = cl.getDiffusionCoefficient(gridIndex);
+
+	auto rate = this->_clusterData->extraData.leftSideRates;
+	double s_m = 0.0;
+	if (rate.size() > 0) {
+		// Look for the correct index
+		auto sMap = this->_clusterData->extraData.sinkMap;
+		auto i = 0;
+		for (i; i < sMap.size(); i++) {
+			if (sMap(i) == this->_reactant)
+				break;
+		}
+		s_m = rate(i) / dc;
+	}
+
+	double strength = this->getSinkBias() * this->getSinkStrength() * dc;
+
+	return strength;
+}
 } // namespace network
 } // namespace core
 } // namespace xolotl

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

@@ -65,48 +65,66 @@ void
 PSIReactionNetwork<TSpeciesEnum>::updateExtraClusterData(
 	const std::vector<double>& gridTemps, const std::vector<double>& gridDepths)
 {
-	if (!this->_enableTrapMutation) {
-		return;
-	}
-
-	// Check which temperature index to use
-	IdType tempId = 0;
-	for (tempId = 0; tempId < gridDepths.size(); tempId++) {
-		if (gridDepths[tempId] > 0.01)
-			break;
-	}
+	if (this->_enableTrapMutation) {
+		// Check which temperature index to use
+		IdType tempId = 0;
+		for (tempId = 0; tempId < gridDepths.size(); tempId++) {
+			if (gridDepths[tempId] > 0.01)
+				break;
+		}
 
-	_tmHandler->updateData(gridTemps[tempId]);
+		_tmHandler->updateData(gridTemps[tempId]);
 
-	auto& tmData = this->_clusterData.h_view().extraData.trapMutationData;
+		auto& tmData = this->_clusterData.h_view().extraData.trapMutationData;
 
-	using Kokkos::HostSpace;
-	using Kokkos::MemoryUnmanaged;
+		using Kokkos::HostSpace;
+		using Kokkos::MemoryUnmanaged;
+
+		auto desorpInit = _tmHandler->getDesorptionInitializer();
+		auto subpaving = this->_subpaving;
+		IndexType desorpId = this->invalidIndex();
+		Kokkos::parallel_reduce(
+			1,
+			KOKKOS_LAMBDA(std::size_t, IndexType & running) {
+				Composition comp{};
+				comp[Species::He] = desorpInit.size;
+				running = static_cast<IndexType>(subpaving.findTileId(comp));
+			},
+			desorpId);
+		auto desorp = create_mirror_view(tmData.desorption);
+		desorp() = detail::Desorption{desorpInit, desorpId};
+		deep_copy(tmData.desorption, desorp);
+
+		auto depths = Kokkos::View<const double[7], HostSpace, MemoryUnmanaged>(
+			_tmHandler->getDepths().data());
+		deep_copy(tmData.tmDepths, depths);
+
+		auto vSizes =
+			Kokkos::View<const AmountType[7], HostSpace, MemoryUnmanaged>(
+				_tmHandler->getVacancySizes().data());
+		deep_copy(tmData.tmVSizes, vSizes);
+
+		this->invalidateDataMirror();
+	}
 
-	auto desorpInit = _tmHandler->getDesorptionInitializer();
-	auto subpaving = this->_subpaving;
-	IndexType desorpId = this->invalidIndex();
-	Kokkos::parallel_reduce(
-		1,
-		KOKKOS_LAMBDA(std::size_t, IndexType & running) {
-			Composition comp{};
-			comp[Species::He] = desorpInit.size;
-			running = static_cast<IndexType>(subpaving.findTileId(comp));
-		},
-		desorpId);
-	auto desorp = create_mirror_view(tmData.desorption);
-	desorp() = detail::Desorption{desorpInit, desorpId};
-	deep_copy(tmData.desorption, desorp);
+	if (this->_enableSink) {
+		using SinkReactionType = typename Superclass::Traits::SinkReactionType;
+		auto sinkReactions =
+			this->_reactions.template getView<SinkReactionType>();
+		this->_clusterData.h_view().extraData.initialize(sinkReactions.size());
 
-	auto depths = Kokkos::View<const double[7], HostSpace, MemoryUnmanaged>(
-		_tmHandler->getDepths().data());
-	deep_copy(tmData.tmDepths, depths);
+		this->copyClusterDataView();
 
-	auto vSizes = Kokkos::View<const AmountType[7], HostSpace, MemoryUnmanaged>(
-		_tmHandler->getVacancySizes().data());
-	deep_copy(tmData.tmVSizes, vSizes);
+		auto& clusterData = this->_clusterData.d_view;
+		Kokkos::parallel_for(
+			"PSIReactionNetwork::updateExtraClusterData", sinkReactions.size(),
+			KOKKOS_LAMBDA(IndexType i) {
+				clusterData().extraData.sinkMap(i) =
+					sinkReactions(i).getReactantId();
+			});
 
-	this->invalidateDataMirror();
+		this->invalidateDataMirror();
+	}
 }
 
 template <typename TSpeciesEnum>
@@ -144,6 +162,24 @@ PSIReactionNetwork<TSpeciesEnum>::computeFluxesPreProcess(
 		updateDesorptionLeftSideRate(concentrations, gridIndex);
 		selectTrapMutationReactions(surfaceDepth, spacing);
 	}
+	if (this->_enableSink) {
+		// Compute the left side rates
+		auto& clusterData = this->_clusterData.d_view;
+		Kokkos::parallel_for(
+			"PSIReactionNetwork::computeFluxesPreProcess",
+			this->_clusterData.h_view().extraData.sinkMap.size(),
+			KOKKOS_LAMBDA(IndexType i) {
+				clusterData().extraData.leftSideRates(i) = 0.0;
+			});
+
+		// Update the sink rate
+		using SinkReactionType = typename Superclass::Traits::SinkReactionType;
+		auto sinkReactions =
+			this->_reactions.template getView<SinkReactionType>();
+		Kokkos::parallel_for(
+			"PSIReactionNetwork::computeFluxesPreProcess", sinkReactions.size(),
+			KOKKOS_LAMBDA(IndexType i) { sinkReactions[i].updateRates(); });
+	}
 }
 
 template <typename TSpeciesEnum>
@@ -156,6 +192,26 @@ PSIReactionNetwork<TSpeciesEnum>::computePartialsPreProcess(
 		updateDesorptionLeftSideRate(concentrations, gridIndex);
 		selectTrapMutationReactions(surfaceDepth, spacing);
 	}
+	if (this->_enableSink) {
+		// Compute the left side rates
+		auto& clusterData = this->_clusterData.d_view;
+		Kokkos::parallel_for(
+			"PSIReactionNetwork::computeFluxesPreProcess",
+			this->_clusterData.h_view().extraData.sinkMap.size(),
+			KOKKOS_LAMBDA(IndexType i) {
+				clusterData().extraData.leftSideRates(i) =
+					this->getLeftSideRate(concentrations,
+						clusterData().extraData.sinkMap(i), gridIndex);
+			});
+
+		// Update the sink rate
+		using SinkReactionType = typename Superclass::Traits::SinkReactionType;
+		auto sinkReactions =
+			this->_reactions.template getView<SinkReactionType>();
+		Kokkos::parallel_for(
+			"PSIReactionNetwork::computeFluxesPreProcess", sinkReactions.size(),
+			KOKKOS_LAMBDA(IndexType i) { sinkReactions[i].updateRates(); });
+	}
 }
 
 template <typename TSpeciesEnum>
@@ -285,7 +341,6 @@ void
 PSIReactionNetwork<TSpeciesEnum>::updateReactionRates(double time)
 {
 	Superclass::updateReactionRates(time);
-
 	using TrapMutationReactionType =
 		typename Superclass::Traits::TrapMutationReactionType;
 	// TODO: is this just the local largest rate? Is it correct?

xolotl/options/include/xolotl/options/IOptions.h

@@ -132,13 +132,13 @@ class IOptions
 	virtual std::string
 	getPerfHandlerName() const = 0;
 
-    /**
-     * Should we write the performance report to a YAML file?
-     *
-     * @return true to enable YAML output
-     */
-    virtual bool
-    usePerfOutputYAML() const = 0;
+	/**
+	 * Should we write the performance report to a YAML file?
+	 *
+	 * @return true to enable YAML output
+	 */
+	virtual bool
+	usePerfOutputYAML() const = 0;
 
 	/**
 	 * Obtain the name of the visualization handler to be used

xolotl/options/include/xolotl/options/Options.h

@@ -68,10 +68,10 @@ class Options : public IOptions
 	 */
 	std::string perfHandlerName;
 
-    /**
-     * Output performance report to YAML file?
-     */
-    bool perfOutputYAMLFlag;
+	/**
+	 * Output performance report to YAML file?
+	 */
+	bool perfOutputYAMLFlag;
 
 	/**
 	 * Name of the viz handler
@@ -406,14 +406,14 @@ class Options : public IOptions
 		return perfHandlerName;
 	}
 
-    /**
-     * \see IOptions.h
-     */
-    bool
-    usePerfOutputYAML() const override
-    {
-        return perfOutputYAMLFlag;
-    }
+	/**
+	 * \see IOptions.h
+	 */
+	bool
+	usePerfOutputYAML() const override
+	{
+		return perfOutputYAMLFlag;
+	}
 
 	/**
 	 * \see IOptions.h