From 2e775c9c97bf87d0794f24488d0a950eab1e2867 Mon Sep 17 00:00:00 2001 From: Sophie Blondel Date: Thu, 2 Mar 2023 16:27:08 -0500 Subject: [PATCH] Fixed the FeCr flux and reaction rates to reproduce SPICES (#100). --- .../xolotl/core/flux/FeCrFitFluxHandler.h | 39 +-- .../xolotl/core/network/ConstantReaction.h | 8 + .../xolotl/core/network/FeCrReaction.h | 10 + .../xolotl/core/network/FeCrReactionNetwork.h | 15 +- .../include/xolotl/core/network/FeCrTraits.h | 55 ++++ .../xolotl/core/network/IReactionNetwork.h | 7 + .../xolotl/core/network/NucleationReaction.h | 8 + .../xolotl/core/network/ReSolutionReaction.h | 8 + .../include/xolotl/core/network/Reaction.h | 25 ++ .../xolotl/core/network/ReactionNetwork.h | 4 + .../xolotl/core/network/SinkReaction.h | 5 + .../xolotl/core/network/TransformReaction.h | 8 + .../core/network/TrapMutationReaction.h | 8 + .../include/xolotl/core/network/ZrReaction.h | 14 + .../xolotl/core/network/impl/FeCrReaction.tpp | 256 ++++++++++++------ .../core/network/impl/FeCrReactionNetwork.tpp | 109 +++++--- .../xolotl/core/network/impl/Reaction.tpp | 30 ++ .../core/network/impl/ReactionNetwork.tpp | 18 ++ .../xolotl/core/network/impl/SinkReaction.tpp | 16 ++ xolotl/viz/src/standard/plot/SeriesPlot.cpp | 3 +- 20 files changed, 507 insertions(+), 139 deletions(-) diff --git a/xolotl/core/include/xolotl/core/flux/FeCrFitFluxHandler.h b/xolotl/core/include/xolotl/core/flux/FeCrFitFluxHandler.h index ce4b65b7e..d1ece277f 100644 --- a/xolotl/core/include/xolotl/core/flux/FeCrFitFluxHandler.h +++ b/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; } diff --git a/xolotl/core/include/xolotl/core/network/ConstantReaction.h b/xolotl/core/include/xolotl/core/network/ConstantReaction.h index 08e0dd523..cc16455c4 100644 --- a/xolotl/core/include/xolotl/core/network/ConstantReaction.h +++ b/xolotl/core/include/xolotl/core/network/ConstantReaction.h @@ -311,6 +311,14 @@ class ConstantReaction : public Reaction return 0.0; } + KOKKOS_INLINE_FUNCTION + double + computeNetSigma(ConcentrationsView concentrations, IndexType clusterId, + IndexType gridIndex) + { + return 0.0; + } + KOKKOS_INLINE_FUNCTION void mapJacobianEntries(Connectivity connectivity) diff --git a/xolotl/core/include/xolotl/core/network/FeCrReaction.h b/xolotl/core/include/xolotl/core/network/FeCrReaction.h index 2a5a9384c..196c4fa31 100644 --- a/xolotl/core/include/xolotl/core/network/FeCrReaction.h +++ b/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 : diff --git a/xolotl/core/include/xolotl/core/network/FeCrReactionNetwork.h b/xolotl/core/include/xolotl/core/network/FeCrReactionNetwork.h index 221f87680..b4c0e48f2 100644 --- a/xolotl/core/include/xolotl/core/network/FeCrReactionNetwork.h +++ b/xolotl/core/include/xolotl/core/network/FeCrReactionNetwork.h @@ -1,8 +1,8 @@ #pragma once +#include #include #include -#include namespace xolotl { @@ -36,6 +36,19 @@ class FeCrReactionNetwork : public ReactionNetwork 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 values, IndexType gridIndex, double surfaceDepth, + double spacing); + private: double checkLatticeParameter(double latticeParameter); diff --git a/xolotl/core/include/xolotl/core/network/FeCrTraits.h b/xolotl/core/include/xolotl/core/network/FeCrTraits.h index 185523ef4..c2062c591 100644 --- a/xolotl/core/include/xolotl/core/network/FeCrTraits.h +++ b/xolotl/core/include/xolotl/core/network/FeCrTraits.h @@ -90,6 +90,61 @@ struct ReactionNetworkTraits using ClusterGenerator = FeCrClusterGenerator; }; + +namespace detail +{ +template +struct ClusterDataExtra +{ + using NetworkType = FeCrReactionNetwork; + + template + using View = ViewType; + + using IndexType = detail::ReactionNetworkIndexType; + + ClusterDataExtra() = default; + + template + KOKKOS_INLINE_FUNCTION + ClusterDataExtra(const ClusterDataExtra& data) : + netSigma(data.netSigma) + { + } + + template + void + deepCopy(const ClusterDataExtra& 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("Net Sigma", numClusters); + } + + View netSigma; +}; +} // namespace detail } // namespace network } // namespace core } // namespace xolotl diff --git a/xolotl/core/include/xolotl/core/network/IReactionNetwork.h b/xolotl/core/include/xolotl/core/network/IReactionNetwork.h index 3ad26b253..fbc9907fb 100644 --- a/xolotl/core/include/xolotl/core/network/IReactionNetwork.h +++ b/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. diff --git a/xolotl/core/include/xolotl/core/network/NucleationReaction.h b/xolotl/core/include/xolotl/core/network/NucleationReaction.h index ffd063be8..11fdbb953 100644 --- a/xolotl/core/include/xolotl/core/network/NucleationReaction.h +++ b/xolotl/core/include/xolotl/core/network/NucleationReaction.h @@ -113,6 +113,14 @@ class NucleationReaction : public Reaction return 0.0; } + KOKKOS_INLINE_FUNCTION + double + computeNetSigma(ConcentrationsView concentrations, IndexType clusterId, + IndexType gridIndex) + { + return 0.0; + } + KOKKOS_INLINE_FUNCTION void mapJacobianEntries(Connectivity connectivity); diff --git a/xolotl/core/include/xolotl/core/network/ReSolutionReaction.h b/xolotl/core/include/xolotl/core/network/ReSolutionReaction.h index e796c82c7..a1d642068 100644 --- a/xolotl/core/include/xolotl/core/network/ReSolutionReaction.h +++ b/xolotl/core/include/xolotl/core/network/ReSolutionReaction.h @@ -108,6 +108,14 @@ class ReSolutionReaction : public Reaction 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); diff --git a/xolotl/core/include/xolotl/core/network/Reaction.h b/xolotl/core/include/xolotl/core/network/Reaction.h index 66cbfcd0e..be643175e 100644 --- a/xolotl/core/include/xolotl/core/network/Reaction.h +++ b/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 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 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); diff --git a/xolotl/core/include/xolotl/core/network/ReactionNetwork.h b/xolotl/core/include/xolotl/core/network/ReactionNetwork.h index db5241570..729afa07c 100644 --- a/xolotl/core/include/xolotl/core/network/ReactionNetwork.h +++ b/xolotl/core/include/xolotl/core/network/ReactionNetwork.h @@ -472,6 +472,10 @@ class ReactionNetwork : public ReactionNetworkInterface::Type getLeftSideRate(ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) override; + double + getNetSigma(ConcentrationsView concentrations, IndexType clusterId, + IndexType gridIndex) override; + IndexType getDiagonalFill(SparseFillMap& fillMap) override; diff --git a/xolotl/core/include/xolotl/core/network/SinkReaction.h b/xolotl/core/include/xolotl/core/network/SinkReaction.h index b0ddf0ef6..95174c3e1 100644 --- a/xolotl/core/include/xolotl/core/network/SinkReaction.h +++ b/xolotl/core/include/xolotl/core/network/SinkReaction.h @@ -139,6 +139,11 @@ class SinkReaction : public Reaction return 0.0; } + KOKKOS_INLINE_FUNCTION + double + computeNetSigma(ConcentrationsView concentrations, IndexType clusterId, + IndexType gridIndex); + KOKKOS_INLINE_FUNCTION void mapJacobianEntries(Connectivity connectivity) diff --git a/xolotl/core/include/xolotl/core/network/TransformReaction.h b/xolotl/core/include/xolotl/core/network/TransformReaction.h index 8c4946221..ace0117d5 100644 --- a/xolotl/core/include/xolotl/core/network/TransformReaction.h +++ b/xolotl/core/include/xolotl/core/network/TransformReaction.h @@ -153,6 +153,14 @@ class TransformReaction : public Reaction return 0.0; } + KOKKOS_INLINE_FUNCTION + double + computeNetSigma(ConcentrationsView concentrations, IndexType clusterId, + IndexType gridIndex) + { + return 0.0; + } + KOKKOS_INLINE_FUNCTION void mapJacobianEntries(Connectivity connectivity) diff --git a/xolotl/core/include/xolotl/core/network/TrapMutationReaction.h b/xolotl/core/include/xolotl/core/network/TrapMutationReaction.h index 975795ffa..0d2bdfeea 100644 --- a/xolotl/core/include/xolotl/core/network/TrapMutationReaction.h +++ b/xolotl/core/include/xolotl/core/network/TrapMutationReaction.h @@ -139,6 +139,14 @@ class TrapMutationReaction : public Reaction 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) diff --git a/xolotl/core/include/xolotl/core/network/ZrReaction.h b/xolotl/core/include/xolotl/core/network/ZrReaction.h index e11d0cf1a..f28a6c6fe 100644 --- a/xolotl/core/include/xolotl/core/network/ZrReaction.h +++ b/xolotl/core/include/xolotl/core/network/ZrReaction.h @@ -51,6 +51,20 @@ class ZrSinkReaction : public SinkReaction 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); diff --git a/xolotl/core/include/xolotl/core/network/impl/FeCrReaction.tpp b/xolotl/core/include/xolotl/core/network/impl/FeCrReaction.tpp index a3bb15c2f..675d0242f 100644 --- a/xolotl/core/include/xolotl/core/network/impl/FeCrReaction.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/FeCrReaction.tpp @@ -14,18 +14,16 @@ template KOKKOS_INLINE_FUNCTION double getRate(const TRegion& pairCl0Reg, const TRegion& pairCl1Reg, const double r0, - const double r1, const double dc0, const double dc1, const double rho) + const double r1, const double dc0, const double dc1, const double rho, + const double sigma0, const double sigma1) { constexpr double pi = ::xolotl::core::pi; - const double zs = 4.0 * pi * (r0 + r1 + ::xolotl::core::fecrCoreRadius); + double zs = 4.0 * pi * (r0 + r1 + ::xolotl::core::fecrCoreRadius); using Species = typename TRegion::EnumIndex; auto lo0 = pairCl0Reg.getOrigin(); auto lo1 = pairCl1Reg.getOrigin(); - if (lo0.isOnAxis(Species::I) and lo1.isOnAxis(Species::I)) - return 1.1 * zs * (dc0 + dc1); - // react_extension if ((lo0.isOnAxis(Species::I) and lo1.isOnAxis(Species::V)) or (lo0.isOnAxis(Species::V) and lo1.isOnAxis(Species::I))) @@ -48,30 +46,88 @@ getRate(const TRegion& pairCl0Reg, const TRegion& pairCl1Reg, const double r0, if (cl0IsTrap || cl1IsTrap) { double rT = cl0IsTrap ? r0 : r1; double rL = cl0IsTrap ? r1 : r0; + double sigmaL = cl0IsTrap ? sigma1 : sigma0; double sigma = pi * (r0 + r1) * (r0 + r1); - if (rT > rL) + if (rT < rL) sigma = 4.0 * pi * r0 * r1; - return 2.0 * (dc0 + dc1) * sigma; + return 2.0 * (dc0 + dc1) * sigma * sigmaL; } - double rS = cl0IsSphere ? r0 : r1; - double rL = cl0IsSphere ? r1 : r0; - double rP = ::xolotl::core::fecrCoreRadius + rS; - if (!cl0IsSphere && !cl1IsSphere) { - rL = r0 + r1; - rP = ::xolotl::core::fecrCoalesceRadius; - } - double ratio = rL / (3.0 * rP); + double rB = (r0 > r1) ? r0 : r1; + double rs = (r0 > r1) ? r1 : r0; + double rint = ::xolotl::core::fecrCoreRadius + rs; + double ratio = rB / (3.0 * rint); double p = 1.0 / (1.0 + ratio * ratio); - double zl = 4.0 * pi * pi * rL / log(1.0 + 8.0 * rL / rP); + double zd = 4.0 * pi * pi * rB / log(1.0 + 8.0 * rB / rint); // Minimum loop size - rL = cl0IsSphere ? r1 : r0; - double zd = -8.0 * pi * pi * rL / - log(pi * rho * (rS + ::xolotl::core::fecrCoreRadius) * - (rS + ::xolotl::core::fecrCoreRadius)); + rB = (r0 > r1) ? r0 : r1; + double zd_alt = 0.0; + if (rho > 0.0) + zd_alt = -8.0 * pi * pi * rB / log(pi * rho * rint * rint); + + // Loop + loop + if (!cl0IsSphere and !cl1IsSphere) { + double sigmaL = dc0 > 0.0 ? sigma0 : sigma1; + double rSessile = dc0 > 0.0 ? r1 : r0; + double rFree = dc0 > 0.0 ? r0 : r1; + double dFree = dc0 > 0.0 ? dc0 : dc1; + + util::Array align; + align[0] = {0.0, 0.0}; + align[1] = {0.0, 0.0}; + align[2] = {0.0, 0.0}; + + // Find what type of cluster is the sessile one + auto loSessile = dc0 > 0.0 ? lo1 : lo0; + if (loSessile[static_cast(Species::Trapped)] > 0) { + align[0] = {1.0, 0.25}; + align[1] = {0.333333, 0.75}; + } + else if (loSessile[static_cast(Species::Junction)] > 0) { + align[0] = {1.0, 0.142857}; + align[1] = {0.333333, 0.4285714}; + align[2] = {0.57735, 0.4285714}; + } + else if (loSessile[static_cast(Species::Loop)] > 0) { + align[0] = {0.57735, 1.0}; + } + + double sigma = 0.0; + for (auto i = 0; i < align.size(); i++) { + double ao = rFree + rSessile + ::xolotl::core::fecrCoalesceRadius; + double bo = rFree + rSessile * align[i][0] + + ::xolotl::core::fecrCoalesceRadius; + double ai = util::max( + rFree + rSessile - ::xolotl::core::fecrCoalesceRadius, 0.0); + double bi = util::max(rFree + rSessile * align[i][0] - + ::xolotl::core::fecrCoalesceRadius, + 0.0); + + sigma += pi * (ao * bo - ai * bi) * align[i][1]; + } + double k_plus = 2.0 * dFree * sigma * sigmaL; + + return k_plus; + } + + double dsphere = cl0IsSphere ? dc0 : dc1; + + double k_plus = dsphere * (p * zs + (1.0 - p) * util::max(zd, zd_alt)); + + // Loop + Sphere point cross section + if (!cl0IsSphere or !cl1IsSphere) { + double sigmaL = cl0IsSphere ? sigma1 : sigma0; + double rsphere = cl0IsSphere ? r0 : r1; + double rloop = cl0IsSphere ? r1 : r0; + double dloop = cl0IsSphere ? dc1 : dc0; + rint = rsphere + ::xolotl::core::fecrCoreRadius; + double sigma = pi * (rloop + rint) * (rloop + rint); + if (rloop > rint) + sigma = 4.0 * pi * rloop * rint; + k_plus += 2.0 * dloop * sigma * sigmaL; + } - double k_plus = (dc0 + dc1) * (p * zs + (1.0 - p) * util::max(zd, zl)); double bias = 1.0; if (lo0.isOnAxis(Species::I) || lo1.isOnAxis(Species::I)) { if (lo0.isOnAxis(Species::Free) || lo1.isOnAxis(Species::Free)) @@ -101,7 +157,75 @@ FeCrProductionReaction::getRateForProduction(IndexType gridIndex) auto rho = this->_clusterData->sinkDensity(); // nm-2 - return getRate(cl0.getRegion(), cl1.getRegion(), r0, r1, dc0, dc1, rho); + return getRate(cl0.getRegion(), cl1.getRegion(), r0, r1, dc0, dc1, rho, + this->_clusterData->extraData.netSigma(_reactants[0]), + this->_clusterData->extraData.netSigma(_reactants[1])); +} + +KOKKOS_INLINE_FUNCTION +double +FeCrProductionReaction::computeNetSigma( + ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) +{ + double preFactor = 0.0; + // Check if our cluster is on the left side of this reaction + if (clusterId == _reactants[0]) { + preFactor = concentrations[_reactants[1]] * this->_coefs(0, 0, 0, 0); + } + if (clusterId == _reactants[1]) { + preFactor = concentrations[_reactants[0]] * this->_coefs(0, 0, 0, 0); + } + + if (clusterId == _reactants[0] or clusterId == _reactants[1]) { + // Compute cross section + auto cl0 = this->_clusterData->getCluster(_reactants[0]); + auto cl1 = this->_clusterData->getCluster(_reactants[1]); + + double r0 = cl0.getReactionRadius(); + double r1 = cl1.getReactionRadius(); + auto lo0 = cl0.getRegion().getOrigin(); + auto lo1 = cl1.getRegion().getOrigin(); + + bool cl0IsSphere = (lo0.isOnAxis(Species::V) || + lo0.isOnAxis(Species::Complex) || lo0.isOnAxis(Species::I)), + cl1IsSphere = (lo1.isOnAxis(Species::V) || + lo1.isOnAxis(Species::Complex) || lo1.isOnAxis(Species::I)); + bool cl0IsTrap = lo0.isOnAxis(Species::Trap), + cl1IsTrap = lo1.isOnAxis(Species::Trap); + // Sphere + sphere + if (cl0IsSphere && cl1IsSphere) { + return 0.0; + } + + // With a trap + if (cl0IsTrap || cl1IsTrap) { + double rT = cl0IsTrap ? r0 : r1; + double rL = cl0IsTrap ? r1 : r0; + double sigma = pi * (r0 + r1) * (r0 + r1); + if (rT < rL) + sigma = 4.0 * pi * r0 * r1; + return preFactor * sigma; + } + + // Loop + loop + if (!cl0IsSphere and !cl1IsSphere) { + return 0.0; + } + + // Loop + sphere + if (!cl0IsSphere or !cl1IsSphere) { + double rsphere = cl0IsSphere ? r0 : r1; + double rloop = cl0IsSphere ? r1 : r0; + double rint = rsphere + ::xolotl::core::fecrCoreRadius; + double sigma = pi * (rloop + rint) * (rloop + rint); + if (rloop > rint) + sigma = 4.0 * pi * rloop * rint; + return preFactor * sigma; + } + } + + // This cluster is not part of the reaction + return 0.0; } KOKKOS_INLINE_FUNCTION @@ -119,7 +243,11 @@ FeCrDissociationReaction::getRateForProduction(IndexType gridIndex) auto rho = this->_clusterData->sinkDensity(); // nm-2 - return getRate(cl0.getRegion(), cl1.getRegion(), r0, r1, dc0, dc1, rho); + double kPlus = getRate(cl0.getRegion(), cl1.getRegion(), r0, r1, dc0, dc1, + rho, this->_clusterData->extraData.netSigma(_products[0]), + this->_clusterData->extraData.netSigma(_products[1])); + + return kPlus; } KOKKOS_INLINE_FUNCTION @@ -139,6 +267,7 @@ FeCrDissociationReaction::computeRate(IndexType gridIndex, double) auto cl0 = this->_clusterData->getCluster(_products[0]); auto cl1 = this->_clusterData->getCluster(_products[1]); + auto lo = this->_clusterData->getCluster(_reactant).getRegion().getOrigin(); auto lo0 = cl0.getRegion().getOrigin(); auto lo1 = cl1.getRegion().getOrigin(); bool cl0IsSphere = (lo0.isOnAxis(Species::V) || @@ -150,74 +279,35 @@ FeCrDissociationReaction::computeRate(IndexType gridIndex, double) double kMinus = kPlus * std::exp(-E_b / (k_B * T)); - // Recoil detrap - if ((lo0.isOnAxis(Species::Trap) && lo1.isOnAxis(Species::Free)) || - (lo0.isOnAxis(Species::Free) && lo1.isOnAxis(Species::Trap))) { - auto cl = this->_clusterData->getCluster(_reactant); - Composition lo = cl.getRegion().getOrigin(); - double kRecoil = exp(-(double)lo[Species::Trapped] / 10000.0); - if (lo[Species::Junction] > 0) - kRecoil = exp(-(double)lo[Species::Junction] / 400.0); - // kRecoil = 0.0; - - return ((kMinus / this->_clusterData->atomicVolume()) + - ::xolotl::core::detrapFrequency) * - kRecoil; - } - // Standard case - if (cl0IsSphere and cl1IsSphere) + if (cl0IsSphere and cl1IsSphere) { return (kMinus / this->_clusterData->atomicVolume()); + } double r0 = cl0.getReactionRadius(); double r1 = cl1.getReactionRadius(); - auto rCoal = ::xolotl::core::fecrCoalesceRadius; - double a = r0 + r1 + rCoal; - double b = util::max(r0 + r1 - rCoal, 0.0); - double sigma = ::xolotl::core::pi * (a * a - b * b); + double rsphere = cl0IsSphere ? r0 : r1; + double rloop = cl0IsSphere ? r1 : r0; - const double jumpDistance = - this->_clusterData->latticeParameter() * sqrt(3.0) / 2.0; + double rint = rsphere + ::xolotl::core::fecrCoreRadius; + double sigma = pi * (rloop + rint) * (rloop + rint); + if (rloop > rint) + sigma = 4.0 * pi * rloop * rint; - return kMinus / (sigma * jumpDistance); - - // Trap case if (cl0IsTrap or cl1IsTrap) { double rT = cl0IsTrap ? r0 : r1; double rL = cl0IsTrap ? r1 : r0; sigma = pi * (r0 + r1) * (r0 + r1); - if (rT > rL) + if (rT < rL) sigma = 4.0 * pi * r0 * r1; } - // Loop and sphere case - else if (cl0IsSphere or cl1IsSphere) { - double rA = cl0IsSphere ? r0 : r1; - double rP = cl0IsSphere ? r1 : r0; - double rO = rP + ::xolotl::core::fecrCoreRadius; - for (auto q : _align4by4) { - double aOut = rA + rO; - double bOut = rA + rO * q; - double aIn = util::max(rA - rO, 0.0); - double bIn = util::max(rA - rO * q, 0.0); - - sigma += _q4by4 * pi * (aOut * bOut - aIn * bIn); - } - } - // Loop and loop - else { - double rO = ::xolotl::core::fecrCoalesceRadius; - for (auto q : _align4by4) { - double aOut = r0 + r1 + rO; - double bOut = r0 + r1 * q + rO; - double aIn = util::max(r0 + r1 - rO, 0.0); - double bIn = util::max(r0 + r1 * q - rO, 0.0); - - sigma += _q4by4 * pi * (aOut * bOut - aIn * bIn); - } - } - return kMinus / (sigma * jumpDistance); + double preFactor = 1.0 / + (sigma * this->_clusterData->latticeParameter() * + ::xolotl::core::fecrBurgers); + + return kMinus * preFactor; } KOKKOS_INLINE_FUNCTION @@ -334,6 +424,20 @@ FeCrSinkReaction::getSinkStrength() return density * Z; } +KOKKOS_INLINE_FUNCTION +double +FeCrSinkReaction::computeNetSigma( + ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) +{ + // Check if our cluster is on the left side of this reaction + if (clusterId == _reactant) { + return this->asDerived()->getSinkStrength(); + } + + // This cluster is not part of the reaction + return 0.0; +} + KOKKOS_INLINE_FUNCTION double FeCrTransformReaction::getSize() diff --git a/xolotl/core/include/xolotl/core/network/impl/FeCrReactionNetwork.tpp b/xolotl/core/include/xolotl/core/network/impl/FeCrReactionNetwork.tpp index d0d5bb5a1..2359a4eda 100644 --- a/xolotl/core/include/xolotl/core/network/impl/FeCrReactionNetwork.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/FeCrReactionNetwork.tpp @@ -12,6 +12,53 @@ namespace core { namespace network { +void +FeCrReactionNetwork::initializeExtraClusterData( + const options::IOptions& options) +{ + this->_clusterData.h_view().extraData.initialize( + this->_clusterData.h_view().numClusters); + this->copyClusterDataView(); + this->invalidateDataMirror(); +} + +void +FeCrReactionNetwork::computeFluxesPreProcess(ConcentrationsView concentrations, + FluxesView fluxes, IndexType gridIndex, double surfaceDepth, double spacing) +{ + auto data = this->_clusterData.h_view(); + auto diffusionFactor = this->getClusterDataMirror().diffusionFactor; + auto sigma = create_mirror_view(data.extraData.netSigma); + for (auto i = 0; i < this->_numClusters; i++) { + if (diffusionFactor(i) == 0.0) { + sigma(i) = 0.0; + continue; + } + sigma(i) = this->getNetSigma(concentrations, i, gridIndex); + } + deep_copy(data.extraData.netSigma, sigma); + this->updateReactionRates(0.0); +} + +void +FeCrReactionNetwork::computePartialsPreProcess( + ConcentrationsView concentrations, Kokkos::View values, + IndexType gridIndex, double surfaceDepth, double spacing) +{ + auto data = this->_clusterData.h_view(); + auto diffusionFactor = this->getClusterDataMirror().diffusionFactor; + auto sigma = create_mirror_view(data.extraData.netSigma); + for (auto i = 0; i < this->_numClusters; i++) { + if (diffusionFactor(i) == 0.0) { + sigma(i) = 0.0; + continue; + } + sigma(i) = this->getNetSigma(concentrations, i, gridIndex); + } + deep_copy(data.extraData.netSigma, sigma); + this->updateReactionRates(0.0); +} + namespace detail { template @@ -494,25 +541,15 @@ FeCrReactionGenerator::operator()(IndexType i, IndexType j, TTag tag) const int maxSize = fReg[Species::Free].end() + tReg[Species::Trapped].end() - 2; - for (auto size = minSize; size <= maxSize; size++) { - // Compute possible ratio - auto maxL = util::min(size - tReg[Species::Trapped].begin(), - fReg[Species::Free].end() - 1); - auto minL = util::max(size - tReg[Species::Trapped].end() + 1, - fReg[Species::Free].begin()); - - double minRatio = - fabs((double)size - (double)(2 * maxL)) / (double)size; - double maxRatio = - fabs((double)size - (double)(2 * minL)) / (double)size; - if (minRatio > maxRatio) { - double temp = minRatio; - minRatio = maxRatio; - maxRatio = temp; - } + double ratio = 2.0 * + fabs((double)fReg[Species::Free].begin() - + (double)tReg[Species::Trapped].begin()) / + ((double)fReg[Species::Free].begin() + + (double)tReg[Species::Trapped].begin()); + for (auto size = minSize; size <= maxSize; size++) { // Junction case - if (minRatio < 0.5) { + if (ratio < 0.5) { Composition comp = Composition::zero(); comp[Species::Junction] = size; comp[Species::Trap] = 1; @@ -525,7 +562,7 @@ FeCrReactionGenerator::operator()(IndexType i, IndexType j, TTag tag) const } } // Trapped case - if (maxRatio >= 0.5) { + if (ratio >= 0.5) { Composition comp = Composition::zero(); comp[Species::Trapped] = size; comp[Species::Trap] = 1; @@ -588,29 +625,15 @@ FeCrReactionGenerator::operator()(IndexType i, IndexType j, TTag tag) const int minSize = fReg[Species::Free].begin() + lReg[Species::Loop].begin(); int maxSize = fReg[Species::Free].end() + lReg[Species::Loop].end() - 2; - for (auto size = minSize; size <= maxSize; size++) { - // Compute possible ratio - auto maxL = util::min(size - lReg[Species::Loop].begin(), - fReg[Species::Free].end() - 1); - auto minL = util::max(size - lReg[Species::Loop].end() + 1, - fReg[Species::Free].begin()); - auto maxK = util::min(size - fReg[Species::Free].begin(), - lReg[Species::Loop].end() - 1); - auto minK = util::max(size - fReg[Species::Free].end() + 1, - lReg[Species::Loop].begin()); - - double minRatio = - fabs((double)size - (double)(2 * maxL)) / (double)size; - double maxRatio = - fabs((double)size - (double)(2 * minL)) / (double)size; - if (minRatio > maxRatio) { - double temp = minRatio; - minRatio = maxRatio; - maxRatio = temp; - } + double ratio = 2.0 * + fabs((double)fReg[Species::Free].begin() - + (double)lReg[Species::Loop].begin()) / + ((double)fReg[Species::Free].begin() + + (double)lReg[Species::Loop].begin()); + for (auto size = minSize; size <= maxSize; size++) { // Junction case - if (minRatio < 0.5) { + if (ratio < 0.5) { Composition comp = Composition::zero(); comp[Species::Junction] = size; comp[Species::Trap] = 1; @@ -622,10 +645,10 @@ FeCrReactionGenerator::operator()(IndexType i, IndexType j, TTag tag) const // No dissociation } } - // Trapped of Loop case - if (maxRatio > 0.5) { + // Trapped or Loop case + if (ratio > 0.5) { // Loop - if (minL < maxK) { + if (fReg[Species::Free].begin() < lReg[Species::Loop].begin()) { Composition comp = Composition::zero(); comp[Species::Loop] = size; comp[Species::Trap] = 1; @@ -638,7 +661,7 @@ FeCrReactionGenerator::operator()(IndexType i, IndexType j, TTag tag) const } } // Trapped - if (minK < maxL) { + else { Composition comp = Composition::zero(); comp[Species::Trapped] = size; comp[Species::Trap] = 1; diff --git a/xolotl/core/include/xolotl/core/network/impl/Reaction.tpp b/xolotl/core/include/xolotl/core/network/impl/Reaction.tpp index 60440fd38..cb1eea28f 100644 --- a/xolotl/core/include/xolotl/core/network/impl/Reaction.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/Reaction.tpp @@ -1835,6 +1835,36 @@ ProductionReaction::computeLeftSideRate( return 0.0; } +template +KOKKOS_INLINE_FUNCTION +double +ProductionReaction::computeNetSigma( + ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) +{ + double preFactor = 0.0; + // Check if our cluster is on the left side of this reaction + if (clusterId == _reactants[0]) { + preFactor = concentrations[_reactants[1]] * this->_coefs(0, 0, 0, 0); + } + if (clusterId == _reactants[1]) { + preFactor = concentrations[_reactants[0]] * this->_coefs(0, 0, 0, 0); + } + + if (clusterId == _reactants[0] or clusterId == _reactants[1]) { + // Compute cross section + auto cl0 = this->_clusterData->getCluster(_reactants[0]); + auto cl1 = this->_clusterData->getCluster(_reactants[1]); + + double r0 = cl0.getReactionRadius(); + double r1 = cl1.getReactionRadius(); + + return preFactor * (r0 + r1) * (r0 + r1) * ::xolotl::core::pi; + } + + // This cluster is not part of the reaction + return 0.0; +} + template KOKKOS_INLINE_FUNCTION void diff --git a/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp b/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp index f7dc60d85..d8331ff12 100644 --- a/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/ReactionNetwork.tpp @@ -658,6 +658,24 @@ ReactionNetwork::getLeftSideRate( return leftSideRate; } +template +double +ReactionNetwork::getNetSigma( + ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) +{ + double sigma = 0.0; + // Loop on all the reactions + _reactions.reduce( + "ReactionNetwork::getNetSigma", + DEVICE_LAMBDA(auto&& reaction, double& lsum) { + lsum += reaction.contributeNetSigma( + concentrations, clusterId, gridIndex); + }, + sigma); + + return sigma; +} + template double ReactionNetwork::getTotalConcentration( diff --git a/xolotl/core/include/xolotl/core/network/impl/SinkReaction.tpp b/xolotl/core/include/xolotl/core/network/impl/SinkReaction.tpp index 986558e0f..deef96e96 100644 --- a/xolotl/core/include/xolotl/core/network/impl/SinkReaction.tpp +++ b/xolotl/core/include/xolotl/core/network/impl/SinkReaction.tpp @@ -24,6 +24,22 @@ SinkReaction::computeRate(IndexType gridIndex, double time) return strength; } + +template +KOKKOS_INLINE_FUNCTION +double +SinkReaction::computeNetSigma( + ConcentrationsView concentrations, IndexType clusterId, IndexType gridIndex) +{ + // Check if our cluster is on the left side of this reaction + if (clusterId == _reactant) { + return this->asDerived()->getSinkBias() * + this->asDerived()->getSinkStrength(); + } + + // This cluster is not part of the reaction + return 0.0; +} } // namespace network } // namespace core } // namespace xolotl diff --git a/xolotl/viz/src/standard/plot/SeriesPlot.cpp b/xolotl/viz/src/standard/plot/SeriesPlot.cpp index f4b43728d..fe56c2c08 100644 --- a/xolotl/viz/src/standard/plot/SeriesPlot.cpp +++ b/xolotl/viz/src/standard/plot/SeriesPlot.cpp @@ -97,7 +97,8 @@ SeriesPlot::render(const std::string& fileName) // Accumulate the bounds of our data to focus camera fieldBounds.X = dataSet.GetCoordinateSystem().GetBounds().X; - dataSet.GetField(plotDataProviders->at(i)->getDataName()).GetRange(&fieldBounds.Y); + dataSet.GetField(plotDataProviders->at(i)->getDataName()) + .GetRange(&fieldBounds.Y); bounds.Include(fieldBounds); // Add Plot to our scene for later rendering