Commit 1dec55c7 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Add a glass making master option

parent 91bb16c2
......@@ -270,6 +270,18 @@ elif test "$gravity_force_checks" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
fi
# Check whether we want to switch on glass making
AC_ARG_ENABLE([glass-making],
[AS_HELP_STRING([--enable-glass-making],
[Activate the glass-making procedure by reversing the sign of gravity @<:@yes/no@:>@]
)],
[gravity_glass_making="$enableval"],
[gravity_glass_making="no"]
)
if test "$gravity_glass_making" == "yes"; then
AC_DEFINE([SWIFT_MAKE_GRAVITY_GLASS], 1, [Make the code run in a way to produce a glass file for gravity/cosmology])
fi
# Check if we want to zero the gravity forces for all particles below some ID.
AC_ARG_ENABLE([no-gravity-below-id],
[AS_HELP_STRING([--enable-no-gravity-below-id],
......@@ -1408,6 +1420,7 @@ AC_MSG_RESULT([
Gravity scheme : $with_gravity
Multipole order : $with_multipole_order
No gravity below ID : $no_gravity_below_id
Make gravity glass : $gravity_glass_making
External potential : $with_potential
Cooling function : $with_cooling
......
......@@ -23,3 +23,4 @@ and keep on your desk.
parameter_file
what_about_mpi
running_on_large_systems
special_modes
......@@ -19,38 +19,51 @@ In the sections ``Snapshots`` and ``Statistics``, you can specify the options ``
The ``output_list_on`` enable or not the output list and ``output_list`` is the filename containing the output times.
With the file header, you can choose between writing redshifts, scale factors or times.
Example of file containing with times (in internal units):
::
# Time
0.5
1.5
3.0
12.5
Example of file with scale factors:
::
# Scale Factor
0.1
0.2
0.3
Example of file with redshift:
::
# Redshift
20
15
10
5
Example of file containing with times (in internal units)::
# Time
0.5
1.5
3.0
12.5
Example of file with scale factors::
# Scale Factor
0.1
0.2
0.3
Example of file with redshift::
# Redshift
20
15
10
5
Output Selection
~~~~~~~~~~~~~~~~
With SWIFT, you can select the particle fields to output in snapshot using the parameter file.
In section ``SelectOutput``, you can remove a field by adding a parameter formatted in the
following way ``field_parttype`` where ``field`` is the name of the field that you
want to remove (e.g. ``Masses``) and ``parttype`` is the type of particles that
contains this field (e.g. ``Gas``, ``DM`` or ``Star``). For a parameter, the only
values accepted are 0 (skip this field when writing) or 1 (default, do not skip
this field when writing).
With SWIFT, you can select the particle fields to output in snapshot
using the parameter file. In section ``SelectOutput``, you can remove
a field by adding a parameter formatted in the following way
``field_parttype`` where ``field`` is the name of the field that you
want to remove (e.g. ``Masses``) and ``parttype`` is the type of
particles that contains this field (e.g. ``Gas``, ``DM`` or ``Star``).
For a parameter, the only values accepted are 0 (skip this field when
writing) or 1 (default, do not skip this field when writing). By
default all fields are written.
This field is mostly used to remove unnecessary output by listing them
with 0's. A classic use-case for this feature is a DM-only simulation
(pure n-body) where all particles have the same mass. Outputing the
mass field in the snapshots results in extra i/o time and unnecessary
waste of disk space. The corresponding section of the ``yaml``
parameter file would look like this::
SelectOutput:
Masses_DM: 0
You can generate a ``yaml`` file containing all the possible fields with ``./swift -o output.yml``. By default, all the fields are written.
You can generate a ``yaml`` file containing all the possible fields
available for a given configuration of SWIFT by running ``./swift -o output.yml``.
.. Special modes
Matthieu Schaller, 20/08/2018
Special modes
=============
SWIFT comes with a few special modes of operating to perform additional tasks.
Static particles
~~~~~~~~~~~~~~~~
For some test problems it is convenient to have a set of particles that do not
perceive any gravitational forces and just act as sources for the force
calculation. This can be achieved by configuring SWIFT with the option
``--enable-no-gravity-below-id=N``. This will zero the *accelerations* of all
particles with ``id`` (strictly) lower than ``N`` at every time-step. Note that
if these particles have an initial velocity they will keep moving at that
speed.
This will also naturally set their time-step to the maximal value
(``TimeIntegration:dt_max``) set in the parameter file.
A typical use-case for this feature is to study the evolution of one particle
orbiting a static halo made of particles. This can be used to assess the
quality of the gravity tree and time integration. As more particles are added
to the halo, the orbits will get closer to the analytic solution as the noise
in the sampling of the halo is reduced.
Note also that this does not affect the hydrodynamic forces. This mode is
purely designed for gravity-only accuracy tests.
Gravity glasses
~~~~~~~~~~~~~~~
For many problems in cosmology, it is important to start a simulation with no
initial noise in the particle distribution. Such a "glass" can be created by
starting from a random distribution of particles and running with the sign of
gravity reversed until all the particles reach a steady state. To run SWIFT in
this mode, configure the code with ``--enable-glass-making``.
Note that this will *not* generate the initial random distribution of
particles. An initial condition file with particles still has to be provided.
This example can be used to generate a glass file for gravity calculation.
The makeIC.py script will generate a uniform Poisson distribution of particles
in a cubic box with zero initial velocities.
By running the code with the SWIFT configuration option --enable-glass-making
the code will run with gravity as a repulsive force and the particles will
moe towards a state of minimal energy. These glass files can be used to then
start simulations with a minimal level of noise.
......@@ -20,7 +20,7 @@
import h5py
import sys
from numpy import *
import numpy as np
# Generates a swift IC file containing a cartesian distribution of DM particles
# with a density of 1
......@@ -30,7 +30,7 @@ periodic= 1 # 1 For periodic box
boxSize = 1.
rho = 1.
L = int(sys.argv[1]) # Number of particles along one axis
fileName = "uniformDMBox_%d.hdf5"%L
fileName = "uniform_DM_box.hdf5"
#---------------------------------------------------
numPart = L**3
......@@ -69,27 +69,16 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType1")
v = zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = zeros(1)
v = np.zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f', data=v)
m = full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
m = zeros(1)
m = np.full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f', data=m)
ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ids = np.linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
coords = zeros((numPart, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
coords = np.random.rand(numPart, 3) * boxSize
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd', data=coords)
file.close()
......@@ -6,33 +6,39 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Let's overwrite G to make this more effective
PhysicalConstants:
G: 1.
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 100. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1. # The maximal time-step size of the simulation (in internal units).
time_begin: 0.
time_end: 100.
dt_min: 1e-6
dt_max: 1.
Scheduler:
max_top_level_cells: 8
cell_split_size: 50
# Parameters governing the snapshots
Snapshots:
basename: uniformDMBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 10. # Time difference between consecutive outputs (in internal units)
basename: uniform_DM_box
time_first: 0.
delta_time: 1.
compression: 4
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.8 # Opening angle (Multipole acceptance criterion)
epsilon: 0.01 # Softening length (in internal units).
eta: 0.025
theta: 0.3
mesh_side_length: 32
comoving_softening: 0.001
max_physical_softening: 0.001
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 5. # Time between statistics output
delta_time: 0.1
# Parameters related to the initial conditions
InitialConditions:
file_name: ./uniformDMBox_16.hdf5 # The file to read
file_name: ./uniform_DM_box.hdf5
#!/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)
......@@ -1757,6 +1757,14 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Finish the force calculation */
gravity_end_force(gp, G_newton, potential_normalisation, periodic);
#ifdef SWIFT_MAKE_GRAVITY_GLASS
/* Negate the gravity forces */
gp->a_grav[0] *= -1.f;
gp->a_grav[1] *= -1.f;
gp->a_grav[2] *= -1.f;
#endif
#ifdef SWIFT_NO_GRAVITY_BELOW_ID
/* Get the ID of the gpart */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment