Skip to content
GitLab
Commit 7b0c442e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Do hydrodynamics in 1D and 2D for the Gadget-SPH scheme.

parent 9bd7e668
Hide whitespace changes
Inline Side-by-side
src/Makefile.am
View file @ 7b0c442e
......@@ -56,7 +56,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
equation_of_state.h \
dimension.h equation_of_state.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
......
src/const.h
View file @ 7b0c442e
......@@ -36,6 +36,11 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
/* Dimensionality of the problem */
//#define HYDRO_DIMENSION_3D
#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
//#define HYDRO_GAMMA_4_3
......
src/dimension.h 0 → 100644
View file @ 7b0c442e
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DIMENSION_H
#define SWIFT_DIMENSION_H
/* Config parameters. */
#include "../config.h"
/* Local headers. */
#include "const.h"
#include "inline.h"
#include "vector.h"
/* First define some constants */
#if defined(HYDRO_DIMENSION_3D)
#define hydro_dimension 3.f
#define hydro_dimension_inv 0.3333333333f
#elif defined(HYDRO_DIMENSION_2D)
#define hydro_dimension 2.f
#define hydro_dimension_inv 0.5f
#elif defined(HYDRO_DIMENSION_1D)
#define hydro_dimension 1.f
#define hydro_dimension_inv 1.f
#else
#error "A problem dimensionality must be chosen in const.h !"
#endif
/**
* @brief Returns the argument to the power given by the dimension
*
* Computes \f$x^d\f$.
*/
__attribute__((always_inline)) INLINE static float pow_dimension(float x) {
#if defined(HYDRO_DIMENSION_3D)
return x * x * x;
#elif defined(HYDRO_DIMENSION_2D)
return x * x;
#elif defined(HYDRO_DIMENSION_1D)
return x;
#else
error("The dimension is not defined !");
return 0.f;
#endif
}
/**
* @brief Returns the argument to the power given by the dimension plus one
*
* Computes \f$x^{d+1}\f$.
*/
__attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
float x) {
#if defined(HYDRO_DIMENSION_3D)
const float x2 = x * x;
return x2 * x2;
#elif defined(HYDRO_DIMENSION_2D)
return x * x * x;
#elif defined(HYDRO_DIMENSION_1D)
return x * x;
#else
error("The dimension is not defined !");
return 0.f;
#endif
}
#ifdef WITH_VECTORIZATION
__attribute__((always_inline)) INLINE static vector pow_dimension_vec(
vector x) {
#if defined(HYDRO_DIMENSION_3D)
return (vector)(x.v * x.v * x.v);
#elif defined(HYDRO_DIMENSION_2D)
return (vector)(x.v * x.v);
#elif defined(HYDRO_DIMENSION_1D)
return x;
#else
error("The dimension is not defined !");
return vec_set(0.f);
#endif
}
__attribute__((always_inline)) INLINE static vector pow_dimension_plus_one_vec(
vector x) {
#if defined(HYDRO_DIMENSION_3D)
const vector x2 = x.v * x.v;
return (vector)(x2.v * x2.v);
#elif defined(HYDRO_DIMENSION_2D)
return (vector)(x.v * x.v * x.v);
#elif defined(HYDRO_DIMENSION_1D)
return (vector)(x.v * x.v);
#else
error("The dimension is not defined !");
return vec_set(0.f);
#endif
}
#endif
#endif /* SWIFT_DIMENSION_H */
src/drift.h
View file @ 7b0c442e
......@@ -25,6 +25,7 @@
/* Local headers. */
#include "const.h"
#include "debug.h"
#include "dimension.h"
#include "hydro.h"
#include "part.h"
......@@ -85,7 +86,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
p->h *= expf(w1);
/* Predict density */
const float w2 = -3.0f * w1;
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f)
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
else
......
src/hydro/Gadget2/hydro.h
View file @ 7b0c442e
......@@ -18,6 +18,7 @@
******************************************************************************/
#include "adiabatic_index.h"
#include "dimension.h"
#include "equation_of_state.h"
/**
......@@ -169,33 +170,33 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Some smoothing length multiples. */
const float h = p->h;
const float ih = 1.0f / h;
const float ih2 = ih * ih;
const float ih4 = ih2 * ih2;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
p->rho_dh -= 3.0f * p->mass * kernel_root;
p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= ih * ih2;
p->rho_dh *= ih4;
p->rho *= h_inv_dim;
p->rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= ih4 * irho;
p->density.rot_v[1] *= ih4 * irho;
p->density.rot_v[2] *= ih4 * irho;
p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
/* Finish calculation of the velocity divergence */
p->density.div_v *= ih4 * irho;
p->density.div_v *= h_inv_dim_plus_one * irho;
}
/**
......@@ -308,7 +309,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
p->force.h_dt *= p->h * 0.333333333f;
p->force.h_dt *= p->h * hydro_dimension_inv;
p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
......
src/hydro/Gadget2/hydro_iact.h
View file @ 7b0c442e
......@@ -58,7 +58,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
pi->rho_dh -= mj * (3.f * wi + ui * wi_dx);
pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
......@@ -71,7 +71,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pj->rho += mi * wj;
pj->rho_dh -= mi * (3.f * wj + uj * wj_dx);
pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
......@@ -191,7 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
......@@ -199,7 +199,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Compute density of pj. */
rhoj.v = mi.v * wj.v;
rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v;
......@@ -253,7 +253,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
pi->rho_dh -= mj * (3.f * wi + u * wi_dx);
pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
......@@ -356,7 +356,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
......@@ -402,17 +402,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hi2_inv = hi_inv * hi_inv;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
const float wi_dr = hi2_inv * hi2_inv * wi_dx;
const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
const float hj2_inv = hj_inv * hj_inv;
const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hj2_inv * hj2_inv * wj_dx;
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
......@@ -485,7 +485,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
vector hi2_inv, hj2_inv;
vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector piPOrho, pjPOrho, pirho, pjrho;
vector mi, mj;
......@@ -580,19 +580,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
hi2_inv.v = hi_inv.v * hi_inv.v;
hid_inv = pow_dimension_plus_one_vec(hi_inv); /* 1/h^(d+1) */
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
hj2_inv.v = hj_inv.v * hj_inv.v;
hjd_inv = pow_dimension_plus_one_vec(hj_inv); /* 1/h^(d+1) */
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
......@@ -677,17 +677,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hi2_inv = hi_inv * hi_inv;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
const float wi_dr = hi2_inv * hi2_inv * wi_dx;
const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
const float hj2_inv = hj_inv * hj_inv;
const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hj2_inv * hj2_inv * wj_dx;
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
......@@ -753,7 +753,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
vector hi2_inv, hj2_inv;
vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector piPOrho, pjPOrho, pirho, pjrho;
vector mj;
......@@ -845,19 +845,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
hi2_inv.v = hi_inv.v * hi_inv.v;
hid_inv = pow_dimension_plus_one_vec(hi_inv);
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
hj2_inv.v = hj_inv.v * hj_inv.v;
hjd_inv = pow_dimension_plus_one_vec(hj_inv);
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
......
src/hydro_properties.c
View file @ 7b0c442e
......@@ -27,6 +27,7 @@
/* Local headers. */
#include "adiabatic_index.h"
#include "common_io.h"
#include "dimension.h"
#include "error.h"
#include "hydro.h"
#include "kernel_hydro.h"
......@@ -39,8 +40,7 @@ void hydro_props_init(struct hydro_props *p,
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
const float eta3 = p->eta_neighbours * p->eta_neighbours * p->eta_neighbours;
p->target_neighbours = eta3 * kernel_norm;
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
/* Ghost stuff */
......@@ -51,21 +51,26 @@ void hydro_props_init(struct hydro_props *p,
p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
const float max_volume_change = parser_get_opt_param_float(
params, "SPH:max_volume_change", hydro_props_default_volume_change);
p->log_max_h_change = logf(powf(max_volume_change, 0.33333333333f));
p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
}
void hydro_props_print(const struct hydro_props *p) {
message("Adiabatic index gamma: %f.", hydro_gamma);
message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
(int)hydro_dimension);
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours,
p->eta_neighbours);
message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
message(
"Hydrodynamic integration: Max change of volume: %.2f "
"(max|dlog(h)/dt|=%f).",
powf(expf(p->log_max_h_change), 3.f), p->log_max_h_change);
powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
......@@ -76,6 +81,7 @@ void hydro_props_print(const struct hydro_props *p) {
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
writeAttribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
writeAttribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
writeAttribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
writeAttribute_s(h_grpsph, "Kernel function", kernel_name);
writeAttribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
......
src/kernel_hydro.h
View file @ 7b0c442e
......@@ -24,6 +24,7 @@
/* Includes. */
#include "const.h"
#include "dimension.h"
#include "error.h"
#include "inline.h"
#include "vector.h"
......@@ -35,8 +36,16 @@
#define kernel_name "Cubic spline (M4)"
#define kernel_degree 3 /* Degree of the polynomial */
#define kernel_ivals 2 /* Number of branches */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(1.825742))
#define kernel_constant ((float)(16. * M_1_PI))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(1.778002))
#define kernel_constant ((float)(80. * M_1_PI / 7.))
#elif defined(HYDRO_DIMENSION_1D)
#define kernel_gamma ((float)(1.732051))
#define kernel_constant ((float)(8. / 3.))
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {3.f, -3.f, 0.f, 0.5f, /* 0 < u < 0.5 */
-1.f, 3.f, -3.f, 1.f, /* 0.5 < u < 1 */
......@@ -47,10 +56,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Quartic spline (M5)"
#define kernel_degree 4
#define kernel_ivals 5
#define kernel_degree 4 /* Degree of the polynomial */
#define kernel_ivals 5 /* Number of branches */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.018932))
#define kernel_constant ((float)(15625. * M_1_PI / 512.))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(1.977173))
#define kernel_constant ((float)(46875. * M_1_PI / 2398.))
#elif defined(HYDRO_DIMENSION_1D)
#define kernel_gamma ((float)(1.936492))
#define kernel_constant ((float)(3125. / 768.))
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
6.f, 0.f, -2.4f, 0.f, 0.368f, /* 0 < u < 0.2 */
......@@ -65,10 +82,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Quintic spline (M6)"
#define kernel_degree 5
#define kernel_ivals 3
#define kernel_degree 5 /* Degree of the polynomial */
#define kernel_ivals 3 /* Number of branches */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.195775))
#define kernel_constant ((float)(2187. * M_1_PI / 40.))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(2.158131))
#define kernel_constant ((float)(15309. * M_1_PI / 478.))
#elif defined(HYDRO_DIMENSION_1D)
#define kernel_gamma ((float)(2.121321))
#define kernel_constant ((float)(243. / 40.))
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
-10.f, 10.f, 0.f,
......@@ -85,10 +110,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C2"
#define kernel_degree 5
#define kernel_ivals 1
#define kernel_degree 5 /* Degree of the polynomial */
#define kernel_ivals 1 /* Number of branches */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(1.936492))
#define kernel_constant ((float)(21. * M_1_PI / 2.))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(1.897367))
#define kernel_constant ((float)(7. * M_1_PI))
#elif defined(HYDRO_DIMENSION_1D)
#error "Wendland C2 kernel not defined in 1D."
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
4.f, -15.f, 20.f, -10.f, 0.f, 1.f, /* 0 < u < 1 */
......@@ -99,10 +131,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C4"
#define kernel_degree 8
#define kernel_degree 8 /* Degree of the polynomial */
#define kernel_ivals 1
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.207940))
#define kernel_constant ((float)(495. * M_1_PI / 32.))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(2.171239))
#define kernel_constant ((float)(9. * M_1_PI))
#elif defined(HYDRO_DIMENSION_1D)
#error "Wendland C4 kernel not defined in 1D."
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
11.666667f, -64.f, 140.f, -149.333333f, 70.f,
......@@ -115,10 +154,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C6"
#define kernel_degree 11
#define kernel_ivals 1
#define kernel_degree 11 /* Degree of the polynomial */
#define kernel_ivals 1 /* Number of branches */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.449490))
#define kernel_constant ((float)(1365. * M_1_PI / 64.))
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_gamma ((float)(2.415230))
#define kernel_constant ((float)(78. * M_1_PI / 7.))
#elif defined(HYDRO_DIMENSION_1D)
#error "Wendland C6 kernel not defined in 1D."
#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
32.f, -231.f, 704.f, -1155.f, 1056.f, -462.f,
......@@ -151,16 +197,32 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
#define kernel_ivals_f ((float)(kernel_ivals))
/* Kernel self contribution (i.e. W(0,h)) */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_root \
((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma3)
#elif defined(HYDRO_DIMENSION_2D)
#define kernel_root \
((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma2)
#elif defined(HYDRO_DIMENSION_1D)
#define kernel_root \
((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma)
#endif
/* Kernel normalisation constant (volume term) */
#if defined(HYDRO_DIMENSION_3D)
#define kernel_norm ((float)(4.0 * M_PI * kernel_gamma3 / 3.0))
#elif defined(HYDRO_DIMENSION_2D)