diff --git a/configure.ac b/configure.ac
index 7c93d6fcd6794bee108c10bf725ca2e084e0718f..6c0a06005778c499c53ae8e596a0685a78886542 100644
--- a/configure.ac
+++ b/configure.ac
@@ -799,6 +799,15 @@ case "$with_potential" in
;;
esac
+# Gravity multipole order
+AC_ARG_WITH([multipole-order],
+ [AS_HELP_STRING([--with-multipole-order=<order>],
+ [order of the multipole and gravitational field expansion @<:@ default: 3@:>@]
+ )],
+ [with_multipole_order="$withval"],
+ [with_multipole_order="3"]
+)
+AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Multipole order])
# Check for git, needed for revision stamps.
@@ -846,6 +855,7 @@ AC_MSG_RESULT([
Riemann solver : $with_riemann
Cooling function : $with_cooling
External potential : $with_potential
+ Multipole order : $with_multipole_order
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index bb5f97f029e1d50d81bbdccae9ac620e9e0e6f08..67ce7c530263f7905a0d5147832dfc39148d753d 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -13,6 +13,9 @@ TimeIntegration:
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+Scheduler:
+ cell_split_size: 50
+
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
@@ -23,6 +26,13 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ epsilon: 0.0001 # Softening length (in internal units).
+ a_smooth: 1000.
+ r_cut: 4.
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 12a413b7e2c45443601c0b9753383b90942298b0..c755768bcfafebf3efe6307080e9e85d3a0a4bf5 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -23,6 +23,13 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ epsilon: 0.0001 # Softening length (in internal units).
+ a_smooth: 1000.
+ r_cut: 4.
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/UniformDMBox/plot_gravity_checks.py b/examples/UniformDMBox/plot_gravity_checks.py
new file mode 100644
index 0000000000000000000000000000000000000000..5efd5847ca9749fffaee48e586c0a1976fbac9d5
--- /dev/null
+++ b/examples/UniformDMBox/plot_gravity_checks.py
@@ -0,0 +1,219 @@
+#!/usr/bin/env python
+
+import sys
+import glob
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+params = {'axes.labelsize': 14,
+'axes.titlesize': 18,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 14,
+'ytick.labelsize': 14,
+'text.usetex': True,
+'figure.figsize': (10, 10),
+'figure.subplot.left' : 0.06,
+'figure.subplot.right' : 0.99 ,
+'figure.subplot.bottom' : 0.06 ,
+'figure.subplot.top' : 0.985 ,
+'figure.subplot.wspace' : 0.14 ,
+'figure.subplot.hspace' : 0.14 ,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+plt.rcParams.update(params)
+plt.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+min_error = 1e-6
+max_error = 1e-1
+num_bins = 51
+
+# Construct the bins
+bin_edges = np.linspace(np.log10(min_error), np.log10(max_error), num_bins + 1)
+bin_size = (np.log10(max_error) - np.log10(min_error)) / num_bins
+bins = 0.5*(bin_edges[1:] + bin_edges[:-1])
+bin_edges = 10**bin_edges
+bins = 10**bins
+
+# Colours
+cols = ['b', 'g', 'r', 'm']
+
+# Time-step to plot
+step = int(sys.argv[1])
+
+# Find the files for the different expansion orders
+order_list = glob.glob("gravity_checks_step%d_order*.dat"%step)
+num_order = len(order_list)
+
+# Get the multipole orders
+order = np.zeros(num_order)
+for i in range(num_order):
+ order[i] = int(order_list[i][26])
+
+# Start the plot
+plt.figure()
+
+# Get the Gadget-2 data if existing
+gadget2_file_list = glob.glob("forcetest_gadget2.txt")
+if len(gadget2_file_list) != 0:
+
+ gadget2_data = np.loadtxt(gadget2_file_list[0])
+ gadget2_ids = gadget2_data[:,0]
+ gadget2_pos = gadget2_data[:,1:4]
+ gadget2_a_exact = gadget2_data[:,4:7]
+ gadget2_a_grav = gadget2_data[:, 7:10]
+
+ # Sort stuff
+ sort_index = np.argsort(gadget2_ids)
+ gadget2_ids = gadget2_ids[sort_index]
+ gadget2_pos = gadget2_pos[sort_index, :]
+ gadget2_a_exact = gadget2_a_exact[sort_index, :]
+ gadget2_a_grav = gadget2_a_grav[sort_index, :]
+
+ # Compute the error norm
+ diff = gadget2_a_exact - gadget2_a_grav
+
+ norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
+ norm_a = np.sqrt(gadget2_a_exact[:,0]**2 + gadget2_a_exact[:,1]**2 + gadget2_a_exact[:,2]**2)
+
+ norm_error = norm_diff / norm_a
+ error_x = abs(diff[:,0]) / norm_a
+ error_y = abs(diff[:,1]) / norm_a
+ error_z = abs(diff[:,2]) / norm_a
+
+ # Bin the error
+ norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+
+ norm_median = np.median(norm_error)
+ median_x = np.median(error_x)
+ median_y = np.median(error_y)
+ median_z = np.median(error_z)
+
+ norm_per95 = np.percentile(norm_error,95)
+ per95_x = np.percentile(error_x,95)
+ per95_y = np.percentile(error_y,95)
+ per95_z = np.percentile(error_z,95)
+
+ plt.subplot(221)
+ plt.semilogx(bins, norm_error_hist, 'k--', label="Gadget-2")
+ plt.plot([norm_median, norm_median], [2.7, 3], 'k-', lw=1)
+ plt.plot([norm_per95, norm_per95], [2.7, 3], 'k:', lw=1)
+ plt.subplot(222)
+ plt.semilogx(bins, error_x_hist, 'k--', label="Gadget-2")
+ plt.plot([median_x, median_x], [1.8, 2], 'k-', lw=1)
+ plt.plot([per95_x, per95_x], [1.8, 2], 'k:', lw=1)
+ plt.subplot(223)
+ plt.semilogx(bins, error_y_hist, 'k--', label="Gadget-2")
+ plt.plot([median_y, median_y], [1.8, 2], 'k-', lw=1)
+ plt.plot([per95_y, per95_y], [1.8, 2], 'k:', lw=1)
+ plt.subplot(224)
+ plt.semilogx(bins, error_z_hist, 'k--', label="Gadget-2")
+ plt.plot([median_z, median_z], [1.8, 2], 'k-', lw=1)
+ plt.plot([per95_z, per95_z], [1.8, 2], 'k:', lw=1)
+
+
+# Plot the different histograms
+for i in range(num_order-1, -1, -1):
+ data = np.loadtxt(order_list[i])
+ ids = data[:,0]
+ pos = data[:,1:4]
+ a_exact = data[:,4:7]
+ a_grav = data[:, 7:10]
+
+ # Sort stuff
+ sort_index = np.argsort(ids)
+ ids = ids[sort_index]
+ pos = pos[sort_index, :]
+ a_exact = a_exact[sort_index, :]
+ a_grav = a_grav[sort_index, :]
+
+ # Cross-checks
+ if not np.array_equal(ids, gadget2_ids):
+ print "Comparing different IDs !"
+
+ if not np.array_equal(pos, gadget2_pos):
+ print "Comparing different positions ! max difference:", np.max(pos - gadget2_pos)
+
+ if not np.array_equal(a_exact, gadget2_a_exact):
+ print "Comparing different exact accelerations ! max difference:", np.max(a_exact - gadget2_a_exact)
+
+
+ # Compute the error norm
+ diff = a_exact - a_grav
+
+ norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
+ norm_a = np.sqrt(a_exact[:,0]**2 + a_exact[:,1]**2 + a_exact[:,2]**2)
+
+ norm_error = norm_diff / norm_a
+ error_x = abs(diff[:,0]) / norm_a
+ error_y = abs(diff[:,1]) / norm_a
+ error_z = abs(diff[:,2]) / norm_a
+
+ # Bin the error
+ norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+
+ norm_median = np.median(norm_error)
+ median_x = np.median(error_x)
+ median_y = np.median(error_y)
+ median_z = np.median(error_z)
+
+ norm_per95 = np.percentile(norm_error,95)
+ per95_x = np.percentile(error_x,95)
+ per95_y = np.percentile(error_y,95)
+ per95_z = np.percentile(error_z,95)
+
+ plt.subplot(221)
+ plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.plot([norm_median, norm_median], [2.7, 3],'-', color=cols[i], lw=1)
+ plt.plot([norm_per95, norm_per95], [2.7, 3],':', color=cols[i], lw=1)
+ plt.subplot(222)
+ plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.plot([median_x, median_x], [1.8, 2],'-', color=cols[i], lw=1)
+ plt.plot([per95_x, per95_x], [1.8, 2],':', color=cols[i], lw=1)
+ plt.subplot(223)
+ plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.plot([median_y, median_y], [1.8, 2],'-', color=cols[i], lw=1)
+ plt.plot([per95_y, per95_y], [1.8, 2],':', color=cols[i], lw=1)
+ plt.subplot(224)
+ plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.plot([median_z, median_z], [1.8, 2],'-', color=cols[i], lw=1)
+ plt.plot([per95_z, per95_z], [1.8, 2],':', color=cols[i], lw=1)
+
+
+plt.subplot(221)
+plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,3)
+plt.legend(loc="upper left")
+plt.subplot(222)
+plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+plt.subplot(223)
+plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+plt.subplot(224)
+plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
+plt.ylabel("Density")
+plt.xlim(min_error, 2*max_error)
+plt.ylim(0,2)
+#plt.legend(loc="center left")
+
+
+
+plt.savefig("gravity_checks_step%d.png"%step)
diff --git a/examples/main.c b/examples/main.c
index 034b800887928c049a610c27ef7c916573c71be6..2b741f75b236de9a3d64d382c442d4f53713e1b2 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -323,14 +323,14 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */
if (myrank == 0) {
- message("sizeof(part) is %4zi bytes.", sizeof(struct part));
- message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
- message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
- message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
- message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
- message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
- message("sizeof(task) is %4zi bytes.", sizeof(struct task));
- message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
+ message("sizeof(part) is %4zi bytes.", sizeof(struct part));
+ message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
+ message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
+ message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
+ message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
+ message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
+ message("sizeof(task) is %4zi bytes.", sizeof(struct task));
+ message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
}
/* Read the parameter file */
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index 978448b3cd049c6ff31a92c7255851390ccc700c..1be59d1c8449970321b8ef9053ddf24b4559dabd 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -57,14 +57,15 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init", "ghost", "extra_ghost", "drift", "kick1", "kick2",
- "timestep", "send", "recv", "grav_gather_m", "grav_fft",
- "grav_mm", "grav_up", "cooling", "sourceterms", "count"]
+ "timestep", "send", "recv", "grav_top_level", "grav_long_range",
+ "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
- "tend", "xv", "rho", "gpart", "count"]
+ "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
# Task/subtypes of interest.
-FULLTYPES = ["self/force", "self/density", "sub_self/force",
- "sub_self/density", "pair/force", "pair/density", "sub_pair/force",
+FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
+ "sub_self/density", "pair/force", "pair/density", "pair/grav",
+ "sub_pair/force",
"sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
"recv/tend", "send/tend"]
diff --git a/src/Makefile.am b/src/Makefile.am
index a29cf5521f2b25cbcc91ef85e277d1340ee69dd5..e95a00125e3c86ec32d3ae632065fd1d071877ac 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
- hydro_space.h
+ vector_power.h hydro_space.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
diff --git a/src/cell.c b/src/cell.c
index 753bdd55061ebd2b2cdd2067691bd8319ca5a623..874f03f0cad3257d866b70ec63fd8a6bbcd4b6a7 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1114,7 +1114,7 @@ void cell_check_multipole(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
struct gravity_tensors ma;
- const double tolerance = 1e-5; /* Relative */
+ const double tolerance = 1e-3; /* Relative */
/* First recurse */
if (c->split)
@@ -1127,8 +1127,7 @@ void cell_check_multipole(struct cell *c, void *data) {
gravity_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */
- if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole,
- tolerance)) {
+ if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth);
message("Correct answer:");
gravity_multipole_print(&ma.m_pole);
@@ -1487,7 +1486,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
} else if (ti_current > ti_old_multipole) {
/* Drift the multipole */
- gravity_multipole_drift(c->multipole, dt);
+ gravity_drift(c->multipole, dt);
}
/* Update the time of the last drift */
@@ -1504,7 +1503,21 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
* @param e The #engine (to get ti_current).
*/
void cell_drift_multipole(struct cell *c, const struct engine *e) {
- error("To be implemented");
+
+ const double timeBase = e->timeBase;
+ const integertime_t ti_old_multipole = c->ti_old_multipole;
+ const integertime_t ti_current = e->ti_current;
+
+ /* Drift from the last time the cell was drifted to the current time */
+ const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+ if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+
+ /* Update the time of the last drift */
+ c->ti_old_multipole = ti_current;
}
/**
diff --git a/src/const.h b/src/const.h
index 5e0fac1b5c8db823a528f095f2ab1a2f1e856d33..6962ee8bca32e92664e3f20cdb23e7cb6fbc4abd 100644
--- a/src/const.h
+++ b/src/const.h
@@ -39,9 +39,6 @@
/* Thermal energy per unit mass used as a constant for the isothermal EoS */
#define const_isothermal_internal_energy 20.2615290634f
-/* Self gravity stuff. */
-#define const_gravity_multipole_order 1
-
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
//#define GRADIENTS_SPH
diff --git a/src/engine.c b/src/engine.c
index 98f5cba933158c4f202377e270a6ac51b6541de3..d87d7bafda1ab5b7fb95fa769cf5991a7404ad6b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3029,6 +3029,13 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Computing initial gas densities.");
+ /* Initialise the softening lengths */
+ if (e->policy & engine_policy_self_gravity) {
+
+ for (size_t i = 0; i < s->nr_gparts; ++i)
+ gravity_init_softening(&s->gparts[i], e->gravity_properties);
+ }
+
/* Construct all cells and tasks to start everything */
engine_rebuild(e);
@@ -3061,16 +3068,17 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
}
#ifdef SWIFT_DEBUG_CHECKS
- /* Let's store the total mass in the system for future checks */
- e->s->total_mass = 0.;
- for (size_t i = 0; i < s->nr_gparts; ++i)
- e->s->total_mass += s->gparts[i].mass;
-#ifdef WITH_MPI
- if (MPI_Allreduce(MPI_IN_PLACE, &e->s->total_mass, 1, MPI_DOUBLE, MPI_SUM,
- MPI_COMM_WORLD) != MPI_SUCCESS)
- error("Failed to all-reduce total mass in the system.");
-#endif
- message("Total mass in the system: %e", e->s->total_mass);
+ /* Check that we have the correct total mass in the top-level multipoles */
+ size_t num_gpart_mpole = 0;
+ if (e->policy & engine_policy_self_gravity) {
+ for (int i = 0; i < e->s->nr_cells; ++i)
+ num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+ if (num_gpart_mpole != e->s->nr_gparts)
+ error(
+ "Multipoles don't contain the total number of gpart s->nr_gpart=%zd, "
+ "m_poles=%zd",
+ e->s->nr_gparts, num_gpart_mpole);
+ }
#endif
/* Now time to get ready for the first time-step */
@@ -3085,9 +3093,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Run the brute-force gravity calculation for some gparts */
+ gravity_exact_force_compute(e->s, e);
+#endif
+
/* Run the 0th time-step */
engine_launch(e, e->nr_threads);
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Check the accuracy of the gravity calculation */
+ gravity_exact_force_check(e->s, e, 1e-1);
+#endif
+
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
@@ -3160,6 +3178,17 @@ void engine_step(struct engine *e) {
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we have the correct total mass in the top-level multipoles */
+ size_t num_gpart_mpole = 0;
+ if (e->policy & engine_policy_self_gravity) {
+ for (int i = 0; i < e->s->nr_cells; ++i)
+ num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+ if (num_gpart_mpole != e->s->nr_gparts)
+ error("Multipoles don't contain the total number of gpart");
+ }
+#endif
+
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Run the brute-force gravity calculation for some gparts */
gravity_exact_force_compute(e->s, e);
diff --git a/src/gravity.c b/src/gravity.c
index 86f9fa82e3eb693cc3c051420fb9c7bff277eb9f..369c6b1b0ab0458f7c6ab5057261a7cade97a64c 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -20,6 +20,9 @@
/* Config parameters. */
#include "../config.h"
+/* Some standard headers. */
+#include <stdio.h>
+
/* This object's header. */
#include "gravity.h"
@@ -53,9 +56,7 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
gpart_is_active(gpi, e)) {
/* Be ready for the calculation */
- gpi->a_grav[0] = 0.f;
- gpi->a_grav[1] = 0.f;
- gpi->a_grav[2] = 0.f;
+ double a_grav[3] = {0., 0., 0.};
/* Interact it with all other particles in the space.*/
for (size_t j = 0; j < s->nr_gparts; ++j) {
@@ -66,36 +67,55 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
struct gpart *gpj = &s->gparts[j];
/* Compute the pairwise distance. */
- float dx[3] = {gpi->x[0] - gpj->x[0], // x
- gpi->x[1] - gpj->x[1], // y
- gpi->x[2] - gpj->x[2]}; // z
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ const double dx[3] = {gpi->x[0] - gpj->x[0], // x
+ gpi->x[1] - gpj->x[1], // y
+ gpi->x[2] - gpj->x[2]}; // z
+ const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- runner_iact_grav_pp_nonsym(0.f, r2, dx, gpi, gpj);
- }
+ const double r = sqrtf(r2);
+ const double ir = 1.f / r;
+ const double mj = gpj->mass;
+ const double hi = gpi->epsilon;
+ double f;
+ const double f_lr = 1.;
- /* Finish the calculation */
- gravity_end_force(gpi, const_G);
+ if (r >= hi) {
- /* Store the exact answer */
- gpi->a_grav_exact[0] = gpi->a_grav[0];
- gpi->a_grav_exact[1] = gpi->a_grav[1];
- gpi->a_grav_exact[2] = gpi->a_grav[2];
+ /* Get Newtonian gravity */
+ f = mj * ir * ir * ir * f_lr;
+
+ } else {
+
+ const double hi_inv = 1.f / hi;
+ const double hi_inv3 = hi_inv * hi_inv * hi_inv;
+ const double ui = r * hi_inv;
+ float W;
+
+ kernel_grav_eval(ui, &W);
- /* Restore everything */
- gpi->a_grav[0] = 0.f;
- gpi->a_grav[1] = 0.f;
- gpi->a_grav[2] = 0.f;
+ /* Get softened gravity */
+ f = mj * hi_inv3 * W * f_lr;
+ }
+
+ const double fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
+
+ a_grav[0] -= fdx[0];
+ a_grav[1] -= fdx[1];
+ a_grav[2] -= fdx[2];
+ }
+
+ /* Store the exact answer */
+ gpi->a_grav_exact[0] = a_grav[0] * const_G;
+ gpi->a_grav_exact[1] = a_grav[1] * const_G;
+ gpi->a_grav_exact[2] = a_grav[2] * const_G;
counter++;
}
}
- message("Computed exact gravity for %d gparts.", counter);
+ message("Computed exact gravity for %d gparts (took %.3f %s). ", counter,
+ clocks_from_ticks(getticks() - tic), clocks_getunit());
- if (e->verbose)
- message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
- clocks_getunit());
#else
error("Gravity checking function called without the corresponding flag.");
#endif
@@ -118,17 +138,24 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
- const double const_G = e->physical_constants->const_newton_G;
-
int counter = 0;
/* Some accumulators */
- float err_rel[3];
- float err_rel_max[3] = {0.f, 0.f, 0.f};
- float err_rel_min[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
- float err_rel_mean[3] = {0.f, 0.f, 0.f};
- float err_rel_mean2[3] = {0.f, 0.f, 0.f};
- float err_rel_std[3] = {0.f, 0.f, 0.f};
+ float err_rel_max = 0.f;
+ float err_rel_min = FLT_MAX;
+ float err_rel_mean = 0.f;
+ float err_rel_mean2 = 0.f;
+ float err_rel_std = 0.f;
+
+ char file_name[100];
+ sprintf(file_name, "gravity_checks_step%d_order%d.dat", e->step,
+ SELF_GRAVITY_MULTIPOLE_ORDER);
+ FILE *file = fopen(file_name, "w");
+ fprintf(file, "# Gravity accuracy test G = %16.8e\n",
+ e->physical_constants->const_newton_G);
+ fprintf(file, "# %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
+ "pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
+ "a_exact[2]", "a_grav[0]", "a_grav[1]", "a_grav[2]");
for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -138,48 +165,62 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
+ fprintf(file,
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e "
+ "%16.8e \n",
+ gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2]);
+
+ const float diff[3] = {gpi->a_grav[0] - gpi->a_grav_exact[0],
+ gpi->a_grav[1] - gpi->a_grav_exact[1],
+ gpi->a_grav[2] - gpi->a_grav_exact[2]};
+
+ const float diff_norm =
+ sqrtf(diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]);
+ const float a_norm = sqrtf(gpi->a_grav_exact[0] * gpi->a_grav_exact[0] +
+ gpi->a_grav_exact[1] * gpi->a_grav_exact[1] +
+ gpi->a_grav_exact[2] * gpi->a_grav_exact[2]);
+
/* Compute relative error */
- for (int k = 0; k < 3; ++k)
- if (fabsf(gpi->a_grav_exact[k]) > FLT_EPSILON * const_G)
- err_rel[k] = (gpi->a_grav[k] - gpi->a_grav_exact[k]) /
- fabsf(gpi->a_grav_exact[k]);
- else
- err_rel[k] = 0.f;
+ const float err_rel = diff_norm / a_norm;
/* Check that we are not above tolerance */
- if (fabsf(err_rel[0]) > rel_tol || fabsf(err_rel[1]) > rel_tol ||
- fabsf(err_rel[2]) > rel_tol)
- error("Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e]",
- gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
- gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
-
- /* Construct some statistics */
- for (int k = 0; k < 3; ++k) {
- err_rel_max[k] = max(err_rel_max[k], fabsf(err_rel[k]));
- err_rel_min[k] = min(err_rel_min[k], fabsf(err_rel[k]));
- err_rel_mean[k] += err_rel[k];
- err_rel_mean2[k] += err_rel[k] * err_rel[k];
+ if (err_rel > rel_tol) {
+ message(
+ "Error too large ! gp->a_grav=[%3.6e %3.6e %3.6e] "
+ "gp->a_exact=[%3.6e %3.6e %3.6e], "
+ "gp->num_interacted=%lld, err=%f",
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
+ gpi->num_interacted, err_rel);
+
+ continue;
}
+ /* Construct some statistics */
+ err_rel_max = max(err_rel_max, fabsf(err_rel));
+ err_rel_min = min(err_rel_min, fabsf(err_rel));
+ err_rel_mean += err_rel;
+ err_rel_mean2 += err_rel * err_rel;
counter++;
}
}
+ /* Be nice */
+ fclose(file);
+
/* Final operation on the stats */
if (counter > 0) {
- for (int k = 0; k < 3; ++k) {
- err_rel_mean[k] /= counter;
- err_rel_mean2[k] /= counter;
- err_rel_std[k] =
- sqrtf(err_rel_mean2[k] - err_rel_mean[k] * err_rel_mean[k]);
- }
+ err_rel_mean /= counter;
+ err_rel_mean2 /= counter;
+ err_rel_std = sqrtf(err_rel_mean2 - err_rel_mean * err_rel_mean);
}
/* Report on the findings */
message("Checked gravity for %d gparts.", counter);
- for (int k = 0; k < 3; ++k)
- message("Error on a_grav[%d]: min=%e max=%e mean=%e std=%e", k,
- err_rel_min[k], err_rel_max[k], err_rel_mean[k], err_rel_std[k]);
+ message("Error on |a_grav|: min=%e max=%e mean=%e std=%e", err_rel_min,
+ err_rel_max, err_rel_mean, err_rel_std);
#else
error("Gravity checking function called without the corresponding flag.");
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 93a65e2f5a70e09ad14280cf9c334753359fb8b5..790aed2550a5a75fd123d08f813f7328491c0fc6 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -27,6 +27,8 @@
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
+ * We use Gadget-2's type 0 time-step criterion.
+ *
* @param gp Pointer to the g-particle data.
* @param grav_props Constants used in the gravity scheme.
*/
@@ -38,9 +40,10 @@ gravity_compute_timestep_self(const struct gpart* const gp,
gp->a_grav[1] * gp->a_grav[1] +
gp->a_grav[2] * gp->a_grav[2];
- const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
+ const float ac_inv = (ac2 > 0.f) ? 1.f/sqrtf(ac2) : FLT_MAX;
- const float dt = sqrtf(2.f * grav_props->eta * gp->epsilon / ac);
+ /* Note that 0.714285714 = 2. (from Gadget) / 2.8 (Plummer softening) */
+ const float dt = sqrtf(0.714285714f * grav_props->eta * gp->epsilon * ac_inv);
return dt;
}
@@ -62,7 +65,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
gp->a_grav[2] = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
- gp->mass_interacted = 0.;
+ gp->num_interacted = 0;
#endif
}
@@ -113,9 +116,22 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
gp->time_bin = 0;
- gp->epsilon = 0.; // MATTHIEU
+ gp->epsilon = 0.f;
gravity_init_gpart(gp);
}
+/**
+ * @brief Initialises the softening of the g-particles
+ *
+ * @param gp The particle to act upon.
+ * @param grav_props The properties of the gravity scheme.
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_softening(
+ struct gpart* gp, const struct gravity_props* grav_props) {
+
+ /* Note 2.8 is the Plummer-equivalent correction */
+ gp->epsilon = 2.8f * grav_props->epsilon;
+}
+
#endif /* SWIFT_DEFAULT_GRAVITY_H */
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 7c7c5b9d244534ef3f6b5f509062019bfcd5f9fb..eca5c2491cbdcf5f0eca01417c8e6b29efc53459 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -44,6 +44,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
const float u = r * rlr_inv;
float f_lr, fi, fj, W;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
/* Get long-range correction */
kernel_long_grav_eval(u, &f_lr);
@@ -107,6 +111,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
const float u = r * rlr_inv;
float f_lr, f, W;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r == 0.f) error("Interacting particles with 0 distance");
+#endif
+
/* Get long-range correction */
kernel_long_grav_eval(u, &f_lr);
@@ -141,51 +149,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
float rlr_inv, float r2, const float *dx, struct gpart *gp,
const struct multipole *multi) {
- /* Apply the gravitational acceleration. */
- const float r = sqrtf(r2);
- const float ir = 1.f / r;
- const float u = r * rlr_inv;
- const float mrinv3 = multi->mass * ir * ir * ir;
-
- /* Get long-range correction */
- float f_lr;
- kernel_long_grav_eval(u, &f_lr);
-
-#if const_gravity_multipole_order < 2
-
- /* 0th and 1st order terms */
- gp->a_grav[0] += mrinv3 * f_lr * dx[0];
- gp->a_grav[1] += mrinv3 * f_lr * dx[1];
- gp->a_grav[2] += mrinv3 * f_lr * dx[2];
-
-#elif const_gravity_multipole_order == 2
- /* Terms up to 2nd order (quadrupole) */
-
- /* Follows the notation in Bonsai */
- const float mrinv5 = mrinv3 * ir * ir;
- const float mrinv7 = mrinv5 * ir * ir;
-
- const float D1 = -mrinv3;
- const float D2 = 3.f * mrinv5;
- const float D3 = -15.f * mrinv7;
-
- const float q = multi->I_xx + multi->I_yy + multi->I_zz;
- const float qRx =
- multi->I_xx * dx[0] + multi->I_xy * dx[1] + multi->I_xz * dx[2];
- const float qRy =
- multi->I_xy * dx[0] + multi->I_yy * dx[1] + multi->I_yz * dx[2];
- const float qRz =
- multi->I_xz * dx[0] + multi->I_yz * dx[1] + multi->I_zz * dx[2];
- const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
- const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
-
- gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx);
- gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy);
- gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz);
-
-#else
-#error "Multipoles of order >2 not yet implemented."
-#endif
+ error("Dead function");
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 00ae0f5b05cd95750c34b60e2353a9fc1d0a5c32..bb1307a1a2fc1dcd71202e3426c99f3e30c0de9a 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -52,8 +52,8 @@ struct gpart {
#ifdef SWIFT_DEBUG_CHECKS
- /* Total mass this gpart interacted with */
- double mass_interacted;
+ /* Numer of gparts this gpart interacted with */
+ long long num_interacted;
/* Time of the last drift */
integertime_t ti_drift;
@@ -66,7 +66,7 @@ struct gpart {
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Brute-force particle acceleration. */
- float a_grav_exact[3];
+ double a_grav_exact[3];
#endif
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 4730d9df5dc573de74fde422e1b7dafc0ee0994a..e5c7722bc3e5ebed50b1062c74e3dad830c9b145 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -19,14 +19,31 @@
#ifndef SWIFT_GRAVITY_DERIVATIVE_H
#define SWIFT_GRAVITY_DERIVATIVE_H
+/**
+ * @file gravity_derivatives.h
+ * @brief Derivatives (up to 3rd order) of the gravitational potential.
+ *
+ * We use the notation of Dehnen, Computational Astrophysics and Cosmology,
+ * 1, 1, pp. 24 (2014), arXiv:1405.2255
+ */
+
/* Some standard headers. */
#include <math.h>
/* Local headers. */
#include "inline.h"
+/*************************/
+/* 0th order derivatives */
+/*************************/
+
/**
* @brief \f$ \phi(r_x, r_y, r_z) \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
*/
__attribute__((always_inline)) INLINE static double D_000(double r_x,
double r_y,
@@ -36,8 +53,17 @@ __attribute__((always_inline)) INLINE static double D_000(double r_x,
return r_inv;
}
+/*************************/
+/* 1st order derivatives */
+/*************************/
+
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
*/
__attribute__((always_inline)) INLINE static double D_100(double r_x,
double r_y,
@@ -49,6 +75,11 @@ __attribute__((always_inline)) INLINE static double D_100(double r_x,
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
*/
__attribute__((always_inline)) INLINE static double D_010(double r_x,
double r_y,
@@ -60,6 +91,11 @@ __attribute__((always_inline)) INLINE static double D_010(double r_x,
/**
* @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
*/
__attribute__((always_inline)) INLINE static double D_001(double r_x,
double r_y,
@@ -69,4 +105,306 @@ __attribute__((always_inline)) INLINE static double D_001(double r_x,
return -r_z * r_inv * r_inv * r_inv;
}
+/*************************/
+/* 2nd order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_200(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_x * r_x * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_020(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_y * r_y * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_z^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_002(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_z * r_z * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_110(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_x * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_101(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_x * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_011(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_y * r_z * r_inv5;
+}
+
+/*************************/
+/* 3rd order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_300(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_x * r_inv7 + 9. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_030(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_y * r_y * r_y * r_inv7 + 9. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_003(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_z * r_z * r_z * r_inv7 + 9. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_y}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_210(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_y * r_inv7 + 3. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_201(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_z * r_inv7 + 3. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_120(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_y * r_y * r_inv7 + 3. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^2\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_021(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_z * r_y * r_y * r_inv7 + 3. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_102(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_z * r_z * r_inv7 + 3. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_012(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_y * r_z * r_z * r_inv7 + 3. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z\partial
+ * r_y\partial r_z} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+__attribute__((always_inline)) INLINE static double D_111(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+ const double r_inv3 = r_inv * r_inv * r_inv;
+ const double r_inv7 = r_inv3 * r_inv3 * r_inv;
+ return -15. * r_x * r_y * r_z * r_inv7;
+}
+
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index a1e382a21d04b7354aaf215069e999627e56ee07..a9425c7dda8f864a1d6b9c2a658b308c7e7e0624 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -25,6 +25,7 @@
/* Includes. */
#include "const.h"
#include "inline.h"
+#include "minmax.h"
#include "vector.h"
/* The gravity kernel is defined as a degree 6 polynomial in the distance
@@ -72,7 +73,7 @@ __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)fminf(u * (float)kernel_grav_ivals, kernel_grav_ivals);
+ 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)];
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 6952681999f833bce7755a72aaee742a7fa0ed22..7b1c5984647c3be232770dc32fc1b112ad8bee94 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -37,14 +37,15 @@
__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
float u, float *const W) {
- const float arg1 = u * 0.5f;
- const float arg2 = u * one_over_sqrt_pi;
- const float arg3 = -arg1 * arg1;
+ /* const float arg1 = u * 0.5f; */
+ /* const float arg2 = u * one_over_sqrt_pi; */
+ /* const float arg3 = -arg1 * arg1; */
- const float term1 = erfcf(arg1);
- const float term2 = arg2 * expf(arg3);
+ /* const float term1 = erfcf(arg1); */
+ /* const float term2 = arg2 * expf(arg3); */
- *W = term1 + term2;
+ /* *W = term1 + term2; */
+ *W = 1.f;
}
#endif // SWIFT_KERNEL_LONG_GRAVITY_H
diff --git a/src/multipole.h b/src/multipole.h
index b5c9335ee8fabf46740cefe310fcfecbea3fd77e..9be734ff5f4eea2f74946c8555911b1b970a7105 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -20,6 +20,9 @@
#ifndef SWIFT_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H
+/* Config parameters. */
+#include "../config.h"
+
/* Some standard headers. */
#include <math.h>
#include <string.h>
@@ -34,25 +37,104 @@
#include "kernel_gravity.h"
#include "minmax.h"
#include "part.h"
+#include "vector_power.h"
#define multipole_align 128
-struct acc_tensor {
+struct grav_tensor {
- double F_000;
-};
+ /* 0th order terms */
+ float F_000;
-struct pot_tensor {
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- double F_000;
+ /* 1st order terms */
+ float F_100, F_010, F_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order terms */
+ float F_200, F_020, F_002;
+ float F_110, F_101, F_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 3rd order terms */
+ float F_300, F_030, F_003;
+ float F_210, F_201;
+ float F_120, F_021;
+ float F_102, F_012;
+ float F_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 4th order terms */
+ float F_400, F_040, F_004;
+ float F_310, F_301;
+ float F_130, F_031;
+ float F_103, F_013;
+ float F_220, F_202, F_022;
+ float F_211, F_121, F_112;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Total number of gpart this field tensor interacted with */
+ long long num_interacted;
+
+#endif
};
struct multipole {
- float mass;
-
- /*! Bulk velocity */
+ /* Bulk velocity */
float vel[3];
+
+ /* 0th order terms */
+ float M_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+ /* 1st order terms */
+ float M_100, M_010, M_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order terms */
+ float M_200, M_020, M_002;
+ float M_110, M_101, M_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 3rd order terms */
+ float M_300, M_030, M_003;
+ float M_210, M_201;
+ float M_120, M_021;
+ float M_102, M_012;
+ float M_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 4th order terms */
+ float M_400, M_040, M_004;
+ float M_310, M_301;
+ float M_130, M_031;
+ float M_103, M_013;
+ float M_220, M_202, M_022;
+ float M_211, M_121, M_112;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Total number of gpart in this multipole */
+ long long num_gpart;
+
+#endif
};
/**
@@ -72,24 +154,8 @@ struct gravity_tensors {
/*! Multipole mass */
struct multipole m_pole;
- /*! Field tensor for acceleration along x */
- struct acc_tensor a_x;
-
- /*! Field tensor for acceleration along y */
- struct acc_tensor a_y;
-
- /*! Field tensor for acceleration along z */
- struct acc_tensor a_z;
-
/*! Field tensor for the potential */
- struct pot_tensor pot;
-
-#ifdef SWIFT_DEBUG_CHECKS
-
- /* Total mass this gpart interacted with */
- double mass_interacted;
-
-#endif
+ struct grav_tensor pot;
};
} SWIFT_STRUCT_ALIGN;
@@ -104,15 +170,141 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
bzero(m, sizeof(struct gravity_tensors));
}
-INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
+/**
+ * @brief Drifts a #multipole forward in time.
+ *
+ * @param m The #multipole.
+ * @param dt The drift time-step.
+ */
+INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
- bzero(&m->a_x, sizeof(struct acc_tensor));
- bzero(&m->a_y, sizeof(struct acc_tensor));
- bzero(&m->a_z, sizeof(struct acc_tensor));
- bzero(&m->pot, sizeof(struct pot_tensor));
+ /* Move the whole thing according to bulk motion */
+ m->CoM[0] += m->m_pole.vel[0];
+ m->CoM[1] += m->m_pole.vel[1];
+ m->CoM[2] += m->m_pole.vel[2];
+}
+
+INLINE static void gravity_field_tensors_init(struct grav_tensor *l) {
+ bzero(l, sizeof(struct grav_tensor));
+}
+
+/**
+ * @brief Adds field tensrs to other ones (i.e. does la += lb).
+ *
+ * @param la The gravity tensors to add to.
+ * @param lb The gravity tensors to add.
+ */
+INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
+ const struct grav_tensor *lb) {
#ifdef SWIFT_DEBUG_CHECKS
- m->mass_interacted = 0.;
+ if (lb->num_interacted == 0) error("Adding tensors that did not interact");
+ la->num_interacted += lb->num_interacted;
+#endif
+
+ /* Add 0th order terms */
+ la->F_000 += lb->F_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* Add 1st order terms */
+ la->F_100 += lb->F_100;
+ la->F_010 += lb->F_010;
+ la->F_001 += lb->F_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ /* Add 2nd order terms */
+ la->F_200 += lb->F_200;
+ la->F_020 += lb->F_020;
+ la->F_002 += lb->F_002;
+ la->F_110 += lb->F_110;
+ la->F_101 += lb->F_101;
+ la->F_011 += lb->F_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ /* Add 3rd order terms */
+ la->F_300 += lb->F_300;
+ la->F_030 += lb->F_030;
+ la->F_003 += lb->F_003;
+ la->F_210 += lb->F_210;
+ la->F_201 += lb->F_201;
+ la->F_120 += lb->F_120;
+ la->F_021 += lb->F_021;
+ la->F_102 += lb->F_102;
+ la->F_012 += lb->F_012;
+ la->F_111 += lb->F_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ /* Add 4th order terms */
+ la->F_400 += lb->F_400;
+ la->F_040 += lb->F_040;
+ la->F_004 += lb->F_004;
+ la->F_310 += lb->F_310;
+ la->F_301 += lb->F_301;
+ la->F_130 += lb->F_130;
+ la->F_031 += lb->F_031;
+ la->F_103 += lb->F_103;
+ la->F_013 += lb->F_013;
+ la->F_220 += lb->F_220;
+ la->F_202 += lb->F_202;
+ la->F_022 += lb->F_022;
+ la->F_211 += lb->F_211;
+ la->F_121 += lb->F_121;
+ la->F_112 += lb->F_112;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+}
+
+/**
+ * @brief Prints the content of a #grav_tensor to stdout.
+ *
+ * Note: Uses directly printf(), not a call to message().
+ *
+ * @param l The #grav_tensor to print.
+ */
+INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {
+
+ printf("-------------------------\n");
+ printf("F_000= %12.5e\n", l->F_000);
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ printf("-------------------------\n");
+ printf("F_100= %12.5e F_010= %12.5e F_001= %12.5e\n", l->F_100, l->F_010,
+ l->F_001);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ printf("-------------------------\n");
+ printf("F_200= %12.5e F_020= %12.5e F_002= %12.5e\n", l->F_200, l->F_020,
+ l->F_002);
+ printf("F_110= %12.5e F_101= %12.5e F_011= %12.5e\n", l->F_110, l->F_101,
+ l->F_011);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ printf("-------------------------\n");
+ printf("F_300= %12.5e F_030= %12.5e F_003= %12.5e\n", l->F_300, l->F_030,
+ l->F_003);
+ printf("F_210= %12.5e F_201= %12.5e F_120= %12.5e\n", l->F_210, l->F_201,
+ l->F_120);
+ printf("F_021= %12.5e F_102= %12.5e F_012= %12.5e\n", l->F_021, l->F_102,
+ l->F_012);
+ printf("F_111= %12.5e\n", l->F_111);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ printf("-------------------------\n");
+ printf("F_400= %12.5e F_040= %12.5e F_004= %12.5e\n", l->F_400, l->F_040,
+ l->F_004);
+ printf("F_310= %12.5e F_301= %12.5e F_130= %12.5e\n", l->F_310, l->F_301,
+ l->F_130);
+ printf("F_031= %12.5e F_103= %12.5e F_013= %12.5e\n", l->F_031, l->F_103,
+ l->F_013);
+ printf("F_220= %12.5e F_202= %12.5e F_022= %12.5e\n", l->F_220, l->F_202,
+ l->F_022);
+ printf("F_211= %12.5e F_121= %12.5e F_112= %12.5e\n", l->F_211, l->F_121,
+ l->F_112);
+#endif
+ printf("-------------------------\n");
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
#endif
}
@@ -125,9 +317,48 @@ INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
*/
INLINE static void gravity_multipole_print(const struct multipole *m) {
- // printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
- printf("Mass= %12.5e\n", m->mass);
- printf("Vel= [%12.5e %12.5e %12.5e\n", m->vel[0], m->vel[1], m->vel[2]);
+ printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
+ printf("-------------------------\n");
+ printf("M_000= %12.5e\n", m->M_000);
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ printf("-------------------------\n");
+ printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
+ m->M_001);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ printf("-------------------------\n");
+ printf("M_200= %12.5e M_020= %12.5e M_002= %12.5e\n", m->M_200, m->M_020,
+ m->M_002);
+ printf("M_110= %12.5e M_101= %12.5e M_011= %12.5e\n", m->M_110, m->M_101,
+ m->M_011);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ printf("-------------------------\n");
+ printf("M_300= %12.5e M_030= %12.5e M_003= %12.5e\n", m->M_300, m->M_030,
+ m->M_003);
+ printf("M_210= %12.5e M_201= %12.5e M_120= %12.5e\n", m->M_210, m->M_201,
+ m->M_120);
+ printf("M_021= %12.5e M_102= %12.5e M_012= %12.5e\n", m->M_021, m->M_102,
+ m->M_012);
+ printf("M_111= %12.5e\n", m->M_111);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ printf("-------------------------\n");
+ printf("M_400= %12.5e M_040= %12.5e M_004= %12.5e\n", m->M_400, m->M_040,
+ m->M_004);
+ printf("M_310= %12.5e M_301= %12.5e M_130= %12.5e\n", m->M_310, m->M_301,
+ m->M_130);
+ printf("M_031= %12.5e M_103= %12.5e M_013= %12.5e\n", m->M_031, m->M_103,
+ m->M_013);
+ printf("M_220= %12.5e M_202= %12.5e M_022= %12.5e\n", m->M_220, m->M_202,
+ m->M_022);
+ printf("M_211= %12.5e M_121= %12.5e M_112= %12.5e\n", m->M_211, m->M_121,
+ m->M_112);
+#endif
+ printf("-------------------------\n");
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
}
/**
@@ -139,88 +370,335 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_add(struct multipole *ma,
const struct multipole *mb) {
- const float mass = ma->mass + mb->mass;
- const float imass = 1.f / mass;
+ const float M_000 = ma->M_000 + mb->M_000;
+ const float inv_M_000 = 1.f / M_000;
/* Add the bulk velocities */
- ma->vel[0] = (ma->vel[0] * ma->mass + mb->vel[0] * mb->mass) * imass;
- ma->vel[1] = (ma->vel[1] * ma->mass + mb->vel[1] * mb->mass) * imass;
- ma->vel[2] = (ma->vel[2] * ma->mass + mb->vel[2] * mb->mass) * imass;
+ ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
+ ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
+ ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
+
+ /* Add 0th order terms */
+ ma->M_000 = M_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* Add 1st order terms */
+ ma->M_100 += mb->M_100;
+ ma->M_010 += mb->M_010;
+ ma->M_001 += mb->M_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ /* Add 2nd order terms */
+ ma->M_200 += mb->M_200;
+ ma->M_020 += mb->M_020;
+ ma->M_002 += mb->M_002;
+ ma->M_110 += mb->M_110;
+ ma->M_101 += mb->M_101;
+ ma->M_011 += mb->M_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ /* Add 3rd order terms */
+ ma->M_300 += mb->M_300;
+ ma->M_030 += mb->M_030;
+ ma->M_003 += mb->M_003;
+ ma->M_210 += mb->M_210;
+ ma->M_201 += mb->M_201;
+ ma->M_120 += mb->M_120;
+ ma->M_021 += mb->M_021;
+ ma->M_102 += mb->M_102;
+ ma->M_012 += mb->M_012;
+ ma->M_111 += mb->M_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ /* Add 4th order terms */
+ ma->M_400 += mb->M_400;
+ ma->M_040 += mb->M_040;
+ ma->M_004 += mb->M_004;
+ ma->M_310 += mb->M_310;
+ ma->M_301 += mb->M_301;
+ ma->M_130 += mb->M_130;
+ ma->M_031 += mb->M_031;
+ ma->M_103 += mb->M_103;
+ ma->M_013 += mb->M_013;
+ ma->M_220 += mb->M_220;
+ ma->M_202 += mb->M_202;
+ ma->M_022 += mb->M_022;
+ ma->M_211 += mb->M_211;
+ ma->M_121 += mb->M_121;
+ ma->M_112 += mb->M_112;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
- /* Add the masses */
- ma->mass = mass;
+ // MATTHIEU
+ ma->M_100 = 0.f;
+ ma->M_010 = 0.f;
+ ma->M_001 = 0.f;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ ma->num_gpart += mb->num_gpart;
+#endif
}
/**
* @brief Verifies whether two #multipole's are equal or not.
*
- * @param ma The first #multipole.
- * @param mb The second #multipole.
- * @param tolerance The maximal allowed difference for the fields.
- * @return 1 if the multipoles are equal 0 otherwise.
+ * @param ga The first #multipole.
+ * @param gb The second #multipole.
+ * @param tolerance The maximal allowed relative difference for the fields.
+ * @return 1 if the multipoles are equal, 0 otherwise
*/
-INLINE static int gravity_multipole_equal(const struct multipole *ma,
- const struct multipole *mb,
+INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
+ const struct gravity_tensors *gb,
double tolerance) {
/* Check CoM */
- /* if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) >
- * tolerance) */
- /* return 0; */
- /* if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) >
- * tolerance) */
- /* return 0; */
- /* if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) >
- * tolerance) */
- /* return 0; */
-
- /* Check bulk velocity (if non-zero)*/
- if (fabsf(ma->vel[0] + mb->vel[0]) > 0.f &&
+ if (fabs(ga->CoM[0] - gb->CoM[0]) / fabs(ga->CoM[0] + gb->CoM[0]) >
+ tolerance) {
+ message("CoM[0] different");
+ return 0;
+ }
+ if (fabs(ga->CoM[1] - gb->CoM[1]) / fabs(ga->CoM[1] + gb->CoM[1]) >
+ tolerance) {
+ message("CoM[1] different");
+ return 0;
+ }
+ if (fabs(ga->CoM[2] - gb->CoM[2]) / fabs(ga->CoM[2] + gb->CoM[2]) >
+ tolerance) {
+ message("CoM[2] different");
+ return 0;
+ }
+
+ /* Helper pointers */
+ const struct multipole *ma = &ga->m_pole;
+ const struct multipole *mb = &gb->m_pole;
+
+ const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
+ ma->vel[2] * ma->vel[2];
+
+ /* Check bulk velocity (if non-zero and component > 1% of norm)*/
+ if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
+ (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
- tolerance)
+ tolerance) {
+ message("v[0] different");
return 0;
- if (fabsf(ma->vel[1] + mb->vel[1]) > 0.f &&
+ }
+ if (fabsf(ma->vel[1] + mb->vel[1]) > 1e-10 &&
+ (ma->vel[1] * ma->vel[1]) > 0.0001 * v2 &&
fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
- tolerance)
+ tolerance) {
+ message("v[1] different");
return 0;
- if (fabsf(ma->vel[2] + mb->vel[2]) > 0.f &&
+ }
+ if (fabsf(ma->vel[2] + mb->vel[2]) > 1e-10 &&
+ (ma->vel[2] * ma->vel[2]) > 0.0001 * v2 &&
fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
- tolerance)
+ tolerance) {
+ message("v[2] different");
return 0;
+ }
- /* Check mass */
- if (fabsf(ma->mass - mb->mass) / fabsf(ma->mass + mb->mass) > tolerance)
+ /* Check 0th order terms */
+ if (fabsf(ma->M_000 - mb->M_000) / fabsf(ma->M_000 + mb->M_000) > tolerance) {
+ message("M_000 term different");
return 0;
+ }
- /* All is good */
- return 1;
-}
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* Check 1st order terms */
+ if (fabsf(ma->M_100 + mb->M_100) > 1e-6 * ma->M_000 &&
+ fabsf(ma->M_100 - mb->M_100) / fabsf(ma->M_100 + mb->M_100) > tolerance) {
+ message("M_100 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_010 + mb->M_010) > 1e-6 * ma->M_000 &&
+ fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
+ message("M_010 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_001 + mb->M_001) > 1e-6 * ma->M_000 &&
+ fabsf(ma->M_001 - mb->M_001) / fabsf(ma->M_001 + mb->M_001) > tolerance) {
+ message("M_001 term different");
+ return 0;
+ }
+#endif
-/**
- * @brief Drifts a #multipole forward in time.
- *
- * @param m The #multipole.
- * @param dt The drift time-step.
- */
-INLINE static void gravity_multipole_drift(struct gravity_tensors *m,
- double dt) {
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ /* Check 2nd order terms */
+ if (fabsf(ma->M_200 + mb->M_200) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_200 - mb->M_200) / fabsf(ma->M_200 + mb->M_200) > tolerance) {
+ message("M_200 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_020 + mb->M_020) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_020 - mb->M_020) / fabsf(ma->M_020 + mb->M_020) > tolerance) {
+ message("M_020 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_002 + mb->M_002) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_002 - mb->M_002) / fabsf(ma->M_002 + mb->M_002) > tolerance) {
+ message("M_002 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_110 + mb->M_110) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_110 - mb->M_110) / fabsf(ma->M_110 + mb->M_110) > tolerance) {
+ message("M_110 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_101 + mb->M_101) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_101 - mb->M_101) / fabsf(ma->M_101 + mb->M_101) > tolerance) {
+ message("M_101 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_011 + mb->M_011) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_011 - mb->M_011) / fabsf(ma->M_011 + mb->M_011) > tolerance) {
+ message("M_011 term different");
+ return 0;
+ }
+#endif
- /* Move the whole thing according to bulk motion */
- m->CoM[0] += m->m_pole.vel[0];
- m->CoM[1] += m->m_pole.vel[1];
- m->CoM[2] += m->m_pole.vel[2];
-}
+ tolerance *= 10.;
-/**
- * @brief Applies the forces due to particles j onto particles i directly.
- *
- * @param gparts_i The #gpart to update.
- * @param gcount_i The number of particles to update.
- * @param gparts_j The #gpart that source the gravity field.
- * @param gcount_j The number of sources.
- */
-INLINE static void gravity_P2P(struct gpart *gparts_i, int gcount_i,
- const struct gpart *gparts_j, int gcount_j) {}
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ /* Check 3rd order terms */
+ if (fabsf(ma->M_300 + mb->M_300) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_300 - mb->M_300) / fabsf(ma->M_300 + mb->M_300) > tolerance) {
+ message("M_300 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_030 + mb->M_030) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_030 - mb->M_030) / fabsf(ma->M_030 + mb->M_030) > tolerance) {
+ message("M_030 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_003 + mb->M_003) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_003 - mb->M_003) / fabsf(ma->M_003 + mb->M_003) > tolerance) {
+ message("M_003 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_210 + mb->M_210) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_210 - mb->M_210) / fabsf(ma->M_210 + mb->M_210) > tolerance) {
+ message("M_210 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_201 + mb->M_201) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_201 - mb->M_201) / fabsf(ma->M_201 + mb->M_201) > tolerance) {
+ message("M_201 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_120 + mb->M_120) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_120 - mb->M_120) / fabsf(ma->M_120 + mb->M_120) > tolerance) {
+ message("M_120 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_021 + mb->M_021) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_021 - mb->M_021) / fabsf(ma->M_021 + mb->M_021) > tolerance) {
+ message("M_021 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_102 + mb->M_102) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_102 - mb->M_102) / fabsf(ma->M_102 + mb->M_102) > tolerance) {
+ message("M_102 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_012 + mb->M_012) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_012 - mb->M_012) / fabsf(ma->M_012 + mb->M_012) > tolerance) {
+ message("M_012 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_111 + mb->M_111) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_111 - mb->M_111) / fabsf(ma->M_111 + mb->M_111) > tolerance) {
+ message("M_111 term different");
+ return 0;
+ }
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ /* Check 4th order terms */
+ if (fabsf(ma->M_400 + mb->M_400) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_400 - mb->M_400) / fabsf(ma->M_400 + mb->M_400) > tolerance) {
+ message("M_400 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_040 + mb->M_040) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_040 - mb->M_040) / fabsf(ma->M_040 + mb->M_040) > tolerance) {
+ message("M_040 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_004 + mb->M_004) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_004 - mb->M_004) / fabsf(ma->M_004 + mb->M_004) > tolerance) {
+ message("M_003 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_310 + mb->M_310) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_310 - mb->M_310) / fabsf(ma->M_310 + mb->M_310) > tolerance) {
+ message("M_310 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_301 + mb->M_301) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_301 - mb->M_301) / fabsf(ma->M_301 + mb->M_301) > tolerance) {
+ message("M_301 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_130 + mb->M_130) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_130 - mb->M_130) / fabsf(ma->M_130 + mb->M_130) > tolerance) {
+ message("M_130 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_031 + mb->M_031) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_031 - mb->M_031) / fabsf(ma->M_031 + mb->M_031) > tolerance) {
+ message("M_031 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_103 + mb->M_103) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_103 - mb->M_103) / fabsf(ma->M_103 + mb->M_103) > tolerance) {
+ message("M_103 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_013 + mb->M_013) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_013 - mb->M_013) / fabsf(ma->M_013 + mb->M_013) > tolerance) {
+ message("M_013 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_220 + mb->M_220) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_220 - mb->M_220) / fabsf(ma->M_220 + mb->M_220) > tolerance) {
+ message("M_220 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_202 + mb->M_202) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_202 - mb->M_202) / fabsf(ma->M_202 + mb->M_202) > tolerance) {
+ message("M_202 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_022 + mb->M_022) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_022 - mb->M_022) / fabsf(ma->M_022 + mb->M_022) > tolerance) {
+ message("M_022 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_211 + mb->M_211) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_211 - mb->M_211) / fabsf(ma->M_211 + mb->M_211) > tolerance) {
+ message("M_211 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_121 + mb->M_121) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_121 - mb->M_121) / fabsf(ma->M_121 + mb->M_121) > tolerance) {
+ message("M_121 term different");
+ return 0;
+ }
+ if (fabsf(ma->M_112 + mb->M_112) > 1e-5 * ma->M_000 &&
+ fabsf(ma->M_112 - mb->M_112) / fabsf(ma->M_112 + mb->M_112) > tolerance) {
+ message("M_112 term different");
+ return 0;
+ }
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+
+ /* All is good */
+ return 1;
+}
/**
* @brief Constructs the #multipole of a bunch of particles around their
@@ -235,10 +713,6 @@ INLINE static void gravity_P2P(struct gpart *gparts_i, int gcount_i,
INLINE static void gravity_P2M(struct gravity_tensors *m,
const struct gpart *gparts, int gcount) {
-#if const_gravity_multipole_order >= 2
-#error "Implementation of P2M kernel missing for this order."
-#endif
-
/* Temporary variables */
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
@@ -257,16 +731,154 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
vel[2] += gparts[k].v_full[2] * m;
}
+ /* Final operation on CoM */
const double imass = 1.0 / mass;
+ com[0] *= imass;
+ com[1] *= imass;
+ com[2] *= imass;
+ vel[0] *= imass;
+ vel[1] *= imass;
+ vel[2] *= imass;
+
+/* Prepare some local counters */
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ float M_100 = 0.f, M_010 = 0.f, M_001 = 0.f;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ float M_200 = 0.f, M_020 = 0.f, M_002 = 0.f;
+ float M_110 = 0.f, M_101 = 0.f, M_011 = 0.f;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ float M_300 = 0.f, M_030 = 0.f, M_003 = 0.f;
+ float M_210 = 0.f, M_201 = 0.f, M_120 = 0.f;
+ float M_021 = 0.f, M_102 = 0.f, M_012 = 0.f;
+ float M_111 = 0.f;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ float M_400 = 0.f, M_040 = 0.f, M_004 = 0.f;
+ float M_310 = 0.f, M_301 = 0.f, M_130 = 0.f;
+ float M_031 = 0.f, M_103 = 0.f, M_013 = 0.f;
+ float M_220 = 0.f, M_202 = 0.f, M_022 = 0.f;
+ float M_211 = 0.f, M_121 = 0.f, M_112 = 0.f;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+
+ /* Construce the higher order terms */
+ for (int k = 0; k < gcount; k++) {
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const float m = gparts[k].mass;
+ const double dx[3] = {gparts[k].x[0] - com[0], gparts[k].x[1] - com[1],
+ gparts[k].x[2] - com[2]};
+
+ M_100 += -m * X_100(dx);
+ M_010 += -m * X_010(dx);
+ M_001 += -m * X_001(dx);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ M_200 += m * X_200(dx);
+ M_020 += m * X_020(dx);
+ M_002 += m * X_002(dx);
+ M_110 += m * X_110(dx);
+ M_101 += m * X_101(dx);
+ M_011 += m * X_011(dx);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ M_300 += -m * X_300(dx);
+ M_030 += -m * X_030(dx);
+ M_003 += -m * X_003(dx);
+ M_210 += -m * X_210(dx);
+ M_201 += -m * X_201(dx);
+ M_120 += -m * X_120(dx);
+ M_021 += -m * X_021(dx);
+ M_102 += -m * X_102(dx);
+ M_012 += -m * X_012(dx);
+ M_111 += -m * X_111(dx);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ M_400 += m * X_400(dx);
+ M_040 += m * X_040(dx);
+ M_004 += m * X_004(dx);
+ M_310 += m * X_310(dx);
+ M_301 += m * X_301(dx);
+ M_130 += m * X_130(dx);
+ M_031 += m * X_031(dx);
+ M_103 += m * X_103(dx);
+ M_013 += m * X_013(dx);
+ M_220 += m * X_220(dx);
+ M_202 += m * X_202(dx);
+ M_022 += m * X_022(dx);
+ M_211 += m * X_211(dx);
+ M_121 += m * X_121(dx);
+ M_112 += m * X_112(dx);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+ }
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ M_100 = M_010 = M_001 = 0.f; /* Matthieu */
+#endif
/* Store the data on the multipole. */
- m->m_pole.mass = mass;
- m->CoM[0] = com[0] * imass;
- m->CoM[1] = com[1] * imass;
- m->CoM[2] = com[2] * imass;
- m->m_pole.vel[0] = vel[0] * imass;
- m->m_pole.vel[1] = vel[1] * imass;
- m->m_pole.vel[2] = vel[2] * imass;
+ m->m_pole.M_000 = mass;
+ m->CoM[0] = com[0];
+ m->CoM[1] = com[1];
+ m->CoM[2] = com[2];
+ m->m_pole.vel[0] = vel[0];
+ m->m_pole.vel[1] = vel[1];
+ m->m_pole.vel[2] = vel[2];
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ m->m_pole.M_100 = M_100;
+ m->m_pole.M_010 = M_010;
+ m->m_pole.M_001 = M_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ m->m_pole.M_200 = M_200;
+ m->m_pole.M_020 = M_020;
+ m->m_pole.M_002 = M_002;
+ m->m_pole.M_110 = M_110;
+ m->m_pole.M_101 = M_101;
+ m->m_pole.M_011 = M_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ m->m_pole.M_300 = M_300;
+ m->m_pole.M_030 = M_030;
+ m->m_pole.M_003 = M_003;
+ m->m_pole.M_210 = M_210;
+ m->m_pole.M_201 = M_201;
+ m->m_pole.M_120 = M_120;
+ m->m_pole.M_021 = M_021;
+ m->m_pole.M_102 = M_102;
+ m->m_pole.M_012 = M_012;
+ m->m_pole.M_111 = M_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ m->m_pole.M_400 = M_400;
+ m->m_pole.M_040 = M_040;
+ m->m_pole.M_004 = M_004;
+ m->m_pole.M_310 = M_310;
+ m->m_pole.M_301 = M_301;
+ m->m_pole.M_130 = M_130;
+ m->m_pole.M_031 = M_031;
+ m->m_pole.M_103 = M_103;
+ m->m_pole.M_013 = M_013;
+ m->m_pole.M_220 = M_220;
+ m->m_pole.M_202 = M_202;
+ m->m_pole.M_022 = M_022;
+ m->m_pole.M_211 = M_211;
+ m->m_pole.M_121 = M_121;
+ m->m_pole.M_112 = M_112;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m->m_pole.num_gpart = gcount;
+#endif
}
/**
@@ -284,184 +896,406 @@ INLINE static void gravity_M2M(struct multipole *m_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
-
- m_a->mass = m_b->mass;
-
+ /* Shift bulk velocity */
m_a->vel[0] = m_b->vel[0];
m_a->vel[1] = m_b->vel[1];
m_a->vel[2] = m_b->vel[2];
+
+ /* Shift 0th order term */
+ m_a->M_000 = m_b->M_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
+ pos_a[2] - pos_b[2]};
+
+ /* Shift 1st order term */
+ m_a->M_100 = m_b->M_100 + X_100(dx) * m_b->M_000;
+ m_a->M_010 = m_b->M_010 + X_010(dx) * m_b->M_000;
+ m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* Shift 2nd order term */
+ m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
+ m_a->M_020 = m_b->M_020 + X_010(dx) * m_b->M_010 + X_020(dx) * m_b->M_000;
+ m_a->M_002 = m_b->M_002 + X_001(dx) * m_b->M_001 + X_002(dx) * m_b->M_000;
+ m_a->M_110 = m_b->M_110 + X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +
+ X_110(dx) * m_b->M_000;
+ m_a->M_101 = m_b->M_101 + X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +
+ X_101(dx) * m_b->M_000;
+ m_a->M_011 = m_b->M_011 + X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +
+ X_011(dx) * m_b->M_000;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* Shift 3rd order term */
+ m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + X_200(dx) * m_b->M_100 +
+ X_300(dx) * m_b->M_000;
+ m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + X_020(dx) * m_b->M_010 +
+ X_030(dx) * m_b->M_000;
+ m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + X_002(dx) * m_b->M_001 +
+ X_003(dx) * m_b->M_000;
+ m_a->M_210 = m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
+ X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 +
+ X_210(dx) * m_b->M_000;
+ m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 +
+ X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +
+ X_201(dx) * m_b->M_000;
+ m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 +
+ X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +
+ X_120(dx) * m_b->M_000;
+ m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 +
+ X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +
+ X_021(dx) * m_b->M_000;
+ m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 +
+ X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +
+ X_102(dx) * m_b->M_000;
+ m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 +
+ X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +
+ X_012(dx) * m_b->M_000;
+ m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 +
+ X_001(dx) * m_b->M_110 + X_110(dx) * m_b->M_001 +
+ X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +
+ X_111(dx) * m_b->M_000;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#error "Missing implementation for order >3"
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m_a->num_gpart = m_b->num_gpart;
+#endif
}
+
/**
* @brief Compute the field tensors due to a multipole.
*
* Corresponds to equation (28b).
*
- * @param l_a The field tensor to compute.
- * @param m_b The multipole creating the field.
- * @param pos_a The position of the field tensor.
- * @param pos_b The position of the multipole.
+ * @param l_b The field tensor to compute.
+ * @param m_a The multipole creating the field.
+ * @param pos_b The position of the field tensor.
+ * @param pos_a The position of the multipole.
* @param periodic Is the calculation periodic ?
*/
-INLINE static void gravity_M2L(struct gravity_tensors *l_a,
- const struct multipole *m_b,
- const double pos_a[3], const double pos_b[3],
+INLINE static void gravity_M2L(struct grav_tensor *l_b,
+ const struct multipole *m_a,
+ const double pos_b[3], const double pos_a[3],
int periodic) {
double dx, dy, dz;
if (periodic) {
- dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.);
- dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.);
- dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.);
+ dx = box_wrap(pos_b[0] - pos_a[0], 0., 1.);
+ dy = box_wrap(pos_b[1] - pos_a[1], 0., 1.);
+ dz = box_wrap(pos_b[2] - pos_a[2], 0., 1.);
} else {
- dx = pos_a[0] - pos_b[0];
- dy = pos_a[1] - pos_b[1];
- dz = pos_a[2] - pos_b[2];
+ dx = pos_b[0] - pos_a[0];
+ dy = pos_b[1] - pos_a[1];
+ dz = pos_b[2] - pos_a[2];
}
const double r2 = dx * dx + dy * dy + dz * dz;
const double r_inv = 1. / sqrt(r2);
- /* 1st order multipole term */
- l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
- l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
- l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
+ /* 0th order term */
+ l_b->F_000 += m_a->M_000 * D_000(dx, dy, dz, r_inv);
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* 1st order multipole term (addition to rank 0)*/
+ l_b->F_000 += m_a->M_100 * D_100(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_010(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_001(dx, dy, dz, r_inv);
+
+ /* 1st order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_000 * D_100(dx, dy, dz, r_inv);
+ l_b->F_010 += m_a->M_000 * D_010(dx, dy, dz, r_inv);
+ l_b->F_001 += m_a->M_000 * D_001(dx, dy, dz, r_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order multipole term (addition to rank 0)*/
+ l_b->F_000 += m_a->M_200 * D_200(dx, dy, dz, r_inv) +
+ m_a->M_020 * D_020(dx, dy, dz, r_inv) +
+ m_a->M_002 * D_002(dx, dy, dz, r_inv);
+ l_b->F_000 += m_a->M_110 * D_110(dx, dy, dz, r_inv) +
+ m_a->M_101 * D_101(dx, dy, dz, r_inv) +
+ m_a->M_011 * D_011(dx, dy, dz, r_inv);
+
+ /* 2nd order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_100 * D_200(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_110(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_101(dx, dy, dz, r_inv);
+ l_b->F_010 += m_a->M_100 * D_110(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_020(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_011(dx, dy, dz, r_inv);
+ l_b->F_001 += m_a->M_100 * D_101(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_011(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_002(dx, dy, dz, r_inv);
+
+ /* 2nd order multipole term (addition to rank 2)*/
+ l_b->F_200 += m_a->M_000 * D_200(dx, dy, dz, r_inv);
+ l_b->F_020 += m_a->M_000 * D_020(dx, dy, dz, r_inv);
+ l_b->F_002 += m_a->M_000 * D_002(dx, dy, dz, r_inv);
+ l_b->F_110 += m_a->M_000 * D_110(dx, dy, dz, r_inv);
+ l_b->F_101 += m_a->M_000 * D_101(dx, dy, dz, r_inv);
+ l_b->F_011 += m_a->M_000 * D_011(dx, dy, dz, r_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 3rd order multipole term (addition to rank 0)*/
+ l_b->F_000 += m_a->M_300 * D_300(dx, dy, dz, r_inv) +
+ m_a->M_030 * D_030(dx, dy, dz, r_inv) +
+ m_a->M_003 * D_003(dx, dy, dz, r_inv);
+ l_b->F_000 += m_a->M_210 * D_210(dx, dy, dz, r_inv) +
+ m_a->M_201 * D_201(dx, dy, dz, r_inv) +
+ m_a->M_120 * D_120(dx, dy, dz, r_inv);
+ l_b->F_000 += m_a->M_021 * D_021(dx, dy, dz, r_inv) +
+ m_a->M_102 * D_102(dx, dy, dz, r_inv) +
+ m_a->M_012 * D_012(dx, dy, dz, r_inv);
+ l_b->F_000 += m_a->M_111 * D_111(dx, dy, dz, r_inv);
+
+ /* 3rd order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_200 * D_300(dx, dy, dz, r_inv) +
+ m_a->M_020 * D_120(dx, dy, dz, r_inv) +
+ m_a->M_002 * D_102(dx, dy, dz, r_inv);
+ l_b->F_100 += m_a->M_110 * D_210(dx, dy, dz, r_inv) +
+ m_a->M_101 * D_201(dx, dy, dz, r_inv) +
+ m_a->M_011 * D_111(dx, dy, dz, r_inv);
+ l_b->F_010 += m_a->M_200 * D_210(dx, dy, dz, r_inv) +
+ m_a->M_020 * D_030(dx, dy, dz, r_inv) +
+ m_a->M_002 * D_012(dx, dy, dz, r_inv);
+ l_b->F_010 += m_a->M_110 * D_120(dx, dy, dz, r_inv) +
+ m_a->M_101 * D_111(dx, dy, dz, r_inv) +
+ m_a->M_011 * D_021(dx, dy, dz, r_inv);
+ l_b->F_001 += m_a->M_200 * D_201(dx, dy, dz, r_inv) +
+ m_a->M_020 * D_021(dx, dy, dz, r_inv) +
+ m_a->M_002 * D_003(dx, dy, dz, r_inv);
+ l_b->F_001 += m_a->M_110 * D_111(dx, dy, dz, r_inv) +
+ m_a->M_101 * D_102(dx, dy, dz, r_inv) +
+ m_a->M_011 * D_012(dx, dy, dz, r_inv);
+
+ /* 3rd order multipole term (addition to rank 2)*/
+ l_b->F_200 += m_a->M_100 * D_300(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_210(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_201(dx, dy, dz, r_inv);
+ l_b->F_020 += m_a->M_100 * D_120(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_030(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_021(dx, dy, dz, r_inv);
+ l_b->F_002 += m_a->M_100 * D_102(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_012(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_003(dx, dy, dz, r_inv);
+ l_b->F_110 += m_a->M_100 * D_210(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_120(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_111(dx, dy, dz, r_inv);
+ l_b->F_101 += m_a->M_100 * D_201(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_111(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_102(dx, dy, dz, r_inv);
+ l_b->F_011 += m_a->M_100 * D_111(dx, dy, dz, r_inv) +
+ m_a->M_010 * D_021(dx, dy, dz, r_inv) +
+ m_a->M_001 * D_012(dx, dy, dz, r_inv);
+
+ /* 3rd order multipole term (addition to rank 2)*/
+ l_b->F_300 += m_a->M_000 * D_300(dx, dy, dz, r_inv);
+ l_b->F_030 += m_a->M_000 * D_030(dx, dy, dz, r_inv);
+ l_b->F_003 += m_a->M_000 * D_003(dx, dy, dz, r_inv);
+ l_b->F_210 += m_a->M_000 * D_210(dx, dy, dz, r_inv);
+ l_b->F_201 += m_a->M_000 * D_201(dx, dy, dz, r_inv);
+ l_b->F_120 += m_a->M_000 * D_120(dx, dy, dz, r_inv);
+ l_b->F_021 += m_a->M_000 * D_021(dx, dy, dz, r_inv);
+ l_b->F_102 += m_a->M_000 * D_102(dx, dy, dz, r_inv);
+ l_b->F_012 += m_a->M_000 * D_012(dx, dy, dz, r_inv);
+ l_b->F_111 += m_a->M_000 * D_111(dx, dy, dz, r_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#error "Missing implementation for order >3"
+#endif
#ifdef SWIFT_DEBUG_CHECKS
- l_a->mass_interacted += m_b->mass;
+
+ l_b->num_interacted += m_a->num_gpart;
#endif
}
/**
- * @brief Creates a copy of #acc_tensor shifted to a new location.
+ * @brief Creates a copy of #grav_tensor shifted to a new location.
*
* Corresponds to equation (28e).
*
- * @param l_a The #acc_tensor copy (content will be overwritten).
- * @param l_b The #acc_tensor to shift.
+ * @param la The #grav_tensor copy (content will be overwritten).
+ * @param lb The #grav_tensor to shift.
* @param pos_a The position to which m_b will be shifted.
* @param pos_b The current postion of the multipole to shift.
* @param periodic Is the calculation periodic ?
*/
-INLINE static void gravity_L2L(struct gravity_tensors *l_a,
- const struct gravity_tensors *l_b,
+INLINE static void gravity_L2L(struct grav_tensor *la,
+ const struct grav_tensor *lb,
const double pos_a[3], const double pos_b[3],
- int periodic) {}
-
-/**
- * @brief Applies the #acc_tensor to a set of #gpart.
- *
- * Corresponds to equation (28a).
- */
-INLINE static void gravity_L2P(const struct gravity_tensors *l,
- struct gpart *gparts, int gcount) {
+ int periodic) {
- for (int i = 0; i < gcount; ++i) {
+ /* Initialise everything to zero */
+ gravity_field_tensors_init(la);
#ifdef SWIFT_DEBUG_CHECKS
- struct gpart *gp = &gparts[i];
-
- // if(gpart_is_active(gp, e)){
-
- gp->mass_interacted += l->mass_interacted;
+ if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
+ la->num_interacted = lb->num_interacted;
#endif
- //}
- }
-}
-#if 0
+ /* Distance to shift by */
+ const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
+ pos_a[2] - pos_b[2]};
-/* Multipole function prototypes. */
-void multipole_add(struct gravity_tensors *m_sum, const struct gravity_tensors *m_term);
-void multipole_init(struct gravity_tensors *m, const struct gpart *gparts,
- int gcount);
-void multipole_reset(struct gravity_tensors *m);
+ /* Shift 0th order term */
+ la->F_000 += X_000(dx) * lb->F_000;
-/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
-/* double *shift); */
-/* void multipole_addpart(struct multipole *m, struct gpart *p); */
-/* void multipole_addparts(struct multipole *m, struct gpart *p, int N); */
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* Shift 1st order multipole term (addition to rank 0)*/
+ la->F_000 +=
+ X_100(dx) * lb->F_100 + X_010(dx) * lb->F_010 + X_001(dx) * lb->F_001;
-/**
- * @brief Compute the pairwise interaction between two multipoles.
- *
- * @param ma The first #multipole.
- * @param mb The second #multipole.
- * @param shift The periodicity correction.
- */
-__attribute__((always_inline)) INLINE static void multipole_iact_mm(
- struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
- /* float dx[3], ir, r, r2 = 0.0f, acc; */
- /* int k; */
-
- /* /\* Compute the multipole distance. *\/ */
- /* for (k = 0; k < 3; k++) { */
- /* dx[k] = ma->x[k] - mb->x[k] - shift[k]; */
- /* r2 += dx[k] * dx[k]; */
- /* } */
-
- /* /\* Compute the normalized distance vector. *\/ */
- /* ir = 1.0f / sqrtf(r2); */
- /* r = r2 * ir; */
-
- /* /\* Evaluate the gravity kernel. *\/ */
- /* kernel_grav_eval(r, &acc); */
-
- /* /\* Scale the acceleration. *\/ */
- /* acc *= const_G * ir * ir * ir; */
-
- /* /\* Compute the forces on both multipoles. *\/ */
- /* #if const_gravity_multipole_order == 1 */
- /* float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
- /* for (k = 0; k < 3; k++) { */
- /* ma->a[k] -= dx[k] * acc * mmb; */
- /* mb->a[k] += dx[k] * acc * mma; */
- /* } */
- /* #else */
- /* #error( "Multipoles of order %i not yet implemented." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
+ /* Shift 1st order multipole term (addition to rank 1)*/
+ la->F_100 += X_000(dx) * lb->F_100;
+ la->F_010 += X_000(dx) * lb->F_010;
+ la->F_001 += X_000(dx) * lb->F_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* Shift 2nd order multipole term (addition to rank 0)*/
+ la->F_000 +=
+ X_200(dx) * lb->F_200 + X_020(dx) * lb->F_020 + X_002(dx) * lb->F_002;
+ la->F_000 +=
+ X_110(dx) * lb->F_110 + X_101(dx) * lb->F_101 + X_011(dx) * lb->F_011;
+
+ /* Shift 2nd order multipole term (addition to rank 1)*/
+ la->F_100 +=
+ X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
+ la->F_010 +=
+ X_100(dx) * lb->F_110 + X_010(dx) * lb->F_020 + X_001(dx) * lb->F_011;
+ la->F_001 +=
+ X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;
+
+ /* Shift 2nd order multipole term (addition to rank 2)*/
+ la->F_200 += X_000(dx) * lb->F_200;
+ la->F_020 += X_000(dx) * lb->F_020;
+ la->F_002 += X_000(dx) * lb->F_002;
+ la->F_110 += X_000(dx) * lb->F_110;
+ la->F_101 += X_000(dx) * lb->F_101;
+ la->F_011 += X_000(dx) * lb->F_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* Shift 3rd order multipole term (addition to rank 0)*/
+ la->F_000 +=
+ X_300(dx) * lb->F_300 + X_030(dx) * lb->F_030 + X_003(dx) * lb->F_003;
+ la->F_000 +=
+ X_210(dx) * lb->F_210 + X_201(dx) * lb->F_201 + X_120(dx) * lb->F_120;
+ la->F_000 +=
+ X_021(dx) * lb->F_021 + X_102(dx) * lb->F_102 + X_012(dx) * lb->F_012;
+ la->F_000 += X_111(dx) * lb->F_111;
+
+ /* Shift 3rd order multipole term (addition to rank 1)*/
+ la->F_100 +=
+ X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
+ la->F_100 +=
+ X_110(dx) * lb->F_210 + X_101(dx) * lb->F_201 + X_011(dx) * lb->F_111;
+ la->F_010 +=
+ X_200(dx) * lb->F_210 + X_020(dx) * lb->F_030 + X_002(dx) * lb->F_012;
+ la->F_010 +=
+ X_110(dx) * lb->F_120 + X_101(dx) * lb->F_111 + X_011(dx) * lb->F_021;
+ la->F_001 +=
+ X_200(dx) * lb->F_201 + X_020(dx) * lb->F_021 + X_002(dx) * lb->F_003;
+ la->F_001 +=
+ X_110(dx) * lb->F_111 + X_101(dx) * lb->F_102 + X_011(dx) * lb->F_012;
+
+ /* Shift 3rd order multipole term (addition to rank 2)*/
+ la->F_200 +=
+ X_100(dx) * lb->F_300 + X_010(dx) * lb->F_210 + X_001(dx) * lb->F_201;
+ la->F_020 +=
+ X_100(dx) * lb->F_120 + X_010(dx) * lb->F_030 + X_001(dx) * lb->F_021;
+ la->F_002 +=
+ X_100(dx) * lb->F_102 + X_010(dx) * lb->F_012 + X_001(dx) * lb->F_003;
+ la->F_110 +=
+ X_100(dx) * lb->F_210 + X_010(dx) * lb->F_120 + X_001(dx) * lb->F_111;
+ la->F_101 +=
+ X_100(dx) * lb->F_201 + X_010(dx) * lb->F_111 + X_001(dx) * lb->F_102;
+ la->F_011 +=
+ X_100(dx) * lb->F_111 + X_010(dx) * lb->F_021 + X_001(dx) * lb->F_012;
+
+ /* Shift 3rd order multipole term (addition to rank 2)*/
+ la->F_300 += X_000(dx) * lb->F_300;
+ la->F_030 += X_000(dx) * lb->F_030;
+ la->F_003 += X_000(dx) * lb->F_003;
+ la->F_210 += X_000(dx) * lb->F_210;
+ la->F_201 += X_000(dx) * lb->F_201;
+ la->F_120 += X_000(dx) * lb->F_120;
+ la->F_021 += X_000(dx) * lb->F_021;
+ la->F_102 += X_000(dx) * lb->F_102;
+ la->F_012 += X_000(dx) * lb->F_012;
+ la->F_111 += X_000(dx) * lb->F_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#error "Missing implementation for order >3"
+#endif
}
/**
- * @brief Compute the interaction of a multipole on a particle.
+ * @brief Applies the #grav_tensor to a #gpart.
*
- * @param m The #multipole.
- * @param p The #gpart.
- * @param shift The periodicity correction.
+ * Corresponds to equation (28a).
+ *
+ * @param lb The gravity field tensor to apply.
+ * @param loc The position of the gravity field tensor.
+ * @param gp The #gpart to update.
*/
-__attribute__((always_inline)) INLINE static void multipole_iact_mp(
- struct gravity_tensors *m, struct gpart *p, double *shift) {
-
- /* float dx[3], ir, r, r2 = 0.0f, acc; */
- /* int k; */
-
- /* /\* Compute the multipole distance. *\/ */
- /* for (k = 0; k < 3; k++) { */
- /* dx[k] = m->x[k] - p->x[k] - shift[k]; */
- /* r2 += dx[k] * dx[k]; */
- /* } */
-
- /* /\* Compute the normalized distance vector. *\/ */
- /* ir = 1.0f / sqrtf(r2); */
- /* r = r2 * ir; */
-
- /* /\* Evaluate the gravity kernel. *\/ */
- /* kernel_grav_eval(r, &acc); */
-
- /* /\* Scale the acceleration. *\/ */
- /* acc *= const_G * ir * ir * ir * m->coeffs[0]; */
-
- /* /\* Compute the forces on both multipoles. *\/ */
- /* #if const_gravity_multipole_order == 1 */
- /* for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
- /* #else */
- /* #error( "Multipoles of order %i not yet implemented." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
-}
+INLINE static void gravity_L2P(const struct grav_tensor *lb,
+ const double loc[3], struct gpart *gp) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (lb->num_interacted == 0) error("Interacting with empty field tensor");
+ gp->num_interacted += lb->num_interacted;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* Distance to the multipole */
+ const double dx[3] = {gp->x[0] - loc[0], gp->x[1] - loc[1],
+ gp->x[2] - loc[2]};
+
+ /* 0th order interaction */
+ gp->a_grav[0] += X_000(dx) * lb->F_100;
+ gp->a_grav[1] += X_000(dx) * lb->F_010;
+ gp->a_grav[2] += X_000(dx) * lb->F_001;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 1st order interaction */
+ gp->a_grav[0] +=
+ X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
+ gp->a_grav[1] +=
+ X_100(dx) * lb->F_110 + X_010(dx) * lb->F_020 + X_001(dx) * lb->F_011;
+ gp->a_grav[2] +=
+ X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;
#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 2nd order interaction */
+ gp->a_grav[0] +=
+ X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
+ gp->a_grav[0] +=
+ X_110(dx) * lb->F_210 + X_101(dx) * lb->F_201 + X_011(dx) * lb->F_111;
+ gp->a_grav[1] +=
+ X_200(dx) * lb->F_210 + X_020(dx) * lb->F_030 + X_002(dx) * lb->F_012;
+ gp->a_grav[1] +=
+ X_110(dx) * lb->F_120 + X_101(dx) * lb->F_111 + X_011(dx) * lb->F_021;
+ gp->a_grav[2] +=
+ X_200(dx) * lb->F_201 + X_020(dx) * lb->F_021 + X_002(dx) * lb->F_003;
+ gp->a_grav[2] +=
+ X_110(dx) * lb->F_111 + X_101(dx) * lb->F_102 + X_011(dx) * lb->F_012;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#error "Missing implementation for order >3"
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+}
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner.c b/src/runner.c
index 8240f77763d73d3ebc1e7a402194054b5c9ea7d5..c4504c22be21df4f4c6154fe570c9ee7d982c5e5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -498,7 +498,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Reset the gravity acceleration tensors */
if (e->policy & engine_policy_self_gravity)
- gravity_field_tensor_init(c->multipole);
+ gravity_field_tensors_init(&c->multipole->pot);
/* Recurse? */
if (c->split) {
@@ -1373,13 +1373,17 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
}
#ifdef SWIFT_DEBUG_CHECKS
- if (e->policy & engine_policy_self_gravity) {
- gp->mass_interacted += gp->mass;
- if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
+ if (e->policy & engine_policy_self_gravity && gpart_is_active(gp, e)) {
+
+ /* Check that this gpart has interacted with all the other particles
+ * (via direct or multipoles) in the box */
+ gp->num_interacted++;
+ if (gp->num_interacted != (long long)e->s->nr_gparts)
error(
- "g-particle did not interact gravitationally with all other "
- "particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
- gp->mass_interacted, e->s->total_mass, gp->mass);
+ "g-particle (id=%lld, type=%d) did not interact gravitationally "
+ "with all other gparts gp->num_interacted=%lld, total_gparts=%zd",
+ gp->id_or_neg_offset, gp->type, gp->num_interacted,
+ e->s->nr_gparts);
}
#endif
}
@@ -1846,7 +1850,7 @@ void *runner_main(void *data) {
// runner_do_grav_mm(r, t->ci, 1);
break;
case task_type_grav_down:
- runner_do_grav_down(r, t->ci);
+ runner_do_grav_down(r, t->ci, 1);
break;
case task_type_grav_top_level:
// runner_do_grav_top_level(r);
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9988b7d553a962ee541f20cab64cc534c991957a..b6c93e72f1c8d31d555a3306bf35098f0125ddf5 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -26,51 +26,54 @@
#include "part.h"
/**
- * @brief Compute the recursive upward sweep, i.e. construct the
- * multipoles in a cell hierarchy.
+ * @brief Recursively propagate the multipoles down the tree by applying the
+ * L2L and L2P kernels.
*
* @param r The #runner.
- * @param c The top-level #cell.
+ * @param c The #cell we are working on.
+ * @param timer Are we timing this ?
*/
-void runner_do_grav_up(struct runner *r, struct cell *c) {
+void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
- /* if (c->split) { /\* Regular node *\/ */
-
- /* /\* Recurse. *\/ */
- /* for (int k = 0; k < 8; k++) */
- /* if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]); */
-
- /* /\* Collect the multipoles from the progeny. *\/ */
- /* multipole_reset(&c->multipole); */
- /* for (int k = 0; k < 8; k++) { */
- /* if (c->progeny[k] != NULL) */
- /* multipole_add(&c->multipole, &c->progeny[k]->multipole); */
- /* } */
-
- /* } else { /\* Leaf node. *\/ */
- /* /\* Just construct the multipole from the gparts. *\/ */
- /* multipole_init(&c->multipole, c->gparts, c->gcount); */
- /* } */
-}
+ const struct engine *e = r->e;
+ const int periodic = e->s->periodic;
-void runner_do_grav_down(struct runner *r, struct cell *c) {
+ TIMER_TIC;
if (c->split) {
for (int k = 0; k < 8; ++k) {
struct cell *cp = c->progeny[k];
- struct gravity_tensors temp;
+ struct grav_tensor temp;
if (cp != NULL) {
- gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
- 1);
+
+ /* Shift the field tensor */
+ gravity_L2L(&temp, &c->multipole->pot, cp->multipole->CoM,
+ c->multipole->CoM, 0 * periodic);
+ /* Add it to this level's tensor */
+ gravity_field_tensors_add(&cp->multipole->pot, &temp);
+
+ /* Recurse */
+ runner_do_grav_down(r, cp, 0);
}
}
- } else {
+ } else { /* Leaf case */
+
+ const struct engine *e = r->e;
+ struct gpart *gparts = c->gparts;
+ const int gcount = c->gcount;
- gravity_L2P(c->multipole, c->gparts, c->gcount);
+ /* Apply accelerations to the particles */
+ for (int i = 0; i < gcount; ++i) {
+ struct gpart *gp = &gparts[i];
+ if (gpart_is_active(gp, e))
+ gravity_L2P(&c->multipole->pot, c->multipole->CoM, gp);
+ }
}
+
+ if (timer) TIMER_TOC(timer_dograv_down);
}
/**
@@ -81,9 +84,8 @@ void runner_do_grav_down(struct runner *r, struct cell *c) {
* @param ci The #cell with field tensor to interact.
* @param cj The #cell with the multipole.
*/
-__attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
- const struct runner *r, const struct cell *restrict ci,
- const struct cell *restrict cj) {
+void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
+ struct cell *restrict cj) {
const struct engine *e = r->e;
const int periodic = e->s->periodic;
@@ -94,14 +96,20 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (multi_j->mass == 0.0) error("Multipole does not seem to have been set.");
+ if (ci == cj) error("Interacting a cell with itself using M2L");
+
+ if (multi_j->M_000 == 0.f) error("Multipole does not seem to have been set.");
#endif
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
- gravity_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
- periodic);
+ /* Do we need to drift the multipole ? */
+ if (cj->ti_old_multipole != e->ti_current) cell_drift_multipole(cj, e);
+
+ /* Let's interact at this level */
+ gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
+ cj->multipole->CoM, periodic * 0);
TIMER_TOC(timer_dopair_grav_mm);
}
@@ -114,65 +122,11 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
* @param ci The #cell with particles to interct.
* @param cj The #cell with the multipole.
*/
-__attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
- const struct runner *r, const struct cell *restrict ci,
- const struct cell *restrict cj) {
-
- const struct engine *e = r->e;
- const int gcount = ci->gcount;
- struct gpart *restrict gparts = ci->gparts;
- const struct gravity_tensors *multi = cj->multipole;
- const float a_smooth = e->gravity_properties->a_smooth;
- const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
-
- TIMER_TIC;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (gcount == 0) error("Empty cell!");
-
- if (multi->m_pole.mass == 0.0)
- error("Multipole does not seem to have been set.");
-#endif
-
- /* Anything to do here? */
- if (!cell_is_active(ci, e)) return;
-
-#if ICHECK > 0
- for (int pid = 0; pid < gcount; pid++) {
-
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts[pid];
+void runner_dopair_grav_pm(const struct runner *r,
+ const struct cell *restrict ci,
+ const struct cell *restrict cj) {
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
- gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2],
- cj->width[0], cj->gcount);
- }
-#endif
-
- /* Loop over every particle in leaf. */
- for (int pid = 0; pid < gcount; pid++) {
-
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts[pid];
-
- if (!gpart_is_active(gp, e)) continue;
-
- /* Compute the pairwise distance. */
- const float dx[3] = {multi->CoM[0] - gp->x[0], // x
- multi->CoM[1] - gp->x[1], // y
- multi->CoM[2] - gp->x[2]}; // z
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
- /* Interact !*/
- runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole);
-
-#ifdef SWIFT_DEBUG_CHECKS
- gp->mass_interacted += multi->m_pole.mass;
-#endif
- }
-
- TIMER_TOC(timer_dopair_grav_pm);
+ error("Function should not be called");
}
/**
@@ -185,8 +139,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
*
* @todo Use a local cache for the particles.
*/
-__attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
- struct runner *r, struct cell *ci, struct cell *cj) {
+void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
const struct engine *e = r->e;
const int gcount_i = ci->gcount;
@@ -199,7 +152,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (ci->width[0] != cj->width[0]) // MATTHIEU sanity check
+ if (ci->width[0] != cj->width[0])
error("Non matching cell sizes !! h_i=%f h_j=%f", ci->width[0],
cj->width[0]);
#endif
@@ -237,49 +190,55 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
/* Get a hold of the ith part in ci. */
struct gpart *restrict gpi = &gparts_i[pid];
+ if (!gpart_is_active(gpi, e)) continue;
+
/* Loop over every particle in the other cell. */
for (int pjd = 0; pjd < gcount_j; pjd++) {
/* Get a hold of the jth part in cj. */
- struct gpart *restrict gpj = &gparts_j[pjd];
+ const struct gpart *restrict gpj = &gparts_j[pjd];
/* Compute the pairwise distance. */
- float dx[3] = {gpi->x[0] - gpj->x[0], // x
- gpi->x[1] - gpj->x[1], // y
- gpi->x[2] - gpj->x[2]}; // z
+ const float dx[3] = {gpi->x[0] - gpj->x[0], // x
+ gpi->x[1] - gpj->x[1], // y
+ gpi->x[2] - gpj->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact ! */
- if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
-
- runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
- gpj->mass_interacted += gpi->mass;
+ gpi->num_interacted++;
#endif
+ }
+ }
- } else {
+ /* Loop over all particles in cj... */
+ for (int pjd = 0; pjd < gcount_j; pjd++) {
- if (gpart_is_active(gpi, e)) {
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gpj = &gparts_j[pjd];
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+ if (!gpart_is_active(gpj, e)) continue;
-#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
-#endif
- } else if (gpart_is_active(gpj, e)) {
+ /* Loop over every particle in the other cell. */
+ for (int pid = 0; pid < gcount_i; pid++) {
- dx[0] = -dx[0];
- dx[1] = -dx[1];
- dx[2] = -dx[2];
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+ /* Get a hold of the ith part in ci. */
+ const struct gpart *restrict gpi = &gparts_i[pid];
+
+ /* Compute the pairwise distance. */
+ const float dx[3] = {gpj->x[0] - gpi->x[0], // x
+ gpj->x[1] - gpi->x[1], // y
+ gpj->x[2] - gpi->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* Interact ! */
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
- gpj->mass_interacted += gpi->mass;
+ gpj->num_interacted++;
#endif
- }
- }
}
}
@@ -294,8 +253,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
*
* @todo Use a local cache for the particles.
*/
-__attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
- struct runner *r, struct cell *c) {
+void runner_doself_grav_pp(struct runner *r, struct cell *c) {
const struct engine *e = r->e;
const int gcount = c->gcount;
@@ -306,8 +264,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (c->gcount == 0) // MATTHIEU sanity check
- error("Empty cell !");
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
/* Anything to do here? */
@@ -350,8 +307,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
- gpj->mass_interacted += gpi->mass;
+ gpi->num_interacted++;
+ gpj->num_interacted++;
#endif
} else {
@@ -361,7 +318,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
+ gpi->num_interacted++;
#endif
} else if (gpart_is_active(gpj, e)) {
@@ -372,7 +329,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
- gpj->mass_interacted += gpi->mass;
+ gpj->num_interacted++;
#endif
}
}
@@ -393,8 +350,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
*
* @todo Use a local cache for the particles.
*/
-static void runner_dopair_grav(struct runner *r, struct cell *ci,
- struct cell *cj, int gettimer) {
+void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
+ int gettimer) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -402,7 +359,8 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
const int gcount_j = cj->gcount;
/* Early abort? */
- if (gcount_i == 0 || gcount_j == 0) error("Empty cell !");
+ if (gcount_i == 0 || gcount_j == 0)
+ error("Doing pair gravity on an empty cell !");
/* Bad stuff will happen if cell sizes are different */
if (ci->width[0] != cj->width[0])
@@ -468,9 +426,9 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
} else {
- /* Ok, here we can go for particle-multipole interactions */
- runner_dopair_grav_pm(r, ci->progeny[j], cj->progeny[k]);
- runner_dopair_grav_pm(r, cj->progeny[k], ci->progeny[j]);
+ /* Ok, here we can go for multipole-multipole interactions */
+ runner_dopair_grav_mm(r, ci->progeny[j], cj->progeny[k]);
+ runner_dopair_grav_mm(r, cj->progeny[k], ci->progeny[j]);
}
}
}
@@ -494,12 +452,12 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
*
* @todo Use a local cache for the particles.
*/
-static void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
+void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
#ifdef SWIFT_DEBUG_CHECKS
/* Early abort? */
- if (c->gcount == 0) error("Empty cell !");
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
TIMER_TIC;
@@ -532,8 +490,8 @@ static void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
if (gettimer) TIMER_TOC(timer_dosub_self_grav);
}
-static void runner_dosub_grav(struct runner *r, struct cell *ci,
- struct cell *cj, int timer) {
+void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
+ int timer) {
/* Is this a single cell? */
if (cj == NULL) {
@@ -551,8 +509,7 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
}
}
-static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
- int timer) {
+void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
#if ICHECK > 0
for (int pid = 0; pid < ci->gcount; pid++) {
@@ -567,6 +524,8 @@ static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
}
#endif
+ TIMER_TIC;
+
/* Recover the list of top-level cells */
const struct engine *e = r->e;
struct cell *cells = e->s->cells_top;
@@ -579,6 +538,9 @@ static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
+ /* Drift our own multipole if need be */
+ if (ci->ti_old_multipole != e->ti_current) cell_drift_multipole(ci, e);
+
/* Loop over all the cells and go for a p-m interaction if far enough but not
* too far */
for (int i = 0; i < nr_cells; ++i) {
@@ -586,16 +548,18 @@ static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
struct cell *cj = &cells[i];
if (ci == cj) continue;
+ if (cj->gcount == 0) continue;
/* const double dx[3] = {cj->loc[0] - pos_i[0], // x */
/* cj->loc[1] - pos_i[1], // y */
/* cj->loc[2] - pos_i[2]}; // z */
/* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
-
- // if (r2 > max_d2) continue;
+ /* if (r2 > max_d2) continue; */
if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}
+
+ if (timer) TIMER_TOC(timer_dograv_long_range);
}
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/space.c b/src/space.c
index 5899cc396b84cd3aa92aab10a687caaa8d6d81e9..64a9ab15c960e7664afdf6be4293bbad3176fc76 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2100,10 +2100,10 @@ void space_split_recursive(struct space *s, struct cell *c,
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct gravity_tensors *m = c->progeny[k]->multipole;
- CoM[0] += m->CoM[0] * m->m_pole.mass;
- CoM[1] += m->CoM[1] * m->m_pole.mass;
- CoM[2] += m->CoM[2] * m->m_pole.mass;
- mass += m->m_pole.mass;
+ CoM[0] += m->CoM[0] * m->m_pole.M_000;
+ CoM[1] += m->CoM[1] * m->m_pole.M_000;
+ CoM[2] += m->CoM[2] * m->m_pole.M_000;
+ mass += m->m_pole.M_000;
}
}
c->multipole->CoM[0] = CoM[0] / mass;
diff --git a/src/space.h b/src/space.h
index 55ebdad91db4de7f7baa422b2e00b7dde438d8ca..73bd50da928c55890a91415f6e07c5100a7b71e7 100644
--- a/src/space.h
+++ b/src/space.h
@@ -76,9 +76,6 @@ struct space {
/*! Are we doing gravity? */
int gravity;
- /*! Total mass in the system */
- double total_mass;
-
/*! Width of the top-level cells. */
double width[3];
diff --git a/src/task.c b/src/task.c
index dfc3dc538c75297cbe7e79b7d64c1b7e13a015dc..0ac7bc0c59eb678383adba9587a8026cf872ee6d 100644
--- a/src/task.c
+++ b/src/task.c
@@ -276,8 +276,8 @@ float task_overlap(const struct task *restrict ta,
*/
void task_unlock(struct task *t) {
- const int type = t->type;
- const int subtype = t->subtype;
+ const enum task_types type = t->type;
+ const enum task_subtypes subtype = t->subtype;
struct cell *ci = t->ci, *cj = t->cj;
/* Act based on task type. */
@@ -300,6 +300,7 @@ void task_unlock(struct task *t) {
case task_type_sub_self:
if (subtype == task_subtype_grav) {
cell_gunlocktree(ci);
+ cell_munlocktree(ci);
} else {
cell_unlocktree(ci);
}
@@ -310,15 +311,25 @@ void task_unlock(struct task *t) {
if (subtype == task_subtype_grav) {
cell_gunlocktree(ci);
cell_gunlocktree(cj);
+ cell_munlocktree(ci);
+ cell_munlocktree(cj);
} else {
cell_unlocktree(ci);
cell_unlocktree(cj);
}
break;
- case task_type_grav_mm:
+ case task_type_grav_down:
cell_gunlocktree(ci);
+ cell_munlocktree(ci);
+ break;
+
+ case task_type_grav_top_level:
+ case task_type_grav_long_range:
+ case task_type_grav_mm:
+ cell_munlocktree(ci);
break;
+
default:
break;
}
@@ -331,8 +342,8 @@ void task_unlock(struct task *t) {
*/
int task_lock(struct task *t) {
- const int type = t->type;
- const int subtype = t->subtype;
+ const enum task_types type = t->type;
+ const enum task_subtypes subtype = t->subtype;
struct cell *ci = t->ci, *cj = t->cj;
#ifdef WITH_MPI
int res = 0, err = 0;
@@ -379,7 +390,14 @@ int task_lock(struct task *t) {
case task_type_self:
case task_type_sub_self:
if (subtype == task_subtype_grav) {
- if (cell_glocktree(ci) != 0) return 0;
+ /* Lock the gparts and the m-pole */
+ if (ci->ghold || ci->mhold) return 0;
+ if (cell_glocktree(ci) != 0)
+ return 0;
+ else if (cell_mlocktree(ci) != 0) {
+ cell_gunlocktree(ci);
+ return 0;
+ }
} else {
if (cell_locktree(ci) != 0) return 0;
}
@@ -388,13 +406,24 @@ int task_lock(struct task *t) {
case task_type_pair:
case task_type_sub_pair:
if (subtype == task_subtype_grav) {
+ /* Lock the gparts and the m-pole in both cells */
if (ci->ghold || cj->ghold) return 0;
if (cell_glocktree(ci) != 0) return 0;
if (cell_glocktree(cj) != 0) {
cell_gunlocktree(ci);
return 0;
+ } else if (cell_mlocktree(ci) != 0) {
+ cell_gunlocktree(ci);
+ cell_gunlocktree(cj);
+ return 0;
+ } else if (cell_mlocktree(cj) != 0) {
+ cell_gunlocktree(ci);
+ cell_gunlocktree(cj);
+ cell_munlocktree(ci);
+ return 0;
}
} else {
+ /* Lock the parts in both cells */
if (ci->hold || cj->hold) return 0;
if (cell_locktree(ci) != 0) return 0;
if (cell_locktree(cj) != 0) {
@@ -404,8 +433,23 @@ int task_lock(struct task *t) {
}
break;
+ case task_type_grav_down:
+ /* Lock the gparts and the m-poles */
+ if (ci->ghold || ci->mhold) return 0;
+ if (cell_glocktree(ci) != 0)
+ return 0;
+ else if (cell_mlocktree(ci) != 0) {
+ cell_gunlocktree(ci);
+ return 0;
+ }
+ break;
+
+ case task_type_grav_top_level:
+ case task_type_grav_long_range:
case task_type_grav_mm:
- cell_glocktree(ci);
+ /* Lock the m-poles */
+ if (ci->mhold) return 0;
+ if (cell_mlocktree(ci) != 0) return 0;
break;
default:
diff --git a/src/timers.h b/src/timers.h
index 4cb4d7e0ba60003ba4caefffe257c929c59a8d9e..39bcf30fba0ec36d9209ddcbf3c71035a5851dbb 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -49,6 +49,8 @@ enum {
timer_dopair_grav_mm,
timer_dopair_grav_pp,
timer_dograv_external,
+ timer_dograv_down,
+ timer_dograv_long_range,
timer_dosource,
timer_dosub_self_density,
timer_dosub_self_gradient,
diff --git a/src/vector_power.h b/src/vector_power.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e7c197f8db74f5865ba4d55a79c0e0abaab9baf
--- /dev/null
+++ b/src/vector_power.h
@@ -0,0 +1,413 @@
+/*******************************************************************************
+ * 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_VECTOR_POWER_H
+#define SWIFT_VECTOR_POWER_H
+
+/**
+ * @file vector_power.h
+ * @brief Powers of 3D vectors to a multi-index with factorial.
+ *
+ * These expressions are to be used in 3D Taylor series.
+ *
+ * We use the notation of Dehnen, Computational Astrophysics and Cosmology,
+ * 1, 1, pp. 24 (2014), arXiv:1405.2255.
+ *
+ * We compute \f$ \frac{1}{\vec{m}!}\vec{v}^{\vec{m}} \f$ for all relevant m.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/***************************/
+/* 0th order vector powers */
+/***************************/
+
+/**
+ * @brief \f$ \frac{1}{(0,0,0)!}\vec{v}^{(0,0,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_000(const double v[3]) {
+
+ return 1.;
+}
+
+/***************************/
+/* 1st order vector powers */
+/***************************/
+
+/**
+ * @brief \f$ \frac{1}{(1,0,0)!}\vec{v}^{(1,0,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_100(const double v[3]) {
+
+ return v[0];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,1,0)!}\vec{v}^{(0,1,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_010(const double v[3]) {
+
+ return v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,0,1)!}\vec{v}^{(0,0,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_001(const double v[3]) {
+
+ return v[2];
+}
+
+/***************************/
+/* 2nd order vector powers */
+/***************************/
+
+/**
+ * @brief \f$ \frac{1}{(2,0,0)!}\vec{v}^{(2,0,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_200(const double v[3]) {
+
+ return 0.5 * v[0] * v[0];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,2,0)!}\vec{v}^{(0,2,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_020(const double v[3]) {
+
+ return 0.5 * v[1] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,0,2)!}\vec{v}^{(0,0,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_002(const double v[3]) {
+
+ return 0.5 * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,1,0)!}\vec{v}^{(1,1,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_110(const double v[3]) {
+
+ return v[0] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,0,1)!}\vec{v}^{(1,0,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_101(const double v[3]) {
+
+ return v[0] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,1,1)!}\vec{v}^{(0,1,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_011(const double v[3]) {
+
+ return v[1] * v[2];
+}
+
+/***************************/
+/* 3rd order vector powers */
+/***************************/
+
+/**
+ * @brief \f$ \frac{1}{(3,0,0)!}\vec{v}^{(3,0,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_300(const double v[3]) {
+
+ return 0.1666666666666667 * v[0] * v[0] * v[0];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,3,0)!}\vec{v}^{(0,3,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_030(const double v[3]) {
+
+ return 0.1666666666666667 * v[1] * v[1] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,0,3)!}\vec{v}^{(0,0,3)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_003(const double v[3]) {
+
+ return 0.1666666666666667 * v[2] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(2,1,0)!}\vec{v}^{(2,1,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_210(const double v[3]) {
+
+ return 0.5 * v[0] * v[0] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(2,0,1)!}\vec{v}^{(2,0,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_201(const double v[3]) {
+
+ return 0.5 * v[0] * v[0] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,2,0)!}\vec{v}^{(1,2,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_120(const double v[3]) {
+
+ return 0.5 * v[0] * v[1] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,2,1)!}\vec{v}^{(0,2,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_021(const double v[3]) {
+
+ return 0.5 * v[1] * v[1] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,0,2)!}\vec{v}^{(1,0,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_102(const double v[3]) {
+
+ return 0.5 * v[0] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,1,2)!}\vec{v}^{(0,1,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_012(const double v[3]) {
+
+ return 0.5 * v[1] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,1,1)!}\vec{v}^{(1,1,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_111(const double v[3]) {
+
+ return v[0] * v[1] * v[2];
+}
+
+/***************************/
+/* 4th order vector powers */
+/***************************/
+
+/**
+ * @brief \f$ \frac{1}{(4,0,0)!}\vec{v}^{(4,0,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_400(const double v[3]) {
+
+ const double vv = v[0] * v[0];
+ return 0.041666666666666667 * vv * vv;
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,4,0)!}\vec{v}^{(0,4,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_040(const double v[3]) {
+
+ const double vv = v[1] * v[1];
+ return 0.041666666666666667 * vv * vv;
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,0,4)!}\vec{v}^{(0,0,4)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_004(const double v[3]) {
+
+ const double vv = v[2] * v[2];
+ return 0.041666666666666667 * vv * vv;
+}
+
+/**
+ * @brief \f$ \frac{1}{(3,1,0)!}\vec{v}^{(3,1,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_310(const double v[3]) {
+
+ return 0.1666666666666667 * v[0] * v[0] * v[0] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(3,0,1)!}\vec{v}^{(3,0,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_301(const double v[3]) {
+
+ return 0.1666666666666667 * v[0] * v[0] * v[0] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,3,0)!}\vec{v}^{(1,3,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_130(const double v[3]) {
+
+ return 0.1666666666666667 * v[0] * v[1] * v[1] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,3,1)!}\vec{v}^{(0,3,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_031(const double v[3]) {
+
+ return 0.1666666666666667 * v[1] * v[1] * v[1] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,0,3)!}\vec{v}^{(1,0,3)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_103(const double v[3]) {
+
+ return 0.1666666666666667 * v[0] * v[2] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,1,3)!}\vec{v}^{(0,1,3)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_013(const double v[3]) {
+
+ return 0.1666666666666667 * v[1] * v[2] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(2,2,0)!}\vec{v}^{(2,2,0)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_220(const double v[3]) {
+
+ return 0.25 * v[0] * v[0] * v[1] * v[1];
+}
+
+/**
+ * @brief \f$ \frac{1}{(2,0,2)!}\vec{v}^{(2,0,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_202(const double v[3]) {
+
+ return 0.25 * v[0] * v[0] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(0,2,2)!}\vec{v}^{(0,2,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_022(const double v[3]) {
+
+ return 0.25 * v[1] * v[1] * v[2] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(2,1,1)!}\vec{v}^{(2,1,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_211(const double v[3]) {
+
+ return 0.5 * v[0] * v[0] * v[1] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,2,1)!}\vec{v}^{(1,2,1)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_121(const double v[3]) {
+
+ return 0.5 * v[0] * v[1] * v[1] * v[2];
+}
+
+/**
+ * @brief \f$ \frac{1}{(1,1,2)!}\vec{v}^{(1,1,2)} \f$.
+ *
+ * @param v vector (\f$ v \f$).
+ */
+__attribute__((always_inline)) INLINE static double X_112(const double v[3]) {
+
+ return 0.5 * v[0] * v[1] * v[2] * v[2];
+}
+
+#endif /* SWIFT_VECTOR_POWER_H */
diff --git a/tests/testKernelGrav.c b/tests/testKernelGrav.c
index 41e085945693efaf658c219d0e5f992ddf023d74..b4a5e4d9f1ff05d8f34840dd19b2a2ccb9ec79b5 100644
--- a/tests/testKernelGrav.c
+++ b/tests/testKernelGrav.c
@@ -31,26 +31,31 @@
/**
* @brief The Gadget-2 gravity kernel function
*
+ * Taken from Gadget-2.0.7's forcetree.c lines 2755-2800
+ *
* @param r The distance between particles
- * @param h The cut-off distance of the kernel
+ * @param epsilon The cut-off distance of the kernel
*/
-float gadget(float r, float h) {
- float fac;
- const float r2 = r * r;
- if (r >= h)
- fac = 1.0f / (r2 * r);
- else {
- const float h_inv = 1. / h;
- const float h_inv3 = h_inv * h_inv * h_inv;
- const float u = r * h_inv;
+float gadget(float r, float epsilon) {
+
+ const float h = epsilon;
+ const float h_inv = 1.f / h;
+
+ const float u = r * h_inv;
+
+ if (u >= 1) {
+ const float r_inv = 1. / r;
+
+ return r_inv * r_inv * r_inv;
+ } else {
if (u < 0.5)
- fac = h_inv3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
+ return h_inv * h_inv * h_inv *
+ (10.666666666667 + u * u * (32.0 * u - 38.4));
else
- fac =
- h_inv3 * (21.333333333333 - 48.0 * u + 38.4 * u * u -
- 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
+ return h_inv * h_inv * h_inv *
+ (21.333333333333 - 48.0 * u + 38.4 * u * u -
+ 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
- return fac;
}
int main() {
@@ -61,20 +66,22 @@ int main() {
for (int k = 1; k < numPoints; ++k) {
const float r = (r_max * k) / numPoints;
-
- const float u = r / h;
-
const float gadget_w = gadget(r, h);
+ const float h_inv = 1.f / h;
+ const float h_inv3 = h_inv * h_inv * h_inv;
+ const float u = r * h_inv;
+
float swift_w;
- if (u < 1.) {
- kernel_grav_eval(u, &swift_w);
- swift_w *= (1 / (h * h * h));
- } else {
+ if (r >= h) {
swift_w = 1 / (r * r * r);
+
+ } else {
+ kernel_grav_eval(u, &swift_w);
+ swift_w *= h_inv3;
}
- if (fabsf(gadget_w - swift_w) > 2e-7) {
+ if (fabsf(gadget_w - swift_w) > 1e-5 * fabsf(gadget_w)) {
printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
gadget_w, swift_w);
@@ -87,28 +94,30 @@ int main() {
printf("\nAll values are consistent\n");
/* Now test the long range function */
- const float a_smooth = 4.5f;
+ /* const float a_smooth = 4.5f; */
- for (int k = 1; k < numPoints; ++k) {
+ /* for (int k = 1; k < numPoints; ++k) { */
- const float r = (r_max * k) / numPoints;
+ /* const float r = (r_max * k) / numPoints; */
- const float u = r / a_smooth;
+ /* const float u = r / a_smooth; */
- float swift_w;
- kernel_long_grav_eval(u, &swift_w);
+ /* float swift_w; */
+ /* kernel_long_grav_eval(u, &swift_w); */
- float gadget_w = erfc(u / 2) + u * exp(-u * u / 4) / sqrt(M_PI);
+ /* float gadget_w = erfcf(u / 2) + u * expf(-u * u / 4) / sqrtf(M_PI); */
- if (fabsf(gadget_w - swift_w) > 2e-7) {
+ /* if (fabsf(gadget_w - swift_w) > 1e-4 * fabsf(gadget_w)) { */
- printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth,
- u, swift_w, gadget_w);
+ /* printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r,
+ * a_smooth, */
+ /* u, swift_w, gadget_w); */
- printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
- return 1;
- }
- }
+ /* printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
+ */
+ /* return 1; */
+ /* } */
+ /* } */
return 0;
}
diff --git a/theory/Multipoles/multipoles.py b/theory/Multipoles/multipoles.py
new file mode 100644
index 0000000000000000000000000000000000000000..06a886e70a7b48bf67d9b88ef28cb486f188a853
--- /dev/null
+++ b/theory/Multipoles/multipoles.py
@@ -0,0 +1,179 @@
+import numpy as np
+import sys
+
+def factorial(x):
+ if x == 0:
+ return 1
+ else:
+ return x * factorial(x-1)
+
+SUFFIXES = {1: 'st', 2: 'nd', 3: 'rd'}
+def ordinal(num):
+ suffix = SUFFIXES.get(num % 10, 'th')
+ return str(num) + suffix
+
+# Get the order
+order = int(sys.argv[1])
+
+print "-------------------------------------------------"
+print "Generating code for multipoles of order", order, "(only)."
+print "-------------------------------------------------\n"
+
+print "-------------------------------------------------"
+print "Multipole structure:"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d\n"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "float M_%d%d%d;"%(i,j,k)
+
+if order > 0:
+ print "#endif"
+
+print ""
+print "-------------------------------------------------"
+
+print "Field tensor structure:"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d\n"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "float F_%d%d%d;"%(i,j,k)
+if order > 0:
+ print "#endif"
+
+print ""
+print "-------------------------------------------------"
+
+print "gravity_field_tensors_add():"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "la->F_%d%d%d += lb->F_%d%d%d;"%(i,j,k,i,j,k)
+if order > 0:
+ print "#endif"
+
+print ""
+print "-------------------------------------------------"
+
+print "gravity_multipole_add():"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "ma->M_%d%d%d += mb->M_%d%d%d;"%(i,j,k,i,j,k)
+if order > 0:
+ print "#endif"
+
+print ""
+print "-------------------------------------------------"
+
+print "gravity_P2M(): (loop)"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ if order % 2 == 0:
+ print "M_%d%d%d += m * X_%d%d%d(dx);"%(i,j,k,i,j,k)
+ else:
+ print "M_%d%d%d += -m * X_%d%d%d(dx);"%(i,j,k,i,j,k)
+if order > 0:
+ print "#endif"
+
+print ""
+print "-------------------------------------------------"
+
+print "gravity_P2M(): (storing)"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "m->m_pole.M_%d%d%d = M_%d%d%d;"%(i,j,k,i,j,k)
+if order > 0:
+ print "#endif"
+
+
+print ""
+print "-------------------------------------------------"
+
+print "gravity_M2M():"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* Shift %s order terms */"%ordinal(order)
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "m_a->M_%d%d%d = m_b->M_%d%d%d"%(i,j,k,i,j,k),
+
+ for ii in range(order+1):
+ for jj in range(order+1):
+ for kk in range(order+1):
+
+ if not(ii == 0 and jj == 0 and kk == 0):
+ for iii in range(order+1):
+ for jjj in range(order+1):
+ for kkk in range(order+1):
+ if ii+iii == i and jj + jjj == j and kk+kkk == k:
+ print "+ X_%d%d%d(dx) * m_b->M_%d%d%d"%(ii, jj, kk, iii, jjj, kkk),
+
+
+ print ";"
+if order > 0:
+ print "#endif"
+
+
diff --git a/theory/Multipoles/vector_powers.py b/theory/Multipoles/vector_powers.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b1cbd4436f00ff5a1720badd1e8c22bd1762511
--- /dev/null
+++ b/theory/Multipoles/vector_powers.py
@@ -0,0 +1,55 @@
+import numpy as np
+import sys
+
+def factorial(x):
+ if x == 0:
+ return 1
+ else:
+ return x * factorial(x-1)
+
+SUFFIXES = {1: 'st', 2: 'nd', 3: 'rd'}
+def ordinal(num):
+ suffix = SUFFIXES.get(num % 10, 'th')
+ return str(num) + suffix
+
+# Get the order
+order = int(sys.argv[1])
+
+print "-------------------------------------------------"
+print "Generating code for vector powers of order", order, "(only)."
+print "-------------------------------------------------\n"
+
+print "/***************************/"
+print "/* %s order vector powers */"%ordinal(order)
+print "/***************************/\n"
+
+# Create all the terms relevent for this order
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ fact = factorial(i) * factorial(j) * factorial(k)
+ print "/**"
+ print "* @brief \\f$ \\frac{1}{(%d,%d,%d)!}\\vec{v}^{(%d,%d,%d)} \\f$."%(i,j,k,i,j,k)
+ print "*"
+ print "* Note \\f$ \\frac{1}{(%d,%d,%d)!} = 1/(%d!*%d!*%d!) = 1/%d! = %e"%(i,j,k,i,j,k, fact, 1./fact)
+ print "*"
+ print "* @param v vector (\\f$ v \\f$)."
+ print "*/"
+ print "__attribute__((always_inline)) INLINE static double X_%d%d%d(const double v[3]) {"%(i,j,k)
+ print ""
+ print " return",
+ if fact != 1:
+ print "%12.15e"%(1./fact),
+ else:
+ print "1.",
+ for ii in range(i):
+ print "* v[0]",
+ for jj in range(j):
+ print "* v[1]",
+ for kk in range(k):
+ print "* v[2]",
+ print ";"
+ print "}"
+
+