Skip to content
GitLab
Commit 1ff74dff authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Use a branch-free Wendland C2 kernel for the gravitational softening

parent 7eb23d60
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 1ff74dff
......@@ -91,6 +91,7 @@ theory/SPH/*.pdf
theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf
theory/Multipoles/potential.pdf
m4/libtool.m4
m4/ltoptions.m4
......
src/kernel_gravity.h
View file @ 1ff74dff
......@@ -28,107 +28,46 @@
#include "minmax.h"
#include "vector.h"
/* The gravity kernel is defined as a degree 6 polynomial in the distance
r. The resulting value should be post-multiplied with r^-3, resulting
in a polynomial with terms ranging from r^-3 to r^3, which are
sufficient to model both the direct potential as well as the splines
near the origin.
As in the hydro case, the 1/h^3 needs to be multiplied in afterwards */
/* Coefficients for the kernel. */
#define kernel_grav_name "Gadget-2 softening kernel"
#define kernel_grav_degree 6 /* Degree of the polynomial */
#define kernel_grav_ivals 2 /* Number of branches */
static const float
kernel_grav_coeffs[(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)]
__attribute__((aligned(16))) = {32.f,
-38.4f,
0.f,
10.66666667f,
0.f,
0.f,
0.f, /* 0 < u < 0.5 */
-10.66666667f,
38.4f,
-48.f,
21.3333333f,
0.f,
0.f,
-0.066666667f, /* 0.5 < u < 1 */
0.f,
0.f,
0.f,
0.f,
0.f,
0.f,
0.f}; /* 1 < u */
/**
* @brief Computes the gravity softening function.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_eval(
float u, float *const W) {
/* Pick the correct branch of the kernel */
const int ind = (int)min(u * (float)kernel_grav_ivals, kernel_grav_ivals);
const float *const coeffs =
&kernel_grav_coeffs[ind * (kernel_grav_degree + 1)];
/* First two terms of the polynomial ... */
float w = coeffs[0] * u + coeffs[1];
/* ... and the rest of them */
for (int k = 2; k <= kernel_grav_degree; k++) w = u * w + coeffs[k];
/* Return everything */
*W = w / (u * u * u);
/* W(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
*W = 21.f * u - 90.f;
*W = *W * u + 140.f;
*W = *W * u - 84.f;
*W = *W * u;
*W = *W * u + 14.f;
*W = *W * u;
}
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/**
* @brief Computes the gravity softening function in double precision.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
double u, double *const W) {
static const double kernel_grav_coeffs_double[(kernel_grav_degree + 1) *
(kernel_grav_ivals + 1)]
__attribute__((aligned(16))) = {32.,
-38.4,
0.,
10.66666667,
0.,
0.,
0., /* 0 < u < 0.5 */
-10.66666667,
38.4,
-48.,
21.3333333,
0.,
0.,
-0.066666667, /* 0.5 < u < 1 */
0.,
0.,
0.,
0.,
0.,
0.,
0.}; /* 1 < u */
/* Pick the correct branch of the kernel */
const int ind = (int)min(u * (double)kernel_grav_ivals, kernel_grav_ivals);
const double *const coeffs =
&kernel_grav_coeffs_double[ind * (kernel_grav_degree + 1)];
/* First two terms of the polynomial ... */
double w = coeffs[0] * u + coeffs[1];
/* ... and the rest of them */
for (int k = 2; k <= kernel_grav_degree; k++) w = u * w + coeffs[k];
/* Return everything */
*W = w / (u * u * u);
/* W(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
*W = 21. * u - 90.;
*W = *W * u + 140.;
*W = *W * u - 84.;
*W = *W * u;
*W = *W * u + 14.;
*W = *W * u;
}
#endif /* SWIFT_GRAVITY_FORCE_CHECKS */
......
theory/Multipoles/potential.py 0 → 100644
View file @ 1ff74dff
###############################################################################
# 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/>.
#
##############################################################################
import matplotlib
matplotlib.use("Agg")
from pylab import *
from scipy import integrate
import distinct_colours as colours
from scipy.optimize import curve_fit
from scipy.optimize import fsolve
from matplotlib.font_manager import FontProperties
import numpy
import math
params = {'axes.labelsize': 9,
'axes.titlesize': 10,
'font.size': 12,
'legend.fontsize': 10,
'xtick.labelsize': 9,
'ytick.labelsize': 9,
'text.usetex': True,
'figure.figsize' : (3.15,3.15),
'figure.subplot.left' : 0.115,
'figure.subplot.right' : 0.99 ,
'figure.subplot.bottom' : 0.08 ,
'figure.subplot.top' : 0.99 ,
'figure.subplot.wspace' : 0. ,
'figure.subplot.hspace' : 0. ,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
# Parameters
epsilon = 2.
epsilon_plummer = epsilon
r_min = 0.
r_max = 4
r_max_plot = 2.6
# Radius
r = linspace(r_min, r_max, 401)
r[0] += 1e-9
u = r / epsilon
# Newtonian solution
phi_newton = 1. / r
F_newton = 1. / r**2
# Softened potential
phi = np.zeros(np.size(r))
W = np.zeros(np.size(r))
F = np.zeros(np.size(r))
phi_plummer = np.zeros(np.size(r))
F_plummer = np.zeros(np.size(r))
for i in range(np.size(r)):
if r[i] > epsilon:
phi[i] = 1. / r[i]
W[i] = 0.
F[i] = 1. / r[i]**2
else:
phi[i] = (-1./epsilon) * (3.*u[i]**7 - 15.*u[i]**6 + 28.*u[i]**5 - 21.*u[i]**4 + 7.*u[i]**2 - 3.)
W[i] = (21. / 2.*math.pi) * (4.*u[i]**5 - 15.*u[i]**4 + 20.*u[i]**3 - 10.*u[i]**2 + 1.) / epsilon**6
F[i] = (1./epsilon**2) * (21.*u[i]**6 - 90*u[i]**5 + 140.*u[i]**4 - 84.*u[i]**3 + 14*u[i])
phi_plummer[i] = (1./epsilon_plummer) * (1 + (u[i])**2)**(-1./2.)
F_plummer[i] = (1./epsilon_plummer**3) * r[i] / (1 + (u[i])**2)**(3./2.)
figure()
# Potential
subplot(311)
plot(r, phi_newton, 'g--', lw=0.8, label="${\\rm Newtonian}$")
plot(r, phi_plummer, 'b-.', lw=0.8, label="${\\rm Plummer}$")
plot(r, phi, 'r-', lw=1, label="${\\rm SWIFT}$")
plot([epsilon, epsilon], [0, 10], 'k:', alpha=0.5, lw=0.5)
legend(loc="upper right", frameon=False, handletextpad=0.1, handlelength=3.2, fontsize=8)
ylim(0, 2.1)
ylabel("$\\phi(r)$", labelpad=0)
xlim(0,r_max_plot)
xticks([0., 0.5, 1., 1.5, 2., 2.5], ["", "", "", "", "", ""])
# Force
subplot(312)
plot(r, F_newton, 'g--', lw=0.8)
plot(r, F_plummer, 'b-.', lw=0.8)
plot(r, F, 'r-', lw=1)
plot([epsilon, epsilon], [0, 10], 'k:', alpha=0.5, lw=0.5)
xlim(0,r_max_plot)
xticks([0., 0.5, 1., 1.5, 2., 2.5], ["", "", "", "", "", ""])
ylim(0, 0.95)
ylabel("$|\\mathbf{\\nabla}\\phi(r)|$", labelpad=0)
# Density
subplot(313)
plot(r, W, 'r-', lw=1)
plot([epsilon, epsilon], [0, 10], 'k:', alpha=0.5, lw=0.5)
xlim(0,r_max_plot)
xticks([0., 0.5, 1., 1.5, 2., 2.5], ["$0$", "$0.5$", "$1$", "$1.5$", "$2$", "$2.5$"])
xlabel("$r$", labelpad=-1)
ylim(0., 0.58)
ylabel("$\\rho(r)$", labelpad=0)
savefig("potential.pdf")
#Construct potential
# phi = np.zeros(np.size(r))
# for i in range(np.size(r)):
# if r[i] > 2*epsilon:
# phi[i] = 1./ r[i]
# elif r[i] > epsilon:
# phi[i] = -(1./epsilon) * ((32./3.)*u[i]**2 - (48./3.)*u[i]**3 + (38.4/4.)*u[i]**4 - (32./15.)*u[i]**5 + (2./30.)*u[i]**(-1) - (9/5.))
# else:
# phi[i] = -(1./epsilon) * ((32./6.)*u[i]**2 - (38.4/4.)*u[i]**4 + (32./5.)*u[i]**4 - (7./5.))
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment