From 36110b1085b095a81379786efe9b7401678c1604 Mon Sep 17 00:00:00 2001 From: Sebastian Schunert Date: Thu, 9 Jan 2020 09:29:53 -0700 Subject: [PATCH] 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..07b5f34c 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 = 1; + // 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