Made the writing of the rates in the HDF5 file work in parallel (ORNL…

…-Fusion#185).

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

@@ -38,6 +38,7 @@ class IReactionNetwork
 	using SubMapView = Kokkos::View<AmountType*, Kokkos::MemoryUnmanaged>;
 	using OwnedSubMapView = Kokkos::View<AmountType*>;
 	using BelongingView = Kokkos::View<bool*>;
+	using RateConstantView = Kokkos::View<double**>;
 	using Connectivity = detail::ClusterConnectivity<>;
 	using SparseFillMap = std::unordered_map<int, std::vector<int>>;
 	using Bounds = std::vector<std::vector<AmountType>>;
@@ -482,10 +483,10 @@ class IReactionNetwork
 	virtual PhaseSpace
 	getPhaseSpace() = 0;
 
-	virtual std::vector<std::vector<double>>
+	virtual RateConstantView
 	getAllProdRates(IndexType gridIndex = 0) = 0;
 
-	virtual std::vector<std::vector<double>>
+	virtual RateConstantView
 	getAllDissoRates(IndexType gridIndex = 0) = 0;
 
 	/**

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

@@ -51,6 +51,7 @@ class Reaction
 	using ConnectivitiesPairView = IReactionNetwork::ConnectivitiesPairView;
 	using BelongingView = IReactionNetwork::BelongingView;
 	using OwnedSubMapView = IReactionNetwork::OwnedSubMapView;
+	using RateConstantView = IReactionNetwork::RateConstantView;
 	using Connectivity = typename IReactionNetwork::Connectivity;
 	using ReactionDataRef = typename Types::ReactionDataRef;
 	using ClusterData = typename Types::ClusterData;
@@ -105,10 +106,11 @@ class Reaction
 	}
 
 	KOKKOS_INLINE_FUNCTION
-	std::vector<double>
-	contributeRateVector(IndexType gridIndex)
+	void
+	contributeRateVector(
+		RateConstantView toReturn, IndexType id, IndexType gridIndex)
 	{
-		return asDerived()->computeRateVector(gridIndex);
+		return asDerived()->computeRateVector(toReturn, id, gridIndex);
 	}
 
 	KOKKOS_INLINE_FUNCTION
@@ -365,6 +367,7 @@ class ProductionReaction : public Reaction<TNetwork, TDerived>
 	using ConnectivitiesPairView = typename Superclass::ConnectivitiesPairView;
 	using BelongingView = typename Superclass::BelongingView;
 	using OwnedSubMapView = typename Superclass::OwnedSubMapView;
+	using RateConstantView = typename Superclass::RateConstantView;
 	using Composition = typename Superclass::Composition;
 	using Region = typename Superclass::Region;
 	using AmountType = typename Superclass::AmountType;
@@ -417,8 +420,9 @@ class ProductionReaction : public Reaction<TNetwork, TDerived>
 	computeReducedConnectivity(const Connectivity& connectivity);
 
 	KOKKOS_INLINE_FUNCTION
-	std::vector<double>
-	computeRateVector(IndexType gridIndex);
+	void
+	computeRateVector(
+		RateConstantView toReturn, IndexType id, IndexType gridIndex);
 
 	KOKKOS_INLINE_FUNCTION
 	void
@@ -544,6 +548,7 @@ class DissociationReaction : public Reaction<TNetwork, TDerived>
 	using ConnectivitiesPairView = typename Superclass::ConnectivitiesPairView;
 	using BelongingView = typename Superclass::BelongingView;
 	using OwnedSubMapView = typename Superclass::OwnedSubMapView;
+	using RateConstantView = typename Superclass::RateConstantView;
 	using AmountType = typename Superclass::AmountType;
 	using ReactionDataRef = typename Superclass::ReactionDataRef;
 	using ClusterData = typename Superclass::ClusterData;
@@ -593,8 +598,9 @@ class DissociationReaction : public Reaction<TNetwork, TDerived>
 	computeReducedConnectivity(const Connectivity& connectivity);
 
 	KOKKOS_INLINE_FUNCTION
-	std::vector<double>
-	computeRateVector(IndexType gridIndex);
+	void
+	computeRateVector(
+		RateConstantView toReturn, IndexType id, IndexType gridIndex);
 
 	KOKKOS_INLINE_FUNCTION
 	void

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

@@ -80,6 +80,7 @@ class ReactionNetwork : public ReactionNetworkInterface<TImpl>::Type
 	using SubMapView = typename IReactionNetwork::SubMapView;
 	using OwnedSubMapView = typename IReactionNetwork::OwnedSubMapView;
 	using BelongingView = typename IReactionNetwork::BelongingView;
+	using RateConstantView = typename IReactionNetwork::RateConstantView;
 	using SparseFillMap = typename IReactionNetwork::SparseFillMap;
 	using ClusterData = typename Types::ClusterData;
 	using ClusterDataMirror = typename Types::ClusterDataMirror;
@@ -510,10 +511,10 @@ class ReactionNetwork : public ReactionNetworkInterface<TImpl>::Type
 	{
 	}
 
-	std::vector<std::vector<double>>
+	RateConstantView
 	getAllProdRates(IndexType gridIndex = 0) final;
 
-	std::vector<std::vector<double>>
+	RateConstantView
 	getAllDissoRates(IndexType gridIndex = 0) final;
 
 	void

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

@@ -682,13 +682,14 @@ ProductionReaction<TNetwork, TDerived>::computeReducedConnectivity(
 
 template <typename TNetwork, typename TDerived>
 KOKKOS_INLINE_FUNCTION
-std::vector<double>
-ProductionReaction<TNetwork, TDerived>::computeRateVector(IndexType gridIndex)
+void
+ProductionReaction<TNetwork, TDerived>::computeRateVector(
+	RateConstantView toReturn, IndexType id, IndexType gridIndex)
 {
-	std::vector<double> toReturn = {_reactants[0], _reactants[1], _products[0],
-		this->_coefs(0, 0, 0, 0) * this->_rate(gridIndex)};
-
-	return toReturn;
+	toReturn(id, 0) = _reactants[0];
+	toReturn(id, 1) = _reactants[1];
+	toReturn(id, 2) = _products[0];
+	toReturn(id, 3) = this->_coefs(0, 0, 0, 0) * this->_rate(gridIndex);
 }
 
 template <typename TNetwork, typename TDerived>
@@ -2557,13 +2558,14 @@ DissociationReaction<TNetwork, TDerived>::computeReducedConnectivity(
 
 template <typename TNetwork, typename TDerived>
 KOKKOS_INLINE_FUNCTION
-std::vector<double>
-DissociationReaction<TNetwork, TDerived>::computeRateVector(IndexType gridIndex)
+void
+DissociationReaction<TNetwork, TDerived>::computeRateVector(
+	RateConstantView toReturn, IndexType id, IndexType gridIndex)
 {
-	std::vector<double> toReturn = {_reactant, _products[0], _products[1],
-		this->_coefs(0, 0, 0, 0) * this->_rate(gridIndex)};
-
-	return toReturn;
+	toReturn(id, 0) = _reactant;
+	toReturn(id, 1) = _products[0];
+	toReturn(id, 2) = _products[1];
+	toReturn(id, 3) = this->_coefs(0, 0, 0, 0) * this->_rate(gridIndex);
 }
 
 template <typename TNetwork, typename TDerived>

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

@@ -544,34 +544,51 @@ ReactionNetwork<TImpl>::generateClusterData(const ClusterGenerator& generator)
 }
 
 template <typename TImpl>
-std::vector<std::vector<double>>
+Kokkos::View<double**>
 ReactionNetwork<TImpl>::getAllProdRates(IndexType gridIndex)
 {
-	std::vector<std::vector<double>> toReturn;
-	auto vecPtr = &toReturn;
-
-	_reactions.template forEachOn<typename Traits::ProductionReactionType>(
-		"ReactionNetwork::getAllRates", DEVICE_LAMBDA(auto&& reaction) {
-			vecPtr->push_back(reaction.contributeRateVector(gridIndex));
+	// Create the views at the correct side and for the device
+	auto prodReactions =
+		_reactions.template getView<typename Traits::ProductionReactionType>();
+	int nReactions = prodReactions.size();
+	RateConstantView toReturn =
+		Kokkos::View<double**>("Production Rates", nReactions, 4);
+	auto dToReturn = create_mirror_view(toReturn);
+
+	Kokkos::parallel_for(
+		nReactions, KOKKOS_LAMBDA(const IndexType i) {
+			prodReactions(i).contributeRateVector(dToReturn, i, gridIndex);
 		});
 	Kokkos::fence();
 
+	// Copy the data back for the host
+	deep_copy(toReturn, dToReturn);
+
 	return toReturn;
 }
 
 template <typename TImpl>
-std::vector<std::vector<double>>
+Kokkos::View<double**>
 ReactionNetwork<TImpl>::getAllDissoRates(IndexType gridIndex)
 {
-	std::vector<std::vector<double>> toReturn;
-	auto vecPtr = &toReturn;
-
-	_reactions.template forEachOn<typename Traits::DissociationReactionType>(
-		"ReactionNetwork::getAllRates", DEVICE_LAMBDA(auto&& reaction) {
-			vecPtr->push_back(reaction.contributeRateVector(gridIndex));
+	// Create the views at the correct side and for the device
+	auto dissoReactions =
+		_reactions
+			.template getView<typename Traits::DissociationReactionType>();
+	int nReactions = dissoReactions.size();
+	RateConstantView toReturn =
+		Kokkos::View<double**>("Dissociation Rates", nReactions, 4);
+	auto dToReturn = create_mirror_view(toReturn);
+
+	Kokkos::parallel_for(
+		nReactions, KOKKOS_LAMBDA(const IndexType i) {
+			dissoReactions(i).contributeRateVector(dToReturn, i, gridIndex);
 		});
 	Kokkos::fence();
 
+	// Copy the data back for the host
+	deep_copy(toReturn, dToReturn);
+
 	return toReturn;
 }
 

xolotl/core/src/flux/FluxHandler.cpp

@@ -133,8 +133,12 @@ FluxHandler::syncFluxIndices()
 void
 FluxHandler::syncIncidentFluxVec()
 {
-	incidentFlux = Kokkos::View<double**>(
-		"Incident Flux Vec", incidentFluxVec.size(), incidentFluxVec[0].size());
+	if (incidentFluxVec.size() == 0) {
+		incidentFlux = Kokkos::View<double**>();
+	}
+	else
+		incidentFlux = Kokkos::View<double**>("Incident Flux Vec",
+			incidentFluxVec.size(), incidentFluxVec[0].size());
 	auto incidentFlux_h = create_mirror_view(incidentFlux);
 	for (std::size_t i = 0; i < incidentFluxVec.size(); ++i) {
 		for (std::size_t j = 0; j < incidentFluxVec[i].size(); ++j) {

xolotl/solver/src/monitor/PetscMonitor.cpp

@@ -260,30 +260,30 @@ PetscMonitor::startStop(TS ts, PetscInt timestep, PetscReal time, Vec solution)
 			timestep, time, previousTime, currentTimeStep);
 
 		// Production reactions
-		auto rateVector = network.getAllProdRates();
-		auto size = rateVector.size();
+		auto rateView = network.getAllProdRates();
+		auto size = rateView.extent(0);
 
 		// Make it an array
 		double rateArray[size][4];
 		for (auto i = 0; i < size; i++) {
-			rateArray[i][0] = rateVector[i][0];
-			rateArray[i][1] = rateVector[i][1];
-			rateArray[i][2] = rateVector[i][2];
-			rateArray[i][3] = rateVector[i][3];
+			rateArray[i][0] = rateView(i, 0);
+			rateArray[i][1] = rateView(i, 1);
+			rateArray[i][2] = rateView(i, 2);
+			rateArray[i][3] = rateView(i, 3);
 		}
 
 		tsGroup->writeReactionDataset(size, rateArray, true, 0);
 
 		// Dissociation reactions
-		rateVector = network.getAllDissoRates();
-		size = rateVector.size();
+		rateView = network.getAllDissoRates();
+		size = rateView.extent(0);
 
 		// Make it an array
 		for (auto i = 0; i < size; i++) {
-			rateArray[i][0] = rateVector[i][0];
-			rateArray[i][1] = rateVector[i][1];
-			rateArray[i][2] = rateVector[i][2];
-			rateArray[i][3] = rateVector[i][3];
+			rateArray[i][0] = rateView(i, 0);
+			rateArray[i][1] = rateView(i, 1);
+			rateArray[i][2] = rateView(i, 2);
+			rateArray[i][3] = rateView(i, 3);
 		}
 
 		tsGroup->writeReactionDataset(size, rateArray, true, 1);