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 "PolyatomicDisplacementFunctionBase.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 :
21 : // gsl includes
22 : #include "gsl/gsl_sf_legendre.h"
23 : #include "gsl/gsl_integration.h"
24 :
25 67 : PolyatomicDisplacementFunctionBase::PolyatomicDisplacementFunctionBase(
26 : std::vector<MyTRIM_NS::Element> polyatomic_material,
27 : nrt_type damage_function_type,
28 67 : std::vector<std::vector<Real>> Ecap)
29 67 : : _damage_function_type(damage_function_type),
30 67 : _quad_order(4),
31 67 : _n_species(polyatomic_material.size()),
32 67 : _problem_size(_damage_function_type == ENERGY ? _n_species
33 63 : : _damage_function_type == NET_DERIVATIVE ? _n_species * _n_species * _n_species
34 : : _n_species * _n_species),
35 67 : _energy_history({0.0}),
36 268 : _displacement_function({std::vector<Real>(_problem_size)})
37 : {
38 134 : _simconf = std::make_unique<MyTRIM_NS::SimconfType>(1);
39 134 : _material = std::make_unique<MyTRIM_NS::MaterialBase>(_simconf.get(), 1);
40 171 : for (auto & elem : polyatomic_material)
41 : {
42 104 : _material->_element.push_back(elem);
43 208 : _ions.push_back(std::make_unique<MyTRIM_NS::IonBase>(elem._Z, elem._m, 0));
44 : }
45 67 : _material->prepare();
46 :
47 : // compute _lambda
48 67 : _lambda.resize(_n_species);
49 171 : for (unsigned int i = 0; i < _n_species; ++i)
50 : {
51 104 : _lambda[i].resize(_n_species);
52 282 : for (unsigned int j = 0; j < _n_species; ++j)
53 : {
54 178 : Real Ai = _material->_element[i]._m;
55 178 : Real Aj = _material->_element[j]._m;
56 178 : _lambda[i][j] = 4.0 * Ai * Aj / (Ai + Aj) / (Ai + Aj);
57 : }
58 : }
59 :
60 67 : if (Ecap.size() == 1 && Ecap[0].size() == 0)
61 : {
62 : // default is Ecap_ij = Ed_j
63 51 : _Ecap.resize(_n_species);
64 123 : for (unsigned int i = 0; i < _n_species; ++i)
65 : {
66 72 : _Ecap[i].resize(_n_species);
67 186 : for (unsigned int j = 0; j < _n_species; ++j)
68 114 : _Ecap[i][j] = _material->_element[j]._Edisp;
69 : }
70 : }
71 : else
72 16 : _Ecap = Ecap;
73 :
74 : // set up the gsl ODE machinery
75 : auto func = &PolyatomicDisplacementFunction::gslInterface;
76 67 : _sys = {func, NULL, _problem_size, this};
77 67 : _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
78 :
79 : // set up integration rule
80 67 : gslQuadRule(_quad_order, _quad_points, _quad_weights);
81 201 : }
82 :
83 67 : PolyatomicDisplacementFunctionBase::~PolyatomicDisplacementFunctionBase()
84 : {
85 67 : gsl_odeiv2_driver_free(_ode_driver);
86 67 : }
87 :
88 : int
89 113653 : PolyatomicDisplacementFunctionBase::gslInterface(Real energy,
90 : const Real disp[],
91 : Real f[],
92 : void * params)
93 : {
94 : (void)disp;
95 : PolyatomicDisplacementFunctionBase * padf = (PolyatomicDisplacementFunctionBase *)params;
96 113653 : std::vector<Real> rhs = padf->getRHS(energy);
97 483771 : for (unsigned int j = 0; j < rhs.size(); ++j)
98 370118 : f[j] = rhs[j];
99 113653 : return GSL_SUCCESS;
100 : }
101 :
102 : void
103 2423 : PolyatomicDisplacementFunctionBase::advanceDisplacements(Real new_energy)
104 : {
105 2423 : if (_energy_history.back() > new_energy)
106 0 : return;
107 :
108 : // add the new energy point to energy history and store old_energy
109 2423 : Real old_energy = _energy_history.back();
110 2423 : _energy_history.push_back(new_energy);
111 :
112 : // add new energy step to _displacement_function; initial guess is the old energy
113 2423 : _displacement_function.push_back(_displacement_function.back());
114 :
115 2423 : int status = gsl_odeiv2_driver_apply(
116 : _ode_driver, &old_energy, new_energy, _displacement_function.back().data());
117 :
118 2423 : if (status != GSL_SUCCESS)
119 0 : std ::cout << "gsl_odeiv2_driver_apply returned error code " << status << std::endl;
120 : }
121 :
122 : void
123 5 : PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral()
124 : {
125 : // clear the _displacement_function_integral array because this might
126 : // have been called before; not the most efficient if the user keeps on calling
127 : // this function but they are at their own peril and should know to call it
128 : // after finishing the computation of disp function
129 5 : _displacement_function_integral.clear();
130 :
131 : // resize the array; note that initial value is zero (entry for
132 : // _displacement_function_integral[0])
133 5 : _displacement_function_integral.resize(nEnergySteps(), std::vector<Real>(_problem_size));
134 :
135 256 : for (unsigned int e = 1; e < nEnergySteps(); ++e)
136 : {
137 : // initialize the integral to the integral of the previous step
138 251 : _displacement_function_integral[e] = _displacement_function_integral[e - 1];
139 :
140 : // set energy points for integration
141 : Real lower = energyPoint(e - 1);
142 : Real upper = energyPoint(e);
143 251 : Real f = 0.5 * (upper - lower);
144 1255 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
145 : {
146 1004 : Real energy = f * (_quad_points[qp] + 1) + lower;
147 1004 : Real w = f * _quad_weights[qp];
148 :
149 5020 : for (unsigned int n = 0; n < _problem_size; ++n)
150 4016 : _displacement_function_integral[e][n] += w * linearInterpolationFromFlatIndex(energy, n);
151 : }
152 : }
153 5 : }
154 :
155 : void
156 218403 : PolyatomicDisplacementFunctionBase::gslQuadRule(unsigned int quad_order,
157 : std::vector<Real> & quad_points,
158 : std::vector<Real> & quad_weights)
159 : {
160 : // set up integration rule
161 218403 : auto * qp_table = gsl_integration_glfixed_table_alloc(quad_order);
162 218403 : quad_points.resize(quad_order);
163 218403 : quad_weights.resize(quad_order);
164 22052271 : for (std::size_t j = 0; j < quad_order; ++j)
165 : {
166 : Real point, weight;
167 21833868 : gsl_integration_glfixed_point(-1.0, 1.0, j, &point, &weight, qp_table);
168 21833868 : quad_points[j] = point;
169 21833868 : quad_weights[j] = weight;
170 : }
171 218403 : gsl_integration_glfixed_table_free(qp_table);
172 218403 : }
173 :
174 : Real
175 189656 : PolyatomicDisplacementFunctionBase::stoppingPower(unsigned int species, Real energy)
176 : {
177 : assert(species < _n_species);
178 189656 : _ions[species]->_E = energy;
179 : // need to divide by atomic density because we need Se for PC equations
180 : // units are eV * Angstrom**2
181 189656 : return _material->getrstop(_ions[species].get()) / _material->_arho;
182 : }
183 :
184 : Real
185 135890 : PolyatomicDisplacementFunctionBase::stoppingPowerDerivative(unsigned int species,
186 : unsigned int changing_species,
187 : Real energy)
188 : {
189 135890 : _ions[species]->_E = energy;
190 135890 : return _material->getDrstopDcomp(_ions[species].get(), _material->_element[changing_species]);
191 : }
192 :
193 : /**
194 : * The scattering cross section is computed using Lindhard's universal formula
195 : * with f(xi) evaluated using the Thomas Fermi potential. Units are Angstrom**2 / eV
196 : */
197 : Real
198 152233696 : PolyatomicDisplacementFunctionBase::scatteringCrossSection(unsigned int i,
199 : unsigned int j,
200 : Real energy,
201 : Real recoil_energy) const
202 : {
203 : assert(i < _n_species);
204 : assert(j < _n_species);
205 :
206 : // get the current A & Z for projectile: i and target: j
207 152233696 : Real Ai = _material->_element[i]._m;
208 152233696 : Real Aj = _material->_element[j]._m;
209 152233696 : Real Zi = _material->_element[i]._Z;
210 152233696 : Real Zj = _material->_element[j]._Z;
211 :
212 : // Z only appears as 1/3 power, so it's better to apply the 1/3 power here to get a sqrt
213 152233696 : Real Z = std::sqrt(std::pow(Zi, 2.0 / 3.0) + std::pow(Zj, 2.0 / 3.0));
214 152233696 : Real a = 0.8853 * _abohr / Z;
215 :
216 : // compute El and Tm
217 152233696 : Real Tm = _lambda[i][j] * energy;
218 152233696 : Real El = 30.7514664 * Zi * Zj * Z * (Ai + Aj) / Aj;
219 :
220 : // compute sqrt(t) and then the cross section using Lindhard's formula
221 152233696 : Real sqrt_t = std::sqrt((energy * energy / El / El) * recoil_energy / Tm);
222 152233696 : return 0.5 * universalF(sqrt_t) / recoil_energy / sqrt_t * M_PI * a * a;
223 : }
224 :
225 : Real
226 152233696 : PolyatomicDisplacementFunctionBase::universalF(Real xi) const
227 : {
228 : Real lp = 1.309;
229 152233696 : return lp * std::pow(xi, 1.0 / 3.0) *
230 152233696 : std::pow(1.0 + std::pow(2.0 * lp * std::pow(xi, 4.0 / 3.0), 2.0 / 3.0), -1.5);
231 : }
232 :
233 : unsigned int
234 8916 : PolyatomicDisplacementFunctionBase::findSpeciesIndex(int atomic_number, Real mass_number) const
235 : {
236 9822 : for (unsigned int j = 0; j < _n_species; ++j)
237 9822 : if (atomic_number == _material->_element[j]._Z &&
238 8916 : std::abs(mass_number - _material->_element[j]._m) < 1.0e-10)
239 8916 : return j;
240 : return libMesh::invalid_uint;
241 : }
242 :
243 : Real
244 105367424 : PolyatomicDisplacementFunctionBase::displacementProbability(unsigned int k, Real energy) const
245 : {
246 105367424 : return energy < _material->_element[k]._Edisp ? 0 : 1;
247 : }
248 :
249 : Real
250 54220888 : PolyatomicDisplacementFunctionBase::nonCaptureProbability(unsigned int i,
251 : unsigned int k,
252 : Real energy,
253 : Real recoil_energy) const
254 : {
255 : // Parkin & Coulter, JNM 101, (1981)
256 : return 1 -
257 54220888 : displacementProbability(k, recoil_energy) * (energy - recoil_energy < _Ecap[i][k] ? 1 : 0);
258 : }
259 :
260 : unsigned int
261 148318925 : PolyatomicDisplacementFunctionBase::energyIndex(Real energy) const
262 : {
263 148318925 : return energy >= _energy_history.back()
264 : ? nEnergySteps() - 1
265 148273255 : : std::distance(
266 : _energy_history.begin(),
267 148318925 : std::upper_bound(_energy_history.begin(), _energy_history.end(), energy));
268 : }
269 :
270 : /**
271 : * weightedScatteringIntegral computes the integral:
272 : * int_0^{energy_limit} d(sigma_ij) / dT T dT
273 : */
274 : Real
275 908364 : PolyatomicDisplacementFunctionBase::weightedScatteringIntegral(Real energy,
276 : Real energy_limit,
277 : unsigned int i,
278 : unsigned int j) const
279 : {
280 : // uses the t -> 0 limit to WSS form of Lindhard's universal cross section with Thomas-Fermi
281 : // potential for recoil_energy < _asymptotic_threshold
282 908460 : Real limit = std::min(_asymptotic_threshold, energy_limit);
283 908364 : Real Mi = _material->_element[i]._m;
284 908364 : Real Mj = _material->_element[j]._m;
285 908364 : Real Zi = _material->_element[i]._Z;
286 908364 : Real Zj = _material->_element[j]._Z;
287 908364 : Real Z = std::pow(std::pow(Zi, 2.0 / 3.0) + std::pow(Zj, 2.0 / 3.0), 1.5);
288 908364 : Real a = _abohr * 0.8853 * std::pow(Z, -1.0 / 3.0);
289 908364 : Real alpha_ij = 0.5 * 1.309 * M_PI * a * a * std::pow(Mi / Mj, 1.0 / 3.0) *
290 908364 : std::pow(Zi * Zj * 61.5038 * std::pow(Z, 1.0 / 3.0), 2.0 / 3.0);
291 908364 : Real integral = 1.5 * std::pow(limit, 2.0 / 3.0) * alpha_ij * std::pow(energy, -1.0 / 3.0);
292 :
293 908364 : if (energy_limit <= _asymptotic_threshold)
294 : return integral;
295 :
296 : // numerical integration from _asymptotic_threshold to energy_limit
297 26784 : Real spacing = std::max(energy_limit / _asymptotic_threshold /
298 13392 : std::ceil(energy_limit / _asymptotic_threshold / 1.05),
299 13392 : 1.02);
300 :
301 13392 : std::vector<Real> energies = {_asymptotic_threshold};
302 13392 : Real current_energy = _asymptotic_threshold;
303 : for (;;)
304 : {
305 1505088 : current_energy *= spacing;
306 1505088 : if (current_energy >= energy_limit)
307 : {
308 13392 : energies.push_back(energy_limit);
309 : break;
310 : }
311 1491696 : energies.push_back(current_energy);
312 : }
313 :
314 1518480 : for (unsigned int l = 1; l < energies.size(); ++l)
315 : {
316 1505088 : Real lower = energies[l - 1];
317 1505088 : Real upper = energies[l];
318 1505088 : Real scale = 0.5 * (upper - lower);
319 7525440 : for (unsigned int qp = 0; qp < _quad_order; ++qp)
320 : {
321 6020352 : Real recoil_energy = scale * (_quad_points[qp] + 1) + lower;
322 6020352 : integral += scale * _quad_weights[qp] * scatteringCrossSection(i, j, energy, recoil_energy) *
323 : recoil_energy;
324 : }
325 : }
326 : return integral;
327 : }
328 :
329 : Real
330 148149258 : PolyatomicDisplacementFunctionBase::linearInterpolation(Real energy,
331 : unsigned int i,
332 : unsigned int j,
333 : unsigned int l) const
334 : {
335 148149258 : unsigned int energy_index = energyIndex(energy);
336 148149258 : unsigned int index = mapIndex(i, j, l);
337 148149258 : if (energy_index == 0)
338 4556883 : return _displacement_function[0][index];
339 :
340 143592375 : return linearInterpolationHelper(energy, energy_index, index);
341 : }
342 :
343 : Real
344 4016 : PolyatomicDisplacementFunctionBase::linearInterpolationFromFlatIndex(Real energy,
345 : unsigned int index) const
346 : {
347 4016 : unsigned int energy_index = energyIndex(energy);
348 4016 : if (energy_index == 0)
349 0 : return _displacement_function[0][index];
350 :
351 4016 : return linearInterpolationHelper(energy, energy_index, index);
352 : }
353 :
354 : Real
355 143596391 : PolyatomicDisplacementFunctionBase::linearInterpolationHelper(Real energy,
356 : unsigned int energy_index,
357 : unsigned int index) const
358 : {
359 : // linear interpolation
360 143596391 : Real e1 = _energy_history[energy_index - 1];
361 143596391 : Real e2 = _energy_history[energy_index];
362 143596391 : Real v1 = _displacement_function[energy_index - 1][index];
363 143596391 : Real v2 = _displacement_function[energy_index][index];
364 :
365 143596391 : return v1 + (energy - e1) / (e2 - e1) * (v2 - v1);
366 : }
367 :
368 : Real
369 29761 : PolyatomicDisplacementFunctionBase::linearInterpolationIntegralDamageFunction(Real energy,
370 : unsigned int i,
371 : unsigned int j,
372 : unsigned int l) const
373 : {
374 29761 : unsigned int energy_index = energyIndex(energy);
375 29761 : unsigned int index = mapIndex(i, j, l);
376 :
377 29761 : if (energy_index == 0)
378 16112 : return _displacement_function_integral[0][index];
379 :
380 : // linear interpolation
381 13649 : Real e1 = _energy_history[energy_index - 1];
382 13649 : Real e2 = _energy_history[energy_index];
383 13649 : Real v1 = _displacement_function_integral[energy_index - 1][index];
384 13649 : Real v2 = _displacement_function_integral[energy_index][index];
385 :
386 13649 : return v1 + (energy - e1) / (e2 - e1) * (v2 - v1);
387 : }
388 :
389 : unsigned int
390 148514246 : PolyatomicDisplacementFunctionBase::mapIndex(unsigned int i, unsigned int j, unsigned int l) const
391 : {
392 : assert(_n_indices > 0);
393 : assert(j == 0 || _n_indices > 1);
394 : assert(l == 0 || _n_indices > 2);
395 148514246 : return i + j * _n_species + l * _n_species * _n_species;
396 : };
397 :
398 : #endif
|
