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 "PolyatomicDisplacementFunction.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 35 : PolyatomicDisplacementFunction::PolyatomicDisplacementFunction(
27 : std::vector<MyTRIM_NS::Element> polyatomic_material,
28 : nrt_type damage_function_type,
29 35 : std::vector<std::vector<Real>> Ecap)
30 : : PolyatomicDisplacementFunctionBase(polyatomic_material, damage_function_type, Ecap),
31 35 : _total_displacement_function(_damage_function_type == TOTAL)
32 : {
33 35 : if (damage_function_type != TOTAL && damage_function_type != NET)
34 0 : throw std::exception();
35 :
36 : // set the number of indices
37 35 : _n_indices = 2;
38 :
39 : Real Edisp_min = std::numeric_limits<Real>::max();
40 93 : for (unsigned int j = 0; j < _n_species; ++j)
41 58 : if (_material->_element[j]._Edisp < Edisp_min)
42 : Edisp_min = _material->_element[j]._Edisp;
43 35 : _energy_history[0] = Edisp_min;
44 :
45 : // set initial conditions for _displacement_function,
46 : // note: for net displacement function, the gii = 1, for total displacement function nij = 0
47 35 : if (_damage_function_type == NET)
48 81 : for (unsigned int i = 0; i < _n_species; ++i)
49 50 : _displacement_function[0][mapIndex(i, i, 0)] = 1;
50 35 : }
51 :
52 : std::vector<Real>
53 65529 : PolyatomicDisplacementFunction::getRHS(Real energy)
54 : {
55 65529 : std::vector<Real> f(_problem_size);
56 173739 : for (unsigned int i = 0; i < nSpecies(); ++i)
57 : {
58 108210 : Real stopping_power = stoppingPower(i, energy);
59 193572 : for (unsigned int j = 0; j < nSpecies(); ++j)
60 : {
61 : // working on the right hand side for nu_ij
62 193572 : unsigned int n = mapIndex(i, j, 0);
63 193572 : f[n] = 0;
64 :
65 557868 : for (unsigned int k = 0; k < nSpecies(); ++k)
66 364296 : f[n] += numberFraction(k) *
67 364296 : (integralTypeI(energy, i, j, k) + integralTypeII(energy, i, j, k)) / stopping_power;
68 : }
69 : }
70 65529 : return f;
71 0 : }
72 :
73 : Real
74 498368 : PolyatomicDisplacementFunction::integralTypeI(Real energy,
75 : unsigned int i,
76 : unsigned int j,
77 : unsigned int k) const
78 : {
79 498368 : Real upper_integration_limit = energy * _lambda[i][k];
80 : Real integral = 0;
81 498368 : Real delta_kj = k == j && _total_displacement_function ? 1 : 0;
82 498368 : Real Ebind = _material->_element[k]._Elbind;
83 :
84 : // the integration follows the already existing energies
85 11851192 : for (unsigned int l = 1; l < _energy_history.size(); ++l)
86 : {
87 : // adjust integration limits to be within E_{l-1} ... min(E_l, upper_integration_limit)
88 11352824 : Real lower = _energy_history[l - 1];
89 11352824 : Real upper = std::min(_energy_history[l], upper_integration_limit);
90 :
91 : // nothing to integrate
92 11352824 : if (lower > upper)
93 1989596 : continue;
94 :
95 : // now integrate from lower to upper
96 9363228 : Real f = 0.5 * (upper - lower);
97 46816140 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
98 : {
99 37452912 : Real recoil_energy = f * (_quad_points[qp] + 1) + lower;
100 37452912 : integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) *
101 37452912 : displacementProbability(k, recoil_energy) *
102 37452912 : (delta_kj + linearInterpolation(recoil_energy - Ebind, k, j));
103 : }
104 : }
105 498368 : return integral;
106 : }
107 :
108 : Real
109 498368 : PolyatomicDisplacementFunction::integralTypeII(Real energy,
110 : unsigned int i,
111 : unsigned int j,
112 : unsigned int k) const
113 : {
114 498368 : Real upper_integration_limit = energy * _lambda[i][k];
115 498368 : Real threshold = std::min(_asymptotic_threshold, energy * _lambda[i][k]);
116 :
117 : // store the current displacement function value
118 498368 : Real current_value = linearInterpolation(energy, i, j);
119 :
120 : // estimate the derivative d(nu_i) / dE:
121 498368 : Real dE = _energy_history.back() - _energy_history[_energy_history.size() - 2];
122 498368 : Real derivative = (current_value - linearInterpolation(energy - dE, i, j)) / dE;
123 :
124 : // integrate up to threshold and multiply by estimate of the derivative
125 498368 : Real integral = -weightedScatteringIntegral(energy, threshold, i, k) * derivative;
126 :
127 498368 : if (energy * _lambda[i][k] <= _asymptotic_threshold)
128 : return integral;
129 :
130 : // the integration follows the already existing energies
131 12349560 : for (unsigned int l = 0; l < _energy_history.size(); ++l)
132 : {
133 : // adjust integration limits to be within
134 : Real lower;
135 11851192 : if (l == 0)
136 498368 : lower = _asymptotic_threshold;
137 : else
138 11352824 : lower = std::max(_energy_history[l - 1], _asymptotic_threshold);
139 :
140 11851192 : Real upper = std::min(_energy_history[l], upper_integration_limit);
141 :
142 : // nothing to integrate
143 11851192 : if (lower > upper)
144 1989596 : continue;
145 :
146 : // now integrate from lower to upper
147 9861596 : Real f = 0.5 * (upper - lower);
148 49307980 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
149 : {
150 39446384 : Real recoil_energy = f * (_quad_points[qp] + 1) + lower;
151 39446384 : integral += f * _quad_weights[qp] * scatteringCrossSection(i, k, energy, recoil_energy) *
152 39446384 : (nonCaptureProbability(i, k, energy, recoil_energy) *
153 39446384 : linearInterpolation(energy - recoil_energy, i, j) -
154 : current_value);
155 : }
156 : }
157 : return integral;
158 : }
159 :
160 : #endif
|
