Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 : #ifdef GSL_ENABLED
9 :
10 : #include "PolyatomicDamageEnergyFunction.h"
11 :
12 : // mytrim includes
13 : #include <mytrim/simconf.h>
14 : #include <mytrim/ion.h>
15 : #include <mytrim/element.h>
16 :
17 : // general includes
18 : #include <assert.h>
19 : #include <limits>
20 : #include <exception>
21 :
22 : // gsl includes
23 : #include "gsl/gsl_sf_legendre.h"
24 : #include "gsl/gsl_integration.h"
25 :
26 4 : PolyatomicDamageEnergyFunction::PolyatomicDamageEnergyFunction(
27 : std::vector<MyTRIM_NS::Element> polyatomic_material,
28 : nrt_type damage_function_type,
29 4 : std::vector<std::vector<Real>> Ecap)
30 8 : : PolyatomicDisplacementFunctionBase(polyatomic_material, damage_function_type, Ecap)
31 : {
32 4 : if (damage_function_type != ENERGY)
33 0 : throw std::exception();
34 :
35 : // set the number of indices
36 4 : _n_indices = 1;
37 :
38 : /*
39 : * The ENERGY mode has no lower cutoff, because is nu_i(0) = 0. However, this will
40 : * lead to issues with evaluating the scattering XS. We use Lindhard's initial condition
41 : * stated as nu(E) / E -> 1 as E -> 0 [Integral Equations Governing Radiation Effects,
42 : * Mat. Fys . Medd . Dan . Vid. Selsk . 33, no . 10 (1963)].
43 : * Threshold is set at 0.01 eV
44 : */
45 4 : _energy_history[0] = 0.01;
46 :
47 12 : for (unsigned int i = 0; i < _n_species; ++i)
48 8 : _displacement_function[0][i] = _energy_history[0];
49 4 : }
50 :
51 : std::vector<Real>
52 17472 : PolyatomicDamageEnergyFunction::getRHS(Real energy)
53 : {
54 17472 : std::vector<Real> f(_problem_size);
55 52416 : for (unsigned int i = 0; i < nSpecies(); ++i)
56 : {
57 34944 : f[i] = 0;
58 34944 : Real denominator = stoppingPower(i, energy);
59 :
60 : // range of energies from the threshold to E
61 104832 : for (unsigned int j = 0; j < nSpecies(); ++j)
62 : {
63 69888 : f[i] += numberFraction(j) * integralTypeI(energy, i, j);
64 :
65 69888 : if (energy <= taylorSeriesThreshold())
66 13440 : denominator +=
67 13440 : numberFraction(j) * weightedScatteringIntegral(energy, energy * lambda(i, j), i, j);
68 : else
69 56448 : f[i] += numberFraction(j) * integralTypeII(energy, i, j);
70 : }
71 34944 : f[i] /= denominator;
72 : }
73 17472 : return f;
74 : }
75 :
76 : Real
77 69888 : PolyatomicDamageEnergyFunction::integralTypeI(Real energy, unsigned int i, unsigned int j)
78 : {
79 69888 : Real upper_integration_limit = energy * _lambda[i][j];
80 69888 : Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][j]);
81 :
82 : // integrate up to threshold, 0, ..., threshold
83 69888 : Real integral = weightedScatteringIntegral(energy, threshold, i, j);
84 :
85 69888 : if (energy * _lambda[i][j] <= _asymptotic_threshold)
86 : return integral;
87 :
88 : // start at asymptotic threshold and integrate up to energy
89 : // the integration follows the already existing energies
90 5530176 : for (unsigned int l = 0; l < _energy_history.size(); ++l)
91 : {
92 : Real lower;
93 5460336 : if (l == 0)
94 69840 : lower = _asymptotic_threshold;
95 : else
96 5460336 : lower = std::max(_energy_history[l - 1], _asymptotic_threshold);
97 :
98 5460336 : Real upper = std::min(_energy_history[l], upper_integration_limit);
99 :
100 : // nothing to integrate
101 5460336 : if (lower > upper)
102 99424 : continue;
103 :
104 : // now integrate from lower to upper
105 5360912 : Real scale = 0.5 * (upper - lower);
106 26804560 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
107 : {
108 21443648 : Real recoil_energy = scale * (_quad_points[qp] + 1) + lower;
109 42887296 : integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) *
110 21443648 : linearInterpolation(recoil_energy, j);
111 : }
112 : }
113 : return integral;
114 : }
115 :
116 : Real
117 56448 : PolyatomicDamageEnergyFunction::integralTypeII(Real energy, unsigned int i, unsigned int j)
118 : {
119 56448 : Real upper_integration_limit = energy * _lambda[i][j];
120 56448 : Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][j]);
121 :
122 : // store the current displacement function value
123 56448 : Real current_value = linearInterpolation(energy, i);
124 :
125 : // estimate the derivative d(nu_i) / dE:
126 56448 : Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2];
127 56448 : Real derivative = (current_value - linearInterpolation(energy - dE, i)) / dE;
128 :
129 : // integrate up to threshold and multiply by estimate of the derivative
130 56448 : Real integral = -weightedScatteringIntegral(energy, threshold, i, j) * derivative;
131 :
132 56448 : if (energy * _lambda[i][j] <= _asymptotic_threshold)
133 : return integral;
134 :
135 : // integrate from threshold up to energy
136 4985024 : for (unsigned int l = 0; l < _energy_history.size(); ++l)
137 : {
138 : Real lower;
139 4928576 : if (l == 0)
140 56448 : lower = _asymptotic_threshold;
141 : else
142 4928576 : lower = std::max(_energy_history[l - 1], _asymptotic_threshold);
143 :
144 4928576 : Real upper = std::min(_energy_history[l], upper_integration_limit);
145 :
146 : // nothing to integrate
147 4928576 : if (lower > upper)
148 78008 : continue;
149 :
150 : // now integrate from lower to upper
151 4850568 : Real scale = 0.5 * (upper - lower);
152 24252840 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
153 : {
154 19402272 : Real recoil_energy = scale * (_quad_points[qp] + 1) + lower;
155 19402272 : integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) *
156 19402272 : (linearInterpolation(energy - recoil_energy, i) - current_value);
157 : }
158 : }
159 : return integral;
160 : }
161 :
162 : #endif
|
