Skip to content

Commit

Permalink
Address comments (idaholab#410)
Browse files Browse the repository at this point in the history
  • Loading branch information
Sebastian Schunert committed Jan 9, 2020
1 parent e005643 commit 36110b1
Show file tree
Hide file tree
Showing 8 changed files with 59 additions and 106 deletions.
13 changes: 0 additions & 13 deletions include/utils/PolyatomicDamageEnergyFunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
12 changes: 0 additions & 12 deletions include/utils/PolyatomicDisplacementDerivativeFunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};
Expand Down
12 changes: 0 additions & 12 deletions include/utils/PolyatomicDisplacementFunction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};
Expand Down
18 changes: 9 additions & 9 deletions include/utils/PolyatomicDisplacementFunctionBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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;

Expand Down
16 changes: 4 additions & 12 deletions src/utils/PolyatomicDamageEnergyFunction.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
15 changes: 3 additions & 12 deletions src/utils/PolyatomicDisplacementDerivativeFunction.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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
14 changes: 3 additions & 11 deletions src/utils/PolyatomicDisplacementFunction.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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
65 changes: 40 additions & 25 deletions src/utils/PolyatomicDisplacementFunctionBase.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
}
Expand Down Expand Up @@ -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);
}
Expand All @@ -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

0 comments on commit 36110b1

Please sign in to comment.