Peano 4
Loading...
Searching...
No Matches
kernel_hydro.h
Go to the documentation of this file.
1#pragma once
2#include "cmath"
3#include "tarch/tarch.h"
4
17
18
19namespace swift2 {
20 namespace kernels {
21 namespace legacy {
22 namespace kernelHydro {
23
24#if defined(QUARTIC_SPLINE_KERNEL)
25 /* ------------------------------------------------------------------------- */
26
27 /* Quartic spline M5 kernel -- SWIFT version */
28
29#define kernel_name "Quartic spline (M5)"
30/* Coefficients for the kernel. */
31#define kernel_degree 4 /* Degree of the polynomial */
32#define kernel_ivals 5 /* Number of branches */
33
34#if HYDRO_DIMENSION == 1
35#define kernel_gamma ((float)(1.936492))
36#define kernel_constant ((float)(3125. / 768.))
37#elif HYDRO_DIMENSION == 2
38#define kernel_gamma ((float)(1.977173))
39#define kernel_constant ((float)(46875. * M_1_PI / 2398.))
40#elif HYDRO_DIMENSION == 3
41#define kernel_gamma ((float)(2.018932))
42#define kernel_constant ((float)(15625. * M_1_PI / 512.))
43#else
44#pragma error "Hydro dimension undefined"
45#endif
46
47 static const double kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)] __attribute__((aligned(AlignmentOnHeap
48 )))
49 = {
50 6.0, 0.0, -2.4, 0.0, 0.368, /* 0 < u < 0.2 */
51 -4.0, 8.0, -4.8, 0.32, 0.352, /* 0.2 < u < 0.4 */
52 -4.0, 8.0, -4.8, 0.32, 0.352, /* 0.4 < u < 0.6 */
53 1.0, -4.0, 6.0, -4.0, 1.0, /* 0.6 < u < 0.8 */
54 1.0, -4.0, 6.0, -4.0, 1.0, /* 0.8 < u < 1 */
55 0.0, 0.0, 0.0, 0.0, 0.0}; /* 1 < u */
56
57/* ------------------------------------------------------------------------- */
58
59// @TODO implement additional kernels (M4)
60#else
61
62#error "An SPH projection kernel function must be chosen !"
63/* ------------------------------------------------------------------------- */
64#endif
65
66/* ------------------------------------------------------------------------- */
67
68/* Define generic values (i.e. common to all kernels) */
69
70/* First some powers of gamma = H/h */
71#define kernel_gamma_inv ((float)(1. / kernel_gamma))
72#define kernel_gamma2 ((float)(kernel_gamma * kernel_gamma))
73
74/* define gamma^d, gamma^(d+1), 1/gamma^d and 1/gamma^(d+1) */
75#if HYDRO_DIMENSION == 1
76#define kernel_gamma_dim ((float)(kernel_gamma))
77#define kernel_gamma_dim_plus_one ((float)(kernel_gamma * kernel_gamma))
78#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma)))
79#define kernel_gamma_inv_dim_plus_one ((float)(1. / (kernel_gamma * kernel_gamma)))
80#elif HYDRO_DIMENSION == 2
81#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma))
82#define kernel_gamma_dim_plus_one ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
83#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma * kernel_gamma)))
84#define kernel_gamma_inv_dim_plus_one ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
85#elif HYDRO_DIMENSION == 3
86#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
87#define kernel_gamma_dim_plus_one ((float)(kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma))
88#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
89#define kernel_gamma_inv_dim_plus_one ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma)))
90#endif
91
92/* The number of branches (floating point conversion) */
93#define kernel_ivals_f ((float)(kernel_ivals))
94
95/* Kernel self contribution (i.e. W(0,h)) */
96#define kernel_root ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_gamma_inv_dim)
97
98 /* ------------------------------------------------------------------------- */
99
111 static void kernel_eval(double u, double& W) InlineMethod {
112
113 /* Go to the range [0,1[ from [0,H[ */
114 const double x = u * kernel_gamma_inv;
115
116 /* Pick the correct branch of the kernel */
117 const int temp = (int)(x * kernel_ivals_f);
118 const int ind = temp > kernel_ivals ? kernel_ivals : temp;
119 const double* const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
120
121 /* First two terms of the polynomial ... */
122 double w = coeffs[0] * x + coeffs[1];
123
124 /* ... and the rest of them */
125 for (int k = 2; k <= kernel_degree; k++)
126 w = x * w + coeffs[k];
127
128 w = std::max(w, 0.0);
129
130 /* Return everything */
131 W = w * kernel_constant * kernel_gamma_inv_dim;
132 }
133
196 static void kernel_deval(double u, double& W, double& dW_dx) InlineMethod {
197 /* Go to the range [0,1[ from [0,H[ */
198 const double x = u * kernel_gamma_inv;
199
200 /* Pick the correct branch of the kernel */
201 const int temp = (int)(x * kernel_ivals_f);
202 const int ind = temp > kernel_ivals ? kernel_ivals : temp;
203 const double* const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
204
205 /* w term in Horner's scheme */
206 double w = coeffs[0];
207 for (int k = 1; k <= kernel_degree; k++) {
208 w = w * x + coeffs[k];
209 }
210
211 /* dw_dx term in Horner's scheme */
212
213 double dw_dx = coeffs[0] * kernel_degree;
214
215 /* ... and the rest of them */
216 for (int k = 1; k <= kernel_degree - 1; k++) {
217 dw_dx = dw_dx * x + coeffs[k] * (kernel_degree - k);
218 }
219
220 w = std::max(w, 0.0);
221 dw_dx = std::min(dw_dx, 0.0);
222
223 /* Return everything */
224 W = w * kernel_constant * kernel_gamma_inv_dim;
225 dW_dx = dw_dx * kernel_constant * kernel_gamma_inv_dim_plus_one;
226 }
227
228
229 } // namespace kernelHydro
230 } // namespace legacy
231 }
232} // namespace swift2
int __attribute__((optimize("O0"))) toolbox
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define AlignmentOnHeap
Definition LinuxAMD.h:23
#define kernel_gamma_inv
#define kernel_ivals_f
This file is part of the ExaHyPE project.
static void kernel_eval(double u, double &W) InlineMethod
Computes the kernel function.
static void kernel_deval(double u, double &W, double &dW_dx) InlineMethod
Computes the kernel function and its derivative.
This file is part of the SWIFT 2 project.
#define InlineMethod
This is the marker that is to be used after the argument list of a function declaration.
Definition tarch.h:58