Commit 0ee1775e authored by James Willis's avatar James Willis
Browse files

Merge branch 'dopair-vectorisation-merge' of...

Merge branch 'dopair-vectorisation-merge' of gitlab.cosma.dur.ac.uk:swift/swiftsim into dopair-vectorisation-merge
parents 5556f231 2ede7d27
...@@ -73,6 +73,9 @@ tests/testRiemannExact ...@@ -73,6 +73,9 @@ tests/testRiemannExact
tests/testRiemannTRRS tests/testRiemannTRRS
tests/testRiemannHLLC tests/testRiemannHLLC
tests/testMatrixInversion tests/testMatrixInversion
tests/testVoronoi1D
tests/testVoronoi2D
tests/testVoronoi3D
tests/testDump tests/testDump
tests/testLogger tests/testLogger
tests/benchmarkInteractions tests/benchmarkInteractions
......
...@@ -587,7 +587,7 @@ fi ...@@ -587,7 +587,7 @@ fi
# Hydro scheme. # Hydro scheme.
AC_ARG_WITH([hydro], AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>], [AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo default: gadget2@:>@] [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax default: gadget2@:>@]
)], )],
[with_hydro="$withval"], [with_hydro="$withval"],
[with_hydro="gadget2"] [with_hydro="gadget2"]
...@@ -608,6 +608,9 @@ case "$with_hydro" in ...@@ -608,6 +608,9 @@ case "$with_hydro" in
gizmo) gizmo)
AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH]) AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH])
;; ;;
shadowfax)
AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
;;
*) *)
AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro]) AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
...@@ -796,6 +799,15 @@ case "$with_potential" in ...@@ -796,6 +799,15 @@ case "$with_potential" in
;; ;;
esac 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. # Check for git, needed for revision stamps.
...@@ -843,6 +855,7 @@ AC_MSG_RESULT([ ...@@ -843,6 +855,7 @@ AC_MSG_RESULT([
Riemann solver : $with_riemann Riemann solver : $with_riemann
Cooling function : $with_cooling Cooling function : $with_cooling
External potential : $with_potential External potential : $with_potential
Multipole order : $with_multipole_order
Task debugging : $enable_task_debugging Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks Debugging checks : $enable_debugging_checks
......
...@@ -13,6 +13,9 @@ TimeIntegration: ...@@ -13,6 +13,9 @@ TimeIntegration:
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). 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). dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
Scheduler:
cell_split_size: 50
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
basename: eagle # Common part of the name of output files basename: eagle # Common part of the name of output files
...@@ -23,6 +26,13 @@ Snapshots: ...@@ -23,6 +26,13 @@ Snapshots:
Statistics: Statistics:
delta_time: 1e-2 # Time between statistics output 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 # Parameters for the hydrodynamics scheme
SPH: SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
......
...@@ -23,6 +23,13 @@ Snapshots: ...@@ -23,6 +23,13 @@ Snapshots:
Statistics: Statistics:
delta_time: 1e-2 # Time between statistics output 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 # Parameters for the hydrodynamics scheme
SPH: SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
......
#!/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)
...@@ -323,14 +323,14 @@ int main(int argc, char *argv[]) { ...@@ -323,14 +323,14 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */ /* How large are the parts? */
if (myrank == 0) { if (myrank == 0) {
message("sizeof(part) is %4zi bytes.", sizeof(struct part)); message("sizeof(part) is %4zi bytes.", sizeof(struct part));
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart)); message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart)); message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart)); message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole)); message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor)); message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
message("sizeof(task) is %4zi bytes.", sizeof(struct task)); message("sizeof(task) is %4zi bytes.", sizeof(struct task));
message("sizeof(cell) is %4zi bytes.", sizeof(struct cell)); message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
} }
/* Read the parameter file */ /* Read the parameter file */
......
...@@ -57,14 +57,15 @@ pl.rcParams.update(PLOT_PARAMS) ...@@ -57,14 +57,15 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h. # Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init", "ghost", "extra_ghost", "drift", "kick1", "kick2", "init", "ghost", "extra_ghost", "drift", "kick1", "kick2",
"timestep", "send", "recv", "grav_gather_m", "grav_fft", "timestep", "send", "recv", "grav_top_level", "grav_long_range",
"grav_mm", "grav_up", "cooling", "sourceterms", "count"] "grav_mm", "grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav", SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "count"] "tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
# Task/subtypes of interest. # Task/subtypes of interest.
FULLTYPES = ["self/force", "self/density", "sub_self/force", FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
"sub_self/density", "pair/force", "pair/density", "sub_pair/force", "sub_self/density", "pair/force", "pair/density", "pair/grav",
"sub_pair/force",
"sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho", "sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
"recv/tend", "send/tend"] "recv/tend", "send/tend"]
......
...@@ -45,7 +45,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ ...@@ -45,7 +45,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \ physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.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 \ sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h sort.h xmf.h gravity_properties.h gravity_derivatives.h dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
vector_power.h hydro_space.h sort_part.h
# Common source files # Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
...@@ -55,7 +56,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ ...@@ -55,7 +56,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \ physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \ runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \ statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c part_type.c xmf.c gravity_properties.c gravity.c \
hydro_space.c
# Include files for distribution, not installation. # Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
...@@ -75,8 +77,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h ...@@ -75,8 +77,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \ hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \ hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \ hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \ hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h \
hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \ hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h \
hydro/Gizmo/hydro_part.h \
hydro/Gizmo/hydro_gradients_gizmo.h \
hydro/Gizmo/hydro_gradients.h \
hydro/Gizmo/hydro_gradients_sph.h \
hydro/Gizmo/hydro_slope_limiters_cell.h \
hydro/Gizmo/hydro_slope_limiters_face.h \
hydro/Gizmo/hydro_slope_limiters.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
hydro/Shadowswift/hydro_io.h \
hydro/Shadowswift/hydro_part.h \
hydro/Shadowswift/hydro_slope_limiters_cell.h \
hydro/Shadowswift/hydro_slope_limiters_face.h \
hydro/Shadowswift/hydro_slope_limiters.h \
hydro/Shadowswift/voronoi1d_algorithm.h \
hydro/Shadowswift/voronoi1d_cell.h \
hydro/Shadowswift/voronoi2d_algorithm.h \
hydro/Shadowswift/voronoi2d_cell.h \
hydro/Shadowswift/voronoi3d_algorithm.h \
hydro/Shadowswift/voronoi3d_cell.h \
hydro/Shadowswift/voronoi_algorithm.h \
hydro/Shadowswift/voronoi_cell.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \ riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \ riemann/riemann_exact.h riemann/riemann_vacuum.h \
stars.h stars_io.h \ stars.h stars_io.h \
......
...@@ -26,7 +26,7 @@ ...@@ -26,7 +26,7 @@
#include "cell.h" #include "cell.h"
#include "error.h" #include "error.h"
#include "part.h" #include "part.h"
#include "sort.h" #include "sort_part.h"
#include "vector.h" #include "vector.h"
#define NUM_VEC_PROC 2 #define NUM_VEC_PROC 2
......
...@@ -1114,7 +1114,7 @@ void cell_check_multipole(struct cell *c, void *data) { ...@@ -1114,7 +1114,7 @@ void cell_check_multipole(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
struct gravity_tensors ma; struct gravity_tensors ma;
const double tolerance = 1e-5; /* Relative */ const double tolerance = 1e-3; /* Relative */
/* First recurse */ /* First recurse */
if (c->split) if (c->split)
...@@ -1127,8 +1127,7 @@ void cell_check_multipole(struct cell *c, void *data) { ...@@ -1127,8 +1127,7 @@ void cell_check_multipole(struct cell *c, void *data) {
gravity_P2M(&ma, c->gparts, c->gcount); gravity_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */ /* Now compare the multipole expansion */
if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole, if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth); message("Multipoles are not equal at depth=%d!", c->depth);
message("Correct answer:"); message("Correct answer:");
gravity_multipole_print(&ma.m_pole); gravity_multipole_print(&ma.m_pole);
...@@ -1487,7 +1486,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) { ...@@ -1487,7 +1486,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
} else if (ti_current > ti_old_multipole) { } else if (ti_current > ti_old_multipole) {
/* Drift the multipole */ /* Drift the multipole */
gravity_multipole_drift(c->multipole, dt); gravity_drift(c->multipole, dt);
} }
/* Update the time of the last drift */ /* Update the time of the last drift */
...@@ -1504,7 +1503,21 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) { ...@@ -1504,7 +1503,21 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
* @param e The #engine (to get ti_current). * @param e The #engine (to get ti_current).
*/ */
void cell_drift_multipole(struct cell *c, const struct engine *e) { 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;
} }
/** /**
......
...@@ -39,9 +39,6 @@ ...@@ -39,9 +39,6 @@
/* Thermal energy per unit mass used as a constant for the isothermal EoS */ /* Thermal energy per unit mass used as a constant for the isothermal EoS */
#define const_isothermal_internal_energy 20.2615290634f #define const_isothermal_internal_energy 20.2615290634f
/* Self gravity stuff. */
#define const_gravity_multipole_order 1
/* Type of gradients to use (GIZMO_SPH only) */ /* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */ /* If no option is chosen, no gradients are used (first order scheme) */
//#define GRADIENTS_SPH //#define GRADIENTS_SPH
...@@ -57,6 +54,23 @@ ...@@ -57,6 +54,23 @@
//#define GIZMO_FIX_PARTICLES //#define GIZMO_FIX_PARTICLES
//#define GIZMO_TOTAL_ENERGY //#define GIZMO_TOTAL_ENERGY
/* Types of gradients to use for SHADOWFAX_SPH */
/* If no option is chosen, no gradients are used (first order scheme) */
#define SHADOWFAX_GRADIENTS
/* SHADOWFAX_SPH slope limiters */
#define SHADOWFAX_SLOPE_LIMITER_PER_FACE
#define SHADOWFAX_SLOPE_LIMITER_CELL_WIDE
/* Options to control SHADOWFAX_SPH */
/* This option disables cell movement */
//#define SHADOWFAX_FIX_CELLS
/* This option enables cell steering, i.e. trying to keep the cells regular by
adding a correction to the cell velocities.*/
#define SHADOWFAX_STEER_CELL_MOTION
/* This option evolves the total energy instead of the thermal energy */
//#define SHADOWFAX_TOTAL_ENERGY
/* Source terms */ /* Source terms */
#define SOURCETERMS_NONE #define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK //#define SOURCETERMS_SN_FEEDBACK
......
...@@ -49,6 +49,8 @@ ...@@ -49,6 +49,8 @@
#include "./hydro/Default/hydro_debug.h" #include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_SPH) #elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_debug.h" #include "./hydro/Gizmo/hydro_debug.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_debug.h"
#else