Skip to content
GitLab
Commit fa3be81e authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge remote-tracking branch 'origin/master' into new_time_line 

Conflicts:
	src/runner.c
parents c6fcd48a fb49f1ee
Hide whitespace changes
Inline Side-by-side
examples/CoolingHaloWithSpin/cooling_halo.yml
View file @ fa3be81e
......@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
......@@ -32,9 +32,6 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
SoftenedIsothermalPotential:
......@@ -43,12 +40,12 @@ SoftenedIsothermalPotential:
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # controls time step
epsilon: 0.1 #softening for the isothermal potential
epsilon: 1.0 #softening for the isothermal potential
# Cooling parameters
LambdaCooling:
lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
cooling_tstep_mult: 0.1 # Dimensionless pre-factor for the time-step condition
examples/HydrostaticHalo/README
View file @ fa3be81e
Hydrostatic halo in equilibrium in an isothermal potential. Running
for 10 dynamical times.
To make the initial conditions we distribute gas particles randomly in
a cube with a side length twice that of the virial radius. The density
profile of the gas is proportional to r^(-2) where r is the distance
from the centre of the cube.
The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
The parameter v_rot (in makeIC.py and hydrostatic.yml) sets the circular
velocity of the halo, and by extension, the viral radius, viral mass,
and the internal energy of the gas such that hydrostatic equilibrium
is achieved.
......@@ -12,10 +14,12 @@ is achieved.
To run this example, make such that the code is compiled with either
the isothermal potential or softened isothermal potential set in
src/const.h. In the latter case, a (small) value of epsilon needs to
be set in cooling.yml. 0.1 kpc should work well.
be set in hydrostatic.yml. ~1 kpc should work well.
The plotting scripts produce a plot of the density, internal energy
and radial velocity profile for each
snapshot. test_energy_conservation.py shows the evolution of energy
and radial velocity profile for each snapshot and divides the profile
by the expected profile.
The script test_energy_conservation.py shows the evolution of energy
with time. These can be used to check if the example has run properly.
examples/HydrostaticHalo/hydrostatic.yml
View file @ fa3be81e
......@@ -9,7 +9,7 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1.0 # The end time of the simulation (in internal units).
time_end: 30. # 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: 1e-2 # The maximal time-step size of the simulation (in internal units).
......@@ -20,21 +20,18 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
basename: Hydrostatic # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: Hydrostatic.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
SoftenedIsothermalPotential:
......@@ -42,6 +39,6 @@ SoftenedIsothermalPotential:
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
epsilon: 0.1
epsilon: 1.0
timestep_mult: 0.03 # controls time step
examples/HydrostaticHalo/makeIC.py
View file @ fa3be81e
......@@ -116,9 +116,9 @@ print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
print np.mean(coords[:,0])
print np.mean(coords[:,1])
print np.mean(coords[:,2])
#print np.mean(coords[:,0])
#print np.mean(coords[:,1])
#print np.mean(coords[:,2])
#now find the particles which are within the box
......
examples/HydrostaticHalo/run.sh
View file @ fa3be81e
......@@ -2,14 +2,23 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 10000
python makeIC.py 100000
../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
# Run for 10 dynamical times
../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
python density_profile.py 2. 200 100
echo "Plotting density profiles"
mkdir plots
mkdir plots/density_profile
python density_profile.py 2. 200 300
python internal_energy_profile.py 2. 200 100
echo "Plotting internal energy profiles"
mkdir plots/internal_energy
python internal_energy_profile.py 2. 200 300
python velocity_profile.py 2. 200 100
echo "Plotting radial velocity profiles"
mkdir plots/radial_velocity_profile
python velocity_profile.py 2. 200 300
python test_energy_conservation.py
echo "Plotting energy as a function of time"
python test_energy_conservation.py 300
examples/HydrostaticHalo/test_energy_conservation.py
View file @ fa3be81e
......@@ -3,7 +3,7 @@ import h5py as h5
import matplotlib.pyplot as plt
import sys
n_snaps = 5
n_snaps = int(sys.argv[1])
#some constants
OMEGA = 0.3 # Cosmological matter fraction at z = 0
......@@ -24,7 +24,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
......
examples/IsothermalPotential/GravityOnly/makeIC.py → examples/IsothermalPotential/makeIC.py
View file @ fa3be81e
......@@ -141,17 +141,17 @@ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = numpy.zeros(1)
m = numpy.full((numPart, ), mass)
m = numpy.full((numPart, ), mass, dtype='f')
ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
h = numpy.full((numPart, ), 1.1255 * boxSize / L)
h = numpy.full((numPart, ), 1.1255 * boxSize / L, dtype='f')
ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
ds[()] = h
h = numpy.zeros(1)
u = numpy.full((numPart, ), internalEnergy)
u = numpy.full((numPart, ), internalEnergy, dtype='f')
ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
ds[()] = u
u = numpy.zeros(1)
......
examples/IsothermalPotential/GravityOnly/run.sh → examples/IsothermalPotential/run.sh
View file @ fa3be81e
......@@ -7,4 +7,4 @@ then
python makeIC.py 1000 1
fi
../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log
../swift -g -t 1 isothermal.yml 2>&1 | tee output.log
examples/Makefile.am
View file @ fa3be81e
......@@ -65,6 +65,9 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/test.pro \
GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py \
IsothermalPotential/README IsothermalPotential/run.sh IsothermalPotential/test.pro IsothermalPotential/isothermal.yml IsothermalPotential/makeIC.py \
KelvinHelmholtz_2D/kelvinHelmholtz.yml KelvinHelmholtz_2D/makeIC.py KelvinHelmholtz_2D/plotSolution.py KelvinHelmholtz_2D/run.sh \
MultiTypes/makeIC.py MultiTypes/multiTypes.yml MultiTypes/run.sh \
PerturbedBox_2D/makeIC.py PerturbedBox_2D/perturbedPlane.yml \
......
src/cooling/const_lambda/cooling.h
View file @ fa3be81e
......@@ -24,6 +24,7 @@
#define SWIFT_COOLING_CONST_LAMBDA_H
/* Some standard headers. */
#include <float.h>
#include <math.h>
/* Local includes. */
......@@ -84,7 +85,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Calculate du_dt */
const float du_dt = cooling_rate(phys_const, us, cooling, p);
/* Intergrate cooling equation, but enforce energy floor */
/* Integrate cooling equation, but enforce energy floor */
float u_new;
if (u_old + du_dt * dt > u_floor) {
u_new = u_old + du_dt * dt;
......@@ -92,6 +93,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
u_new = u_floor;
}
/* Don't allow particle to cool too much in one timestep */
if (u_new < 0.5f * u_old) u_new = 0.5f * u_old;
/* Update the internal energy */
hydro_set_internal_energy(p, u_new);
......@@ -112,13 +116,16 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
const struct phys_const* restrict phys_const,
const struct UnitSystem* restrict us, const struct part* restrict p) {
/* Get du_dt */
const float du_dt = cooling_rate(phys_const, us, cooling, p);
/* Get current internal energy (dt=0) */
const float u = hydro_get_internal_energy(p, 0.f);
const float du_dt = cooling_rate(phys_const, us, cooling, p);
return u / fabsf(du_dt);
/* If we are close to (or below) the energy floor, we ignore cooling timestep
*/
if (u < 1.01f * cooling->min_energy)
return FLT_MAX;
else
return cooling->cooling_tstep_mult * u / fabsf(du_dt);
}
/**
......
src/engine.c
View file @ fa3be81e
......@@ -2246,7 +2246,7 @@ void engine_prepare(struct engine *e, int nodrift) {
TIMER_TOC(timer_prepare);
if (e->verbose)
message("took %.3f %s (including marktask, rebuild and reweight).",
message("took %.3f %s (including drift all, rebuild and reweight).",
clocks_from_ticks(getticks() - tic), clocks_getunit());
}
......
src/runner.c
View file @ fa3be81e
......@@ -133,16 +133,16 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
return;
}
} else {
if (count > 0) {
if (count > 0) {
/* do sourceterms in this cell? */
const int incell =
sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
if (incell == 1) {
sourceterms_apply(r, sourceterms, c);
/* do sourceterms in this cell? */
const int incell =
sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
if (incell == 1) {
sourceterms_apply(r, sourceterms, c);
}
}
}
......@@ -174,23 +174,18 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
return;
}
#ifdef TASK_VERBOSE
OUT;
#endif
/* Loop over the gparts in this cell. */
for (int i = 0; i < gcount; i++) {
} else {
/* Get a direct pointer on the part. */
struct gpart *restrict gp = &gparts[i];
/* Loop over the gparts in this cell. */
for (int i = 0; i < gcount; i++) {
/* Is this part within the time step? */
if (gpart_is_active(gp, e)) {
/* Get a direct pointer on the part. */
struct gpart *restrict gp = &gparts[i];
external_gravity_acceleration(time, potential, constants, gp);
/* Is this part within the time step? */
if (gpart_is_active(gp, e)) {
external_gravity_acceleration(time, potential, constants, gp);
}
}
}
......@@ -222,26 +217,22 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
return;
}
#ifdef TASK_VERBOSE
OUT;
#endif
} else {
/* Loop over the parts in this cell. */
for (int i = 0; i < count; i++) {
/* Loop over the parts in this cell. */
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[i];
struct xpart *restrict xp = &xparts[i];
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[i];
struct xpart *restrict xp = &xparts[i];
/* Kick has already updated ti_end, so need to check ti_begin */
if (p->ti_begin == ti_current) {
/* Kick has already updated ti_end, so need to check ti_begin */
if (p->ti_begin == ti_current) {
const double dt = (p->ti_end - p->ti_begin) * timeBase;
const double dt = (p->ti_end - p->ti_begin) * timeBase;
cooling_cool_part(constants, us, cooling_func, p, xp, dt);
cooling_cool_part(constants, us, cooling_func, p, xp, dt);
}
}
}
......@@ -504,7 +495,6 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
return;
} else {
/* Loop over the parts in this cell. */
......@@ -543,8 +533,9 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_extra_ghost(struct runner *r, struct cell *c) {
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
#ifdef EXTRA_HYDRO_LOOP
......@@ -552,14 +543,15 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
const int count = c->count;
const struct engine *e = r->e;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k]);
return;
if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
} else {
/* Loop over the parts in this cell. */
......@@ -576,6 +568,8 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
}
}
if (timer) TIMER_TOC(timer_do_extra_ghost);
#else
error("SWIFT was not compiled with the extra hydro loop activated.");
#endif
......@@ -587,8 +581,9 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_ghost(struct runner *r, struct cell *c) {
void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
......@@ -611,140 +606,141 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
return;
}
if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
} else {
/* Init the IDs that have to be updated. */
int *pid = NULL;
if ((pid = malloc(sizeof(int) * count)) == NULL)
error("Can't allocate memory for pid.");
for (int k = 0; k < count; k++) pid[k] = k;
/* Init the IDs that have to be updated. */
int *pid = NULL;
if ((pid = malloc(sizeof(int) * count)) == NULL)
error("Can't allocate memory for pid.");
for (int k = 0; k < count; k++) pid[k] = k;
/* While there are particles that need to be updated... */
for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
/* While there are particles that need to be updated... */
for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
/* Reset the redo-count. */
redo = 0;
/* Reset the redo-count. */
redo = 0;
/* Loop over the parts in this cell. */
for (int i = 0; i < count; i++) {
/* Loop over the parts in this cell. */
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[pid[i]];
struct xpart *restrict xp = &xparts[pid[i]];
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[pid[i]];
struct xpart *restrict xp = &xparts[pid[i]];
/* Is this part within the timestep? */
if (part_is_active(p, e)) {
/* Is this part within the timestep? */
if (part_is_active(p, e)) {
/* Finish the density calculation */
hydro_end_density(p);
/* Finish the density calculation */
hydro_end_density(p);
float h_corr = 0.f;
float h_corr = 0.f;
/* If no derivative, double the smoothing length. */
if (p->density.wcount_dh == 0.0f) h_corr = p->h;
/* If no derivative, double the smoothing length. */
if (p->density.wcount_dh == 0.0f) h_corr = p->h;
/* Otherwise, compute the smoothing length update (Newton step). */
else {
h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
/* Otherwise, compute the smoothing length update (Newton step). */
else {
h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
/* Truncate to the range [ -p->h/2 , p->h ]. */
h_corr = (h_corr < p->h) ? h_corr : p->h;
h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
}
/* Truncate to the range [ -p->h/2 , p->h ]. */
h_corr = (h_corr < p->h) ? h_corr : p->h;
h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
}
/* Did we get the right number density? */
if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
/* Did we get the right number density? */
if (p->density.wcount > max_wcount ||
p->density.wcount < min_wcount) {
/* Ok, correct then */
p->h += h_corr;
/* Ok, correct then */
p->h += h_corr;
/* Flag for another round of fun */
pid[redo] = pid[i];
redo += 1;
/* Flag for another round of fun */
pid[redo] = pid[i];
redo += 1;
/* Re-initialise everything */
hydro_init_part(p);
/* Re-initialise everything */
hydro_init_part(p);
/* Off we go ! */
continue;
}
/* Off we go ! */
continue;
}
/* We now have a particle whose smoothing length has converged */
/* We now have a particle whose smoothing length has converged */
/* As of here, particle force variables will be set. */
/* As of here, particle force variables will be set. */
/* Compute variables required for the force loop */
hydro_prepare_force(p, xp, ti_current, timeBase);
/* Compute variables required for the force loop */
hydro_prepare_force(p, xp, ti_current, timeBase);
/* The particle force values are now set. Do _NOT_
try to read any particle density variables! */
/* The particle force values are now set. Do _NOT_
try to read any particle density variables! */
/* Prepare the particle for the force loop over neighbours */
hydro_reset_acceleration(p);
/* Prepare the particle for the force loop over neighbours */
hydro_reset_acceleration(p);
}
}
}
/* We now need to treat the particles whose smoothing length had not
* converged again */
/* We now need to treat the particles whose smoothing length had not
* converged again */
/* Re-set the counter for the next loop (potentially). */
count = redo;
if (count > 0) {
/* Re-set the counter for the next loop (potentially). */
count = redo;
if (count > 0) {
/* Climb up the cell hierarchy. */
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
/* Climb up the cell hierarchy. */
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
/* Run through this cell's density interactions. */
for (struct link *l = finger->density; l != NULL; l = l->next) {
/* Run through this cell's density interactions. */
for (struct link *l = finger->density; l != NULL; l = l->next) {
/* Self-interaction? */
if (l->t->type == task_type_self)
runner_doself_subset_density(r, finger, parts, pid, count);
/* Self-interaction? */
if (l->t->type == task_type_self)
runner_doself_subset_density(r, finger, parts, pid, count);
/* Otherwise, pair interaction? */
else if (l->t->type == task_type_pair) {
/* Otherwise, pair interaction? */
else if (l->t->type == task_type_pair) {
/* Left or right? */
if (l->t->ci == finger)
runner_dopair_subset_density(r, finger, parts, pid, count,
l->t->cj);
else
runner_dopair_subset_density(r, finger, parts, pid, count,
l->t->ci);
/* Left or right? */
if (l->t->ci == finger)
runner_dopair_subset_density(r, finger, parts, pid, count,
l->t->cj);
else
runner_dopair_subset_density(r, finger, parts, pid, count,
l->t->ci);
}
}
/* Otherwise, sub-self interaction? */
else if (l->t->type == task_type_sub_self)
runner_dosub_subset_density(r, finger, parts, pid, count, NULL, -1,
1);
/* Otherwise, sub-pair interaction? */
else if (l->t->type == task_type_sub_pair) {