Merge branch 'master' into mpi_periodic_gravity

.gitignore

@@ -165,6 +165,8 @@ m4/lt~obsolete.m4
 /stamp-h1
 /test-driver
 
+src/equation_of_state/planetary/*.txt
+
 # Intel compiler optimization reports
 *.optrpt
 
@@ -316,3 +318,6 @@ sympy-plots-for-*.tex/
 
 #ctags
 *tags
+
+# vim
+*.swp

configure.ac

@@ -576,7 +576,7 @@ if test "x$with_profiler" != "xno"; then
       proflibs="-lprofiler"
    fi
    AC_CHECK_LIB([profiler],[ProfilerFlush],
-    [have_profiler="yes" 
+    [have_profiler="yes"
       AC_DEFINE([WITH_PROFILER],1,[Link against the gperftools profiling library.])],
     [have_profiler="no"], $proflibs)
 
@@ -973,7 +973,7 @@ esac
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -1012,10 +1012,11 @@ case "$with_hydro" in
    shadowfax)
       AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
    ;;
-   minimal-multi-mat)
-      AC_DEFINE([MINIMAL_MULTI_MAT_SPH], [1], [Minimal Multiple Material SPH])
+   planetary)
+      AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
    ;;
 
+
    *)
       AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
    ;;

examples/ConstantCosmoVolume/constant_volume.yml

+# 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
+
+Cosmology:
+  Omega_m: 1.
+  Omega_lambda: 0.
+  Omega_b: 1.
+  h: 1.
+  a_begin: 0.00990099
+  a_end: 1.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).
+
+# 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)
+  scale_factor_first:  0.00991
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          2. # Time between statistics output
+
+# 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.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./constantBox.hdf5       # The file to read
+
+Scheduler:
+  max_top_level_cells: 8
+  cell_split_size:     50
+  
+Gravity:
+  mesh_side_length:   32
+  eta: 0.025
+  theta: 0.3
+  r_cut_max: 5.
+  comoving_softening: 0.05
+  max_physical_softening: 0.05

examples/ConstantCosmoVolume/makeIC.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+import h5py
+from numpy import *
+
+# Parameters
+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"
+fileName = "constantBox.hdf5"
+
+
+# Some units
+Mpc_in_m = 3.08567758e22
+Msol_in_kg = 1.98848e30
+Gyr_in_s = 3.08567758e19
+mH_in_kg = 1.6737236e-27
+
+# Some constants
+kB_in_SI = 1.38064852e-23
+G_in_SI = 6.67408e-11
+
+# Some useful variables in h-full units
+H_0 = 1. / Mpc_in_m * 10**5 # h s^-1
+rho_0 = 3. * H_0**2 / (8* math.pi * G_in_SI) # h^2 kg m^-3
+lambda_i = 64. / H_0 * 10**5 # h^-1 m (= 64 h^-1 Mpc)
+x_min = -0.5 * lambda_i
+x_max = 0.5 * lambda_i
+
+# SI system of units
+unit_l_in_si = Mpc_in_m
+unit_m_in_si = Msol_in_kg * 1.e10
+unit_t_in_si = Gyr_in_s
+unit_v_in_si = unit_l_in_si / unit_t_in_si
+unit_u_in_si = unit_v_in_si**2
+
+#---------------------------------------------------
+
+# Read the glass file
+#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()
+
+# Total number of particles
+#numPart = size(h)
+#if numPart != numPart_1D**3:
+#  print "Non-matching glass file"
+numPart = numPart_1D**3
+
+# Set box size and interparticle distance
+boxSize = x_max - x_min
+delta_x = boxSize / numPart_1D
+
+# Get the particle mass
+a_i = 1. / (1. + z_i)
+m_i = boxSize**3 * rho_0 / numPart
+
+# Build the arrays
+#pos *= boxSize
+#h *= boxSize
+coords = zeros((numPart, 3))
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+h = zeros(numPart)
+u = zeros(numPart)
+
+# Set the particles on the left
+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
+      u[index] = kB_in_SI * T_i / (gamma - 1.) / mH_in_kg
+      h[index] = 1.2348 * delta_x
+      m[index] = m_i
+      v[index,0] = 0.
+      v[index,1] = 0.
+      v[index,2] = 0.
+
+# Unit conversion
+coords /= unit_l_in_si
+v /= unit_v_in_si
+m /= unit_m_in_si
+h /= unit_l_in_si
+u /= unit_u_in_si
+
+boxSize /= unit_l_in_si
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, boxSize, boxSize]
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 100. * unit_l_in_si
+grp.attrs["Unit mass in cgs (U_M)"] = 1000. * unit_m_in_si
+grp.attrs["Unit time in cgs (U_t)"] = 1. * unit_t_in_si
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+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('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)
+grp.create_dataset('InternalEnergy', data=u, dtype='f', compression="gzip", shuffle=True)
+grp.create_dataset('ParticleIDs', data=ids, dtype='L', compression="gzip", shuffle=True)
+
+file.close()

examples/ConstantCosmoVolume/plotSolution.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+# Computes the analytical solution of the Zeldovich pancake and compares with
+# the simulation result
+
+# Parameters
+T_i = 100.           # Initial temperature of the gas (in K)
+z_c = 1.             # Redshift of caustic formation (non-linear collapse)
+z_i = 100.           # Initial redshift
+gas_gamma = 5./3.    # Gas adiabatic index
+
+# Physical constants needed for internal energy to temperature conversion
+k_in_J_K = 1.38064852e-23
+mH_in_kg = 1.6737236e-27
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left'    : 0.08,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.06,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.2,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("box_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+redshift = sim["/Header"].attrs["Redshift"][0]
+a = sim["/Header"].attrs["Scale-factor"][0]
+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"]
+git = sim["Code"].attrs["Git Revision"]
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+sim.close()
+
+# Mean quantities over time
+z = np.zeros(120)
+a = np.zeros(120)
+S_mean = np.zeros(120)
+S_std = np.zeros(120)
+u_mean = np.zeros(120)
+u_std = np.zeros(120)
+P_mean = np.zeros(120)
+P_std = np.zeros(120)
+rho_mean = np.zeros(120)
+rho_std = np.zeros(120)
+
+vx_mean = np.zeros(120)
+vy_mean = np.zeros(120)
+vz_mean = np.zeros(120)
+vx_std = np.zeros(120)
+vy_std = np.zeros(120)
+vz_std = np.zeros(120)
+
+for i in range(120):
+    sim = h5py.File("box_%04d.hdf5"%i, "r")
+
+    z[i] = sim["/Cosmology"].attrs["Redshift"][0]
+    a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
+    
+    S = sim["/PartType0/Entropy"][:]
+    S_mean[i] = np.mean(S)
+    S_std[i] = np.std(S)
+    
+    u = sim["/PartType0/InternalEnergy"][:]
+    u_mean[i] = np.mean(u)
+    u_std[i] = np.std(u)
+
+    P = sim["/PartType0/Pressure"][:]
+    P_mean[i] = np.mean(P)
+    P_std[i] = np.std(P)
+
+    rho = sim["/PartType0/Density"][:]
+    rho_mean[i] = np.mean(rho)
+    rho_std[i] = np.std(rho)
+
+    v = sim["/PartType0/Velocities"][:,:]
+    vx_mean[i] = np.mean(v[:,0])
+    vy_mean[i] = np.mean(v[:,1])
+    vz_mean[i] = np.mean(v[:,2])
+    vx_std[i] = np.std(v[:,0])
+    vy_std[i] = np.std(v[:,1])
+    vz_std[i] = np.std(v[:,2])
+    
+# Move to physical quantities
+rho_mean_phys = rho_mean / a**3
+u_mean_phys = u_mean / a**(3*(gas_gamma - 1.))
+S_mean_phys = S_mean
+
+# Solution in physical coordinates
+#T_solution = np.ones(T) / a
+
+figure()
+
+# Density evolution --------------------------------
+subplot(231)#, yscale="log")
+semilogx(a, rho_mean, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$", labelpad=0.)
+ylabel("${\\rm Comoving~density}$", labelpad=0.)
+
+# Thermal energy evolution --------------------------------
+subplot(232)#, yscale="log")
+semilogx(a, u_mean, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$", labelpad=0.)
+ylabel("${\\rm Comoving~internal~energy}$", labelpad=0.)
+
+# Entropy evolution --------------------------------
+subplot(233)#, yscale="log")
+semilogx(a, S_mean, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$", labelpad=0.)
+ylabel("${\\rm Comoving~entropy}$", labelpad=0.)
+
+# Peculiar velocity evolution ---------------------
+subplot(234)
+semilogx(a, vx_mean, '-', color='r', lw=1)
+semilogx(a, vy_mean, '-', color='g', lw=1)
+semilogx(a, vz_mean, '-', color='b', lw=1)
+xlabel("${\\rm Scale-factor}$", labelpad=0.)
+ylabel("${\\rm Peculiar~velocity~mean}$", labelpad=0.)
+
+# Peculiar velocity evolution ---------------------
+subplot(235)
+semilogx(a, vx_std, '--', color='r', lw=1)
+semilogx(a, vy_std, '--', color='g', lw=1)
+semilogx(a, vz_std, '--', color='b', lw=1)
+xlabel("${\\rm Scale-factor}$", labelpad=0.)
+ylabel("${\\rm Peculiar~velocity~std-dev}$", labelpad=0.)
+
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+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)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+savefig("ConstantBox_comoving.png", dpi=200)
+
+
+
+figure()
+
+# Density evolution --------------------------------
+subplot(231)#, yscale="log")
+loglog(a, rho_mean_phys, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$")
+ylabel("${\\rm Physical~density}$")
+
+# Thermal energy evolution --------------------------------
+subplot(232)#, yscale="log")
+loglog(a, u_mean_phys, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$")
+ylabel("${\\rm Physical~internal~energy}$")
+
+# Entropy evolution --------------------------------
+subplot(233)#, yscale="log")
+semilogx(a, S_mean_phys, '-', color='r', lw=1)
+xlabel("${\\rm Scale-factor}$")
+ylabel("${\\rm Physical~entropy}$")
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+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)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+savefig("ConstantBox_physical.png", dpi=200)
+
+

examples/ConstantCosmoVolume/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e constantBox.hdf5 ]
+then
+    echo "Generating initial conditions for the uniform cosmo box example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -c -G -t 8 constant_volume.yml 2>&1 | tee output.log
+
+# Plot the result
+python plotSolution.py $i

examples/MoonFormingImpact/README.md

-Canonical Moon-Forming Giant Impact
-===================================
-
-NOTE: This doesn't really work because the EOS are different to Canup (2004) so
-the impactor just glances then flies away!
-
-A version of the canonical moon-forming giant impact of Theia onto the early
-Earth (Canup 2004; Barr 2016). Both bodies are here made of a (Tillotson) iron
-core and granite mantle. Only ~10,000 particles are used for a quick and crude
-simulation.
-
-Setup
------
-
-In `swiftsim/`:
-
-`$ ./configure --with-hydro=minimal-multi-mat --with-equation-of-state=planetary`
-`$ make`
-
-In `swiftsim/examples/MoonFormingImpact/`:
-
-`$ ./get_init_cond.sh`
-
-Run
----
-
-`$ ./run.sh`
-
-Output
-------
-
-`$ python plot.py`
-`$ mplayer anim.mpg`
-

examples/MoonFormingImpact/get_init_cond.sh

-#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/moon_forming_impact.hdf5

examples/MoonFormingImpact/moon_forming_impact.yml

-# Define the system of units to use internally.
-InternalUnitSystem:
-    UnitMass_in_cgs:        5.9724e27   # Grams
-    UnitLength_in_cgs:      6.371e8     # Centimeters
-    UnitVelocity_in_cgs:    6.371e8     # Centimeters per second
-    UnitCurrent_in_cgs:     1           # Amperes
-    UnitTemp_in_cgs:        1           # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
-    time_begin:     0                   # The starting time of the simulation (in internal units).
-    time_end:       100000              # The end time of the simulation (in internal units).
-    dt_min:         0.001               # The minimal time-step size of the simulation (in internal units).
-    dt_max:         100                 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
-                                        # Common part of the name of output files
-    basename:       snapshots/moon_forming_impact
-    time_first:     0                   # Time of the first output (in internal units)
-    delta_time:     100                 # Time difference between consecutive outputs (in internal units)
-    label_delta:    100                 # Integer increment between snapshot output labels
-
-# Parameters governing the conserved quantities statistics
-Statistics:
-    delta_time:     500                 # Time between statistics output
-
-# Parameters for the hydrodynamics scheme
-SPH:
-    resolution_eta:     1.2348          # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-    delta_neighbours:   0.1             # The tolerance for the targetted number of neighbours.
-    CFL_condition:      0.2             # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters for the self-gravity scheme
-Gravity:
-    eta:                    0.025       # Constant dimensionless multiplier for time integration.
-    theta:                  0.7         # Opening angle (Multipole acceptance criterion)
-    comoving_softening:     0.005       # Comoving softening length (in internal units).
-    max_physical_softening: 0.005       # Physical softening length (in internal units).
-
-# Parameters related to the initial conditions
-InitialConditions:
-                                        # The initial conditions file to read
-    file_name:  moon_forming_impact.hdf5
-
-# Parameters related to the equation of state
-EoS:
-    planetary_use_Til:    1                       # Whether to prepare the Tillotson EOS

examples/MoonFormingImpact/plot.py

-"""
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-###############################################################################
-
-Plotting script for the Canonical Moon-Forming Giant Impact example.
-
-Save a figure for each snapshot in `./plots/` then make them into a simple
-animation with ffmpeg in `./`.
-
-Usage:
-    `$ python  plot.py  time_end  delta_time`
-
-    Sys args:
-        + `time_end` | (opt) int | The time of the last snapshot to plot.
-            Default = 100000
-        + `delta_time` | (opt) int | The time between successive snapshots.
-            Default = 100
-"""
-
-from __future__ import division
-import numpy as np
-import matplotlib
-import matplotlib.pyplot as plt
-import h5py
-import sys
-import subprocess
-
-# Particle array fields
-dtype_picle = [
-    ('m', float), ('x', float), ('y', float), ('z', float), ('v_x', float),
-    ('v_y', float), ('v_z', float), ('ID', int), ('rho', float), ('u', float),
-    ('phi', float), ('P', float), ('h', float), ('mat_ID', int), ('r', float)
-    ]
-
-s_to_hour   = 1 / 60**2
-
-# Snapshot info
-file_snap   = "./snapshots/moon_forming_impact_"
-file_plot   = "./plots/moon_forming_impact_"
-# Number of particles in the target body
-num_target  = 9496
-
-# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
-type_factor = 100
-Di_type = {
-    'Til'       : 1,
-    'HM80'      : 2,
-    'ANEOS'     : 3,
-    'SESAME'    : 4,
-}
-Di_material = {
-    # Tillotson
-    'Til_iron'      : Di_type['Til']*type_factor,
-    'Til_granite'   : Di_type['Til']*type_factor + 1,
-    'Til_water'     : Di_type['Til']*type_factor + 2,
-    # Hubbard & MacFarlane (1980) Uranus/Neptune
-    'HM80_HHe'      : Di_type['HM80']*type_factor,      # Hydrogen-helium atmosphere
-    'HM80_ice'      : Di_type['HM80']*type_factor + 1,  # H20-CH4-NH3 ice mix
-    'HM80_rock'     : Di_type['HM80']*type_factor + 2,  # SiO2-MgO-FeS-FeO rock mix
-    # ANEOS
-    'ANEOS_iron'        : Di_type['ANEOS']*type_factor,
-    'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
-    # SESAME
-    'SESAME_iron'   : Di_type['SESAME']*type_factor,
-}
-
-# Material offset for impactor particles
-ID_imp  = 10000
-# Material colours
-Di_mat_colour = {
-    # Target
-    Di_material['Til_iron']             : 'orange',
-    Di_material['Til_granite']          : '#FFF0E0',
-    # Impactor
-    Di_material['Til_iron'] + ID_imp    : 'dodgerblue',
-    Di_material['Til_granite'] + ID_imp : '#A080D0',
-    }
-
-
-def load_snapshot(filename):
-    """ Load the hdf5 snapshot file and return the structured particle array.
-    """
-    # Add extension if needed
-    if (filename[-5:] != ".hdf5"):
-        filename += ".hdf5"
-
-    # Load the hdf5 file
-    with h5py.File(filename, 'r') as f:
-        header      = f['Header'].attrs
-        A2_pos      = f['PartType0/Coordinates'].value
-        A2_vel      = f['PartType0/Velocities'].value
-
-        # Structured array of all particle data
-        A2_picle    = np.empty(header['NumPart_Total'][0],
-                               dtype=dtype_picle)
-
-        A2_picle['x']       = A2_pos[:, 0]
-        A2_picle['y']       = A2_pos[:, 1]
-        A2_picle['z']       = A2_pos[:, 2]
-        A2_picle['v_x']     = A2_vel[:, 0]
-        A2_picle['v_y']     = A2_vel[:, 1]
-        A2_picle['v_z']     = A2_vel[:, 2]
-        A2_picle['m']       = f['PartType0/Masses'].value
-        A2_picle['ID']      = f['PartType0/ParticleIDs'].value
-        A2_picle['rho']     = f['PartType0/Density'].value
-        A2_picle['u']       = f['PartType0/InternalEnergy'].value
-        A2_picle['phi']     = f['PartType0/Potential'].value
-        A2_picle['P']       = f['PartType0/Pressure'].value
-        A2_picle['h']       = f['PartType0/SmoothingLength'].value
-        A2_picle['mat_ID']  = f['PartType0/MaterialID'].value
-
-    return A2_picle
-
-
-def process_particles(A2_picle, num_target):
-    """ Modify things like particle units, material IDs, and coordinate origins.
-    """
-    # Offset material IDs for impactor particles
-    A2_picle['mat_ID'][A2_picle['ID'] >= num_target] += ID_imp
-
-    # Shift coordinates to the centre of the target's core's mass and momentum
-    sel_tar  = np.where(A2_picle['mat_ID'] == Di_material['Til_iron'])[0]
-
-    # Centre of mass
-    m_tot   = np.sum(A2_picle[sel_tar]['m'])
-    x_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['x']) / m_tot
-    y_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['y']) / m_tot
-    z_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['z']) / m_tot
-
-    # Change origin to the centre-of-mass
-    A2_picle['x']   -= x_com
-    A2_picle['y']   -= y_com
-    A2_picle['z']   -= z_com
-    A2_picle['r']   = np.sqrt(
-        A2_picle['x']**2 + A2_picle['y']**2 + A2_picle['z']**2
-        )
-
-    # Centre of momentum
-    v_x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_x']) / m_tot
-    v_y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_y']) / m_tot
-    v_z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_z']) / m_tot
-
-    # Change to the centre-of-momentum frame of reference
-    A2_picle['v_x'] -= v_x_com
-    A2_picle['v_y'] -= v_y_com
-    A2_picle['v_z'] -= v_z_com
-
-    return A2_picle
-
-
-def plot_snapshot(A2_picle, filename, time, ax_lim=100, dz=0.1):
-    """ Plot the snapshot particles and save the figure.
-    """
-    # Add extension if needed
-    if (filename[-5:] != ".png"):
-        filename += ".png"
-
-    fig = plt.figure(figsize=(9, 9))
-    ax  = fig.add_subplot(111, aspect='equal')
-
-    # Plot slices in z below zero
-    for z in np.arange(-ax_lim, 0, dz):
-        sel_z       = np.where((z < A2_picle['z']) & (A2_picle['z'] < z+dz))[0]
-        A2_picle_z  = A2_picle[sel_z]
-
-        # Plot each material
-        for mat_ID, colour in Di_mat_colour.iteritems():
-            sel_col = np.where(A2_picle_z['mat_ID'] == mat_ID)[0]
-
-            ax.scatter(
-                A2_picle_z[sel_col]['x'], A2_picle_z[sel_col]['y'],
-                c=colour, edgecolors='none', marker='.', s=50, alpha=0.7
-                )
-
-    # Axes etc.
-    ax.set_axis_bgcolor('k')
-
-    ax.set_xlabel("x Position ($R_\oplus$)")
-    ax.set_ylabel("y Position ($R_\oplus$)")
-
-    ax.set_xlim(-ax_lim, ax_lim)
-    ax.set_ylim(-ax_lim, ax_lim)
-
-    plt.text(
-        -0.92*ax_lim, 0.85*ax_lim, "%.1f h" % (time*s_to_hour), fontsize=20,
-        color='w'
-        )
-
-    # Font sizes
-    for item in (
-        [ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() +
-        ax.get_yticklabels()
-        ):
-        item.set_fontsize(20)
-
-    plt.tight_layout()
-
-    plt.savefig(filename)
-    plt.close()
-
-
-if __name__ == '__main__':
-    # Sys args
-    try:
-        time_end    = int(sys.argv[1])
-
-        try:
-            delta_time  = int(sys.argv[2])
-        except IndexError:
-            delta_time  = 100
-    except IndexError:
-        time_end    = 100000
-        delta_time  = 100
-
-    # Load and plot each snapshot
-    for i_snap in range(int(time_end/delta_time) + 1):
-        snap_time   = i_snap * delta_time
-        print "\rPlotting snapshot %06d (%d of %d)" % (
-            snap_time, i_snap+1, int(time_end/delta_time)
-            ),
-        sys.stdout.flush()
-
-        # Load particle data
-        filename    = "%s%06d" % (file_snap, snap_time)
-        A2_picle    = load_snapshot(filename)
-
-        # Process particle data
-        A2_picle    = process_particles(A2_picle, num_target)
-
-        # Plot particles
-        filename    = "%s%06d" % (file_plot, snap_time)
-        plot_snapshot(A2_picle, filename, snap_time)
-
-    # Animation
-    command = (
-        "ffmpeg -framerate 12 -i plots/moon_forming_impact_%*.png -r 25 "
-        "anim.mpg -y"
-        )
-    print "\n%s\n" % command
-    subprocess.check_output(command, shell=True)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

examples/MoonFormingImpact/run.sh

-#!/bin/bash
-../swift -G -s -t 8 moon_forming_impact.yml

examples/SineWavePotential_1D/makeIC.py

@@ -31,12 +31,12 @@ cs2 = 2.*uconst/3.
 A = 10.
 
 fileName = "sineWavePotential.hdf5"
-numPart = 100
+numPart = 1000
 boxSize = 1.
 
 coords = np.zeros((numPart, 3))
 v = np.zeros((numPart, 3))
-m = np.zeros(numPart) + 1.
+m = np.zeros(numPart) + 1000. / numPart
 h = np.zeros(numPart) + 2./numPart
 u = np.zeros(numPart) + uconst
 ids = np.arange(numPart, dtype = 'L')

examples/SineWavePotential_1D/plotSolution.py

@@ -23,8 +23,13 @@
 import numpy as np
 import h5py
 import sys
+import matplotlib
+matplotlib.use("Agg")
 import pylab as pl
 
+pl.rcParams["figure.figsize"] = (12, 10)
+pl.rcParams["text.usetex"] = True
+
 # these should be the same as in makeIC.py
 uconst = 20.2615290634
 cs2 = 2.*uconst/3.
@@ -39,15 +44,20 @@ fileName = sys.argv[1]
 file = h5py.File(fileName, 'r')
 coords = np.array(file["/PartType0/Coordinates"])
 rho = np.array(file["/PartType0/Density"])
+P = np.array(file["/PartType0/Pressure"])
 u = np.array(file["/PartType0/InternalEnergy"])
-agrav = np.array(file["/PartType0/GravAcceleration"])
 m = np.array(file["/PartType0/Masses"])
+vs = np.array(file["/PartType0/Velocities"])
 ids = np.array(file["/PartType0/ParticleIDs"])
 
 x = np.linspace(0., 1., 1000)
 rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
 
-P = cs2*rho
+a = A * np.sin(2. * np.pi * x)
+
+apart = A * np.sin(2. * np.pi * coords[:,0])
+tkin = -0.5 * np.dot(apart, coords[:,0])
+print tkin, 0.5 * np.dot(m, vs[:,0]**2)
 
 ids_reverse = np.argsort(ids)
 
@@ -65,13 +75,38 @@ for i in range(len(P)):
     corr = 1.
   idxp1 = ids_reverse[ip1]
   idxm1 = ids_reverse[im1]
-  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
+  gradP[i] = (P[idxp1] - P[idxm1])/(coords[idxp1,0] - coords[idxm1,0])
+
+fig, ax = pl.subplots(2, 2, sharex = True)
+
+ax[0][0].plot(coords[:,0], rho, ".", label = "gizmo-mfm")
+ax[0][0].plot(x, rho_x, "-", label = "stable solution")
+ax[0][0].set_ylabel("$\\rho{}$ (code units)")
+ax[0][0].legend(loc = "best")
+
+ax[0][1].plot(x, a, "-", label = "$\\nabla{}\\Phi{}$ external")
+ax[0][1].plot(coords[:,0], gradP/rho, ".",
+              label = "$\\nabla{}P/\\rho{}$ gizmo-mfm")
+ax[0][1].set_ylabel("force (code units)")
+ax[0][1].legend(loc = "best")
+
+ax[1][0].axhline(y = uconst, label = "isothermal $u$", color = "k",
+                 linestyle = "--")
+ax[1][0].plot(coords[:,0], u, ".", label = "gizmo-mfm")
+ax[1][0].set_ylabel("$u$ (code units)")
+ax[1][0].set_xlabel("$x$ (code units)")
+ax[1][0].legend(loc = "best")
+
+#ax[1][1].plot(coords[:,0], m, "y.")
+#ax[1][1].set_ylabel("$m_i$ (code units)")
+#ax[1][1].set_xlabel("$x$ (code units)")
 
-fig, ax = pl.subplots(2, 2)
+ax[1][1].axhline(y = 0., label = "target", color = "k",
+                 linestyle = "--")
+ax[1][1].plot(coords[:,0], vs[:,0], ".", label = "gizmo-mfm")
+ax[1][1].set_ylabel("$v_x$ (code units)")
+ax[1][1].set_xlabel("$x$ (code units)")
+ax[1][1].legend(loc = "best")
 
-ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
-ax[0][0].plot(x, rho_x, "g-")
-ax[0][1].plot(coords[:,0], gradP/rho, "b.")
-ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
-ax[1][1].plot(coords[:,0], m, "y.")
+pl.tight_layout()
 pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))

examples/SineWavePotential_1D/sineWavePotential.yml

@@ -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-6  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-8  # 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).
 
 # Parameters governing the conserved quantities statistics
@@ -36,3 +36,6 @@ InitialConditions:
 SineWavePotential:
   amplitude: 10.
   timestep_limit: 1.
+
+EoS:
+  isothermal_internal_energy: 20.2615290634

examples/UranusImpact/README.md

-Uranus Giant Impact
-===================
-
-A simple version of the low angular momentum impact onto the early Uranus shown
-in Kegerreis et al. (2018), Fig. 2; with only ~10,000 particles for a quick and
-crude simulation.
-
-The collision of a 2 Earth mass impactor onto a proto-Uranus that can explain
-the spin of the present-day planet, with an angular momentum of 2e36 kg m^2 s^-1
-and velocity at infinity of 5 km s^-1 for a relatively head-on impact.
-
-Both bodies have a rocky core and icy mantle, with a hydrogen-helium atmosphere
-on the target as well. Although with this low number of particles it cannot be
-modelled in any detail.
-
-Setup
------
-
-In `swiftsim/`:
-
-`$ ./configure --with-hydro=minimal-multi-mat --with-equation-of-state=planetary`
-
-`$ make`
-
-In `swiftsim/examples/UranusImpact/`:
-
-`$ ./get_init_cond.sh`
-
-Run
----
-
-`$ ./run.sh`
-
-Analysis
---------
-
-`$ python plot.py`
-
-`$ mplayer anim.mpg`
-

examples/UranusImpact/get_init_cond.sh

-#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/uranus_impact.hdf5

examples/UranusImpact/plot.py

-"""
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-###############################################################################
-
-Plotting script for the Uranus Giant Impact example.
-
-Save a figure for each snapshot in `./plots/` then make them into a simple
-animation with ffmpeg in `./`.
-
-The snapshot plots show all particles with z < 0, coloured by their material.
-
-Usage:
-    `$ python  plot.py  time_end  delta_time`
-
-    Sys args:
-        + `time_end` | (opt) int | The time of the last snapshot to plot.
-            Default = 100000
-        + `delta_time` | (opt) int | The time between successive snapshots.
-            Default = 500
-"""
-
-from __future__ import division
-import numpy as np
-import matplotlib
-import matplotlib.pyplot as plt
-import h5py
-import sys
-import subprocess
-
-# Particle array fields
-dtype_picle = [
-    ('m', float), ('x', float), ('y', float), ('z', float), ('v_x', float),
-    ('v_y', float), ('v_z', float), ('ID', int), ('rho', float), ('u', float),
-    ('phi', float), ('P', float), ('h', float), ('mat_ID', int), ('r', float)
-    ]
-
-s_to_hour   = 1 / 60**2
-R_Ea        = 6.371e6
-
-# Default sys args
-time_end_default    = 100000
-delta_time_default  = 500
-
-# Snapshot info
-file_snap   = "./snapshots/uranus_impact_"
-file_plot   = "./plots/uranus_impact_"
-
-# Number of particles in the target body
-num_target  = 8992
-
-# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
-type_factor = 100
-Di_type = {
-    'Til'       : 1,
-    'HM80'      : 2,
-    'ANEOS'     : 3,
-    'SESAME'    : 4,
-}
-Di_material = {
-    # Tillotson
-    'Til_iron'      : Di_type['Til']*type_factor,
-    'Til_granite'   : Di_type['Til']*type_factor + 1,
-    'Til_water'     : Di_type['Til']*type_factor + 2,
-    # Hubbard & MacFarlane (1980) Uranus/Neptune
-    'HM80_HHe'      : Di_type['HM80']*type_factor,      # Hydrogen-helium atmosphere
-    'HM80_ice'      : Di_type['HM80']*type_factor + 1,  # H20-CH4-NH3 ice mix
-    'HM80_rock'     : Di_type['HM80']*type_factor + 2,  # SiO2-MgO-FeS-FeO rock mix
-    # ANEOS
-    'ANEOS_iron'        : Di_type['ANEOS']*type_factor,
-    'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
-    # SESAME
-    'SESAME_iron'   : Di_type['SESAME']*type_factor,
-}
-
-# Material offset for impactor particles
-ID_imp  = 10000
-# Material colours
-Di_mat_colour = {
-    # Target
-    Di_material['HM80_HHe']     : '#33DDFF',
-    Di_material['HM80_ice']     : 'lightsteelblue',
-    Di_material['HM80_rock']    : 'slategrey',
-    # Impactor
-    Di_material['HM80_ice'] + ID_imp    : '#A080D0',
-    Di_material['HM80_rock'] + ID_imp   : '#706050',
-    }
-
-
-def load_snapshot(filename):
-    """ Load the hdf5 snapshot file and return the structured particle array.
-    """
-    # Add extension if needed
-    if (filename[-5:] != ".hdf5"):
-        filename += ".hdf5"
-
-    # Load the hdf5 file
-    with h5py.File(filename, 'r') as f:
-        header      = f['Header'].attrs
-        A2_pos      = f['PartType0/Coordinates'].value
-        A2_vel      = f['PartType0/Velocities'].value
-
-        # Structured array of all particle data
-        A2_picle    = np.empty(header['NumPart_Total'][0],
-                               dtype=dtype_picle)
-
-        A2_picle['x']       = A2_pos[:, 0]
-        A2_picle['y']       = A2_pos[:, 1]
-        A2_picle['z']       = A2_pos[:, 2]
-        A2_picle['v_x']     = A2_vel[:, 0]
-        A2_picle['v_y']     = A2_vel[:, 1]
-        A2_picle['v_z']     = A2_vel[:, 2]
-        A2_picle['m']       = f['PartType0/Masses'].value
-        A2_picle['ID']      = f['PartType0/ParticleIDs'].value
-        A2_picle['rho']     = f['PartType0/Density'].value
-        A2_picle['u']       = f['PartType0/InternalEnergy'].value
-        A2_picle['phi']     = f['PartType0/Potential'].value
-        A2_picle['P']       = f['PartType0/Pressure'].value
-        A2_picle['h']       = f['PartType0/SmoothingLength'].value
-        A2_picle['mat_ID']  = f['PartType0/MaterialID'].value
-
-    return A2_picle
-
-
-def process_particles(A2_picle, num_target):
-    """ Modify things like particle units, material IDs, and coordinate origins.
-    """
-    # Offset material IDs for impactor particles
-    A2_picle['mat_ID'][A2_picle['ID'] >= num_target] += ID_imp
-
-    # Shift coordinates to the centre of the target's core's mass and momentum
-    sel_tar  = np.where(A2_picle['mat_ID'] == Di_material['HM80_rock'])[0]
-
-    # Centre of mass
-    m_tot   = np.sum(A2_picle[sel_tar]['m'])
-    x_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['x']) / m_tot
-    y_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['y']) / m_tot
-    z_com   = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['z']) / m_tot
-
-    # Change origin to the centre-of-mass
-    A2_picle['x']   -= x_com
-    A2_picle['y']   -= y_com
-    A2_picle['z']   -= z_com
-    A2_picle['r']   = np.sqrt(
-        A2_picle['x']**2 + A2_picle['y']**2 + A2_picle['z']**2
-        )
-
-    # Centre of momentum
-    v_x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_x']) / m_tot
-    v_y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_y']) / m_tot
-    v_z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_z']) / m_tot
-
-    # Change to the centre-of-momentum frame of reference
-    A2_picle['v_x'] -= v_x_com
-    A2_picle['v_y'] -= v_y_com
-    A2_picle['v_z'] -= v_z_com
-
-    return A2_picle
-
-
-def plot_snapshot(A2_picle, filename, time, ax_lim=13, dz=0.1):
-    """ Plot the snapshot particles and save the figure.
-    """
-    # Add extension if needed
-    if (filename[-5:] != ".png"):
-        filename += ".png"
-
-    fig = plt.figure(figsize=(9, 9))
-    ax  = fig.add_subplot(111, aspect='equal')
-
-    # Plot slices in z below zero
-    for z in np.arange(-ax_lim, 0, dz):
-        sel_z       = np.where((z < A2_picle['z']) & (A2_picle['z'] < z+dz))[0]
-        A2_picle_z  = A2_picle[sel_z]
-
-        # Plot each material
-        for mat_ID, colour in Di_mat_colour.iteritems():
-            sel_col = np.where(A2_picle_z['mat_ID'] == mat_ID)[0]
-
-            ax.scatter(
-                A2_picle_z[sel_col]['x'], A2_picle_z[sel_col]['y'],
-                c=colour, edgecolors='none', marker='.', s=50, alpha=0.7
-                )
-
-    # Axes etc.
-    ax.set_axis_bgcolor('k')
-
-    ax.set_xlabel("x Position ($R_\oplus$)")
-    ax.set_ylabel("y Position ($R_\oplus$)")
-
-    ax.set_xlim(-ax_lim, ax_lim)
-    ax.set_ylim(-ax_lim, ax_lim)
-
-    plt.text(
-        -0.92*ax_lim, 0.85*ax_lim, "%.1f h" % (time*s_to_hour), fontsize=20,
-        color='w'
-        )
-
-    # Font sizes
-    for item in (
-        [ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() +
-        ax.get_yticklabels()
-        ):
-        item.set_fontsize(20)
-
-    plt.tight_layout()
-
-    plt.savefig(filename)
-    plt.close()
-
-
-if __name__ == '__main__':
-    # Sys args
-    try:
-        time_end    = int(sys.argv[1])
-        try:
-            delta_time  = int(sys.argv[2])
-        except IndexError:
-            delta_time  = delta_time_default
-    except IndexError:
-        time_end    = time_end_default
-        delta_time  = delta_time_default
-
-    # Load and plot each snapshot
-    for i_snap in range(int(time_end/delta_time) + 1):
-        snap_time   = i_snap * delta_time
-        print "\rPlotting snapshot %06d (%d of %d)" % (
-            snap_time, i_snap+1, int(time_end/delta_time)
-            ),
-        sys.stdout.flush()
-
-        # Load particle data
-        filename    = "%s%06d" % (file_snap, snap_time)
-        A2_picle    = load_snapshot(filename)
-
-        # Process particle data
-        A2_picle    = process_particles(A2_picle, num_target)
-
-        # Plot particles
-        filename    = "%s%06d" % (file_plot, snap_time)
-        plot_snapshot(A2_picle, filename, snap_time)
-
-    # Animation
-    command = (
-        "ffmpeg -framerate 10 -i plots/uranus_impact_%*.png -r 25 anim.mpg -y"
-        )
-    print "\n$ %s\n" % command
-    subprocess.call(command, shell=True)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

examples/UranusImpact/run.sh

-#!/bin/bash
-../swift -G -s -t 8 uranus_impact.yml

examples/UranusImpact/uranus_impact.yml

-# Define the system of units to use internally.
-InternalUnitSystem:
-    UnitMass_in_cgs:        5.9724e27   # Grams
-    UnitLength_in_cgs:      6.371e8     # Centimeters
-    UnitVelocity_in_cgs:    6.371e8     # Centimeters per second
-    UnitCurrent_in_cgs:     1           # Amperes
-    UnitTemp_in_cgs:        1           # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
-    time_begin:     0                   # The starting time of the simulation (in internal units).
-    time_end:       100000              # The end time of the simulation (in internal units).
-    dt_min:         0.001               # The minimal time-step size of the simulation (in internal units).
-    dt_max:         100                 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
-                                        # Common part of the name of output files
-    basename:       snapshots/uranus_impact
-    time_first:     0                   # Time of the first output (in internal units)
-    delta_time:     500                 # Time difference between consecutive outputs (in internal units)
-    label_delta:    500                 # Integer increment between snapshot output labels
-
-# Parameters governing the conserved quantities statistics
-Statistics:
-    delta_time:     1000                # Time between statistics output
-
-# Parameters for the hydrodynamics scheme
-SPH:
-    resolution_eta:     1.2348          # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-    delta_neighbours:   0.1             # The tolerance for the targetted number of neighbours.
-    CFL_condition:      0.2             # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters for the self-gravity scheme
-Gravity:
-    eta:                    0.025       # Constant dimensionless multiplier for time integration.
-    theta:                  0.7         # Opening angle (Multipole acceptance criterion)
-    comoving_softening:     0.01        # Comoving softening length (in internal units).
-    max_physical_softening: 0.01        # Physical softening length (in internal units).
-
-# Parameters related to the initial conditions
-InitialConditions:
-    file_name:      uranus_impact.hdf5  # The initial conditions file to read
-
-# Parameters related to the equation of state
-EoS:
-    planetary_use_HM80:   1                       # Whether to prepare the Hubbard & MacFarlane (1980) EOS
-                                        # Table file paths
-    planetary_HM80_HHe_table_file:    /gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt
-    planetary_HM80_ice_table_file:    /gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt
-    planetary_HM80_rock_table_file:   /gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt

examples/ZeldovichPancake_3D/makeIC.py

@@ -28,29 +28,34 @@ z_c = 1.             # Redshift of caustic formation (non-linear collapse)
 z_i = 100.           # Initial redshift
 gamma = 5./3.        # Gas adiabatic index
 numPart_1D = 32      # Number of particles along each dimension
+fileName = "zeldovichPancake.hdf5"
 
 
 # Some units
-Mpc_in_m = 3.085678e22
-Msol_in_kg = 1.989e30
-Gyr_in_s = 3.085678e19
+Mpc_in_m = 3.08567758e22
+Msol_in_kg = 1.98848e30
+Gyr_in_s = 3.08567758e19
 mH_in_kg = 1.6737236e-27
-k_in_J_K = 1.38064852e-23
 
-# Parameters
-rho_0 = 1.8788e-26 # h^2 kg m^-3
-H_0 = 1. / Mpc_in_m * 10**5 # s^-1
+# Some constants
+kB_in_SI = 1.38064852e-23
+G_in_SI = 6.67408e-11
+
+# Some useful variables in h-full units
+H_0 = 1. / Mpc_in_m * 10**5 # h s^-1
+rho_0 = 3. * H_0**2 / (8* math.pi * G_in_SI) # h^2 kg m^-3
 lambda_i = 64. / H_0 * 10**5 # h^-1 m (= 64 h^-1 Mpc)
 x_min = -0.5 * lambda_i
 x_max = 0.5 * lambda_i
-fileName = "zeldovichPancake.hdf5"
 
+# SI system of units
 unit_l_in_si = Mpc_in_m
 unit_m_in_si = Msol_in_kg * 1.e10
 unit_t_in_si = Gyr_in_s
 unit_v_in_si = unit_l_in_si / unit_t_in_si
 unit_u_in_si = unit_v_in_si**2
 
+# Total number of particles
 numPart = numPart_1D**3
 
 #---------------------------------------------------
@@ -87,7 +92,7 @@ for i in range(numPart_1D):
       coords[index,1] = (j + 0.5) * delta_x
       coords[index,2] = (k + 0.5) * delta_x
       T = T_i * (1. / (1. - zfac * cos(k_i * q)))**(2. / 3.)
-      u[index] = k_in_J_K * T / (gamma - 1.) / mH_in_kg
+      u[index] = kB_in_SI * T / (gamma - 1.) / mH_in_kg
       h[index] = 1.2348 * delta_x
       m[index] = m_i
       v[index,0] = -H_0 * (1. + z_c) / sqrt(1. + z_i) * sin(k_i * q) / k_i
@@ -141,7 +146,7 @@ grp.create_dataset('ParticleIDs', data=ids, dtype='L')
 
 file.close()
 
-import pylab as pl
+#import pylab as pl
 
-pl.plot(coords[:,0], v[:,0], "k.")
-pl.show()
+#pl.plot(coords[:,0], v[:,0], "k.")
+#pl.show()