Merge pull request idaholab#412 from snschune/refactor_poly_atomic_di…

…sp_func_410

Refactor polyatomic displacement function objects

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

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;
@@ -45,12 +39,6 @@ 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
-  unsigned int mapIndex(unsigned int i, unsigned int j, unsigned int l) const
-  {
-    return i + j * _n_species + l * _n_species * _n_species;
-  };
-
   /// the source term in the NRT equatons for the derivative requires the solution of the equations themselves
   const PolyatomicDisplacementFunction * _net_displacement_function;
 };

include/utils/PolyatomicDisplacementFunction.h

@@ -25,21 +25,13 @@ 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;
 
   /// computes terms 2 & 3  in Parkin-Coulter expression (nu_i(E-T) & nu_i(E) terms)
   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
-  unsigned int mapIndex(unsigned int i, unsigned int j) const { return i + j * _n_species; };
-
   /// is the total damage function computed
   bool _total_displacement_function;
 };

include/utils/PolyatomicDisplacementFunctionBase.h

@@ -41,8 +41,12 @@ class PolyatomicDisplacementFunctionBase
 
   virtual ~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; }
@@ -55,6 +59,18 @@ 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;
+  Real linearInterpolationFromFlatIndex(Real energy, unsigned int index) 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);
 
@@ -65,6 +81,9 @@ class PolyatomicDisplacementFunctionBase
   unsigned int energyIndex(Real energy) const;
 
 protected:
+  /// 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
   weightedScatteringIntegral(Real energy, Real energy_limit, unsigned int i, unsigned int j) const;
@@ -83,6 +102,9 @@ 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 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;
 
@@ -92,17 +114,21 @@ 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;
 
   /// 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,8 @@ 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

@@ -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) *
-                  linearInterpolation(j, recoil_energy, l);
+                  linearInterpolation(recoil_energy, j);
     }
   }
   return integral;
@@ -123,11 +126,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 +159,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

@@ -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};
@@ -104,7 +107,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 +121,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 +156,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 +177,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