Refactor mapIndex and linearInterpolation (idaholab#410)

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

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;
   };

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;

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<Real> _energy_history;
 
-  /**
-   * The displacement values nu_ij flattened into 1D array
-   */
+  /// The displacement values nu_ij flattened into 1D array
   std::vector<std::vector<Real>> _displacement_function;
 
+  /// The integral of the displacement function over energy
+  std::vector<std::vector<Real>> _displacement_function_integral;
+
   /// Ecap_ik average residual energy of type i atom to not be trapped at k-site
   std::vector<std::vector<Real>> _Ecap;
 

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

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<MooseEnum>("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;
         }
       }

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

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,41 +174,13 @@ 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) +
          _net_displacement_function->integralTypeII(energy, i, j, l) -
          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

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