diff --git a/include/utils/PolyatomicDamageEnergyFunction.h b/include/utils/PolyatomicDamageEnergyFunction.h index f3da8fb3..005f4ebf 100644 --- 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> Ecap = {{}}); - static int odeRHS(Real energy, const Real disp[], Real f[], void * params); + virtual std::vector 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 index 4bad85d6..c4db7f86 100644 --- a/include/utils/PolyatomicDisplacementDerivativeFunction.h +++ b/include/utils/PolyatomicDisplacementDerivativeFunction.h @@ -25,7 +25,7 @@ class PolyatomicDisplacementDerivativeFunction : public PolyatomicDisplacementFu const PolyatomicDisplacementFunction * net_disp, std::vector> Ecap = {{}}); - static int odeRHS(Real energy, const Real disp[], Real f[], void * params); + virtual std::vector 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 index 2a207343..ca712355 100644 --- 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> Ecap = {{}}); - static int odeRHS(Real energy, const Real disp[], Real f[], void * params); + virtual std::vector 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 index d640da27..413015d2 100644 --- 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 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 index 8b6521b3..036d3fb0 100644 --- 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 +PolyatomicDamageEnergyFunction::getRHS(Real energy) { - (void)disp; - PolyatomicDamageEnergyFunction * padf = (PolyatomicDamageEnergyFunction *)params; - for (unsigned int i = 0; i < padf->nSpecies(); ++i) + std::vector 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 index dd053224..2714d842 100644 --- 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::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 +PolyatomicDisplacementDerivativeFunction::getRHS(Real energy) { - (void)disp; - PolyatomicDisplacementDerivativeFunction * padf = - (PolyatomicDisplacementDerivativeFunction *)params; - for (unsigned int i = 0; i < padf->nSpecies(); ++i) + std::vector 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 index ec507495..34205466 100644 --- 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::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 +PolyatomicDisplacementFunction::getRHS(Real energy) { - (void)disp; - PolyatomicDisplacementFunction * padf = (PolyatomicDisplacementFunction *)params; - for (unsigned int i = 0; i < padf->nSpecies(); ++i) + std::vector 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 index 98a3d86b..c108fb33 100644 --- a/src/utils/PolyatomicDisplacementFunctionBase.C +++ b/src/utils/PolyatomicDisplacementFunctionBase.C @@ -84,6 +84,11 @@ 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() @@ -91,6 +96,20 @@ 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 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) {