Added a 1D Sod shock test-case to the example suite. Includes a plotting script.

examples/SodShock_1D/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 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
+import random
+from numpy import *
+
+# Generates a swift IC file for the Sod Shock in a periodic box
+
+# Parameters
+gamma = 5./3.          # Gas adiabatic index
+numPart_L = 800        # Number of particles in the left state
+x_min = -1.
+x_max = 1.
+rho_L = 1.             # Density left state
+rho_R = 0.125          # Density right state
+v_L = 0.               # Velocity left state
+v_R = 0.               # Velocity right state
+P_L = 1.               # Pressure left state
+P_R = 0.1              # Pressure right state
+fileName = "sodShock.hdf5" 
+
+
+#---------------------------------------------------
+
+# Find how many particles we actually have
+boxSize = x_max - x_min
+numPart_R = int(numPart_L * (rho_R / rho_L))
+numPart = numPart_L + numPart_R
+
+# Now get the distances
+delta_L = (boxSize/2)  / numPart_L
+delta_R = (boxSize/2)  / numPart_R
+offset_L = delta_L / 2
+offset_R = delta_R / 2
+
+# Build the arrays
+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_L):
+    coords[i,0] = x_min + offset_L + i * delta_L
+    u[i] = P_L / (rho_L * (gamma - 1.))
+    h[i] = 1.2348 * delta_L
+    m[i] = boxSize * rho_L / (2. * numPart_L)
+    
+# Set the particles on the right
+for j in range(numPart_R):
+    i = numPart_L + j
+    coords[i,0] = offset_R + j * delta_R
+    u[i] = P_R / (rho_R * (gamma - 1.))
+    h[i] = 1.2348 * delta_R
+    m[i] = boxSize * rho_R / (2. * numPart_R)
+
+# Shift particles
+coords[:,0] -= x_min
+    
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["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
+
+#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)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+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')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+
+
+file.close()
+
+

examples/SodShock_1D/plotSolution.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 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 Sod shock and plots the SPH answer
+ 
+
+# Generates the analytical  solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
+# entropic function on N points between x_min and x_max.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3.      # Polytropic index
+rho_L = 1.             # Density left state
+rho_R = 0.125          # Density right state
+v_L = 0.               # Velocity left state
+v_R = 0.               # Velocity right state
+P_L = 1.               # Pressure left state
+P_R = 0.1              # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# 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.045,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.05,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'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']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][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"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+N = 1000  # Number of points
+x_min = -1.
+x_max = 1.
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+c_L = sqrt(gas_gamma * P_L / rho_L)   # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R)   # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+    u = P_3 / P
+    if u > 1:
+        term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+        term2 = sqrt(2./term1)
+        fp = (u - 1.)*c*term2
+        dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+    else:
+        fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+        dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+    return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009) 
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+    P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+    P_3 = 0.5*(P_R + P_L)
+    f_L = 1.
+    while fabs(P_3 - P_new) > 1e-6:
+        P_3 = P_new
+        (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+        (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+        f = f_L + f_R + (v_R - v_L)
+        df = dfdp_L + dfdp_R
+        dp =  -f/df
+        prnew = P_3 + dp
+    v_3 = v_L - f_L
+    return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+    v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+    v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+    v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+    v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+    x_12 = -fabs(v_left) * time
+else:
+    x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+    x_56 = fabs(v_right) * time
+else:
+    x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+    if x_s[i] <= x_12:
+        rho_s[i] = rho_L
+        P_s[i] = P_L
+        v_s[i] = v_L
+    if x_s[i] >= x_12 and x_s[i] < x_23:
+        if shock_L:
+            rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+            P_s[i] = P_3
+            v_s[i] = v_3
+        else:
+            rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+            P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+            v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+    if x_s[i] >= x_23 and x_s[i] < x_34:
+        if shock_L:
+            rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+        else:
+            rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+        P_s[i] = P_3
+        v_s[i] = v_3
+    if x_s[i] >= x_34 and x_s[i] < x_45:
+        if shock_R:
+            rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+        else:
+            rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+        P_s[i] = P_3
+        v_s[i] = v_3
+    if x_s[i] >= x_45 and x_s[i] < x_56:
+        if shock_R:
+            rho_s[i] = rho_R
+            P_s[i] = P_R
+            v_s[i] = v_R
+        else:
+            rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+            P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+            v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+    if x_s[i] >= x_56:
+        rho_s[i] = rho_R
+        P_s[i] = P_R
+        v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+        
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r')
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("$x$", labelpad=0)
+ylabel("$v$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r')
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("$x$", labelpad=0)
+ylabel("$\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r')
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("$x$", labelpad=0)
+ylabel("$P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r')
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("$x$", labelpad=0)
+ylabel("$u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r')
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("$x$", labelpad=0)
+ylabel("$S$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sod shock with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, scheme, fontsize=10)
+text(-0.49, 0.4, kernel, fontsize=10)
+text(-0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("SodShock.png", dpi=200)

examples/SodShock_1D/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+    echo "Generating initial conditions for the 1D SodShock example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 sodShock.yml
+
+# Plot the result
+python plotSolution.py 1

examples/SodShock_1D/sodShock.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # 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:   0.2   # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # 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 snapshots
+Snapshots:
+  basename:            sodShock # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          0.2     # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-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 (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  max_smoothing_length:  0.4      # Maximal smoothing length allowed (in internal units).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./sodShock.hdf5       # The file to read
+