GADGET-4
kernel.h
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#ifndef KERNEL_H
13#define KERNEL_H
14
16{
17 double dpos[NUMDIMS];
18 double r;
19 double dv[NUMDIMS];
20 double wk, dwk;
21 double hinv, hinv3, hinv4;
22 double mj_wk, mj_dwk_r;
23};
24
26{
27 double dx, dy, dz;
28 double r, vsig, sound_i, sound_j;
29 double dvx, dvy, dvz, vdotr2;
30 double wk_i, wk_j, dwk_i, dwk_j;
32};
33
34#if !defined(CUBIC_SPLINE_KERNEL) && !defined(WENDLAND_C2_KERNEL) && !defined(WENDLAND_C4_KERNEL) && !defined(WENDLAND_C6_KERNEL)
35#define CUBIC_SPLINE_KERNEL /* fall back to cubic spline kernel */
36#endif
37
38/* fall back to three dimensions */
39#if !defined(TWODIMS) && !defined(ONEDIMS)
40#define THREEDIMS
41#endif
42
43/* Norms */
44#ifdef CUBIC_SPLINE_KERNEL
45
46#ifdef THREEDIMS
47#define NORM (8.0 / M_PI)
48#endif
49
50#ifdef TWODIMS
51#define NORM (40.0 / (7.0 * M_PI))
52#endif
53
54#ifdef ONEDIMS
55#define NORM (4.0 / 3.0)
56#endif
57
58#endif /* CUBIC_SPLINE_KERNEL */
59
60#ifdef WENDLAND_C2_KERNEL
61
62#ifdef THREEDIMS
63#define NORM (21.0 / (2.0 * M_PI))
64#endif
65
66#ifdef TWODIMS
67#define NORM (7.0 / M_PI)
68#endif
69
70#ifdef ONEDIMS
71#define NORM (5.0 / 4.0)
72#endif
73
74#endif /* WENDLAND_C2_KERNEL */
75
76#ifdef WENDLAND_C4_KERNEL
77
78#ifdef THREEDIMS
79#define NORM (495.0 / (32.0 * M_PI))
80#endif
81
82#ifdef TWODIMS
83#define NORM (9.0 / M_PI)
84#endif
85
86#ifdef ONEDIMS
87#define NORM (3.0 / 2.0)
88#endif
89
90#endif /* WENDLAND_C4_KERNEL */
91
92#ifdef WENDLAND_C6_KERNEL
93
94#ifdef THREEDIMS
95#define NORM (1365.0 / (64.0 * M_PI))
96#endif
97
98#ifdef TWODIMS
99#define NORM (78.0 / (7.0 * M_PI))
100#endif
101
102#ifdef ONEDIMS
103#define NORM (55.0 / 32.0)
104#endif
105
106#endif /* WENDLAND_C6_KERNEL */
107
108#define COMPUTE_WK -1
109#define COMPUTE_WK_AND_DWK 0
110#define COMPUTE_DWK 1
111
112static inline void kernel_hinv(double h, double *hinv, double *hinv3, double *hinv4)
113{
114 *hinv = 1.0 / h;
115
116#ifdef THREEDIMS
117 *hinv3 = *hinv * *hinv * *hinv;
118#endif
119
120#ifdef TWODIMS
121 *hinv3 = *hinv * *hinv;
122#endif
123
124#ifdef ONEDIMS
125 *hinv3 = *hinv;
126#endif
127
128 *hinv4 = *hinv3 * *hinv;
129}
130
131/* Attention: Here we assume that kernel is only called
132 with range 0..1 for u as done in hydra or density !!
133 Call with mode COMPUTE_WK_AND_DWK to calculate dwk and wk
134 Call with mode COMPUTE_WK to calculate only wk
135 Call with mode COMPUTE_DWK to calculate only dwk */
136static inline void kernel_main(double u, double hinv3, double hinv4, double *wk, double *dwk, int mode)
137{
138#ifdef CUBIC_SPLINE_KERNEL
139#if defined(WENDLAND_C2_KERNEL) || defined(WENDLAND_C4_KERNEL) || defined(WENDLAND_C6_KERNEL)
140#error "Only one SPH kernel can be used"
141#endif
142 if(u < 0.5)
143 {
144 if(mode >= COMPUTE_WK_AND_DWK)
145 *dwk = u * (18.0 * u - 12.0);
146 if(mode <= COMPUTE_WK_AND_DWK)
147 *wk = (1.0 + 6.0 * (u - 1.0) * u * u);
148 }
149 else
150 {
151 double t1 = (1.0 - u);
152 double t2 = t1 * t1;
153 if(mode >= COMPUTE_WK_AND_DWK)
154 *dwk = -6.0 * t2;
155 if(mode <= COMPUTE_WK_AND_DWK)
156 *wk = 2.0 * t2 * t1;
157 }
158#endif
159
160#ifdef WENDLAND_C2_KERNEL /* Dehnen & Aly 2012 */
161#ifdef ONEDIMS
162 double t1 = (1.0 - u);
163 double t2 = (t1 * t1);
164
165 if(mode >= COMPUTE_WK_AND_DWK)
166 *dwk = -12.0 * u * t2;
167 if(mode <= COMPUTE_WK_AND_DWK)
168 *wk = t2 * t1 * (1.0 + u * 3.0);
169
170#else /* 2d or 3d */
171 double t1 = (1.0 - u);
172 double t2 = (t1 * t1);
173 double t4 = t2 * t2;
174 if(mode >= COMPUTE_WK_AND_DWK)
175 *dwk = -20.0 * u * t2 * t1;
176 if(mode <= COMPUTE_WK_AND_DWK)
177 *wk = t4 * (1.0 + u * 4.0);
178
179#endif
180#endif /* WENDLAND_C2_KERNEL */
181
182#ifdef WENDLAND_C4_KERNEL /* Dehnen & Aly 2012 */
183#ifdef ONEDIMS
184 double t1 = (1.0 - u);
185 double t2 = t1 * t1;
186 double t4 = t2 * t2;
187 double t5 = t4 * t1;
188
189 if(mode >= COMPUTE_WK_AND_DWK)
190 *dwk = -14.0 * t4 * (4.0 * u + 1) * u;
191 if(mode <= COMPUTE_WK_AND_DWK)
192 *wk = t5 * (1.0 + u * (5.0 + 8.0 * u));
193
194#else /* 2d or 3d */
195 double t1 = (1.0 - u);
196 double t2 = (t1 * t1);
197 double t4 = t2 * t2;
198 double t6 = t2 * t2 * t2;
199 if(mode >= COMPUTE_WK_AND_DWK)
200 *dwk = -56.0 / 3.0 * u * t4 * t1 * (5.0 * u + 1);
201 if(mode <= COMPUTE_WK_AND_DWK)
202 *wk = t6 * (1.0 + u * (6.0 + 35.0 / 3.0 * u));
203
204#endif
205#endif /* WENDLAND_C4_KERNEL */
206
207#ifdef WENDLAND_C6_KERNEL /* Dehnen & Aly 2012 */
208#ifdef ONEDIMS
209 double t1 = (1.0 - u);
210 double t2 = (t1 * t1);
211 double t4 = t2 * t2;
212 double t6 = t4 * t2;
213 double t7 = t4 * t2 * t1;
214 if(mode >= COMPUTE_WK_AND_DWK)
215 *dwk = -6.0 * u * t6 * (3.0 + u * (18.0 + 35.0 * u));
216 if(mode <= COMPUTE_WK_AND_DWK)
217 *wk = t7 * (1.0 + u * (7.0 + u * (19.0 + 21.0 * u)));
218
219#else /* 2d or 3d */
220 double t1 = (1.0 - u);
221 double t2 = (t1 * t1);
222 double t4 = t2 * t2;
223 double t7 = t4 * t2 * t1;
224 double t8 = t4 * t4;
225 if(mode >= COMPUTE_WK_AND_DWK)
226 *dwk = -22.0 * u * (1.0 + u * (7.0 + 16.0 * u)) * t7;
227 if(mode <= COMPUTE_WK_AND_DWK)
228 *wk = t8 * (1.0 + u * (8.0 + u * (25.0 + 32.0 * u)));
229
230#endif
231#endif /* WENDLAND_C6_KERNEL */
232 if(mode >= COMPUTE_WK_AND_DWK)
233 *dwk *= NORM * hinv4;
234 if(mode <= COMPUTE_WK_AND_DWK)
235 *wk *= NORM * hinv3;
236}
237
238#if defined(WENDLAND_BIAS_CORRECTION) && (!(defined(WENDLAND_C2_KERNEL) || defined(WENDLAND_C4_KERNEL) || defined(WENDLAND_C6_KERNEL)))
239#error "WENDLAND_BIAS_CORRECTION only works with a Wendland kernel"
240#endif
241
242#if defined(WENDLAND_BIAS_CORRECTION) && (defined(WENDLAND_C2_KERNEL) || defined(WENDLAND_C4_KERNEL) || defined(WENDLAND_C6_KERNEL))
243
244#if defined(ONEDIMS) || defined(TWODIMS)
245#error "WENDLAND_BIAS_CORRECTION is only implemented for 3D"
246#endif
247
248static inline void get_bias_correction_parameters(double *alpha, double *eps100)
249{
250#ifdef WENDLAND_C2_KERNEL
251 *eps100 = 0.0294;
252 *alpha = 0.977;
253#endif
254#ifdef WENDLAND_C4_KERNEL
255 *eps100 = 0.01342;
256 *alpha = 1.579;
257#endif
258#ifdef WENDLAND_C6_KERNEL
259 *eps100 = 0.0116;
260 *alpha = 2.236;
261#endif
262}
263static inline double get_density_bias(double hsml, double mass, int DesNumNgb)
264{
265 kernel_density kernel;
266 kernel_hinv(hsml, &kernel.hinv, &kernel.hinv3, &kernel.hinv4);
267 kernel_main(0, kernel.hinv3, kernel.hinv4, &kernel.wk, &kernel.dwk, COMPUTE_WK);
268 double alpha = 0;
269 double eps100 = 0;
270 get_bias_correction_parameters(&alpha, &eps100);
271 double wc_correction = eps100 * pow(DesNumNgb * 0.01, -alpha) * mass * kernel.wk;
272 return wc_correction;
273}
274#endif
275
276#ifdef EXPLICIT_VECTORIZATION
277static inline void kernel_main_vector(Vec4d u, Vec4d hinv3, Vec4d hinv4, Vec4d *wk, Vec4d *dwk)
278{
279#ifdef CUBIC_SPLINE_KERNEL
280 Vec4d ucompl = u - 1.0;
281 Vec4db decision = (u < 0.5);
282 Vec4d wksub = 6.0 * u;
283
284 Vec4d ucompsq = ucompl * ucompl;
285 Vec4d uucomp = u * ucompl;
286 Vec4d wk1 = 1.0 + uucomp * wksub;
287 Vec4d dwk1 = wksub + uucomp * 18.0;
288 Vec4d wk2 = ucompsq * -2.0 * ucompl;
289 Vec4d dwk2 = ucompsq * -6.0;
290 *wk = select(decision, wk1, wk2);
291 *dwk = select(decision, dwk1, dwk2);
292#endif
293
294#ifdef WENDLAND_C2_KERNEL /* Dehnen & Aly 2012 */
295#ifdef ONEDIMS
296 Vec4d t1 = 1.0 - u;
297
298 Vec4d t2 = t1 * t1;
299
300 *dwk = -12.0 * u * t2;
301 *wk = t2 * t1 * (1.0 + u * 3.0);
302
303#else /* 2d or 3d */
304 Vec4d t1 = (1.0 - u);
305 Vec4d t2 = (t1 * t1);
306
307 *dwk = -20.0 * u * t2 * t1;
308 *wk = t2 * t2 * (1.0 + u * 4.0);
309#endif
310#endif /* WENDLAND_C2_KERNEL */
311
312#ifdef WENDLAND_C4_KERNEL /* Dehnen & Aly 2012 */
313#ifdef ONEDIMS
314 Vec4d t1 = (1.0 - u);
315 Vec4d t2 = t1 * t1;
316 Vec4d t4 = t2 * t2;
317
318 *dwk = -14.0 * t4 * (4.0 * u + 1) * u;
319 *wk = t1 * t4 * (1.0 + u * (5.0 + 8.0 * u));
320
321#else /* 2d or 3d */
322 Vec4d t1 = (1.0 - u);
323 Vec4d t2 = (t1 * t1);
324 Vec4d t4 = t2 * t2;
325 Vec4d t6 = t2 * t2 * t2;
326 *dwk = -56.0 / 3.0 * u * t4 * t1 * (5.0 * u + 1);
327
328 *wk = t6 * (1.0 + u * (6.0 + 35.0 / 3.0 * u));
329
330#endif
331#endif /* WENDLAND_C4_KERNEL */
332
333#ifdef WENDLAND_C6_KERNEL /* Dehnen & Aly 2012 */
334#ifdef ONEDIMS
335 Vec4d t1 = (1.0 - u);
336 Vec4d t2 = (t1 * t1);
337 Vec4d t4 = t2 * t2;
338 Vec4d t6 = t4 * t2;
339 Vec4d t7 = t4 * t2 * t1;
340 *dwk = -6.0 * u * t6 * (3.0 + u * (18.0 + 35.0 * u));
341 *wk = t7 * (1.0 + u * (7.0 + u * (19.0 + 21.0 * u)));
342
343#else /* 2d or 3d */
344 Vec4d t1 = (1.0 - u);
345 Vec4d t2 = (t1 * t1);
346 Vec4d t4 = t2 * t2;
347 Vec4d t7 = t4 * t2 * t1;
348 Vec4d t8 = t4 * t4;
349 *dwk = -22.0 * u * (1.0 + u * (7.0 + 16.0 * u)) * t7;
350 *wk = t8 * (1.0 + u * (8.0 + u * (25.0 + 32.0 * u)));
351
352#endif
353#endif /* WENDLAND_C6_KERNEL */
354 *dwk *= hinv4;
355 *wk *= hinv3;
356}
357#endif /* EXPLICIT_VECTORIZATION */
358#endif
#define NUMDIMS
Definition: constants.h:369
#define COMPUTE_WK_AND_DWK
Definition: kernel.h:109
#define COMPUTE_WK
Definition: kernel.h:108
#define NORM
Definition: kernel.h:47
expr pow(half base, half exp)
Definition: half.hpp:2803
double mj_wk
Definition: kernel.h:22
double dpos[NUMDIMS]
Definition: kernel.h:17
double wk
Definition: kernel.h:20
double hinv
Definition: kernel.h:21
double dv[NUMDIMS]
Definition: kernel.h:19
double r
Definition: kernel.h:18
double mj_dwk_r
Definition: kernel.h:22
double hinv4
Definition: kernel.h:21
double hinv3
Definition: kernel.h:21
double dwk
Definition: kernel.h:20
double h_i
Definition: kernel.h:31
double dx
Definition: kernel.h:27
double dvx
Definition: kernel.h:29
double dwk_ij
Definition: kernel.h:31
double dwk_j
Definition: kernel.h:30
double vdotr2
Definition: kernel.h:29
double vsig
Definition: kernel.h:28
double dz
Definition: kernel.h:27
double sound_j
Definition: kernel.h:28
double dvy
Definition: kernel.h:29
double h_j
Definition: kernel.h:31
double r
Definition: kernel.h:28
double dy
Definition: kernel.h:27
double wk_i
Definition: kernel.h:30
double wk_j
Definition: kernel.h:30
double sound_i
Definition: kernel.h:28
double dvz
Definition: kernel.h:29
double rho_ij_inv
Definition: kernel.h:31
double dwk_i
Definition: kernel.h:30