diff --git a/.gitignore b/.gitignore
index d2f69ad4e8759839082a38d4558a5ab7bac0e18d..a2fd8d3d6e932eff75dcd0d80e22a4073f7c2bd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -132,6 +132,7 @@ tests/testEOS*.txt
tests/testEOS*.png
tests/testUtilities
tests/testCbrt
+tests/testFormat.sh
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/examples/ConstantCosmoVolume/Gadget2/README b/examples/ConstantCosmoVolume/Gadget2/README
new file mode 100644
index 0000000000000000000000000000000000000000..8063a5da1e68b608759d35373e6006d17bf5047e
--- /dev/null
+++ b/examples/ConstantCosmoVolume/Gadget2/README
@@ -0,0 +1,6 @@
+This parameter file can be used to run the exact same example
+with the Gadget-2 code.
+
+The Gadget code has to be compiled with at least the following options:
+ - PERIODIC
+ - HAVE_HDF5
diff --git a/examples/ConstantCosmoVolume/Gadget2/constant_volume.param b/examples/ConstantCosmoVolume/Gadget2/constant_volume.param
new file mode 100644
index 0000000000000000000000000000000000000000..a57e3293ae9dce92743737d42605615d3e365f7a
--- /dev/null
+++ b/examples/ConstantCosmoVolume/Gadget2/constant_volume.param
@@ -0,0 +1,138 @@
+
+% System of units
+
+UnitLength_in_cm 3.08567758e24 % 1.0 Mpc
+UnitMass_in_g 1.98848e43 % 1.0e10 solar masses
+UnitVelocity_in_cm_per_s 1e5 % 1 km/sec
+GravityConstantInternal 4.300927e+01 % Same value as SWIFT
+
+% Relevant files
+InitCondFile constantBox
+OutputDir data/
+
+EnergyFile energy.txt
+InfoFile info.txt
+TimingsFile timings.txt
+CpuFile cpu.txt
+
+RestartFile restart
+SnapshotFileBase box
+
+OutputListFilename dummy
+
+% CPU time -limit
+
+TimeLimitCPU 360000 % = 10 hours
+ResubmitOn 0
+ResubmitCommand my-scriptfile
+
+
+% Code options
+
+ICFormat 3
+SnapFormat 3
+ComovingIntegrationOn 1
+
+TypeOfTimestepCriterion 0
+OutputListOn 0
+PeriodicBoundariesOn 1
+
+% Caracteristics of run
+
+TimeBegin 0.00990099 % z = 100
+TimeMax 1. % z = 0.
+
+Omega0 1.0
+OmegaLambda 0.0
+OmegaBaryon 1.0
+HubbleParam 1.0
+BoxSize 64.
+
+% Output frequency
+
+TimeBetSnapshot 1.04
+TimeOfFirstSnapshot 0.00991
+
+CpuTimeBetRestartFile 36000.0 ; here in seconds
+TimeBetStatistics 0.05
+
+NumFilesPerSnapshot 1
+NumFilesWrittenInParallel 1
+
+% Accuracy of time integration
+
+ErrTolIntAccuracy 0.025
+MaxRMSDisplacementFac 0.25
+CourantFac 0.1
+MaxSizeTimestep 0.002
+MinSizeTimestep 1e-7
+
+
+% Tree algorithm, force accuracy, domain update frequency
+
+ErrTolTheta 0.3
+TypeOfOpeningCriterion 0
+ErrTolForceAcc 0.005
+
+TreeDomainUpdateFrequency 0.01
+
+% Further parameters of SPH
+
+DesNumNgb 48
+MaxNumNgbDeviation 1.
+ArtBulkViscConst 0.8
+InitGasTemp 0.
+MinGasTemp 0.
+
+% Memory allocation
+
+PartAllocFactor 1.6
+TreeAllocFactor 0.8
+BufferSize 30
+
+% Softening lengths
+
+MinGasHsmlFractional 0.001
+
+SofteningGas 0.08 # 80 kpc / h = 1/25 of mean inter-particle separation
+SofteningHalo 0
+SofteningDisk 0
+SofteningBulge 0
+SofteningStars 0
+SofteningBndry 0
+
+SofteningGasMaxPhys 0.08 # 80 kpc / h = 1/25 of mean inter-particle separation
+SofteningHaloMaxPhys 0
+SofteningDiskMaxPhys 0
+SofteningBulgeMaxPhys 0
+SofteningStarsMaxPhys 0
+SofteningBndryMaxPhys 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/ConstantCosmoVolume/README b/examples/ConstantCosmoVolume/README
new file mode 100644
index 0000000000000000000000000000000000000000..de84f6909a7c9086603f5d717232d60ff5e312e3
--- /dev/null
+++ b/examples/ConstantCosmoVolume/README
@@ -0,0 +1,7 @@
+This test is a small cosmological volume with constant density and internal energy.
+The ICs are generated from a glass file to minimize the build-up of peculiar velocities
+over time.
+
+The cosmology model is very simple by design. We use Omega_m = 1, Omega_b = 1, h = 1.
+
+The solution script plots the expected solution both in comoving and physical frames.
diff --git a/examples/ConstantCosmoVolume/constant_volume.yml b/examples/ConstantCosmoVolume/constant_volume.yml
index ad31fd1972565b0d7683711a20db78e854c3dc5f..db5b3de536b5b71229a207a3fdcfab8480718ef3 100644
--- a/examples/ConstantCosmoVolume/constant_volume.yml
+++ b/examples/ConstantCosmoVolume/constant_volume.yml
@@ -1,44 +1,44 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
+ UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
+ UnitLength_in_cgs: 3.08567758e24 # 1 Mpc
+ UnitVelocity_in_cgs: 1e5 # 1 km/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
Cosmology:
Omega_m: 1.
Omega_lambda: 0.
Omega_b: 1.
h: 1.
- a_begin: 0.00990099
- a_end: 1.0
+ a_begin: 0.00990099 # z_ini = 100.
+ a_end: 1.0 # z_end = 0.
# Parameters governing the time integration
TimeIntegration:
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 5e-3 # The maximal time-step size of the simulation (in internal units).
+ dt_min: 1e-7
+ dt_max: 2e-3
# Parameters governing the snapshots
Snapshots:
- basename: box # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 1.04 # Time difference between consecutive outputs (in internal units)
+ basename: box
+ delta_time: 1.04
scale_factor_first: 0.00991
compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 2. # Time between statistics output
+ scale_factor_first: 0.00991
+ delta_time: 1.1
# Parameters for the hydrodynamics scheme
SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ resolution_eta: 1.2348 # "48 ngb" for the 3D cubic spline
+ CFL_condition: 0.1
# Parameters related to the initial conditions
InitialConditions:
- file_name: ./constantBox.hdf5 # The file to read
+ file_name: ./constantBox.hdf5
Scheduler:
max_top_level_cells: 8
@@ -48,6 +48,6 @@ Gravity:
mesh_side_length: 32
eta: 0.025
theta: 0.3
- r_cut_max: 5.
- comoving_softening: 0.05
- max_physical_softening: 0.05
+ comoving_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+ max_physical_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+
diff --git a/examples/ConstantCosmoVolume/getGlass.sh b/examples/ConstantCosmoVolume/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..01b4474ac21666c843b7abedfa39a76948934911
--- /dev/null
+++ b/examples/ConstantCosmoVolume/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5
diff --git a/examples/ConstantCosmoVolume/makeIC.py b/examples/ConstantCosmoVolume/makeIC.py
index 970f197400129d2ca3f3a7b6ff2cfdd5a7f53f3f..84ba2a55bdc5eb07a7c041e8255f7a5fec0cff07 100644
--- a/examples/ConstantCosmoVolume/makeIC.py
+++ b/examples/ConstantCosmoVolume/makeIC.py
@@ -25,7 +25,7 @@ T_i = 100. # Initial temperature of the gas (in K)
z_i = 100. # Initial redshift
gamma = 5./3. # Gas adiabatic index
numPart_1D = 32
-#glassFile = "glassCube_32.hdf5"
+glassFile = "gravity_glassCube_32.hdf5"
fileName = "constantBox.hdf5"
@@ -56,17 +56,16 @@ unit_u_in_si = unit_v_in_si**2
#---------------------------------------------------
# Read the glass file
-#glass = h5py.File(glassFile, "r" )
+glass = h5py.File(glassFile, "r" )
# Read particle positions and h from the glass
-#pos = glass["/PartType0/Coordinates"][:,:]
-#h = glass["/PartType0/SmoothingLength"][:] * 0.3
-#glass.close()
+pos = glass["/PartType1/Coordinates"][:,:]
+glass.close()
# Total number of particles
-#numPart = size(h)
-#if numPart != numPart_1D**3:
-# print "Non-matching glass file"
+numPart = size(pos)/3
+if numPart != numPart_1D**3:
+ print "Non-matching glass file"
numPart = numPart_1D**3
# Set box size and interparticle distance
@@ -78,9 +77,7 @@ a_i = 1. / (1. + z_i)
m_i = boxSize**3 * rho_0 / numPart
# Build the arrays
-#pos *= boxSize
-#h *= boxSize
-coords = zeros((numPart, 3))
+pos *= boxSize
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = zeros(numPart)
@@ -92,9 +89,9 @@ for i in range(numPart_1D):
for j in range(numPart_1D):
for k in range(numPart_1D):
index = i * numPart_1D**2 + j * numPart_1D + k
- coords[index,0] = (i + 0.5) * delta_x
- coords[index,1] = (j + 0.5) * delta_x
- coords[index,2] = (k + 0.5) * delta_x
+ #coords[index,0] = (i + 0.5) * delta_x
+ #coords[index,1] = (j + 0.5) * delta_x
+ #coords[index,2] = (k + 0.5) * delta_x
u[index] = kB_in_SI * T_i / (gamma - 1.) / mH_in_kg
h[index] = 1.2348 * delta_x
m[index] = m_i
@@ -103,7 +100,7 @@ for i in range(numPart_1D):
v[index,2] = 0.
# Unit conversion
-coords /= unit_l_in_si
+pos /= unit_l_in_si
v /= unit_v_in_si
m /= unit_m_in_si
h /= unit_l_in_si
@@ -140,7 +137,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=coords, dtype='d', compression="gzip", shuffle=True)
+grp.create_dataset('Coordinates', data=pos, dtype='d', compression="gzip", shuffle=True)
grp.create_dataset('Velocities', data=v, dtype='f',compression="gzip", shuffle=True)
grp.create_dataset('Masses', data=m, dtype='f', compression="gzip", shuffle=True)
grp.create_dataset('SmoothingLength', data=h, dtype='f', compression="gzip", shuffle=True)
diff --git a/examples/ConstantCosmoVolume/run.sh b/examples/ConstantCosmoVolume/run.sh
index 521659b26d6e4d3c07a8322ba92fa3d52f0ba2cf..a2d9cdc63fbef6208caa44e05f91e94a27c4ff20 100755
--- a/examples/ConstantCosmoVolume/run.sh
+++ b/examples/ConstantCosmoVolume/run.sh
@@ -1,6 +1,11 @@
#!/bin/bash
# Generate the initial conditions if they are not present.
+if [ ! -e gravity_glassCube_32.hdf5 ]
+then
+ echo "Fetching initial grvity glass file for the constant cosmological box example..."
+ ./getGlass.sh
+fi
if [ ! -e constantBox.hdf5 ]
then
echo "Generating initial conditions for the uniform cosmo box example..."
diff --git a/examples/SantaBarbara/README b/examples/SantaBarbara/README
new file mode 100644
index 0000000000000000000000000000000000000000..973d09a438e8ff39d7c4a91c5a08d4dfc4a51a60
--- /dev/null
+++ b/examples/SantaBarbara/README
@@ -0,0 +1,16 @@
+Initital conditions for the Santa-Barbara cluster comparison project.
+These have been regenerated from the orinigal Frenk et al. 1999 paper.
+
+The cosmology is Omega_m = 1, Omega_b = 0.1, h = 0.5 and sigma_8 = 0.9.
+
+The ICs are 256^3 particles in a 64^3 Mpc^3 volume. This is about 10x
+higher resolution than in the original paper. The ICs have been
+created for Gadget and the positions and box size are hence expressed
+in h-full units (e.g. box size of 32 / h Mpc). Similarly, the peculiar
+velocitites contain an extra sqrt(a) factor.
+
+We will use SWIFT to cancel the h- and a-factors from the ICs. Gas
+particles will be generated at startup.
+
+MD5 check-sum of the ICS:
+ba9ab4f00a70d39fa601a4a59984b343 SantaBarbara.hdf5
diff --git a/examples/SantaBarbara/getIC.sh b/examples/SantaBarbara/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a3073631ceedea47c8ac218a5e62529efee6fc56
--- /dev/null
+++ b/examples/SantaBarbara/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/SantaBarbara.hdf5
diff --git a/examples/SantaBarbara/santa_barbara.yml b/examples/SantaBarbara/santa_barbara.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b01321c521ec28b78c97cd082f0fd520d7024c3d
--- /dev/null
+++ b/examples/SantaBarbara/santa_barbara.yml
@@ -0,0 +1,59 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98848e43 # 10^10 Msun
+ UnitLength_in_cgs: 3.08567758e24 # 1 Mpc
+ UnitVelocity_in_cgs: 1e5 # 1 km/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.5
+ a_begin: 0.047619048 # z_ini = 20
+ a_end: 1.0 # z_end = 0
+ Omega_m: 1.0
+ Omega_lambda: 0.0
+ Omega_b: 0.1
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_max: 0.01
+ dt_min: 1e-10
+
+Scheduler:
+ max_top_level_cells: 16
+ cell_split_size: 100
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: santabarbara
+ scale_factor_first: 0.05
+ delta_time: 1.02
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02
+ scale_factor_first: 0.05
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025
+ theta: 0.5
+ comoving_softening: 0.01 # 10 kpc = 1/25 mean inter-particle separation
+ max_physical_softening: 0.00263 # 10 ckpc = 2.63 pkpc at z=2.8 (EAGLE-like evolution of softening).
+ mesh_side_length: 128
+
+# Parameters of the hydro scheme
+SPH:
+ resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel
+ CFL_condition: 0.1
+ initial_temperature: 1200. # (1 + z_ini)^2 * 2.72K
+ minimal_temperature: 100.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./SantaBarbara.hdf5
+ cleanup_h_factors: 1 # ICs were generated for Gadget, we need to get rid of h-factors
+ cleanup_velocity_factors: 1 # ICs were generated for Gadget, we need to get rid of sqrt(a) factors in the velocity
+ generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
+ cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
\ No newline at end of file
diff --git a/examples/SmallCosmoVolume/Gadget2/README b/examples/SmallCosmoVolume/Gadget2/README
new file mode 100644
index 0000000000000000000000000000000000000000..8063a5da1e68b608759d35373e6006d17bf5047e
--- /dev/null
+++ b/examples/SmallCosmoVolume/Gadget2/README
@@ -0,0 +1,6 @@
+This parameter file can be used to run the exact same example
+with the Gadget-2 code.
+
+The Gadget code has to be compiled with at least the following options:
+ - PERIODIC
+ - HAVE_HDF5
diff --git a/examples/SmallCosmoVolume/Gadget2/small_cosmo_volume.param b/examples/SmallCosmoVolume/Gadget2/small_cosmo_volume.param
new file mode 100644
index 0000000000000000000000000000000000000000..4eaaab4cb124db898928c75e7a7a03bb850c5a9f
--- /dev/null
+++ b/examples/SmallCosmoVolume/Gadget2/small_cosmo_volume.param
@@ -0,0 +1,137 @@
+
+% System of units
+
+UnitLength_in_cm 3.08567758e24 % 1.0 Mpc
+UnitMass_in_g 1.98848e43 % 1.0e10 solar masses
+UnitVelocity_in_cm_per_s 1e5 % 1 km/sec
+GravityConstantInternal 4.300927e+01 % Same value as SWIFT
+
+% Relevant files
+InitCondFile small_cosmo_volume
+OutputDir data/
+
+EnergyFile energy.txt
+InfoFile info.txt
+TimingsFile timings.txt
+CpuFile cpu.txt
+
+RestartFile restart
+SnapshotFileBase box
+
+OutputListFilename dummy
+
+% CPU time -limit
+
+TimeLimitCPU 360000 % = 10 hours
+ResubmitOn 0
+ResubmitCommand my-scriptfile
+
+
+% Code options
+
+ICFormat 3
+SnapFormat 3
+ComovingIntegrationOn 1
+
+TypeOfTimestepCriterion 0
+OutputListOn 0
+PeriodicBoundariesOn 1
+
+% Caracteristics of run
+
+TimeBegin 0.019607843 % z = 50.
+TimeMax 1. % z = 0.
+
+Omega0 0.276
+OmegaLambda 0.724
+OmegaBaryon 0.0455
+HubbleParam 0.703
+BoxSize 100. % Mpc / h
+
+% Output frequency
+
+TimeBetSnapshot 1.02
+TimeOfFirstSnapshot 0.02
+
+CpuTimeBetRestartFile 36000.0 ; here in seconds
+TimeBetStatistics 0.02
+
+NumFilesPerSnapshot 1
+NumFilesWrittenInParallel 1
+
+% Accuracy of time integration
+
+ErrTolIntAccuracy 0.025
+MaxRMSDisplacementFac 0.25
+CourantFac 0.1
+MaxSizeTimestep 0.01
+MinSizeTimestep 1e-6
+
+
+% Tree algorithm, force accuracy, domain update frequency
+
+ErrTolTheta 0.3
+TypeOfOpeningCriterion 1
+ErrTolForceAcc 0.005
+TreeDomainUpdateFrequency 0.01
+
+% Further parameters of SPH
+
+DesNumNgb 48
+MaxNumNgbDeviation 1.
+ArtBulkViscConst 0.8
+InitGasTemp 0.
+MinGasTemp 0.
+
+% Memory allocation
+
+PartAllocFactor 1.6
+TreeAllocFactor 0.8
+BufferSize 30
+
+% Softening lengths
+
+MinGasHsmlFractional 0.001
+
+SofteningGas 0
+SofteningHalo 0.0625 # 62.5 kpc / h = 1/25 of mean inter-particle separation
+SofteningDisk 0
+SofteningBulge 0
+SofteningStars 0
+SofteningBndry 0
+
+SofteningGasMaxPhys 0
+SofteningHaloMaxPhys 0.0625 # 62.5 kpc / h = 1/25 of mean inter-particle separation
+SofteningDiskMaxPhys 0
+SofteningBulgeMaxPhys 0
+SofteningStarsMaxPhys 0
+SofteningBndryMaxPhys 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/SmallCosmoVolume/README b/examples/SmallCosmoVolume/README
index 68c137aee30c08bb476b760c75dceaa5e1ede87e..14a289cf4a1d638c18f421f23ca8bcf0ced68d1b 100644
--- a/examples/SmallCosmoVolume/README
+++ b/examples/SmallCosmoVolume/README
@@ -6,4 +6,4 @@ The ICs have been generated to run with Gadget-2 so we need to switch
on the options to cancel the h-factors and a-factors at reading time.
MD5 checksum of the ICs:
-2a9c603ffb1f6d29f3d98a3ecb9d3238 small_cosmo_volume.hdf5
+08736c3101fd738e22f5159f78e6022b small_cosmo_volume.hdf5
diff --git a/examples/SmallCosmoVolume/small_cosmo_volume.yml b/examples/SmallCosmoVolume/small_cosmo_volume.yml
index 32ec15db6be35fed4eb0c0168f52f0ba919158ea..b8e907db6cd534d2262778fa929644f87d0ad1a4 100644
--- a/examples/SmallCosmoVolume/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/small_cosmo_volume.yml
@@ -1,29 +1,26 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
+ UnitLength_in_cgs: 3.08567758e24 # 1 Mpc
+ UnitVelocity_in_cgs: 1e5 # 1 km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Structure finding options
StructureFinding:
- config_file_name: stf_input_6dfof_dmonly_sub.cfg # Name of the STF config file.
- basename: ./stf # Common part of the name of output files.
- output_time_format: 0 # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
- scale_factor_first: 0.92 # Scale-factor of the first snaphot (cosmological run)
- time_first: 0.01 # Time of the first structure finding output (in internal units).
- delta_step: 1000 # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
- delta_time: 1.02 # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+ config_file_name: stf_input_6dfof_dmonly_sub.cfg
+ basename: ./stf
+ output_time_format: 1
+ scale_factor_first: 0.02
+ delta_time: 1.02
-# WMAP9 cosmology
-Cosmology:
+Cosmology: # WMAP9 cosmology
Omega_m: 0.276
Omega_lambda: 0.724
Omega_b: 0.0455
h: 0.703
- a_begin: 0.0196078
- a_end: 1.0
+ a_begin: 0.019607843 # z_ini = 50.
+ a_end: 1.0 # z_end = 0.
# Parameters governing the time integration
TimeIntegration:
@@ -34,9 +31,9 @@ TimeIntegration:
Gravity:
eta: 0.025
theta: 0.3
- comoving_softening: 0.08
- max_physical_softening: 0.08
- mesh_side_length: 32
+ comoving_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
+ max_physical_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
+ mesh_side_length: 64
# Parameters governing the snapshots
Snapshots:
diff --git a/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
index 872e0ad6f44d8092ce1da6ac030a949dc4dba5d5..7368e5654204ad600192eff3defdd5f96e986ce5 100644
--- a/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
+++ b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
@@ -104,7 +104,7 @@ Allowed_kinetic_potential_ratio=0.2
#run unbinding of field structures, aka halos
Bound_halos=0
#simple Plummer softening length when calculating gravitational energy. If cosmological simulation with period, is fraction of interparticle spacing
-Softening_length=0.
+Softening_length=0.04
#don't keep background potential when unbinding
Keep_background_potential=0
diff --git a/examples/ZeldovichPancake_3D/plotSolution.py b/examples/ZeldovichPancake_3D/plotSolution.py
index 2a175e346e041a142c6921052ccf13978afa8a38..b7e8fd93813e8906055262ed6deba9a40ec3d0a6 100644
--- a/examples/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/ZeldovichPancake_3D/plotSolution.py
@@ -69,6 +69,7 @@ scheme = sim["/HydroScheme"].attrs["Scheme"]
kernel = sim["/HydroScheme"].attrs["Kernel function"]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"]
git = sim["Code"].attrs["Git Revision"]
# Cosmological parameters
@@ -178,8 +179,8 @@ ylabel("${\\rm{Temperature}}~T$", labelpad=0)
# Information -------------------------------------
subplot(236, frameon=False)
-text(-0.49, 0.9, "Zeldovich pancake with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "$z={0:.2f}$".format(redshift))
+text(-0.49, 0.9, "Zeldovich pancake at z=%.2f "%(redshift), fontsize=10)
+text(-0.49, 0.8, "adiabatic index $\\gamma=%.2f$, viscosity $\\alpha=%.2f$"%(gas_gamma, alpha), fontsize=10)
plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
text(-0.49, 0.4, scheme, fontsize=10)
diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index 5cfa01ff954a959e06076035ae22240bb3c5a120..4d83d805cfebb837f37058167a4e3c974a936317 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -24,7 +24,8 @@ Snapshots:
basename: zeldovichPancake # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1.04 # Time difference between consecutive outputs (in internal units)
- scale_factor_first: 0.00991
+ scale_factor_first: 0.00991
+ compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/src/engine.c b/src/engine.c
index 5502dafe56b7344f2d3195697eeb2098cf1c7219..f039d619c769f6132acadc8b456fd839f3bb5ea7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2159,8 +2159,8 @@ void engine_exchange_proxy_multipoles(struct engine *e) {
const ticks tic = getticks();
/* Start by counting the number of cells to send and receive */
- int count_send_cells = 0;
- int count_recv_cells = 0;
+ int count_send = 0;
+ int count_recv = 0;
int count_send_requests = 0;
int count_recv_requests = 0;
@@ -5439,8 +5439,8 @@ void engine_makeproxies(struct engine *e) {
/* Apply BC */
if (periodic) {
dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[0]);
- dz = nearest(dz, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
}
/* Add to it for the case where the future CoMs are in the
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 4511b2d655b0e7b3293633a466c76757a0237874..37b5628a32761847959fe8c7cb7fa49fde291ba5 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -515,8 +515,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
- p->entropy_dt = 0.5f * cosmo->a2_inv *
- gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
+ p->entropy_dt =
+ 0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
}
/**
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index b2af8909bed1780586a5130370222c9b8157d724..b2d77a261e91a4f8460189cfce99951a820ef8f1 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -479,14 +479,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
- const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
@@ -599,14 +602,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
/* Are the particles moving towards each others ? */
- const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Signal velocity */
@@ -671,7 +677,7 @@ runner_iact_nonsym_1_vec_force(
vector dvx, dvy, dvz;
vector xi, xj;
vector hid_inv, hjd_inv;
- vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
+ vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr, dvdr_Hubble;
vector piax, piay, piaz;
vector pih_dt;
vector v_sig;
@@ -723,14 +729,14 @@ runner_iact_nonsym_1_vec_force(
dvz.v = vec_sub(viz.v, vjz.v);
/* Compute dv dot r. */
- dvdr.v =
- vec_fma(dvx.v, dx->v,
- vec_fma(dvy.v, dy->v,
- vec_fma(dvz.v, dz->v, vec_mul(v_a2_Hubble.v, r2->v))));
+ dvdr.v = vec_fma(dvx.v, dx->v, vec_fma(dvy.v, dy->v, vec_mul(dvz.v, dz->v)));
+
+ /* Add Hubble flow */
+ dvdr_Hubble.v = vec_add(dvdr.v, vec_mul(v_a2_Hubble.v, r2->v));
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
- omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
+ omega_ij.v = vec_fmin(dvdr_Hubble.v, vec_setzero());
mu_ij.v = vec_mul(v_fac_mu.v,
vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
@@ -806,7 +812,7 @@ runner_iact_nonsym_2_vec_force(
vector dvx, dvy, dvz;
vector ui, uj;
vector hid_inv, hjd_inv;
- vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
+ vector wi_dx, wj_dx, wi_dr, wj_dr, dvdr, dvdr_Hubble;
vector piax, piay, piaz;
vector pih_dt;
vector v_sig;
@@ -817,7 +823,7 @@ runner_iact_nonsym_2_vec_force(
vector dvx_2, dvy_2, dvz_2;
vector ui_2, uj_2;
vector hjd_inv_2;
- vector wi_dx_2, wj_dx_2, wi_dr_2, wj_dr_2, dvdr_2;
+ vector wi_dx_2, wj_dx_2, wi_dr_2, wj_dr_2, dvdr_2, dvdr_Hubble_2;
vector piax_2, piay_2, piaz_2;
vector pih_dt_2;
vector v_sig_2;
@@ -903,18 +909,18 @@ runner_iact_nonsym_2_vec_force(
dvz_2.v = vec_sub(viz.v, vjz_2.v);
/* Compute dv dot r. */
- dvdr.v = vec_fma(
- dvx.v, dx.v,
- vec_fma(dvy.v, dy.v, vec_fma(dvz.v, dz.v, vec_mul(v_a2_Hubble.v, r2.v))));
- dvdr_2.v = vec_fma(
- dvx_2.v, dx_2.v,
- vec_fma(dvy_2.v, dy_2.v,
- vec_fma(dvz_2.v, dz_2.v, vec_mul(v_a2_Hubble.v, r2_2.v))));
+ dvdr.v = vec_fma(dvx.v, dx.v, vec_fma(dvy.v, dy.v, vec_mul(dvz.v, dz.v)));
+ dvdr_2.v = vec_fma(dvx_2.v, dx_2.v,
+ vec_fma(dvy_2.v, dy_2.v, vec_mul(dvz_2.v, dz_2.v)));
+
+ /* Add the Hubble flow */
+ dvdr_Hubble.v = vec_add(dvdr.v, vec_mul(v_a2_Hubble.v, r2.v));
+ dvdr_Hubble_2.v = vec_add(dvdr_2.v, vec_mul(v_a2_Hubble.v, r2_2.v));
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
- omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
- omega_ij_2.v = vec_fmin(dvdr_2.v, vec_setzero());
+ omega_ij.v = vec_fmin(dvdr_Hubble.v, vec_setzero());
+ omega_ij_2.v = vec_fmin(dvdr_Hubble_2.v, vec_setzero());
mu_ij.v = vec_mul(v_fac_mu.v,
vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
mu_ij_2.v = vec_mul(
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 01abe55e7267ca04a7f3e9740c10c681f86f29ea..32d6efa3347b0467d0643c66a260a633a24ac2d0 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -522,6 +522,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon.
* @param xp The particle extended data to act upon.
* @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param dt_grav The time-step for this kick (for gravity quantities).
+ * @param dt_hydro The time-step for this kick (for hydro quantities).
+ * @param dt_kick_corr The time-step for this kick (for gravity corrections).
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
@@ -570,6 +573,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
+ /* Convert the physcial internal energy to the comoving one. */
+ /* u' = a^(3(g-1)) u */
+ const float factor = 1.f / cosmo->a_factor_internal_energy;
+ p->u *= factor;
+ xp->u_full *= factor;
+
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 42fd93d6062cbfcea5cf5297eeda0bb6525f3cad..0193886f57413540d110c513ae3119ba1315f384 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -172,10 +172,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Are the particles moving towards each others ? */
- const float omega_ij = min(dvdr, 0.f);
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
@@ -285,10 +288,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
/* Are the particles moving towards each others ? */
- const float omega_ij = min(dvdr, 0.f);
+ const float omega_ij = min(dvdr_Hubble, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 9520781be45f7d9b59534c57e542e0802759aaec..65d009dd39556d9f9ff25c2c4a6d2092a0a58099 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -215,6 +215,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
p->initial_internal_energy);
io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
p->hydrogen_mass_fraction);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha);
}
#endif
diff --git a/theory/Cosmology/coordinates.tex b/theory/Cosmology/coordinates.tex
index 38a571aefea68fbe1bc7a8ebc3867109f1c4736e..a1dbff71c13cbd62acde83c14e9e81f0fbc41214 100644
--- a/theory/Cosmology/coordinates.tex
+++ b/theory/Cosmology/coordinates.tex
@@ -88,13 +88,13 @@ gravitational terms. SPH flavours that evolve the internal energy $u$ instead of
entropy require the additional equation of motion describing the evolution of
$u'$:
\begin{equation}
- \dot{u}_i' = \frac{P_i'}{\rho_i'^2}\left[3H\rho_i' + \frac{1}{a^2}f_i'\sum_jm_j\left(\mathbf{v}_i' -
- \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i)\right],
+ \dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
+ \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
\label{eq:cosmo_eom_u}
\end{equation}
-where the first term in the brackets accounts for the change in energy
-due to the expansion of the Universe. The scale-factors appearing in
-the equations are later absorbed in the time-integration operators
+
+In all these cases, the scale-factors appearing in the equations are
+later absorbed in the time-integration operators
(Sec.~\ref{ssec:operators}) such that the RHS of the equations of
motions is identical for the primed quantities to the ones obtained in
the non-cosmological case for the physical quantities.