Peano
Loading...
Searching...
No Matches
adiabatic_index.h
Go to the documentation of this file.
1/*******************************************************************************
2 * This file is part of SWIFT.
3 * Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl).
4 * Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 ******************************************************************************/
20#ifndef SWIFT_ADIABATIC_INDEX_H
21#define SWIFT_ADIABATIC_INDEX_H
22
29/* Config parameters. */
30#include <config.h>
31
32/* Some standard headers. */
33#include <math.h>
34
35/* Local headers. */
36#include "cbrt.h"
37#include "error.h"
38#include "inline.h"
39
40/* First define some constants */
41#if defined(HYDRO_GAMMA_5_3)
42
43#define hydro_gamma 1.66666666666666667f
44#define hydro_gamma_minus_one 0.66666666666666667f
45#define hydro_gamma_plus_one 2.66666666666666667f
46#define hydro_one_over_gamma_minus_one 1.5f
47#define hydro_gamma_plus_one_over_two_gamma 0.8f
48#define hydro_gamma_minus_one_over_two_gamma 0.2f
49#define hydro_gamma_minus_one_over_gamma_plus_one 0.25f
50#define hydro_two_over_gamma_plus_one 0.75f
51#define hydro_two_over_gamma_minus_one 3.f
52#define hydro_gamma_minus_one_over_two 0.33333333333333333f
53#define hydro_two_gamma_over_gamma_minus_one 5.f
54#define hydro_one_over_gamma 0.6f
55
56#elif defined(HYDRO_GAMMA_7_5)
57
58#define hydro_gamma 1.4f
59#define hydro_gamma_minus_one 0.4f
60#define hydro_gamma_plus_one 2.4f
61#define hydro_one_over_gamma_minus_one 2.5f
62#define hydro_gamma_plus_one_over_two_gamma 0.857142857f
63#define hydro_gamma_minus_one_over_two_gamma 0.142857143f
64#define hydro_gamma_minus_one_over_gamma_plus_one 0.166666667f
65#define hydro_two_over_gamma_plus_one 0.833333333
66#define hydro_two_over_gamma_minus_one 5.f
67#define hydro_gamma_minus_one_over_two 0.2f
68#define hydro_two_gamma_over_gamma_minus_one 7.f
69#define hydro_one_over_gamma 0.714285714f
70
71#elif defined(HYDRO_GAMMA_4_3)
72
73#define hydro_gamma 1.33333333333333333f
74#define hydro_gamma_minus_one 0.33333333333333333f
75#define hydro_gamma_plus_one 2.33333333333333333f
76#define hydro_one_over_gamma_minus_one 3.f
77#define hydro_gamma_plus_one_over_two_gamma 0.875f
78#define hydro_gamma_minus_one_over_two_gamma 0.125f
79#define hydro_gamma_minus_one_over_gamma_plus_one 0.142857143f
80#define hydro_two_over_gamma_plus_one 0.857142857f
81#define hydro_two_over_gamma_minus_one 6.f
82#define hydro_gamma_minus_one_over_two 0.166666666666666666f
83#define hydro_two_gamma_over_gamma_minus_one 8.f
84#define hydro_one_over_gamma 0.75f
85
86#elif defined(HYDRO_GAMMA_2_1)
87
88#define hydro_gamma 2.f
89#define hydro_gamma_minus_one 1.f
90#define hydro_gamma_plus_one 3.f
91#define hydro_one_over_gamma_minus_one 1.f
92#define hydro_gamma_plus_one_over_two_gamma 0.75f
93#define hydro_gamma_minus_one_over_two_gamma 0.25f
94#define hydro_gamma_minus_one_over_gamma_plus_one 0.33333333333333333f
95#define hydro_two_over_gamma_plus_one 0.66666666666666666f
96#define hydro_two_over_gamma_minus_one 2.f
97#define hydro_gamma_minus_one_over_two 0.5f
98#define hydro_two_gamma_over_gamma_minus_one 4.f
99#define hydro_one_over_gamma 0.5f
100
101#else
102
103#error "An adiabatic index needs to be chosen in const.h !"
104
105#endif
106
112__attribute__((always_inline, const)) INLINE static float pow_gamma(float x) {
113
114#if defined(HYDRO_GAMMA_5_3)
115
116#ifdef WITH_ICBRTF
117 const float icbrt = icbrtf(x); /* x^(-1/3) */
118 return icbrt * x * x; /* x^(5/3) */
119#else
120 const float cbrt = cbrtf(x); /* x^(1/3) */
121 return cbrt * cbrt * x; /* x^(5/3) */
122#endif // WITH_ICBRTF
123
124#elif defined(HYDRO_GAMMA_7_5)
125
126 return powf(x, 1.4f); /* x^(7/5) */
127
128#elif defined(HYDRO_GAMMA_4_3)
129
130#ifdef WITH_ICBRTF
131 const float icbrt = icbrtf(x); /* x^(-1/3) */
132 return icbrt * icbrt * x * x; /* x^(4/3) */
133#else
134 return cbrtf(x) * x; /* x^(4/3) */
135#endif // WITH_ICBRTF
136
137#elif defined(HYDRO_GAMMA_2_1)
138
139 return x * x;
140
141#else
142
143 error("The adiabatic index is not defined !");
144 return 0.f;
145
146#endif
147}
148
155__attribute__((always_inline, const)) INLINE static float pow_gamma_minus_one(
156 float x) {
157
158#if defined(HYDRO_GAMMA_5_3)
159
160#ifdef WITH_ICBRTF
161 const float icbrt = icbrtf(x); /* x^(-1/3) */
162 return x * icbrt; /* x^(2/3) */
163#else
164 const float cbrt = cbrtf(x); /* x^(1/3) */
165 return cbrt * cbrt; /* x^(2/3) */
166#endif // WITH_ICBRTF
167
168#elif defined(HYDRO_GAMMA_7_5)
169
170 return powf(x, 0.4f); /* x^(2/5) */
171
172#elif defined(HYDRO_GAMMA_4_3)
173
174#ifdef WITH_ICBRTF
175 const float icbrt = icbrtf(x); /* x^(-1/3) */
176 return x * icbrt * icbrt; /* x^(1/3) */
177#else
178 return cbrtf(x); /* x^(1/3) */
179#endif // WITH_ICBRTF
180
181#elif defined(HYDRO_GAMMA_2_1)
182
183 return x;
184
185#else
186
187 error("The adiabatic index is not defined !");
188 return 0.f;
189
190#endif
191}
192
199__attribute__((always_inline, const)) INLINE static float
200pow_minus_gamma_minus_one(float x) {
201
202#if defined(HYDRO_GAMMA_5_3)
203
204#ifdef WITH_ICBRTF
205 const float icbrt = icbrtf(x); /* x^(-1/3) */
206 return icbrt * icbrt; /* x^(-2/3) */
207#else
208 const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
209 return cbrt_inv * cbrt_inv; /* x^(-2/3) */
210#endif // WITH_ICBRTF
211
212#elif defined(HYDRO_GAMMA_7_5)
213
214 return powf(x, -0.4f); /* x^(-2/5) */
215
216#elif defined(HYDRO_GAMMA_4_3)
217
218#ifdef WITH_ICBRTF
219 return icbrtf(x); /* x^(-1/3) */
220#else
221 return 1.f / cbrtf(x); /* x^(-1/3) */
222#endif // WITH_ICBRTF
223
224#elif defined(HYDRO_GAMMA_2_1)
225
226 return 1.f / x;
227
228#else
229
230 error("The adiabatic index is not defined !");
231 return 0.f;
232
233#endif
234}
235
245__attribute__((always_inline, const)) INLINE static float pow_minus_gamma(
246 float x) {
247
248#if defined(HYDRO_GAMMA_5_3)
249
250#ifdef WITH_ICBRTF
251 const float icbrt = icbrtf(x); /* x^(-1/3) */
252 const float icbrt2 = icbrt * icbrt; /* x^(-2/3) */
253 return icbrt * icbrt2 * icbrt2; /* x^(-5/3) */
254#else
255 const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
256 const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
257 return cbrt_inv * cbrt_inv2 * cbrt_inv2; /* x^(-5/3) */
258#endif // WITH_ICBRTF
259
260#elif defined(HYDRO_GAMMA_7_5)
261
262 return powf(x, -1.4f); /* x^(-7/5) */
263
264#elif defined(HYDRO_GAMMA_4_3)
265
266#ifdef WITH_ICBRTF
267 const float cbrt_inv = icbrtf(x); /* x^(-1/3) */
268#else
269 const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
270#endif // WITH_ICBRTF
271 const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
272 return cbrt_inv2 * cbrt_inv2; /* x^(-4/3) */
273
274#elif defined(HYDRO_GAMMA_2_1)
275
276 const float inv = 1.f / x;
277 return inv * inv;
278
279#else
280
281 error("The adiabatic index is not defined !");
282 return 0.f;
283
284#endif
285}
286
296__attribute__((always_inline, const)) INLINE static float
297pow_two_over_gamma_minus_one(float x) {
298
299#if defined(HYDRO_GAMMA_5_3)
300
301 return x * x * x; /* x^3 */
302
303#elif defined(HYDRO_GAMMA_7_5)
304
305 const float x2 = x * x;
306 const float x3 = x2 * x;
307 return x2 * x3;
308
309#elif defined(HYDRO_GAMMA_4_3)
310
311 const float x3 = x * x * x; /* x^3 */
312 return x3 * x3; /* x^6 */
313
314#elif defined(HYDRO_GAMMA_2_1)
315
316 return x * x; /* x^2 */
317
318#else
319
320 error("The adiabatic index is not defined !");
321 return 0.f;
322
323#endif
324}
325
336__attribute__((always_inline, const)) INLINE static float
337pow_two_gamma_over_gamma_minus_one(float x) {
338
339#if defined(HYDRO_GAMMA_5_3)
340
341 const float x2 = x * x;
342 const float x3 = x2 * x;
343 return x2 * x3;
344
345#elif defined(HYDRO_GAMMA_7_5)
346
347 const float x2 = x * x;
348 const float x4 = x2 * x2;
349 return x4 * x2 * x;
350
351#elif defined(HYDRO_GAMMA_4_3)
352
353 const float x2 = x * x;
354 const float x4 = x2 * x2;
355 return x4 * x4; /* x^8 */
356
357#elif defined(HYDRO_GAMMA_2_1)
358
359 const float x2 = x * x;
360 return x2 * x2; /* x^4 */
361
362#else
363
364 error("The adiabatic index is not defined !");
365 return 0.f;
366
367#endif
368}
369
380__attribute__((always_inline, const)) INLINE static float
381pow_gamma_minus_one_over_two_gamma(float x) {
382
383#if defined(HYDRO_GAMMA_5_3)
384
385 return powf(x, 0.2f); /* x^0.2 */
386
387#elif defined(HYDRO_GAMMA_7_5)
388
389 return powf(x, hydro_gamma_minus_one_over_two_gamma);
390
391#elif defined(HYDRO_GAMMA_4_3)
392
393 return powf(x, 0.125f); /* x^0.125 */
394
395#elif defined(HYDRO_GAMMA_2_1)
396
397 return powf(x, 0.25f); /* x^0.25 */
398
399#else
400
401 error("The adiabatic index is not defined !");
402 return 0.f;
403
404#endif
405}
406
417__attribute__((always_inline, const)) INLINE static float
418pow_minus_gamma_plus_one_over_two_gamma(float x) {
419
420#if defined(HYDRO_GAMMA_5_3)
421
422 return powf(x, -0.8f); /* x^-0.8 */
423
424#elif defined(HYDRO_GAMMA_7_5)
425
426 return powf(x, -hydro_gamma_plus_one_over_two_gamma);
427
428#elif defined(HYDRO_GAMMA_4_3)
429
430 return powf(x, -0.875f); /* x^-0.875 */
431
432#elif defined(HYDRO_GAMMA_2_1)
433
434 return powf(x, -0.75f); /* x^-0.75 */
435
436#else
437
438 error("The adiabatic index is not defined !");
439 return 0.f;
440
441#endif
442}
443
452__attribute__((always_inline, const)) INLINE static float pow_one_over_gamma(
453 float x) {
454
455#if defined(HYDRO_GAMMA_5_3)
456
457 return powf(x, hydro_one_over_gamma); /* x^(3/5) */
458
459#elif defined(HYDRO_GAMMA_7_5)
460
461 return powf(x, hydro_one_over_gamma);
462
463#elif defined(HYDRO_GAMMA_4_3)
464
465 return powf(x, hydro_one_over_gamma); /* x^(3/4) */
466
467#elif defined(HYDRO_GAMMA_2_1)
468
469 return sqrtf(x); /* x^(1/2) */
470
471#else
472
473 error("The adiabatic index is not defined !");
474 return 0.f;
475
476#endif
477}
478
486__attribute__((always_inline, const)) INLINE static float
487pow_three_gamma_minus_two(float x) {
488
489#if defined(HYDRO_GAMMA_5_3)
490
491 return x * x * x; /* x^(3) */
492
493#elif defined(HYDRO_GAMMA_7_5)
494
495 return powf(x, 2.2f); /* x^(11/5) */
496
497#elif defined(HYDRO_GAMMA_4_3)
498
499 return x * x; /* x^(2) */
500
501#elif defined(HYDRO_GAMMA_2_1)
502
503 return x * x * x * x; /* x^(4) */
504
505#else
506
507 error("The adiabatic index is not defined !");
508 return 0.f;
509
510#endif
511}
512
521__attribute__((always_inline, const)) INLINE static float
522pow_three_gamma_minus_five_over_two(float x) {
523
524#if defined(HYDRO_GAMMA_5_3)
525
526 return 1.f; /* x^(0) */
527
528#elif defined(HYDRO_GAMMA_7_5)
529
530 return powf(x, -0.4f); /* x^(-2/5) */
531
532#elif defined(HYDRO_GAMMA_4_3)
533
534 return 1.f / sqrtf(x); /* x^(-1/2) */
535
536#elif defined(HYDRO_GAMMA_2_1)
537
538 return sqrtf(x); /* x^(1/2) */
539
540#else
541
542 error("The adiabatic index is not defined !");
543 return 0.f;
544
545#endif
546}
547
555__attribute__((always_inline, const)) INLINE static float
556pow_three_gamma_minus_one(float x) {
557
558#if defined(HYDRO_GAMMA_5_3)
559
560 return x * x; /* x^(2) */
561
562#elif defined(HYDRO_GAMMA_7_5)
563
564 return powf(x, 1.2f); /* x^(6/5) */
565
566#elif defined(HYDRO_GAMMA_4_3)
567
568 return x; /* x^(1) */
569
570#elif defined(HYDRO_GAMMA_2_1)
571
572 return x * x * x; /* x^(3) */
573
574#else
575
576 error("The adiabatic index is not defined !");
577 return 0.f;
578
579#endif
580}
581
582#endif /* SWIFT_ADIABATIC_INDEX_H */
__attribute__
__attribute__((always_inline, const)) INLINE static float pow_gamma(float x)
Returns the argument to the power given by the adiabatic index.
Definition adiabatic_index.h:112
INLINE
#define INLINE
Defines inline.
Definition inline.h:36