From 3bad2aa37819042f9f9029c5017b6542a2e0bced Mon Sep 17 00:00:00 2001
From: Varun Shah <varun@Varuns-MacBook-Pro.local>
Date: Tue, 8 Aug 2023 19:56:00 +0200
Subject: [PATCH] Addition of DOF for dislocation density (#112)

---
 .../xolotl/core/network/detail/ClusterData.h  | 19 +++++++++++++++--
 .../core/network/impl/PSIReactionNetwork.tpp  |  5 ++++-
 xolotl/solver/src/monitor/PetscMonitor1D.cpp  | 21 +++++++++++--------
 3 files changed, 33 insertions(+), 12 deletions(-)

diff --git a/xolotl/core/include/xolotl/core/network/detail/ClusterData.h b/xolotl/core/include/xolotl/core/network/detail/ClusterData.h
index ec7223fea..1cbc4b60c 100644
--- a/xolotl/core/include/xolotl/core/network/detail/ClusterData.h
+++ b/xolotl/core/include/xolotl/core/network/detail/ClusterData.h
@@ -144,8 +144,9 @@ struct ClusterDataCommon
 	};
 
 	enum IntValsIndex : int
-	{
-		HELIUM_DISLO_ID = 0,
+	{ 
+	    DISLO_ID = 0,
+		HELIUM_DISLO_ID,
 		HELIUM_GB_ID,
 		NUM_INT_VALS
 	};
@@ -223,6 +224,20 @@ struct ClusterDataCommon
 	{
 		setVal(_floatVals, ZETA, val);
 	}
+	
+	KOKKOS_INLINE_FUNCTION
+	IndexType
+	DisloId() const
+	{
+		return _intVals[DISLO_ID];
+	}
+
+	void
+	setDisloId(IndexType val)
+	{
+		setVal(_intVals, DISLO_ID, val);
+	}
+
 
 	KOKKOS_INLINE_FUNCTION
 	IndexType
diff --git a/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp b/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp
index 589e2bcc8..79c0b9486 100644
--- a/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp
+++ b/xolotl/core/include/xolotl/core/network/impl/PSIReactionNetwork.tpp
@@ -65,10 +65,13 @@ void
 PSIReactionNetwork<TSpeciesEnum>::initializeExtraDOFs(
 	const options::IOptions& options)
 {   
+    this->_clusterData.h_view().setDisloId(this->_numDOFs);
+    this->_numDOFs += 1;
+    
 	if (!this->_enableSink) {
 		return;
 	}
-
+	
 	this->_clusterData.h_view().setHeliumDisloId(this->_numDOFs);
 	this->_clusterData.h_view().setHeliumGBId(this->_numDOFs + 1);
 	this->_numDOFs += 2;
diff --git a/xolotl/solver/src/monitor/PetscMonitor1D.cpp b/xolotl/solver/src/monitor/PetscMonitor1D.cpp
index 859244dc6..16ac98ff0 100644
--- a/xolotl/solver/src/monitor/PetscMonitor1D.cpp
+++ b/xolotl/solver/src/monitor/PetscMonitor1D.cpp
@@ -531,7 +531,7 @@ PetscMonitor1D::setup(int loop)
 					outputFile << speciesName << "_surface ";
 				}
 			}
-			outputFile << "Helium_burst Deuterium_burst Tritium_burst ";
+			outputFile << "Helium_burst Deuterium_burst Tritium_burst Disl_dens ";
 			if (_solverHandler->isSink())
 				outputFile << "He_dislo He_GB";
 			outputFile << std::endl;
@@ -981,7 +981,7 @@ PetscMonitor1D::computeHeliumRetention(
 	// Store the concentration over the grid
 	auto numSpecies = network.getSpeciesListSize();
 	auto specIdI = network.getInterstitialSpeciesId();
-	auto myConcData = std::vector<double>(numSpecies + 2 * isSink, 0.0);
+	auto myConcData = std::vector<double>(numSpecies + 1 + 2 * isSink, 0.0);
 
 	// Declare the pointer for the concentrations at a specific grid point
 	PetscReal* gridPointSolution;
@@ -1016,10 +1016,13 @@ PetscMonitor1D::computeHeliumRetention(
 		for (auto id = core::network::SpeciesId(numSpecies); id; ++id) {
 			myConcData[id()] += totals[id()] * hx;
 		}
+		// Summed over 1d profile as of now
+		myConcData[numSpecies] +=
+		    gridPointSolution[dof - 3] * hx;     // Dislocation density
 		if (isSink) {
-			myConcData[numSpecies] +=
-				gridPointSolution[dof - 2] * hx; // He dislo
 			myConcData[numSpecies + 1] +=
+				gridPointSolution[dof - 2] * hx; // He dislo
+			myConcData[numSpecies + 2] +=
 				gridPointSolution[dof - 1] * hx; // He GB
 		}
 	}
@@ -1030,9 +1033,9 @@ PetscMonitor1D::computeHeliumRetention(
 	MPI_Comm_rank(xolotlComm, &procId);
 
 	// Determine total concentrations for He, D, T.
-	auto totalConcData = std::vector<double>(numSpecies + 2 * isSink, 0.0);
+	auto totalConcData = std::vector<double>(numSpecies + 1 + 2 * isSink, 0.0);
 
-	MPI_Reduce(myConcData.data(), totalConcData.data(), numSpecies + 2 * isSink,
+	MPI_Reduce(myConcData.data(), totalConcData.data(), numSpecies + 1 + 2 * isSink,
 		MPI_DOUBLE, MPI_SUM, 0, xolotlComm);
 
 	// Get the delta time from the previous timestep to this timestep
@@ -1249,10 +1252,10 @@ PetscMonitor1D::computeHeliumRetention(
 		}
 		auto tempHandler = _solverHandler->getTemperatureHandler();
 		outputFile << _nHeliumBurst << " " << _nDeuteriumBurst << " "
-				   << _nTritiumBurst << " ";
+				   << _nTritiumBurst << " " << totalConcData[numSpecies] << " ";
 		if (isSink) {
-			outputFile << totalConcData[numSpecies] << " "
-					   << totalConcData[numSpecies + 1];
+			outputFile << totalConcData[numSpecies + 1] << " "
+					   << totalConcData[numSpecies + 2];
 		}
 		outputFile << std::endl;
 		outputFile.close();