An open source method of characteristics neutron transport code.
ExpEvaluator.h
Go to the documentation of this file.
1 
8 #ifndef EXPEVALUATOR_H_
9 #define EXPEVALUATOR_H_
10 
11 #ifdef __cplusplus
12 #define _USE_MATH_DEFINES
13 #include "log.h"
14 #include "Quadrature.h"
15 #include <malloc.h>
16 #include <math.h>
17 #include "exponentials.h"
18 #endif
19 
20 
29 class ExpEvaluator {
30 
31 private:
32 
34  bool _interpolate;
35 
37  bool _linear_source;
38 
40  FP_PRECISION _exp_table_spacing;
41 
43  FP_PRECISION _inverse_exp_table_spacing;
44 
46  FP_PRECISION _sin_theta_no_offset;
47 
49  FP_PRECISION _inverse_sin_theta_no_offset;
50 
52  int _table_size;
53 
55  FP_PRECISION* _exp_table;
56 
58  Quadrature* _quadrature;
59 
61  FP_PRECISION _max_optical_length;
62 
64  FP_PRECISION _exp_precision;
65 
67  int _azim_index;
68 
70  int _polar_index;
71 
73  int _num_exp_terms;
74 
76  int _num_polar_terms;
77 
78 
79 public:
80 
81  ExpEvaluator();
82  virtual ~ExpEvaluator();
83 
84  void setQuadrature(Quadrature* quadrature);
85  void setMaxOpticalLength(FP_PRECISION max_optical_length);
86  void setExpPrecision(FP_PRECISION exp_precision);
87  void useInterpolation();
88  void useIntrinsic();
89  void useLinearSource();
90 
91  FP_PRECISION getMaxOpticalLength();
92  FP_PRECISION getExpPrecision();
93  bool isUsingInterpolation();
94  FP_PRECISION getTableSpacing();
95  int getTableSize();
96  FP_PRECISION* getExpTable();
97  int getExponentialIndex(FP_PRECISION tau);
98  FP_PRECISION getDifference(int index, FP_PRECISION tau);
99  FP_PRECISION convertDistance3Dto2D(FP_PRECISION length);
100  FP_PRECISION getInverseSinTheta();
101 
102  void initialize(int azim_index, int polar_index, bool solve_3D);
103  FP_PRECISION computeExponential(FP_PRECISION tau, int polar_offset);
104  FP_PRECISION computeExponentialF1(int index, int polar_offset,
105  FP_PRECISION dt, FP_PRECISION dt2);
106  FP_PRECISION computeExponentialF2(int index, int polar_offset,
107  FP_PRECISION dt, FP_PRECISION dt2);
108  FP_PRECISION computeExponentialH(int index, int polar_offset,
109  FP_PRECISION dt, FP_PRECISION dt2);
110  void retrieveExponentialComponents(FP_PRECISION tau, int polar_offset,
111  FP_PRECISION* exp_F1,
112  FP_PRECISION* exp_F2,
113  FP_PRECISION* exp_H);
114 
115  FP_PRECISION computeExponentialG2(FP_PRECISION tau);
117 };
118 
119 
126 inline int ExpEvaluator::getExponentialIndex(FP_PRECISION tau) {
127  return int(tau * _inverse_exp_table_spacing);
128 }
129 
130 
138 inline FP_PRECISION ExpEvaluator::getDifference(int index, FP_PRECISION tau) {
139  return tau - index * _exp_table_spacing;
140 }
141 
142 
148 inline FP_PRECISION ExpEvaluator::convertDistance3Dto2D(FP_PRECISION length) {
149  return length * _sin_theta_no_offset;
150 }
151 
152 
157 inline FP_PRECISION ExpEvaluator::getInverseSinTheta() {
158  return _inverse_sin_theta_no_offset;
159 }
160 
161 
168 inline FP_PRECISION ExpEvaluator::computeExponential(FP_PRECISION tau,
169  int polar_offset) {
170 
171 #ifndef THREED
172  FP_PRECISION inv_sin_theta = 1.f / _quadrature->getSinThetaInline(_azim_index,
173  _polar_index + polar_offset);
174 #else
175  FP_PRECISION inv_sin_theta = _inverse_sin_theta_no_offset;
176 #endif
177  FP_PRECISION exp_F1;
178  expF1_fractional(tau * inv_sin_theta, &exp_F1);
179 
180  return inv_sin_theta * exp_F1;
181 }
182 
183 
203 inline FP_PRECISION ExpEvaluator::computeExponentialF1(int index,
204  int polar_offset,
205  FP_PRECISION dt,
206  FP_PRECISION dt2) {
207 
208  /* Calculate full index */
209  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;
210 
211  if (_interpolate) {
212  return _exp_table[full_index] + _exp_table[full_index + 1] * dt +
213  _exp_table[full_index + 2] * dt2;
214  }
215  else {
216  int polar_index = _polar_index + polar_offset;
217  FP_PRECISION tau = index * _exp_table_spacing + dt;
218  FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
219  polar_index);
220  return (1.0 - exp(- tau * inv_sin_theta)) / tau;
221  }
222 }
223 
224 
244 inline FP_PRECISION ExpEvaluator::computeExponentialF2(int index,
245  int polar_offset,
246  FP_PRECISION dt,
247  FP_PRECISION dt2) {
248  /* Calculate full index */
249  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;
250 
251  if (_interpolate)
252  return _exp_table[full_index + 3] + _exp_table[full_index + 4] * dt +
253  _exp_table[full_index + 5] * dt2;
254  else {
255 
256  int polar_index = _polar_index + polar_offset;
257  FP_PRECISION tau = index * _exp_table_spacing + dt;
258  FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
259  polar_index);
260  FP_PRECISION tau_m = tau * inv_sin_theta;
261  FP_PRECISION F1 = (1.0 - exp(- tau_m)) / tau;
262  return 2.0 / tau * (inv_sin_theta - F1) - inv_sin_theta * F1;
263  }
264 }
265 
266 
286 inline FP_PRECISION ExpEvaluator::computeExponentialH(int index,
287  int polar_offset,
288  FP_PRECISION dt,
289  FP_PRECISION dt2) {
290  /* Calculate full index */
291  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;
292 
293  if (_interpolate)
294  return _exp_table[full_index + 6] + _exp_table[full_index + 7] * dt +
295  _exp_table[full_index + 8] * dt2;
296  else {
297  int polar_index = _polar_index + polar_offset;
298  FP_PRECISION tau = index * _exp_table_spacing + dt;
299  FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
300  polar_index);
301  FP_PRECISION tau_m = tau * inv_sin_theta;
302  FP_PRECISION F1 = (1.0 - exp(- tau_m)) / tau;
303  FP_PRECISION G1 = 1.0 / tau + 0.5 * inv_sin_theta - (1.0 + 1.0 / tau_m) * F1;
304  return 0.5 * inv_sin_theta - G1;
305  }
306 }
307 
308 
322 inline FP_PRECISION ExpEvaluator::computeExponentialG2(FP_PRECISION tau) {
323 
324  FP_PRECISION exp_G2;
325  expG2_fractional(tau, &exp_G2);
326 
327  return exp_G2;
328 }
329 
330 
347 inline void ExpEvaluator::retrieveExponentialComponents(FP_PRECISION tau,
348  int polar_offset,
349 #ifdef SWIG //FIXME Find out how to use restrict with SWIG
350  FP_PRECISION* exp_F1,
351  FP_PRECISION* exp_F2,
352  FP_PRECISION* exp_H) {
353 #else
354  FP_PRECISION* __restrict__ exp_F1,
355  FP_PRECISION* __restrict__ exp_F2,
356  FP_PRECISION* __restrict__ exp_H) {
357 #endif
358 
359 #ifndef THREED
360  FP_PRECISION inv_sin_theta = 1.f / _quadrature->getSinThetaInline(_azim_index,
361  _polar_index + polar_offset);
362 #else
363  FP_PRECISION inv_sin_theta = _inverse_sin_theta_no_offset;
364 #endif
365 
366  /* Limit range of tau to avoid numerical errors */
367  tau = std::max(FP_PRECISION(1e-8), tau * inv_sin_theta);
368 
369  /* Compute exponentials from a common exponential */
370  FP_PRECISION exp_G;
371  expG_fractional(tau, &exp_G);
372  *exp_F1 = 1.f - tau*exp_G;
373  *exp_F1 *= inv_sin_theta;
374 
375  exp_G *= inv_sin_theta;
376  *exp_F2 = 2.f*exp_G - *exp_F1;
377  *exp_F2 *= inv_sin_theta;
378 
379  *exp_H = *exp_F1 - exp_G;
380 
381  /* Quadratic exponential interpolation tables */
382 // __builtin_assume_aligned(_exp_table, VEC_ALIGNMENT);
383 //
384 // tau /= inv_sin_theta;
385 // int exp_index = getExponentialIndex(tau);
386 // FP_PRECISION dt = getDifference(exp_index, tau);
387 // FP_PRECISION dt2 = dt * dt;
388 // int full_index = (exp_index * _num_polar_terms + polar_offset)
389 // * _num_exp_terms;
390 // *exp_F1 = _exp_table[full_index] + _exp_table[full_index + 1] * dt +
391 // _exp_table[full_index + 2] * dt2;
392 // *exp_F2 = _exp_table[full_index + 3] + _exp_table[full_index + 4] * dt +
393 // _exp_table[full_index + 5] * dt2;
394 // *exp_H = _exp_table[full_index + 6] + _exp_table[full_index + 7] * dt +
395 // _exp_table[full_index + 8] * dt2;
396 }
397 
398 
399 #endif /* EXPEVALUATOR_H_ */
void expF1_fractional(FP_PRECISION x, FP_PRECISION *expF1)
Computes the F1 exponential term.
Definition: exponentials.h:156
double getSinTheta(size_t azim, size_t polar) const
Returns the value for a particular polar angle.
Definition: Quadrature.cpp:41
void setExpPrecision(FP_PRECISION exp_precision)
Sets the maximum acceptable approximation error for exponentials.
Definition: ExpEvaluator.cpp:68
int getTableSize()
Get the number of entries in the exponential interpolation table.
Definition: ExpEvaluator.cpp:149
void setMaxOpticalLength(FP_PRECISION max_optical_length)
Sets the maximum optical length covered in the exponential interpolation table.
Definition: ExpEvaluator.cpp:50
double getSinThetaInline(size_t azim, size_t polar) const
Returns the value for a particular polar angle.
Definition: Quadrature.h:296
FP_PRECISION getTableSpacing()
Returns the exponential table spacing.
Definition: ExpEvaluator.cpp:135
void retrieveExponentialComponents(FP_PRECISION tau, int polar_offset, FP_PRECISION *exp_F1, FP_PRECISION *exp_F2, FP_PRECISION *exp_H)
Computes the F1, F2, H exponential term.
Definition: ExpEvaluator.h:347
void initialize(int azim_index, int polar_index, bool solve_3D)
If using linear interpolation, builds the table for each polar angle.
Definition: ExpEvaluator.cpp:179
Math approximations for computing exponentials.
void useInterpolation()
Use linear interpolation to compute exponentials.
Definition: ExpEvaluator.cpp:81
FP_PRECISION computeExponentialF1(int index, int polar_offset, FP_PRECISION dt, FP_PRECISION dt2)
Computes the F1 exponential term.
Definition: ExpEvaluator.h:203
The arbitrary quadrature parent class.
Definition: Quadrature.h:76
FP_PRECISION getExpPrecision()
Gets the maximum acceptable approximation error for exponentials.
Definition: ExpEvaluator.cpp:117
int getExponentialIndex(FP_PRECISION tau)
Get the index on the exponential interpolation grid of the value right beneath tau.
Definition: ExpEvaluator.h:126
ExpEvaluator()
Constructor initializes array pointers to NULL.
Definition: ExpEvaluator.cpp:9
ExpEvaluator * deepCopy()
Deep copies an ExpEvaluator, for developing purposes.
Definition: ExpEvaluator.cpp:351
virtual ~ExpEvaluator()
Destructor deletes table for linear interpolation of exponentials.
Definition: ExpEvaluator.cpp:30
FP_PRECISION * getExpTable()
Returns a pointer to the exponential interpolation table.
Definition: ExpEvaluator.cpp:163
This is a class for evaluating exponentials.
Definition: ExpEvaluator.h:29
void expG2_fractional(FP_PRECISION x, FP_PRECISION *expG2)
Computes the G2 exponential term.
Definition: exponentials.h:293
FP_PRECISION computeExponential(FP_PRECISION tau, int polar_offset)
Computes the F1 exponential term.
Definition: ExpEvaluator.h:168
void useIntrinsic()
Use the exponential intrinsic exp(...) to compute exponentials.
Definition: ExpEvaluator.cpp:89
The Quadrature abstract class and subclasses.
FP_PRECISION convertDistance3Dto2D(FP_PRECISION length)
Convert a 3D distance to a 2D based on the evaluator&#39;s polar angle.
Definition: ExpEvaluator.h:148
FP_PRECISION getDifference(int index, FP_PRECISION tau)
Compute the difference between an optical path and an indexed value in the exponential interpolation ...
Definition: ExpEvaluator.h:138
FP_PRECISION computeExponentialF2(int index, int polar_offset, FP_PRECISION dt, FP_PRECISION dt2)
Computes the F2 exponential term.
Definition: ExpEvaluator.h:244
FP_PRECISION computeExponentialH(int index, int polar_offset, FP_PRECISION dt, FP_PRECISION dt2)
Computes the H exponential term.
Definition: ExpEvaluator.h:286
void useLinearSource()
Use linear source exponentials.
Definition: ExpEvaluator.cpp:97
FP_PRECISION getMaxOpticalLength()
Gets the maximum optical length covered with the exponential interpolation table. ...
Definition: ExpEvaluator.cpp:108
void expG_fractional(FP_PRECISION x, FP_PRECISION *expG)
Computes an exponential G term, used to reconstruct F1, F2 and H.
Definition: exponentials.h:110
FP_PRECISION computeExponentialG2(FP_PRECISION tau)
Computes the G2 exponential term for a optical length and polar angle.
Definition: ExpEvaluator.h:322
Utility functions for writing log messages to the screen.
void setQuadrature(Quadrature *quadrature)
Set the Quadrature to use when computing exponentials.
Definition: ExpEvaluator.cpp:40
FP_PRECISION getInverseSinTheta()
Get the inverse of sin theta from the ExpEvaluator.
Definition: ExpEvaluator.h:157
bool isUsingInterpolation()
Returns true if using linear interpolation to compute exponentials.
Definition: ExpEvaluator.cpp:126