diff --git a/examples/InteractingBlastWaves_1D/interactingBlastWaves.yml b/examples/InteractingBlastWaves_1D/interactingBlastWaves.yml new file mode 100644 index 0000000000000000000000000000000000000000..e845599730828fd7b9880ae9aca11420ba50026c --- /dev/null +++ b/examples/InteractingBlastWaves_1D/interactingBlastWaves.yml @@ -0,0 +1,33 @@ +# 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: 3.8e-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-5 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: interactingBlastWaves # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-5 # 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). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./interactingBlastWaves.hdf5 # The file to read diff --git a/examples/InteractingBlastWaves_1D/makeIC.py b/examples/InteractingBlastWaves_1D/makeIC.py new file mode 100644 index 0000000000000000000000000000000000000000..bed0e20c833ccbe54ed571b954cad03ab93f4c0c --- /dev/null +++ b/examples/InteractingBlastWaves_1D/makeIC.py @@ -0,0 +1,87 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com) +# +# 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/>. +# +############################################################################## + +# Generates a regular 1D grid that contains the Woodward & Colella (1984) +# interacting blast waves test. +# Note that the original test requires reflective boundary conditions. Since we +# do not have those in SWIFT, we double the computational domain and mirror the +# setup in the second half. + +import numpy as np +import h5py + +fileName = "interactingBlastWaves.hdf5" +numPart = 800 +boxSize = 2. + +coords = np.zeros((numPart, 3)) +v = np.zeros((numPart, 3)) +m = np.zeros(numPart) + boxSize / numPart +h = np.zeros(numPart) + 2. * boxSize / numPart +u = np.zeros(numPart) +ids = np.arange(numPart, dtype = 'L') +rho = np.ones(numPart) + +for i in range(numPart): + coords[i,0] = (i + 0.5) * boxSize / numPart + if coords[i,0] < 0.1 or coords[i,0] > 1.9: + u[i] = 2500. + elif coords[i,0] > 0.9 and coords[i,0] < 1.1: + u[i] = 250. + else: + u[i] = 0.025 + +#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 +grp.attrs["Dimension"] = 1 + +#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') +grp.create_dataset('Density', data=rho, dtype='f') + +file.close() diff --git a/examples/InteractingBlastWaves_1D/plotSolution.py b/examples/InteractingBlastWaves_1D/plotSolution.py new file mode 100644 index 0000000000000000000000000000000000000000..342b67eae1061a9a94675a59f53ebf9e4ba7717f --- /dev/null +++ b/examples/InteractingBlastWaves_1D/plotSolution.py @@ -0,0 +1,37 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com) +# +# 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 matplotlib +matplotlib.use("Agg") +import pylab as pl +import sys + +snap = int(sys.argv[1]) + +sim = h5py.File("interactingBlastWaves_%04d.hdf5" % snap, "r") +coords = sim["/PartType0/Coordinates"] +rho = sim["/PartType0/Density"] + +pl.xlabel("$x$") +pl.ylabel("$\\rho{}$") +pl.xlim(0.4, 1.) +pl.ylim(0., 7.) +pl.plot(coords[:,0], rho, "k.") +pl.savefig("InteractingBlastWaves.png") diff --git a/examples/InteractingBlastWaves_1D/run.sh b/examples/InteractingBlastWaves_1D/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..65a1a109857da6e32663451430b3284aa7207d5a --- /dev/null +++ b/examples/InteractingBlastWaves_1D/run.sh @@ -0,0 +1,14 @@ +#!/bin/bash + + # Generate the initial conditions if they are not present. +if [ ! -e interactingBlastWaves.hdf5 ] +then + echo "Generating initial conditions for the Sedov blast example..." + python makeIC.py +fi + +# Run SWIFT +../swift -s -t 1 interactingBlastWaves.yml 2>&1 | tee output.log + +# Plot the solution +python plotSolution.py 4