Tracking the defect flux to sink (#100).

xolotl/solver/include/xolotl/solver/handler/ISolverHandler.h

@@ -267,6 +267,22 @@ class ISolverHandler
 	virtual double
 	getHeVRatio() const = 0;
 
+	/**
+	 * Get the sink density.
+	 *
+	 * @return The density
+	 */
+	virtual double
+	getSinkDensity() const = 0;
+
+	/**
+	 * Get the sink portion.
+	 *
+	 * @return The portion
+	 */
+	virtual double
+	getSinkPortion() const = 0;
+
 	/**
 	 * Get the grid left offset.
 	 *

xolotl/solver/include/xolotl/solver/handler/SolverHandler.h

@@ -154,6 +154,12 @@ class SolverHandler : public ISolverHandler
 	//! The ratio of He per V in a bubble.
 	double heVRatio;
 
+	//! The sink density.
+	double sinkDensity;
+
+	//! The sink portion.
+	double sinkPortion;
+
 	//! The value to use to seed the random number generator.
 	unsigned int rngSeed;
 
@@ -397,6 +403,24 @@ class SolverHandler : public ISolverHandler
 		return burstingFactor;
 	}
 
+	/**
+	 * \see ISolverHandler.h
+	 */
+	double
+	getSinkDensity() const override
+	{
+		return sinkDensity;
+	}
+
+	/**
+	 * \see ISolverHandler.h
+	 */
+	double
+	getSinkPortion() const override
+	{
+		return sinkPortion;
+	}
+
 	/**
 	 * \see ISolverHandler.h
 	 */

xolotl/solver/include/xolotl/solver/monitor/PetscMonitor0D.h

@@ -49,6 +49,9 @@ class PetscMonitor0D : public PetscMonitor
 	std::shared_ptr<viz::IPlot> _scatterPlot;
 
 	std::vector<IdType> _clusterOrder;
+
+	std::vector<double> _nSink;
+	std::vector<double> _previousSinkRate;
 };
 } // namespace monitor
 } // namespace solver

xolotl/solver/src/handler/SolverHandler.cpp

@@ -51,6 +51,8 @@ SolverHandler::SolverHandler(
 	burstingFactor(0.1),
 	rngSeed(0),
 	heVRatio(4.0),
+	sinkDensity(0.0),
+	sinkPortion(0.0),
 	previousTime(0.0),
 	nXeGB(0.0)
 {
@@ -448,6 +450,12 @@ SolverHandler::initializeHandlers(core::material::IMaterialHandler* material,
 	// Set the HeV ratio
 	heVRatio = opts.getHeVRatio();
 
+	// Set the sink density
+	sinkDensity = opts.getSinkDensity();
+
+	// Set the sink portion
+	sinkPortion = opts.getSinkPortion();
+
 	// Which type of temperature grid to use
 	if (opts.getTempHandlerName() == "heat")
 		sameTemperatureGrid = false;

xolotl/solver/src/monitor/PetscMonitor0D.cpp

@@ -234,9 +234,16 @@ PetscMonitor0D::setup()
 					   << "_partial_density " << speciesName << "_diameter "
 					   << speciesName << "_partial_diameter ";
 		}
-		outputFile << "black_dot_density" << std::endl;
+		outputFile
+			<< "111_partial_density 111_partial_diameter black_dot_density "
+			   "V_edge V_screw I_edge I_screw"
+			<< std::endl;
 		outputFile.close();
 
+		// Create data to track sink losses
+		_nSink = std::vector<double>(4, 0.0);
+		_previousSinkRate = std::vector<double>(4, 0.0);
+
 		// computeFeCrAl0D will be called at each timestep
 		ierr = TSMonitorSet(_ts, monitor::computeFeCrAl, this, nullptr);
 		checkPetscError(
@@ -765,7 +772,8 @@ PetscMonitor0D::computeFeCrAl(
 	auto& network = dynamic_cast<NetworkType&>(_solverHandler->getNetwork());
 	const auto dof = network.getDOF();
 	auto numSpecies = network.getSpeciesListSize();
-	auto myData = std::vector<double>(numSpecies * 4 + 1, 0.0);
+	auto networkSize = network.getNumClusters();
+	auto myData = std::vector<double>(numSpecies * 4 + 3, 0.0);
 
 	// Get the minimum size for the loop densities and diameters
 	auto minSizes = _solverHandler->getMinSizes();
@@ -776,6 +784,80 @@ PetscMonitor0D::computeFeCrAl(
 	// Get the pointer to the beginning of the solution data for this grid point
 	gridPointSolution = solutionArray[0];
 
+	// Get the delta time from the previous timestep to this timestep
+	double previousTime = _solverHandler->getPreviousTime();
+	double dt = time - previousTime;
+
+	// Compute the total number of impurities that were absorbed by sinks
+	if (timestep > 0) {
+		for (auto i = 0; i < 4; ++i) {
+			_nSink[i] += _previousSinkRate[i] * dt;
+			_previousSinkRate[i] = 0.0;
+		}
+	}
+
+	// Compute the sink strength for the previous rates
+	auto density = _solverHandler->getSinkDensity(); // nm-2
+	auto portion = _solverHandler->getSinkPortion();
+	auto r = 1.0 / sqrt(::xolotl::core::pi * density); // nm
+	auto rCore = ::xolotl::core::fecrCoreRadius;
+	auto temperature = gridPointSolution[dof];
+	constexpr double K = 170.0e9; // GPa
+	constexpr double nu = 0.29;
+	constexpr double b = 0.25; // nm
+	double deltaV = 1.67 * 0.5 * ::xolotl::core::ironLatticeConstant *
+		::xolotl::core::ironLatticeConstant *
+		::xolotl::core::ironLatticeConstant * 1.0e-27; // m3
+	//	constexpr double a0 = 0.91, a1 = -2.16, a2 = -0.92; // Random dipole
+	constexpr double a0 = 0.87, a1 = -5.12, a2 = -0.77; // Full network
+	constexpr double k_B = 1.380649e-23; // J K-1.
+	double L = (K * b * deltaV * (1.0 - 2.0 * nu)) /
+		(2.0 * ::xolotl::core::pi * k_B * temperature * (1.0 - nu));
+	double delta = sqrt(rCore * rCore + (L * L) / 4.0);
+	double edge = portion * density * 2.0 * ::xolotl::core::pi *
+		(a0 + a1 * (delta / r) + a2 * ((delta - rCore) / r)) /
+		std::log(r / delta);
+	double screw = (1.0 - portion) * density * 2.0 * ::xolotl::core::pi *
+		(a0 + a1 * (rCore / r)) / std::log(r / rCore);
+	// Loop on all the clusters
+	for (auto i = 0; i < networkSize; ++i) {
+		auto cluster = network.getCluster(i, plsm::HostMemSpace{});
+		auto diffCoef = cluster.getDiffusionCoefficient(0);
+		if (util::equal(diffCoef, 0.0))
+			continue;
+
+		auto clReg = cluster.getRegion();
+		Composition clLo = clReg.getOrigin();
+		// V case
+		if (clLo[Spec::V] > 0) {
+			double size =
+				clLo[Spec::V] + (double)(clReg[Spec::V].length() - 1) / 2.0;
+			// edge
+			_previousSinkRate[0] +=
+				gridPointSolution[i] * diffCoef * edge * size;
+			// screw
+			_previousSinkRate[1] +=
+				gridPointSolution[i] * diffCoef * screw * size;
+		}
+		else {
+			// I and Free case
+			double bias = 1.0;
+			double size = clLo[Spec::Free] +
+				(double)(clReg[Spec::Free].length() - 1) / 2.0;
+			if (clLo[Spec::I] > 0) {
+				bias = 1.05;
+				size =
+					clLo[Spec::I] + (double)(clReg[Spec::I].length() - 1) / 2.0;
+			}
+			// edge
+			_previousSinkRate[2] +=
+				gridPointSolution[i] * diffCoef * edge * bias * size;
+			// screw
+			_previousSinkRate[3] +=
+				gridPointSolution[i] * diffCoef * screw * bias * size;
+		}
+	}
+
 	using HostUnmanaged =
 		Kokkos::View<double*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>;
 	auto hConcs = HostUnmanaged(gridPointSolution, dof);
@@ -787,14 +869,17 @@ PetscMonitor0D::computeFeCrAl(
 		myData[4 * id()] = network.getTotalConcentration(dConcs, id, 1);
 		myData[(4 * id()) + 1] =
 			network.getTotalConcentration(dConcs, id, minSizes[id()]);
-		myData[(4 * id()) + 2] = 2.0 *
-			network.getTotalRadiusConcentration(dConcs, id, 1) /
-			myData[4 * id()];
+		myData[(4 * id()) + 2] =
+			2.0 * network.getTotalRadiusConcentration(dConcs, id, 1);
 		myData[(4 * id()) + 3] = 2.0 *
-			network.getTotalRadiusConcentration(dConcs, id, minSizes[id()]) /
-			myData[4 * id() + 1];
+			network.getTotalRadiusConcentration(dConcs, id, minSizes[id()]);
 		myData[myData.size() - 1] +=
 			network.getSmallConcentration(dConcs, id, minSizes[id()]);
+		if (network.getSpeciesLabel(id) == "Free" or
+			network.getSpeciesLabel(id) == "Trapped") {
+			myData[myData.size() - 3] += myData[(4 * id()) + 1];
+			myData[myData.size() - 2] += myData[(4 * id()) + 3];
+		}
 	}
 
 	// Set the output precision
@@ -809,9 +894,16 @@ PetscMonitor0D::computeFeCrAl(
 	outputFile << timestep << " " << time << " ";
 	for (auto i = 0; i < numSpecies; ++i) {
 		outputFile << myData[i * 4] << " " << myData[(i * 4) + 1] << " "
-				   << myData[(i * 4) + 2] << " " << myData[(i * 4) + 3] << " ";
+				   << myData[(i * 4) + 2] / myData[i * 4] << " "
+				   << myData[(i * 4) + 3] / myData[(i * 4) + 1] << " ";
+	}
+	outputFile << myData[myData.size() - 3] << " "
+			   << myData[myData.size() - 2] / myData[myData.size() - 3] << " "
+			   << myData[myData.size() - 1];
+	for (auto i = 0; i < 4; ++i) {
+		outputFile << " " << _nSink[i];
 	}
-	outputFile << myData[myData.size() - 1] << std::endl;
+	outputFile << std::endl;
 
 	// Close the output file
 	outputFile.close();