Move GSL interface to base class and have static interface function c…

…all callback (idaholab#410)
snschune · Jan 10, 2020 · 6dbccdd · 6dbccdd
1 parent 259f517
commit 6dbccdd
Show file tree

Hide file tree

Showing 8 changed files with 63 additions and 64 deletions.
diff --git a/include/utils/PolyatomicDamageEnergyFunction.h b/include/utils/PolyatomicDamageEnergyFunction.h
@@ -23,7 +23,7 @@ class PolyatomicDamageEnergyFunction : public PolyatomicDisplacementFunctionBase
                                  nrt_type damage_function_type,
                                  std::vector<std::vector<Real>> Ecap = {{}});
 
-  static int odeRHS(Real energy, const Real disp[], Real f[], void * params);
+  virtual std::vector<Real> getRHS(Real energy) override;
 
   /// a getter needed for accessing this pointer in odeRHS
   Real taylorSeriesThreshold() const { return _taylor_series_threshold; }

diff --git a/include/utils/PolyatomicDisplacementDerivativeFunction.h b/include/utils/PolyatomicDisplacementDerivativeFunction.h
@@ -25,7 +25,7 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu
                                            const PolyatomicDisplacementFunction * net_disp,
                                            std::vector<std::vector<Real>> Ecap = {{}});
 
-  static int odeRHS(Real energy, const Real disp[], Real f[], void * params);
+  virtual std::vector<Real> getRHS(Real energy) override;
 
   /// computes term 1 in Parkin-Coulter expression nu_k(T - Eb)
   Real

diff --git a/include/utils/PolyatomicDisplacementFunction.h b/include/utils/PolyatomicDisplacementFunction.h
@@ -23,7 +23,7 @@ class PolyatomicDisplacementFunction : public PolyatomicDisplacementFunctionBase
                                  nrt_type damage_function_type,
                                  std::vector<std::vector<Real>> Ecap = {{}});
 
-  static int odeRHS(Real energy, const Real disp[], Real f[], void * params);
+  virtual std::vector<Real> getRHS(Real energy) override;
 
   /// 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;

diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h
@@ -41,6 +41,10 @@ class PolyatomicDisplacementFunctionBase
 
   virtual ~PolyatomicDisplacementFunctionBase();
 
+  static int gslInterface(Real energy, const Real disp[], Real f[], void * params);
+
+  virtual std::vector<Real> getRHS(Real energy) = 0;
+
   /// computes the displacement function from current last energy to Emax
   void advanceDisplacements(Real Emax);
 

diff --git a/src/utils/PolyatomicDamageEnergyFunction.C b/src/utils/PolyatomicDamageEnergyFunction.C
@@ -35,11 +35,6 @@ PolyatomicDamageEnergyFunction::PolyatomicDamageEnergyFunction(
   // set the number of indices
   _n_indices = 1;
 
-  // set up the gsl ODE machinery
-  auto func = &PolyatomicDamageEnergyFunction::odeRHS;
-  _sys = {func, NULL, _problem_size, this};
-  _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
-
   /*
    * The ENERGY mode has no lower cutoff, because is nu_i(0) = 0. However, this will
    * lead to issues with evaluating the scattering XS. We use Lindhard's initial condition
@@ -53,30 +48,29 @@ PolyatomicDamageEnergyFunction::PolyatomicDamageEnergyFunction(
     _displacement_function[0][i] = _energy_history[0];
 }
 
-int
-PolyatomicDamageEnergyFunction::odeRHS(Real energy, const Real disp[], Real f[], void * params)
+std::vector<Real>
+PolyatomicDamageEnergyFunction::getRHS(Real energy)
 {
-  (void)disp;
-  PolyatomicDamageEnergyFunction * padf = (PolyatomicDamageEnergyFunction *)params;
-  for (unsigned int i = 0; i < padf->nSpecies(); ++i)
+  std::vector<Real> f(_problem_size);
+  for (unsigned int i = 0; i < nSpecies(); ++i)
   {
     f[i] = 0;
-    Real denominator = padf->stoppingPower(i, energy);
+    Real denominator = stoppingPower(i, energy);
 
     // range of energies from the threshold to E
-    for (unsigned int j = 0; j < padf->nSpecies(); ++j)
+    for (unsigned int j = 0; j < nSpecies(); ++j)
     {
-      f[i] += padf->numberFraction(j) * padf->integralTypeI(energy, i, j);
+      f[i] += numberFraction(j) * integralTypeI(energy, i, j);
 
-      if (energy <= padf->taylorSeriesThreshold())
-        denominator += padf->numberFraction(j) *
-                       padf->weightedScatteringIntegral(energy, energy * padf->lambda(i, j), i, j);
+      if (energy <= taylorSeriesThreshold())
+        denominator +=
+            numberFraction(j) * weightedScatteringIntegral(energy, energy * lambda(i, j), i, j);
       else
-        f[i] += padf->numberFraction(j) * padf->integralTypeII(energy, i, j);
+        f[i] += numberFraction(j) * integralTypeII(energy, i, j);
     }
     f[i] /= denominator;
   }
-  return GSL_SUCCESS;
+  return f;
 }
 
 Real

diff --git a/src/utils/PolyatomicDisplacementDerivativeFunction.C b/src/utils/PolyatomicDisplacementDerivativeFunction.C
@@ -37,11 +37,6 @@ PolyatomicDisplacementDerivativeFunction::PolyatomicDisplacementDerivativeFuncti
   // set the number of indices
   _n_indices = 3;
 
-  // set up the gsl ODE machinery
-  auto func = &PolyatomicDisplacementDerivativeFunction::odeRHS;
-  _sys = {func, NULL, _problem_size, this};
-  _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
-
   Real Edisp_min = std::numeric_limits<Real>::max();
   for (unsigned int j = 0; j < _n_species; ++j)
     if (_material->_element[j]._Edisp < Edisp_min)
@@ -51,34 +46,28 @@ PolyatomicDisplacementDerivativeFunction::PolyatomicDisplacementDerivativeFuncti
   // note that initial conditions are 0s for theta_ijl
 }
 
-int
-PolyatomicDisplacementDerivativeFunction::odeRHS(Real energy,
-                                                 const Real disp[],
-                                                 Real f[],
-                                                 void * params)
+std::vector<Real>
+PolyatomicDisplacementDerivativeFunction::getRHS(Real energy)
 {
-  (void)disp;
-  PolyatomicDisplacementDerivativeFunction * padf =
-      (PolyatomicDisplacementDerivativeFunction *)params;
-  for (unsigned int i = 0; i < padf->nSpecies(); ++i)
+  std::vector<Real> f(_problem_size);
+  for (unsigned int i = 0; i < nSpecies(); ++i)
   {
-    Real stopping_power = padf->stoppingPower(i, energy);
-    for (unsigned int j = 0; j < padf->nSpecies(); ++j)
-      for (unsigned int l = 0; l < padf->nSpecies(); ++l)
+    Real stopping_power = stoppingPower(i, energy);
+    for (unsigned int j = 0; j < nSpecies(); ++j)
+      for (unsigned int l = 0; l < nSpecies(); ++l)
       {
         // working on the right hand side for theta_ijl
-        unsigned int n = padf->mapIndex(i, j, l);
+        unsigned int n = mapIndex(i, j, l);
         f[n] = 0;
 
-        for (unsigned int k = 0; k < padf->nSpecies(); ++k)
-          f[n] +=
-              padf->numberFraction(k) *
-              (padf->integralTypeI(energy, i, j, l, k) + padf->integralTypeII(energy, i, j, l, k)) /
-              stopping_power;
-        f[n] += padf->source(energy, i, j, l) / stopping_power;
+        for (unsigned int k = 0; k < nSpecies(); ++k)
+          f[n] += numberFraction(k) *
+                  (integralTypeI(energy, i, j, l, k) + integralTypeII(energy, i, j, l, k)) /
+                  stopping_power;
+        f[n] += source(energy, i, j, l) / stopping_power;
       }
   }
-  return GSL_SUCCESS;
+  return f;
 }
 
 Real

diff --git a/src/utils/PolyatomicDisplacementFunction.C b/src/utils/PolyatomicDisplacementFunction.C
@@ -36,11 +36,6 @@ PolyatomicDisplacementFunction::PolyatomicDisplacementFunction(
   // set the number of indices
   _n_indices = 2;
 
-  // set up the gsl ODE machinery
-  auto func = &PolyatomicDisplacementFunction::odeRHS;
-  _sys = {func, NULL, _problem_size, this};
-  _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
-
   Real Edisp_min = std::numeric_limits<Real>::max();
   for (unsigned int j = 0; j < _n_species; ++j)
     if (_material->_element[j]._Edisp < Edisp_min)
@@ -54,27 +49,25 @@ PolyatomicDisplacementFunction::PolyatomicDisplacementFunction(
       _displacement_function[0][mapIndex(i, i, 0)] = 1;
 }
 
-int
-PolyatomicDisplacementFunction::odeRHS(Real energy, const Real disp[], Real f[], void * params)
+std::vector<Real>
+PolyatomicDisplacementFunction::getRHS(Real energy)
 {
-  (void)disp;
-  PolyatomicDisplacementFunction * padf = (PolyatomicDisplacementFunction *)params;
-  for (unsigned int i = 0; i < padf->nSpecies(); ++i)
+  std::vector<Real> f(_problem_size);
+  for (unsigned int i = 0; i < nSpecies(); ++i)
   {
-    Real stopping_power = padf->stoppingPower(i, energy);
-    for (unsigned int j = 0; j < padf->nSpecies(); ++j)
+    Real stopping_power = stoppingPower(i, energy);
+    for (unsigned int j = 0; j < nSpecies(); ++j)
     {
       // working on the right hand side for nu_ij
-      unsigned int n = padf->mapIndex(i, j, 0);
+      unsigned int n = mapIndex(i, j, 0);
       f[n] = 0;
 
-      for (unsigned int k = 0; k < padf->nSpecies(); ++k)
-        f[n] += padf->numberFraction(k) *
-                (padf->integralTypeI(energy, i, j, k) + padf->integralTypeII(energy, i, j, k)) /
-                stopping_power;
+      for (unsigned int k = 0; k < nSpecies(); ++k)
+        f[n] += numberFraction(k) *
+                (integralTypeI(energy, i, j, k) + integralTypeII(energy, i, j, k)) / stopping_power;
     }
   }
-  return GSL_SUCCESS;
+  return f;
 }
 
 Real

diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C
@@ -84,13 +84,32 @@ PolyatomicDisplacementFunctionBase::PolyatomicDisplacementFunctionBase(
     _quad_weights[j] = weight;
   }
   gsl_integration_glfixed_table_free(qp_table);
+
+  // set up the gsl ODE machinery
+  auto func = &PolyatomicDisplacementFunction::gslInterface;
+  _sys = {func, NULL, _problem_size, this};
+  _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
 }
 
 PolyatomicDisplacementFunctionBase::~PolyatomicDisplacementFunctionBase()
 {
   gsl_odeiv2_driver_free(_ode_driver);
 }
 
+int
+PolyatomicDisplacementFunctionBase::gslInterface(Real energy,
+                                                 const Real disp[],
+                                                 Real f[],
+                                                 void * params)
+{
+  (void)disp;
+  PolyatomicDisplacementFunctionBase * padf = (PolyatomicDisplacementFunctionBase *)params;
+  std::vector<Real> rhs = padf->getRHS(energy);
+  for (unsigned int j = 0; j < rhs.size(); ++j)
+    f[j] = rhs[j];
+  return GSL_SUCCESS;
+}
+
 void
 PolyatomicDisplacementFunctionBase::advanceDisplacements(Real new_energy)
 {