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 "PolyatomicDisplacementDerivativeFunction.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 25 : PolyatomicDisplacementDerivativeFunction::PolyatomicDisplacementDerivativeFunction(
27 : std::vector<MyTRIM_NS::Element> polyatomic_material,
28 : nrt_type damage_function_type,
29 : const PolyatomicDisplacementFunction * net_disp,
30 25 : std::vector<std::vector<Real>> Ecap)
31 : : PolyatomicDisplacementFunctionBase(polyatomic_material, damage_function_type, Ecap),
32 50 : _net_displacement_function(net_disp)
33 : {
34 25 : if (damage_function_type != NET_DERIVATIVE)
35 0 : throw std::exception();
36 :
37 : // set the number of indices
38 25 : _n_indices = 3;
39 :
40 : Real Edisp_min = std::numeric_limits<Real>::max();
41 60 : for (unsigned int j = 0; j < _n_species; ++j)
42 35 : if (_material->_element[j]._Edisp < Edisp_min)
43 : Edisp_min = _material->_element[j]._Edisp;
44 25 : _energy_history[0] = Edisp_min;
45 :
46 : // note that initial conditions are 0s for theta_ijl
47 25 : }
48 :
49 : std::vector<Real>
50 24940 : PolyatomicDisplacementDerivativeFunction::getRHS(Real energy)
51 : {
52 24940 : std::vector<Real> f(_problem_size);
53 65730 : for (unsigned int i = 0; i < nSpecies(); ++i)
54 : {
55 40790 : Real stopping_power = stoppingPower(i, energy);
56 113280 : for (unsigned int j = 0; j < nSpecies(); ++j)
57 208380 : for (unsigned int l = 0; l < nSpecies(); ++l)
58 : {
59 : // working on the right hand side for theta_ijl
60 135890 : unsigned int n = mapIndex(i, j, l);
61 135890 : f[n] = 0;
62 :
63 398580 : for (unsigned int k = 0; k < nSpecies(); ++k)
64 262690 : f[n] += numberFraction(k) *
65 262690 : (integralTypeI(energy, i, j, l, k) + integralTypeII(energy, i, j, l, k)) /
66 : stopping_power;
67 135890 : f[n] += source(energy, i, j, l) / stopping_power;
68 : }
69 : }
70 24940 : return f;
71 : }
72 :
73 : Real
74 262690 : PolyatomicDisplacementDerivativeFunction::integralTypeI(
75 : Real energy, unsigned int i, unsigned int j, unsigned int l, unsigned int k) const
76 : {
77 262690 : Real upper_integration_limit = energy * _lambda[i][k];
78 : Real integral = 0;
79 262690 : Real Ebind = _material->_element[k]._Elbind;
80 :
81 : // the integration follows the already existing energies
82 4320963 : for (unsigned int n = 1; n < _energy_history.size(); ++n)
83 : {
84 : // adjust integration limits to be within E_{l-1} ... min(E_l, upper_integration_limit)
85 4058273 : Real lower = _energy_history[n - 1];
86 4058273 : Real upper = std::min(_energy_history[n], upper_integration_limit);
87 :
88 : // nothing to integrate
89 4058273 : if (lower > upper)
90 770224 : continue;
91 :
92 : // now integrate from lower to upper
93 3288049 : Real f = 0.5 * (upper - lower);
94 16440245 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
95 : {
96 13152196 : Real recoil_energy = f * (_quad_points[qp] + 1) + lower;
97 13152196 : integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) *
98 26304392 : displacementProbability(k, recoil_energy) *
99 13152196 : linearInterpolation(recoil_energy - Ebind, k, j, l);
100 : }
101 : }
102 262690 : return integral;
103 : }
104 :
105 : Real
106 262690 : PolyatomicDisplacementDerivativeFunction::integralTypeII(
107 : Real energy, unsigned int i, unsigned int j, unsigned int l, unsigned int k) const
108 : {
109 262690 : Real upper_integration_limit = energy * _lambda[i][k];
110 262690 : Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][k]);
111 :
112 : // store the current displacement function value
113 262690 : Real current_value = linearInterpolation(energy, i, j, l);
114 :
115 : // estimate the derivative d(nu_i) / dE:
116 262690 : Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2];
117 262690 : Real derivative = (current_value - linearInterpolation(energy - dE, i, j, l)) / dE;
118 :
119 : // integrate up to threshold and multiply by estimate of the derivative
120 262690 : Real integral = -weightedScatteringIntegral(energy, threshold, i, k) * derivative;
121 :
122 262690 : if (energy * _lambda[i][k] <= _asymptotic_threshold)
123 : return integral;
124 :
125 : // the integration follows the already existing energies
126 4583653 : for (unsigned int n = 0; n < _energy_history.size(); ++n)
127 : {
128 : // adjust integration limits to be within
129 : Real lower;
130 4320963 : if (n == 0)
131 262690 : lower = _asymptotic_threshold;
132 : else
133 4058273 : lower = std::max(_energy_history[n - 1], _asymptotic_threshold);
134 :
135 4320963 : Real upper = std::min(_energy_history[n], upper_integration_limit);
136 :
137 : // nothing to integrate
138 4320963 : if (lower > upper)
139 770224 : continue;
140 :
141 : // now integrate from lower to upper
142 3550739 : Real f = 0.5 * (upper - lower);
143 17753695 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
144 : {
145 14202956 : Real recoil_energy = f * (_quad_points[qp] + 1) + lower;
146 14202956 : integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) *
147 14202956 : (nonCaptureProbability(i, k, energy, recoil_energy) *
148 14202956 : linearInterpolation(energy - recoil_energy, i, j, l) -
149 : current_value);
150 : }
151 : }
152 : return integral;
153 : }
154 :
155 : Real
156 135890 : PolyatomicDisplacementDerivativeFunction::source(Real energy,
157 : unsigned int i,
158 : unsigned int j,
159 : unsigned int l)
160 : {
161 : // compute the derivative of the net displacement function at E
162 135890 : unsigned int index = _net_displacement_function->energyIndex(energy);
163 : Real d_net_dE = 0;
164 : // if energy == Emin => index = 0, the net displacement rate is constant below that value
165 135890 : if (index != 0)
166 : {
167 135890 : Real e1 = _net_displacement_function->energyPoint(index - 1);
168 : Real e2 = _net_displacement_function->energyPoint(index);
169 135890 : Real v1 = _net_displacement_function->linearInterpolation(e1, i, j);
170 135890 : Real v2 = _net_displacement_function->linearInterpolation(e2, i, j);
171 135890 : d_net_dE = (v2 - v1) / (e2 - e1);
172 : }
173 135890 : return _net_displacement_function->integralTypeI(energy, i, j, l) +
174 135890 : _net_displacement_function->integralTypeII(energy, i, j, l) -
175 135890 : d_net_dE * stoppingPowerDerivative(i, l, energy);
176 : }
177 :
178 : #endif
|
