From f80851168d8c8c919b856fd25eb74b0e0fa88337 Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Tue, 7 Jan 2020 16:18:48 -0700 Subject: [PATCH 1/6] Refactor mapIndex and linearInterpolation (#410) --- .../utils/PolyatomicDamageEnergyFunction.h | 11 +++-- ...PolyatomicDisplacementDerivativeFunction.h | 8 +--- .../utils/PolyatomicDisplacementFunction.h | 10 ++--- .../PolyatomicDisplacementFunctionBase.h | 19 +++++++-- src/other/ThreadedRecoilLoopBase.C | 4 +- src/userobjects/PolyatomicRecoil.C | 9 ++-- src/utils/PolyatomicDamageEnergyFunction.C | 34 ++------------- ...PolyatomicDisplacementDerivativeFunction.C | 40 +++--------------- src/utils/PolyatomicDisplacementFunction.C | 41 +++---------------- .../PolyatomicDisplacementFunctionBase.C | 28 +++++++++++++ 10 files changed, 78 insertions(+), 126 deletions(-) diff --git a/include/utils/PolyatomicDamageEnergyFunction.h b/include/utils/PolyatomicDamageEnergyFunction.h index 653b5829..9586305e 100644 --- a/include/utils/PolyatomicDamageEnergyFunction.h +++ b/include/utils/PolyatomicDamageEnergyFunction.h @@ -25,11 +25,8 @@ class PolyatomicDamageEnergyFunction : public PolyatomicDisplacementFunctionBase static int odeRHS(Real energy, const Real disp[], Real f[], void * params); - ///@{ some getters needed for accessing this pointer in odeRHS - Real linearInterpolation(unsigned int i, Real energy) const; - Real linearInterpolation(unsigned int i, Real energy, unsigned int index) const; + /// a getter needed for accessing this pointer in odeRHS Real taylorSeriesThreshold() const { return _taylor_series_threshold; } - ///@} Real integralTypeI(Real energy, unsigned int i, unsigned int j); @@ -40,6 +37,12 @@ class PolyatomicDamageEnergyFunction : public PolyatomicDisplacementFunctionBase * by nu_i(E-T) - nu_i(E) \approx -T d(nu_i) / dE. */ Real _taylor_series_threshold = 1e2; + +protected: + unsigned int mapIndex(unsigned int i, unsigned int /*j*/, unsigned int /*l*/) const override + { + return i; + }; }; #endif // GSL_ENABLED diff --git a/include/utils/PolyatomicDisplacementDerivativeFunction.h b/include/utils/PolyatomicDisplacementDerivativeFunction.h index 4c4a623c..87351a2b 100644 --- a/include/utils/PolyatomicDisplacementDerivativeFunction.h +++ b/include/utils/PolyatomicDisplacementDerivativeFunction.h @@ -27,12 +27,6 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu static int odeRHS(Real energy, const Real disp[], Real f[], void * params); - ///@{ some getters needed for accessing this pointer in odeRHS - Real linearInterpolation(unsigned int i, unsigned int j, unsigned int l, Real energy) const; - Real linearInterpolation( - unsigned int i, unsigned int j, unsigned int l, Real energy, unsigned int index) const; - ///@} - /// computes term 1 in Parkin-Coulter expression nu_k(T - Eb) Real integralTypeI(Real energy, unsigned int i, unsigned int j, unsigned int l, unsigned int k) const; @@ -46,7 +40,7 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu protected: /// maps triple index theta_ijl to single theta_n; l runs faster than j runs faster than i - unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const + unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const override { return i + j * _n_species + l * _n_species * _n_species; }; diff --git a/include/utils/PolyatomicDisplacementFunction.h b/include/utils/PolyatomicDisplacementFunction.h index 303cdc27..770c627c 100644 --- a/include/utils/PolyatomicDisplacementFunction.h +++ b/include/utils/PolyatomicDisplacementFunction.h @@ -25,11 +25,6 @@ class PolyatomicDisplacementFunction : public PolyatomicDisplacementFunctionBase static int odeRHS(Real energy, const Real disp[], Real f[], void * params); - ///@{ some getters needed for accessing this pointer in odeRHS - Real linearInterpolation(unsigned int i, unsigned int j, Real energy) const; - Real linearInterpolation(unsigned int i, unsigned int j, Real energy, unsigned int index) const; - ///@} - /// computes term 1 in Parkin-Coulter expression nu_k(T - Eb) Real integralTypeI(Real energy, unsigned int i, unsigned int j, unsigned int k) const; @@ -38,7 +33,10 @@ class PolyatomicDisplacementFunction : public PolyatomicDisplacementFunctionBase protected: /// maps double index nu_ij to single index nu_n; j runs faster than i - unsigned int mapIndex(unsigned int i, unsigned int j) const { return i + j * _n_species; }; + unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int /*l*/) const override + { + return i + j * _n_species; + }; /// is the total damage function computed bool _total_displacement_function; diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h index a248fe13..e247e317 100644 --- a/include/utils/PolyatomicDisplacementFunctionBase.h +++ b/include/utils/PolyatomicDisplacementFunctionBase.h @@ -41,6 +41,7 @@ class PolyatomicDisplacementFunctionBase virtual ~PolyatomicDisplacementFunctionBase(); + /// computes the displacement function from current last energy to Emax void advanceDisplacements(Real Emax); ///@{ some getters needed for accessing this pointer in odeRHS @@ -55,6 +56,9 @@ class PolyatomicDisplacementFunctionBase Real numberDensity(unsigned int i) const { return _material->_element[i]._t * _material->_arho; } ///@} + /// linear interpolation of the damage function + Real linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + /// gets stopping power for a given species and energy; non-const because it uses _ions so no need to construct ion Real stoppingPower(unsigned int species, Real energy); @@ -65,6 +69,9 @@ class PolyatomicDisplacementFunctionBase unsigned int energyIndex(Real energy) const; protected: + /// this function flattens the arrays i, j, l to a single index + virtual unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const = 0; + /// computes the integral int_0^t dT T * d(sigma_ij) / dT for species combination i, j and small t Real weightedScatteringIntegral(Real energy, Real energy_limit, unsigned int i, unsigned int j) const; @@ -83,6 +90,11 @@ class PolyatomicDisplacementFunctionBase virtual Real nonCaptureProbability(unsigned int i, unsigned int k, Real energy, Real recoil_energy) const; + /// a helper function called from linearInterpolation + Real linearInterpolationHelper( + Real energy, unsigned int index, unsigned int i, unsigned int j, unsigned int l) const; + + /// damage function type [nij and gij, respectively in PK JNM 101, 1981; or nu_i JNM 88, (1980)] nrt_type _damage_function_type; @@ -98,11 +110,12 @@ class PolyatomicDisplacementFunctionBase /// the current maximum energy to which PC equations are solved std::vector _energy_history; - /** - * The displacement values nu_ij flattened into 1D array - */ + /// The displacement values nu_ij flattened into 1D array std::vector> _displacement_function; + /// The integral of the displacement function over energy + std::vector> _displacement_function_integral; + /// Ecap_ik average residual energy of type i atom to not be trapped at k-site std::vector> _Ecap; diff --git a/src/other/ThreadedRecoilLoopBase.C b/src/other/ThreadedRecoilLoopBase.C index 505476cd..23b39b89 100644 --- a/src/other/ThreadedRecoilLoopBase.C +++ b/src/other/ThreadedRecoilLoopBase.C @@ -167,10 +167,10 @@ ThreadedRecoilLoopBase::operator()(const PKARange & pka_list) _trim_parameters.element_prototypes[target_var]._m); // using linear perturbation theory estimate of g_ij(number_fractions) Real w = _pa_nrt[index]->linearInterpolation( - species_index, target_species_index, recoil->_E); + recoil->_E, species_index, target_species_index); for (unsigned int l = 0; l < _pa_derivative_nrt[index]->nSpecies(); ++l) w += _pa_derivative_nrt[index]->linearInterpolation( - species_index, target_species_index, l, recoil->_E) * + recoil->_E, species_index, target_species_index, l) * (number_fractions[l] - _pa_nrt[index]->numberFraction(l)); // increment energy, vacancy and interstitial buffers diff --git a/src/userobjects/PolyatomicRecoil.C b/src/userobjects/PolyatomicRecoil.C index bdf1ce56..fb82aa58 100644 --- a/src/userobjects/PolyatomicRecoil.C +++ b/src/userobjects/PolyatomicRecoil.C @@ -210,10 +210,10 @@ PolyatomicRecoil::finalize() ++derivative) displacement_file << "," << _padf_derivative->linearInterpolation( + _padf_derivative->energyPoint(n), projectile, target, - derivative, - _padf_derivative->energyPoint(n)); + derivative); displacement_file << std::endl; } else @@ -230,7 +230,7 @@ PolyatomicRecoil::finalize() getParam("damage_type") == "TOTAL" && target == projectile ? 1 : 0; displacement_file << "," << delta_ij + displacement_function->linearInterpolation( - projectile, target, _padf->energyPoint(n)); + _padf->energyPoint(n), projectile, target); } displacement_file << std::endl; } @@ -242,8 +242,7 @@ PolyatomicRecoil::finalize() displacement_file << _padf->energyPoint(n); for (unsigned int projectile = 0; projectile < _padf->nSpecies(); ++projectile) displacement_file << "," - << energy_function->linearInterpolation(projectile, - _padf->energyPoint(n)); + << energy_function->linearInterpolation(_padf->energyPoint(n), projectile); displacement_file << std::endl; } } diff --git a/src/utils/PolyatomicDamageEnergyFunction.C b/src/utils/PolyatomicDamageEnergyFunction.C index 9c552c99..10421426 100644 --- a/src/utils/PolyatomicDamageEnergyFunction.C +++ b/src/utils/PolyatomicDamageEnergyFunction.C @@ -110,7 +110,7 @@ PolyatomicDamageEnergyFunction::integralTypeI(Real energy, unsigned int i, unsig { Real recoil_energy = scale * (_quad_points[qp] + 1) + lower; integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) * - linearInterpolation(j, recoil_energy, l); + linearInterpolationHelper(recoil_energy, l, j, 0, 0); } } return integral; @@ -123,11 +123,11 @@ PolyatomicDamageEnergyFunction::integralTypeII(Real energy, unsigned int i, unsi Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][j]); // store the current displacement function value - Real current_value = linearInterpolation(i, energy); + Real current_value = linearInterpolation(energy, i); // estimate the derivative d(nu_i) / dE: Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2]; - Real derivative = (current_value - linearInterpolation(i, energy - dE)) / dE; + Real derivative = (current_value - linearInterpolation(energy - dE, i)) / dE; // integrate up to threshold and multiply by estimate of the derivative Real integral = -weightedScatteringIntegral(energy, threshold, i, j) * derivative; @@ -156,36 +156,10 @@ PolyatomicDamageEnergyFunction::integralTypeII(Real energy, unsigned int i, unsi { Real recoil_energy = scale * (_quad_points[qp] + 1) + lower; integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) * - (linearInterpolation(i, energy - recoil_energy) - current_value); + (linearInterpolation( energy - recoil_energy, i) - current_value); } } return integral; } -Real -PolyatomicDamageEnergyFunction::linearInterpolation(unsigned int i, Real energy) const -{ - unsigned int index = energyIndex(energy); - if (index == 0) - return _displacement_function[0][i]; - - return linearInterpolation(i, energy, index); -} - -Real -PolyatomicDamageEnergyFunction::linearInterpolation(unsigned int i, - Real energy, - unsigned int index) const -{ - // interpolation: power law interpolation df = a * E^b unless one value is zero => linear - // interpolation - Real e1 = _energy_history[index - 1]; - Real e2 = _energy_history[index]; - Real v1 = _displacement_function[index - 1][i]; - Real v2 = _displacement_function[index][i]; - - // do linear interpolation instead of power law - return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); -} - #endif diff --git a/src/utils/PolyatomicDisplacementDerivativeFunction.C b/src/utils/PolyatomicDisplacementDerivativeFunction.C index ac50f93f..fb74011b 100644 --- a/src/utils/PolyatomicDisplacementDerivativeFunction.C +++ b/src/utils/PolyatomicDisplacementDerivativeFunction.C @@ -104,7 +104,7 @@ PolyatomicDisplacementDerivativeFunction::integralTypeI( Real recoil_energy = f * (_quad_points[qp] + 1) + lower; integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) * displacementProbability(k, recoil_energy) * - linearInterpolation(k, j, l, recoil_energy - Ebind); + linearInterpolation(recoil_energy - Ebind, k, j, l); } } return integral; @@ -118,11 +118,11 @@ PolyatomicDisplacementDerivativeFunction::integralTypeII( Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][k]); // store the current displacement function value - Real current_value = linearInterpolation(i, j, l, energy); + Real current_value = linearInterpolation(energy, i, j, l); // estimate the derivative d(nu_i) / dE: Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2]; - Real derivative = (current_value - linearInterpolation(i, j, l, energy - dE)) / dE; + Real derivative = (current_value - linearInterpolation(energy - dE, i, j, l)) / dE; // integrate up to threshold and multiply by estimate of the derivative Real integral = -weightedScatteringIntegral(energy, threshold, i, k) * derivative; @@ -153,7 +153,7 @@ PolyatomicDisplacementDerivativeFunction::integralTypeII( Real recoil_energy = f * (_quad_points[qp] + 1) + lower; integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) * (nonCaptureProbability(i, k, energy, recoil_energy) * - linearInterpolation(i, j, l, energy - recoil_energy) - + linearInterpolation(energy - recoil_energy, i, j, l) - current_value); } } @@ -174,8 +174,8 @@ PolyatomicDisplacementDerivativeFunction::source(Real energy, { Real e1 = _net_displacement_function->energyPoint(index - 1); Real e2 = _net_displacement_function->energyPoint(index); - Real v1 = _net_displacement_function->linearInterpolation(i, j, e1); - Real v2 = _net_displacement_function->linearInterpolation(i, j, e2); + Real v1 = _net_displacement_function->linearInterpolation(e1, i, j); + Real v2 = _net_displacement_function->linearInterpolation(e2, i, j); d_net_dE = (v2 - v1) / (e2 - e1); } return _net_displacement_function->integralTypeI(energy, i, j, l) + @@ -183,32 +183,4 @@ PolyatomicDisplacementDerivativeFunction::source(Real energy, d_net_dE * stoppingPowerDerivative(i, l, energy); } -Real -PolyatomicDisplacementDerivativeFunction::linearInterpolation(unsigned int i, - unsigned int j, - unsigned int l, - Real energy) const -{ - unsigned int index = energyIndex(energy); - if (index == 0) - return _displacement_function[0][mapIndex(i, j, l)]; - - return linearInterpolation(i, j, l, energy, index); -} - -Real -PolyatomicDisplacementDerivativeFunction::linearInterpolation( - unsigned int i, unsigned int j, unsigned int l, Real energy, unsigned int index) const -{ - unsigned int k = mapIndex(i, j, l); - - // linear interpolation - Real e1 = _energy_history[index - 1]; - Real e2 = _energy_history[index]; - Real v1 = _displacement_function[index - 1][k]; - Real v2 = _displacement_function[index][k]; - - return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); -} - #endif diff --git a/src/utils/PolyatomicDisplacementFunction.C b/src/utils/PolyatomicDisplacementFunction.C index 69ff32a7..0bafb05e 100644 --- a/src/utils/PolyatomicDisplacementFunction.C +++ b/src/utils/PolyatomicDisplacementFunction.C @@ -48,7 +48,7 @@ PolyatomicDisplacementFunction::PolyatomicDisplacementFunction( // note: for net displacement function, the gii = 1, for total displacement function nij = 0 if (_damage_function_type == NET) for (unsigned int i = 0; i < _n_species; ++i) - _displacement_function[0][mapIndex(i, i)] = 1; + _displacement_function[0][mapIndex(i, i, 0)] = 1; } int @@ -62,7 +62,7 @@ PolyatomicDisplacementFunction::odeRHS(Real energy, const Real disp[], Real f[], for (unsigned int j = 0; j < padf->nSpecies(); ++j) { // working on the right hand side for nu_ij - unsigned int n = padf->mapIndex(i, j); + unsigned int n = padf->mapIndex(i, j, 0); f[n] = 0; for (unsigned int k = 0; k < padf->nSpecies(); ++k) @@ -103,7 +103,7 @@ PolyatomicDisplacementFunction::integralTypeI(Real energy, Real recoil_energy = f * (_quad_points[qp] + 1) + lower; integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) * displacementProbability(k, recoil_energy) * - (delta_kj + linearInterpolation(k, j, recoil_energy - Ebind)); + (delta_kj + linearInterpolation(recoil_energy - Ebind, k, j)); } } return integral; @@ -119,11 +119,11 @@ PolyatomicDisplacementFunction::integralTypeII(Real energy, Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][k]); // store the current displacement function value - Real current_value = linearInterpolation(i, j, energy); + Real current_value = linearInterpolation(energy, i, j); // estimate the derivative d(nu_i) / dE: Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2]; - Real derivative = (current_value - linearInterpolation(i, j, energy - dE)) / dE; + Real derivative = (current_value - linearInterpolation(energy - dE, i, j)) / dE; // integrate up to threshold and multiply by estimate of the derivative Real integral = -weightedScatteringIntegral(energy, threshold, i, k) * derivative; @@ -154,40 +154,11 @@ PolyatomicDisplacementFunction::integralTypeII(Real energy, Real recoil_energy = f * (_quad_points[qp] + 1) + lower; integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) * (nonCaptureProbability(i, k, energy, recoil_energy) * - linearInterpolation(i, j, energy - recoil_energy) - + linearInterpolation(energy - recoil_energy, i, j) - current_value); } } return integral; } -Real -PolyatomicDisplacementFunction::linearInterpolation(unsigned int i, - unsigned int j, - Real energy) const -{ - unsigned int index = energyIndex(energy); - if (index == 0) - return _displacement_function[0][mapIndex(i, j)]; - - return linearInterpolation(i, j, energy, index); -} - -Real -PolyatomicDisplacementFunction::linearInterpolation(unsigned int i, - unsigned int j, - Real energy, - unsigned int index) const -{ - unsigned int k = mapIndex(i, j); - - // linear interpolation - Real e1 = _energy_history[index - 1]; - Real e2 = _energy_history[index]; - Real v1 = _displacement_function[index - 1][k]; - Real v2 = _displacement_function[index][k]; - - return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); -} - #endif diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C index 007f0142..eb37898b 100644 --- a/src/utils/PolyatomicDisplacementFunctionBase.C +++ b/src/utils/PolyatomicDisplacementFunctionBase.C @@ -267,4 +267,32 @@ PolyatomicDisplacementFunctionBase::weightedScatteringIntegral(Real energy, return integral; } +Real +PolyatomicDisplacementFunctionBase::linearInterpolation(Real energy, + unsigned int i, + unsigned int j, + unsigned int l) const +{ + unsigned int index = energyIndex(energy); + if (index == 0) + return _displacement_function[0][mapIndex(i, j, l)]; + + return linearInterpolationHelper(energy, index, i, j, l); +} + +Real +PolyatomicDisplacementFunctionBase::linearInterpolationHelper(Real energy, unsigned int index, + unsigned int i, unsigned int j, unsigned int l) const +{ + unsigned int k = mapIndex(i, j, l); + + // linear interpolation + Real e1 = _energy_history[index - 1]; + Real e2 = _energy_history[index]; + Real v1 = _displacement_function[index - 1][k]; + Real v2 = _displacement_function[index][k]; + + return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); +} + #endif From dd815ec94bdf86a97e8385515a5e18486a6653f4 Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Tue, 7 Jan 2020 16:46:43 -0700 Subject: [PATCH 2/6] Add inverse mapping function (#410) --- include/utils/PolyatomicDamageEnergyFunction.h | 7 +++++++ .../PolyatomicDisplacementDerivativeFunction.h | 8 +++++++- include/utils/PolyatomicDisplacementFunction.h | 8 +++++++- include/utils/PolyatomicDisplacementFunctionBase.h | 8 ++++++-- src/utils/PolyatomicDamageEnergyFunction.C | 13 ++++++++++++- .../PolyatomicDisplacementDerivativeFunction.C | 12 ++++++++++++ src/utils/PolyatomicDisplacementFunction.C | 11 +++++++++++ 7 files changed, 62 insertions(+), 5 deletions(-) diff --git a/include/utils/PolyatomicDamageEnergyFunction.h b/include/utils/PolyatomicDamageEnergyFunction.h index 9586305e..c471afd4 100644 --- a/include/utils/PolyatomicDamageEnergyFunction.h +++ b/include/utils/PolyatomicDamageEnergyFunction.h @@ -39,10 +39,17 @@ class PolyatomicDamageEnergyFunction : public PolyatomicDisplacementFunctionBase Real _taylor_series_threshold = 1e2; protected: + /// override the mapIndex function: flattens ijl to single index n unsigned int mapIndex(unsigned int i, unsigned int /*j*/, unsigned int /*l*/) const override { return i; }; + + /// override the inverseMapIndex function: retrieves n given ijk + void inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const override; }; #endif // GSL_ENABLED diff --git a/include/utils/PolyatomicDisplacementDerivativeFunction.h b/include/utils/PolyatomicDisplacementDerivativeFunction.h index 87351a2b..8054b6cd 100644 --- a/include/utils/PolyatomicDisplacementDerivativeFunction.h +++ b/include/utils/PolyatomicDisplacementDerivativeFunction.h @@ -39,12 +39,18 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu Real source(Real energy, unsigned int i, unsigned int j, unsigned int l); protected: - /// maps triple index theta_ijl to single theta_n; l runs faster than j runs faster than i + /// override the mapIndex function: flattens ijl to single index n unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const override { return i + j * _n_species + l * _n_species * _n_species; }; + /// override the inverseMapIndex function: retrieves n given ijk + void inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const override; + /// the source term in the NRT equatons for the derivative requires the solution of the equations themselves const PolyatomicDisplacementFunction * _net_displacement_function; }; diff --git a/include/utils/PolyatomicDisplacementFunction.h b/include/utils/PolyatomicDisplacementFunction.h index 770c627c..24b43350 100644 --- a/include/utils/PolyatomicDisplacementFunction.h +++ b/include/utils/PolyatomicDisplacementFunction.h @@ -32,12 +32,18 @@ class PolyatomicDisplacementFunction : public PolyatomicDisplacementFunctionBase Real integralTypeII(Real energy, unsigned int i, unsigned int j, unsigned int k) const; protected: - /// maps double index nu_ij to single index nu_n; j runs faster than i + /// override the mapIndex function: flattens ijl to single index n unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int /*l*/) const override { return i + j * _n_species; }; + /// override the inverseMapIndex function: retrieves n given ijk + void inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const override; + /// is the total damage function computed bool _total_displacement_function; }; diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h index e247e317..13548b68 100644 --- a/include/utils/PolyatomicDisplacementFunctionBase.h +++ b/include/utils/PolyatomicDisplacementFunctionBase.h @@ -57,7 +57,8 @@ class PolyatomicDisplacementFunctionBase ///@} /// linear interpolation of the damage function - Real linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + Real + linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; /// gets stopping power for a given species and energy; non-const because it uses _ions so no need to construct ion Real stoppingPower(unsigned int species, Real energy); @@ -72,6 +73,10 @@ class PolyatomicDisplacementFunctionBase /// this function flattens the arrays i, j, l to a single index virtual unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const = 0; + /// this function unflattens n to indices i, j, l + virtual void + inverseMapIndex(unsigned int n, unsigned int & i, unsigned int & j, unsigned int & l) const = 0; + /// computes the integral int_0^t dT T * d(sigma_ij) / dT for species combination i, j and small t Real weightedScatteringIntegral(Real energy, Real energy_limit, unsigned int i, unsigned int j) const; @@ -94,7 +99,6 @@ class PolyatomicDisplacementFunctionBase Real linearInterpolationHelper( Real energy, unsigned int index, unsigned int i, unsigned int j, unsigned int l) const; - /// damage function type [nij and gij, respectively in PK JNM 101, 1981; or nu_i JNM 88, (1980)] nrt_type _damage_function_type; diff --git a/src/utils/PolyatomicDamageEnergyFunction.C b/src/utils/PolyatomicDamageEnergyFunction.C index 10421426..1aa7d827 100644 --- a/src/utils/PolyatomicDamageEnergyFunction.C +++ b/src/utils/PolyatomicDamageEnergyFunction.C @@ -156,10 +156,21 @@ PolyatomicDamageEnergyFunction::integralTypeII(Real energy, unsigned int i, unsi { Real recoil_energy = scale * (_quad_points[qp] + 1) + lower; integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) * - (linearInterpolation( energy - recoil_energy, i) - current_value); + (linearInterpolation(energy - recoil_energy, i) - current_value); } } return integral; } +void +PolyatomicDamageEnergyFunction::inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const +{ + i = n; + j = 0; + l = 0; +} + #endif diff --git a/src/utils/PolyatomicDisplacementDerivativeFunction.C b/src/utils/PolyatomicDisplacementDerivativeFunction.C index fb74011b..5b161dc2 100644 --- a/src/utils/PolyatomicDisplacementDerivativeFunction.C +++ b/src/utils/PolyatomicDisplacementDerivativeFunction.C @@ -183,4 +183,16 @@ PolyatomicDisplacementDerivativeFunction::source(Real energy, d_net_dE * stoppingPowerDerivative(i, l, energy); } +void +PolyatomicDisplacementDerivativeFunction::inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const +{ + unsigned int t = n % (_n_species * _n_species); + i = t % _n_species; + j = (t - i) / _n_species; + l = (n - t) / (_n_species * _n_species); +} + #endif diff --git a/src/utils/PolyatomicDisplacementFunction.C b/src/utils/PolyatomicDisplacementFunction.C index 0bafb05e..eec82368 100644 --- a/src/utils/PolyatomicDisplacementFunction.C +++ b/src/utils/PolyatomicDisplacementFunction.C @@ -161,4 +161,15 @@ PolyatomicDisplacementFunction::integralTypeII(Real energy, return integral; } +void +PolyatomicDisplacementFunction::inverseMapIndex(unsigned int n, + unsigned int & i, + unsigned int & j, + unsigned int & l) const +{ + i = n % _n_species; + j = (n - i) / _n_species; + l = 0; +} + #endif From 75ebc128224e9197115bf9bd211eec708014b99e Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Tue, 7 Jan 2020 17:10:59 -0700 Subject: [PATCH 3/6] Add capability to compute integral of damage function (#410) --- .../PolyatomicDisplacementFunctionBase.h | 7 +++ .../PolyatomicDisplacementFunctionBase.C | 61 ++++++++++++++++++- 2 files changed, 66 insertions(+), 2 deletions(-) diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h index 13548b68..5eac62ed 100644 --- a/include/utils/PolyatomicDisplacementFunctionBase.h +++ b/include/utils/PolyatomicDisplacementFunctionBase.h @@ -44,6 +44,9 @@ class PolyatomicDisplacementFunctionBase /// computes the displacement function from current last energy to Emax void advanceDisplacements(Real Emax); + /// computes the integral of the displacement function + void computeDisplacementFunctionIntegral(); + ///@{ some getters needed for accessing this pointer in odeRHS unsigned int nSpecies() const { return _n_species; } unsigned int problemSize() const { return _problem_size; } @@ -60,6 +63,10 @@ class PolyatomicDisplacementFunctionBase Real linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + /// linear interpolation of the integral of the damage function + Real + linearInterpolationIntegralDamageFunction(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + /// gets stopping power for a given species and energy; non-const because it uses _ions so no need to construct ion Real stoppingPower(unsigned int species, Real energy); diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C index eb37898b..f6e6d86d 100644 --- a/src/utils/PolyatomicDisplacementFunctionBase.C +++ b/src/utils/PolyatomicDisplacementFunctionBase.C @@ -111,6 +111,42 @@ PolyatomicDisplacementFunctionBase::advanceDisplacements(Real new_energy) std ::cout << "gsl_odeiv2_driver_apply returned error code " << status << std::endl; } +void +PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral() +{ + // clear the _displacement_function_integral array because this might + // have been called before; not the most efficient if the user keeps on calling + // this function but they are at their own peril and should know to call it + // after finishing the computation of disp function + _displacement_function_integral.clear(); + + // resize the array; note that initial value is zero (entry for + // _displacement_function_integral[0]) + _displacement_function_integral.resize(nEnergySteps(), std::vector(_problem_size)); + + for (unsigned int e = 1; e < nEnergySteps(); ++e) + { + _displacement_function_integral.push_back(std::vector(_problem_size)); + + // set energy points for integration + Real lower = energyPoint(e - 1); + Real upper = energyPoint(e); + Real f = 0.5 * (upper - lower); + for (unsigned int qp = 0; qp < _quad_order; ++qp) + { + Real energy = f * (_quad_points[qp] + 1) + lower; + Real w = f * _quad_weights[qp]; + + for (unsigned int n = 0; n < _problem_size; ++n) + { + unsigned int i, j, l; + inverseMapIndex(n, i, j, l); + _displacement_function_integral[e][n] += w * linearInterpolation(energy, i, j, l); + } + } + } +} + Real PolyatomicDisplacementFunctionBase::stoppingPower(unsigned int species, Real energy) { @@ -281,8 +317,8 @@ PolyatomicDisplacementFunctionBase::linearInterpolation(Real energy, } Real -PolyatomicDisplacementFunctionBase::linearInterpolationHelper(Real energy, unsigned int index, - unsigned int i, unsigned int j, unsigned int l) const +PolyatomicDisplacementFunctionBase::linearInterpolationHelper( + Real energy, unsigned int index, unsigned int i, unsigned int j, unsigned int l) const { unsigned int k = mapIndex(i, j, l); @@ -295,4 +331,25 @@ PolyatomicDisplacementFunctionBase::linearInterpolationHelper(Real energy, unsig return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); } +Real +PolyatomicDisplacementFunctionBase::linearInterpolationIntegralDamageFunction(Real energy, + unsigned int i, + unsigned int j, + unsigned int l) const +{ + unsigned int index = energyIndex(energy); + unsigned int k = mapIndex(i, j, l); + + if (index == 0) + return _displacement_function_integral[0][mapIndex(i, j, l)]; + + // linear interpolation + Real e1 = _energy_history[index - 1]; + Real e2 = _energy_history[index]; + Real v1 = _displacement_function_integral[index - 1][k]; + Real v2 = _displacement_function_integral[index][k]; + + return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); +} + #endif From d5d69778799cf2938d269731b46fbd0654dcb12e Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Tue, 7 Jan 2020 17:21:34 -0700 Subject: [PATCH 4/6] Style changes (#410) --- include/utils/PolyatomicDisplacementFunctionBase.h | 6 ++++-- src/userobjects/PolyatomicRecoil.C | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h index 5eac62ed..abaa3821 100644 --- a/include/utils/PolyatomicDisplacementFunctionBase.h +++ b/include/utils/PolyatomicDisplacementFunctionBase.h @@ -64,8 +64,10 @@ class PolyatomicDisplacementFunctionBase linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; /// linear interpolation of the integral of the damage function - Real - linearInterpolationIntegralDamageFunction(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + Real linearInterpolationIntegralDamageFunction(Real energy, + unsigned int i, + unsigned int j = 0, + unsigned int l = 0) const; /// gets stopping power for a given species and energy; non-const because it uses _ions so no need to construct ion Real stoppingPower(unsigned int species, Real energy); diff --git a/src/userobjects/PolyatomicRecoil.C b/src/userobjects/PolyatomicRecoil.C index fb82aa58..b9eae0fd 100644 --- a/src/userobjects/PolyatomicRecoil.C +++ b/src/userobjects/PolyatomicRecoil.C @@ -210,7 +210,7 @@ PolyatomicRecoil::finalize() ++derivative) displacement_file << "," << _padf_derivative->linearInterpolation( - _padf_derivative->energyPoint(n), + _padf_derivative->energyPoint(n), projectile, target, derivative); @@ -242,7 +242,8 @@ PolyatomicRecoil::finalize() displacement_file << _padf->energyPoint(n); for (unsigned int projectile = 0; projectile < _padf->nSpecies(); ++projectile) displacement_file << "," - << energy_function->linearInterpolation(_padf->energyPoint(n), projectile); + << energy_function->linearInterpolation(_padf->energyPoint(n), + projectile); displacement_file << std::endl; } } From e0056431e6be413a0fc91d095ac5f6651f126900 Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Tue, 7 Jan 2020 18:12:44 -0700 Subject: [PATCH 5/6] Bug fix and test (#410) --- .../PolyatomicDisplacementFunctionBase.C | 4 +- unit/src/DisplacementFunctionTest.C | 59 +++++++++++++++++++ 2 files changed, 62 insertions(+), 1 deletion(-) create mode 100644 unit/src/DisplacementFunctionTest.C diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C index f6e6d86d..672a0934 100644 --- a/src/utils/PolyatomicDisplacementFunctionBase.C +++ b/src/utils/PolyatomicDisplacementFunctionBase.C @@ -126,7 +126,8 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral() for (unsigned int e = 1; e < nEnergySteps(); ++e) { - _displacement_function_integral.push_back(std::vector(_problem_size)); + // initialize the integral to the integral of the previous step + _displacement_function_integral[e] = _displacement_function_integral[e - 1]; // set energy points for integration Real lower = energyPoint(e - 1); @@ -141,6 +142,7 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral() { unsigned int i, j, l; inverseMapIndex(n, i, j, l); + _displacement_function_integral[e][n] += w * linearInterpolation(energy, i, j, l); } } diff --git a/unit/src/DisplacementFunctionTest.C b/unit/src/DisplacementFunctionTest.C new file mode 100644 index 00000000..416db1a4 --- /dev/null +++ b/unit/src/DisplacementFunctionTest.C @@ -0,0 +1,59 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MOOSE - Multiphysics Object Oriented Simulation Environment */ +/* */ +/* (c) 2010 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#include "PolyatomicDisplacementFunction.h" +#include +#include + +TEST(DisplacementFunctionTest, integralDispFunction) +{ + std::vector Z = {6, 14}; + std::vector A = {12.0, 28.0}; + std::vector N = {0.5, 0.5}; + std::vector threshold = {16.3, 92.6}; + std::vector poly_mat; + for (unsigned int j = 0; j < Z.size(); ++j) + { + MyTRIM_NS::Element element; + element._Z = Z[j]; + element._m = A[j]; + element._t = N[j]; + element._Edisp = threshold[j]; + element._Elbind = 0.0; + poly_mat.push_back(element); + } + + PolyatomicDisplacementFunction padf(poly_mat, NET); + + Real energy = padf.minEnergy(); + for (unsigned int j = 0; j < 51; ++j) + { + energy += 1; + padf.advanceDisplacements(energy); + } + padf.computeDisplacementFunctionIntegral(); + unsigned int final = padf.nEnergySteps() - 1; + + Real integral = 0; + for (unsigned int j = 0; j <= final; ++j) + { + Real mp = (j == 0 || j == final) ? 0.5 : 1; + Real energy = padf.energyPoint(j); + integral += mp * padf.linearInterpolation(energy, 0, 0, 0); + } + EXPECT_NEAR(integral, + padf.linearInterpolationIntegralDamageFunction(padf.energyPoint(final), 0, 0, 0), + 1e-4) + << "Displacement function integration result is wrong"; +} From 607ba56a20677c024fc279673d18b4873f359a9f Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Thu, 9 Jan 2020 09:29:53 -0700 Subject: [PATCH 6/6] Address comments (#410) --- .../utils/PolyatomicDamageEnergyFunction.h | 13 ---- ...PolyatomicDisplacementDerivativeFunction.h | 12 ---- .../utils/PolyatomicDisplacementFunction.h | 12 ---- .../PolyatomicDisplacementFunctionBase.h | 18 ++--- src/utils/PolyatomicDamageEnergyFunction.C | 16 ++--- ...PolyatomicDisplacementDerivativeFunction.C | 15 +---- src/utils/PolyatomicDisplacementFunction.C | 14 +--- .../PolyatomicDisplacementFunctionBase.C | 65 ++++++++++++------- 8 files changed, 59 insertions(+), 106 deletions(-) diff --git a/include/utils/PolyatomicDamageEnergyFunction.h b/include/utils/PolyatomicDamageEnergyFunction.h index c471afd4..f3da8fb3 100644 --- a/include/utils/PolyatomicDamageEnergyFunction.h +++ b/include/utils/PolyatomicDamageEnergyFunction.h @@ -37,19 +37,6 @@ class PolyatomicDamageEnergyFunction : public PolyatomicDisplacementFunctionBase * by nu_i(E-T) - nu_i(E) \approx -T d(nu_i) / dE. */ Real _taylor_series_threshold = 1e2; - -protected: - /// override the mapIndex function: flattens ijl to single index n - unsigned int mapIndex(unsigned int i, unsigned int /*j*/, unsigned int /*l*/) const override - { - return i; - }; - - /// override the inverseMapIndex function: retrieves n given ijk - void inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const override; }; #endif // GSL_ENABLED diff --git a/include/utils/PolyatomicDisplacementDerivativeFunction.h b/include/utils/PolyatomicDisplacementDerivativeFunction.h index 8054b6cd..4bad85d6 100644 --- a/include/utils/PolyatomicDisplacementDerivativeFunction.h +++ b/include/utils/PolyatomicDisplacementDerivativeFunction.h @@ -39,18 +39,6 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu Real source(Real energy, unsigned int i, unsigned int j, unsigned int l); protected: - /// override the mapIndex function: flattens ijl to single index n - unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const override - { - return i + j * _n_species + l * _n_species * _n_species; - }; - - /// override the inverseMapIndex function: retrieves n given ijk - void inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const override; - /// the source term in the NRT equatons for the derivative requires the solution of the equations themselves const PolyatomicDisplacementFunction * _net_displacement_function; }; diff --git a/include/utils/PolyatomicDisplacementFunction.h b/include/utils/PolyatomicDisplacementFunction.h index 24b43350..2a207343 100644 --- a/include/utils/PolyatomicDisplacementFunction.h +++ b/include/utils/PolyatomicDisplacementFunction.h @@ -32,18 +32,6 @@ class PolyatomicDisplacementFunction : public PolyatomicDisplacementFunctionBase Real integralTypeII(Real energy, unsigned int i, unsigned int j, unsigned int k) const; protected: - /// override the mapIndex function: flattens ijl to single index n - unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int /*l*/) const override - { - return i + j * _n_species; - }; - - /// override the inverseMapIndex function: retrieves n given ijk - void inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const override; - /// is the total damage function computed bool _total_displacement_function; }; diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h index abaa3821..d640da27 100644 --- a/include/utils/PolyatomicDisplacementFunctionBase.h +++ b/include/utils/PolyatomicDisplacementFunctionBase.h @@ -59,9 +59,11 @@ class PolyatomicDisplacementFunctionBase Real numberDensity(unsigned int i) const { return _material->_element[i]._t * _material->_arho; } ///@} - /// linear interpolation of the damage function + ///@{ linear interpolation of the damage function Real linearInterpolation(Real energy, unsigned int i, unsigned int j = 0, unsigned int l = 0) const; + Real linearInterpolationFromFlatIndex(Real energy, unsigned int index) const; + ///@} /// linear interpolation of the integral of the damage function Real linearInterpolationIntegralDamageFunction(Real energy, @@ -79,12 +81,8 @@ class PolyatomicDisplacementFunctionBase unsigned int energyIndex(Real energy) const; protected: - /// this function flattens the arrays i, j, l to a single index - virtual unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const = 0; - - /// this function unflattens n to indices i, j, l - virtual void - inverseMapIndex(unsigned int n, unsigned int & i, unsigned int & j, unsigned int & l) const = 0; + /// override the mapIndex function: flattens ijl to single index n + unsigned int mapIndex(unsigned int i, unsigned int j = 0, unsigned int l = 0) const; /// computes the integral int_0^t dT T * d(sigma_ij) / dT for species combination i, j and small t Real @@ -105,8 +103,7 @@ class PolyatomicDisplacementFunctionBase nonCaptureProbability(unsigned int i, unsigned int k, Real energy, Real recoil_energy) const; /// a helper function called from linearInterpolation - Real linearInterpolationHelper( - Real energy, unsigned int index, unsigned int i, unsigned int j, unsigned int l) const; + Real linearInterpolationHelper(Real energy, unsigned int energy_index, unsigned int index) const; /// damage function type [nij and gij, respectively in PK JNM 101, 1981; or nu_i JNM 88, (1980)] nrt_type _damage_function_type; @@ -117,6 +114,9 @@ class PolyatomicDisplacementFunctionBase /// the number of different species in the material unsigned int _n_species; + /// the number of species indices of this object + unsigned int _n_indices = 0; + /// the size of the problem = _n_species**2 for TOTAL & NET, _n_species for ENERGY unsigned int _problem_size; diff --git a/src/utils/PolyatomicDamageEnergyFunction.C b/src/utils/PolyatomicDamageEnergyFunction.C index 1aa7d827..8b6521b3 100644 --- a/src/utils/PolyatomicDamageEnergyFunction.C +++ b/src/utils/PolyatomicDamageEnergyFunction.C @@ -32,6 +32,9 @@ PolyatomicDamageEnergyFunction::PolyatomicDamageEnergyFunction( if (damage_function_type != ENERGY) throw std::exception(); + // set the number of indices + _n_indices = 1; + // set up the gsl ODE machinery auto func = &PolyatomicDamageEnergyFunction::odeRHS; _sys = {func, NULL, _problem_size, this}; @@ -110,7 +113,7 @@ PolyatomicDamageEnergyFunction::integralTypeI(Real energy, unsigned int i, unsig { Real recoil_energy = scale * (_quad_points[qp] + 1) + lower; integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) * - linearInterpolationHelper(recoil_energy, l, j, 0, 0); + linearInterpolation(recoil_energy, j); } } return integral; @@ -162,15 +165,4 @@ PolyatomicDamageEnergyFunction::integralTypeII(Real energy, unsigned int i, unsi return integral; } -void -PolyatomicDamageEnergyFunction::inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const -{ - i = n; - j = 0; - l = 0; -} - #endif diff --git a/src/utils/PolyatomicDisplacementDerivativeFunction.C b/src/utils/PolyatomicDisplacementDerivativeFunction.C index 5b161dc2..dd053224 100644 --- a/src/utils/PolyatomicDisplacementDerivativeFunction.C +++ b/src/utils/PolyatomicDisplacementDerivativeFunction.C @@ -34,6 +34,9 @@ PolyatomicDisplacementDerivativeFunction::PolyatomicDisplacementDerivativeFuncti if (damage_function_type != NET_DERIVATIVE) throw std::exception(); + // set the number of indices + _n_indices = 3; + // set up the gsl ODE machinery auto func = &PolyatomicDisplacementDerivativeFunction::odeRHS; _sys = {func, NULL, _problem_size, this}; @@ -183,16 +186,4 @@ PolyatomicDisplacementDerivativeFunction::source(Real energy, d_net_dE * stoppingPowerDerivative(i, l, energy); } -void -PolyatomicDisplacementDerivativeFunction::inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const -{ - unsigned int t = n % (_n_species * _n_species); - i = t % _n_species; - j = (t - i) / _n_species; - l = (n - t) / (_n_species * _n_species); -} - #endif diff --git a/src/utils/PolyatomicDisplacementFunction.C b/src/utils/PolyatomicDisplacementFunction.C index eec82368..ec507495 100644 --- a/src/utils/PolyatomicDisplacementFunction.C +++ b/src/utils/PolyatomicDisplacementFunction.C @@ -33,6 +33,9 @@ PolyatomicDisplacementFunction::PolyatomicDisplacementFunction( if (damage_function_type != TOTAL && damage_function_type != NET) throw std::exception(); + // set the number of indices + _n_indices = 2; + // set up the gsl ODE machinery auto func = &PolyatomicDisplacementFunction::odeRHS; _sys = {func, NULL, _problem_size, this}; @@ -161,15 +164,4 @@ PolyatomicDisplacementFunction::integralTypeII(Real energy, return integral; } -void -PolyatomicDisplacementFunction::inverseMapIndex(unsigned int n, - unsigned int & i, - unsigned int & j, - unsigned int & l) const -{ - i = n % _n_species; - j = (n - i) / _n_species; - l = 0; -} - #endif diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C index 672a0934..98a3d86b 100644 --- a/src/utils/PolyatomicDisplacementFunctionBase.C +++ b/src/utils/PolyatomicDisplacementFunctionBase.C @@ -139,12 +139,7 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral() Real w = f * _quad_weights[qp]; for (unsigned int n = 0; n < _problem_size; ++n) - { - unsigned int i, j, l; - inverseMapIndex(n, i, j, l); - - _displacement_function_integral[e][n] += w * linearInterpolation(energy, i, j, l); - } + _displacement_function_integral[e][n] += w * linearInterpolationFromFlatIndex(energy, n); } } } @@ -311,24 +306,35 @@ PolyatomicDisplacementFunctionBase::linearInterpolation(Real energy, unsigned int j, unsigned int l) const { - unsigned int index = energyIndex(energy); - if (index == 0) - return _displacement_function[0][mapIndex(i, j, l)]; + unsigned int energy_index = energyIndex(energy); + unsigned int index = mapIndex(i, j, l); + if (energy_index == 0) + return _displacement_function[0][index]; - return linearInterpolationHelper(energy, index, i, j, l); + return linearInterpolationHelper(energy, energy_index, index); } Real -PolyatomicDisplacementFunctionBase::linearInterpolationHelper( - Real energy, unsigned int index, unsigned int i, unsigned int j, unsigned int l) const +PolyatomicDisplacementFunctionBase::linearInterpolationFromFlatIndex(Real energy, + unsigned int index) const { - unsigned int k = mapIndex(i, j, l); + unsigned int energy_index = energyIndex(energy); + if (energy_index == 0) + return _displacement_function[0][index]; + + return linearInterpolationHelper(energy, energy_index, index); +} +Real +PolyatomicDisplacementFunctionBase::linearInterpolationHelper(Real energy, + unsigned int energy_index, + unsigned int index) const +{ // linear interpolation - Real e1 = _energy_history[index - 1]; - Real e2 = _energy_history[index]; - Real v1 = _displacement_function[index - 1][k]; - Real v2 = _displacement_function[index][k]; + Real e1 = _energy_history[energy_index - 1]; + Real e2 = _energy_history[energy_index]; + Real v1 = _displacement_function[energy_index - 1][index]; + Real v2 = _displacement_function[energy_index][index]; return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); } @@ -339,19 +345,28 @@ PolyatomicDisplacementFunctionBase::linearInterpolationIntegralDamageFunction(Re unsigned int j, unsigned int l) const { - unsigned int index = energyIndex(energy); - unsigned int k = mapIndex(i, j, l); + unsigned int energy_index = energyIndex(energy); + unsigned int index = mapIndex(i, j, l); - if (index == 0) - return _displacement_function_integral[0][mapIndex(i, j, l)]; + if (energy_index == 0) + return _displacement_function_integral[0][index]; // linear interpolation - Real e1 = _energy_history[index - 1]; - Real e2 = _energy_history[index]; - Real v1 = _displacement_function_integral[index - 1][k]; - Real v2 = _displacement_function_integral[index][k]; + Real e1 = _energy_history[energy_index - 1]; + Real e2 = _energy_history[energy_index]; + Real v1 = _displacement_function_integral[energy_index - 1][index]; + Real v2 = _displacement_function_integral[energy_index][index]; return v1 + (energy - e1) / (e2 - e1) * (v2 - v1); } +unsigned int +PolyatomicDisplacementFunctionBase::mapIndex(unsigned int i, unsigned int j, unsigned int l) const +{ + assert(_n_indices > 0); + assert(j == 0 || _n_indices > 1); + assert(l == 0 || _n_indices > 2); + return i + j * _n_species + l * _n_species * _n_species; +}; + #endif