Commit 6775e2f3 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'disc_patch_x' into 'master'

Disc patch x

Replaced disc patch potential in z with a disc patch potential in x. Added a 1D version of the disc patch test. Changed the GIZMO time step criterion (fixed a bug in the original criterion, and tried to make the criterion stricter to cope with the disc patch potential). Removed some gravity variables that were no longer used. Put the flux limiter in a separate file and made it configurable.

See merge request !382
parents b9389d4e 96c1763a
......@@ -39,8 +39,8 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 400.
z_trunc: 300.
z_max: 350.
x_disc: 400.
x_trunc: 300.
x_max: 350.
timestep_mult: 0.03
growth_time: 5.
......@@ -39,7 +39,7 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 400.
z_trunc: 300.
z_max: 380.
x_disc: 400.
x_trunc: 300.
x_max: 380.
timestep_mult: 0.03
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
# Tom Theuns (tom.theuns@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/>.
#
##############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
# Tom Theuns (tom.theuns@durham.ac.uk)
# 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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 sys
......@@ -37,16 +39,16 @@ import random
# Size of the patch -- side_length
# Parameters of the gas disc
surface_density = 10.
surface_density = 10.
scale_height = 100.
gas_gamma = 5./3.
# Parameters of the problem
z_factor = 2
x_factor = 2
side_length = 400.
# File
fileName = "Disc-Patch.hdf5"
fileName = "Disc-Patch.hdf5"
####################################################################
......@@ -64,30 +66,33 @@ unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
unit_velocity_in_cgs = (1e5)
unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
print "UnitLength_in_cgs: %.5e"%unit_length_in_cgs
print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
print "UnitTime_in_cgs: %.5e"%unit_time_in_cgs
print ""
# Derived units
const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * unit_length_in_cgs**-3
const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
unit_length_in_cgs**-3
const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * unit_time_in_cgs**2
const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
unit_time_in_cgs**2
print "--- Some constants [internal units] ---"
print "G_Newton: %.5e"%const_G
print "m_proton: %.5e"%const_mp
print "k_boltzmann: %.5e"%const_kb
print "G_Newton: %.5e"%const_G
print "m_proton: %.5e"%const_mp
print "k_boltzmann: %.5e"%const_kb
print ""
# derived quantities
temp = math.pi * const_G * surface_density * scale_height * const_mp / const_kb
u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
v_disp = math.sqrt(2 * u_therm)
soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
t_dyn = math.sqrt(scale_height / (const_G * surface_density))
t_cross = scale_height / soundspeed
temp = math.pi * const_G * surface_density * scale_height * const_mp / \
const_kb
u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
v_disp = math.sqrt(2 * u_therm)
soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
t_dyn = math.sqrt(scale_height / (const_G * surface_density))
t_cross = scale_height / soundspeed
print "--- Properties of the gas [internal units] ---"
print "Gas temperature: %.5e"%temp
......@@ -101,18 +106,20 @@ print ""
# Problem properties
boxSize_x = side_length
boxSize_y = boxSize_x
boxSize_z = boxSize_x * z_factor
boxSize_z = boxSize_x
boxSize_x *= x_factor
volume = boxSize_x * boxSize_y * boxSize_z
M_tot = boxSize_x * boxSize_y * surface_density * math.tanh(boxSize_z / (2. * scale_height))
M_tot = boxSize_y * boxSize_z * surface_density * \
math.tanh(boxSize_x / (2. * scale_height))
density = M_tot / volume
entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
print "--- Problem properties [internal units] ---"
print "Box: [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
print "Volume: %.5e"%volume
print "Total mass: %.5e"%M_tot
print "Density: %.5e"%density
print "Entropy: %.5e"%entropy
print "Box: [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
print "Volume: %.5e"%volume
print "Total mass: %.5e"%M_tot
print "Density: %.5e"%density
print "Entropy: %.5e"%entropy
print ""
####################################################################
......@@ -123,34 +130,24 @@ one_glass_pos = infile["/PartType0/Coordinates"][:,:]
one_glass_h = infile["/PartType0/SmoothingLength"][:]
# Rescale to the problem size
one_glass_pos *= boxSize_x
one_glass_h *= boxSize_x
#print min(one_glass_p[:,0]), max(one_glass_p[:,0])
#print min(one_glass_p[:,1]), max(one_glass_p[:,1])
#print min(one_glass_p[:,2]), max(one_glass_p[:,2])
one_glass_pos *= side_length
one_glass_h *= side_length
# Now create enough copies to fill the volume in z
# Now create enough copies to fill the volume in x
pos = np.copy(one_glass_pos)
h = np.copy(one_glass_h)
for i in range(1, z_factor):
one_glass_pos[:,2] += boxSize_x
for i in range(1, x_factor):
one_glass_pos[:,0] += side_length
pos = np.append(pos, one_glass_pos, axis=0)
h = np.append(h, one_glass_h, axis=0)
#print min(pos[:,0]), max(pos[:,0])
#print min(pos[:,1]), max(pos[:,1])
#print min(pos[:,2]), max(pos[:,2])
# Compute further properties of ICs
numPart = np.size(h)
mass = M_tot / numPart
print "--- Particle properties [internal units] ---"
print "Number part.: ", numPart
print "Part. mass: %.5e"%mass
print "Part. mass: %.5e"%mass
print ""
# Create additional arrays
......@@ -159,7 +156,6 @@ mass = np.ones(numPart) * mass
vel = np.zeros((numPart, 3))
ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
####################################################################
# Create and write output file
......@@ -169,7 +165,7 @@ file = h5py.File(fileName, 'w')
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
......@@ -200,13 +196,12 @@ ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
####################################################################
print "--- Runtime parameters (YAML file): ---"
print "DiscPatchPotential:surface_density: ", surface_density
print "DiscPatchPotential:scale_height: ", scale_height
print "DiscPatchPotential:z_disc: ", boxSize_z / 2.
print "DiscPatchPotential:x_disc: ", 0.5 * boxSize_x
print ""
print "--- Constant parameters: ---"
......
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
# 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
......@@ -34,9 +35,9 @@ import sys
# Parameters
surface_density = 10.
scale_height = 100.
z_disc = 400.
z_trunc = 300.
z_max = 350.
x_disc = 400.
x_trunc = 300.
x_max = 350.
utherm = 20.2678457288
gamma = 5. / 3.
......@@ -50,14 +51,14 @@ if len(sys.argv) > 2:
# Get the analytic solution for the density
def get_analytic_density(x):
return 0.5 * surface_density / scale_height / \
np.cosh( (x - z_disc) / scale_height )**2
np.cosh( (x - x_disc) / scale_height )**2
# Get the analytic solution for the (isothermal) pressure
def get_analytic_pressure(x):
return (gamma - 1.) * utherm * get_analytic_density(x)
# Get the data fields to plot from the snapshot file with the given name:
# snapshot time, z-coord, density, pressure, velocity norm
# snapshot time, x-coord, density, pressure, velocity norm
def get_data(name):
file = h5py.File(name, "r")
coords = np.array(file["/PartType0/Coordinates"])
......@@ -69,7 +70,7 @@ def get_data(name):
vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
return float(file["/Header"].attrs["Time"]), coords[:,2], rho, P, vtot
return float(file["/Header"].attrs["Time"]), coords[:,0], rho, P, vtot
# scan the folder for snapshot files and plot all of them (within the requested
# range)
......@@ -80,38 +81,38 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
print "processing", f, "..."
zrange = np.linspace(0., 2. * z_disc, 1000)
time, z, rho, P, v = get_data(f)
xrange = np.linspace(0., 2. * x_disc, 1000)
time, x, rho, P, v = get_data(f)
fig, ax = pl.subplots(3, 1, sharex = True)
ax[0].plot(z, rho, "r.")
ax[0].plot(zrange, get_analytic_density(zrange), "k-")
ax[0].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[0].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[0].set_ylim(0., 1.2 * get_analytic_density(z_disc))
ax[0].plot(x, rho, "r.")
ax[0].plot(xrange, get_analytic_density(xrange), "k-")
ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
ax[0].set_ylim(0., 1.2 * get_analytic_density(x_disc))
ax[0].set_ylabel("density")
ax[1].plot(z, v, "r.")
ax[1].plot(zrange, np.zeros(len(zrange)), "k-")
ax[1].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[1].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[1].plot(x, v, "r.")
ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
ax[1].set_ylim(-0.5, 10.)
ax[1].set_ylabel("velocity norm")
ax[2].plot(z, P, "r.")
ax[2].plot(zrange, get_analytic_pressure(zrange), "k-")
ax[2].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
ax[2].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
ax[2].set_xlim(0., 2. * z_disc)
ax[2].set_ylim(0., 1.2 * get_analytic_pressure(z_disc))
ax[2].set_xlabel("z")
ax[2].plot(x, P, "r.")
ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
ax[2].set_xlim(0., 2. * x_disc)
ax[2].set_ylim(0., 1.2 * get_analytic_pressure(x_disc))
ax[2].set_xlabel("x")
ax[2].set_ylabel("pressure")
pl.suptitle("t = {0:.2f}".format(time))
......
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e glassCube_32.hdf5 ]
then
echo "Fetching initial glass file for the disc patch example..."
./getGlass.sh
fi
if [ ! -e Disc-Patch.hdf5 ]
then
echo "Generating initial conditions for the disc patch example..."
python makeIC.py
fi
# Run SWIFT
../../swift -g -s -t 4 disc-patch-icc.yml 2>&1 | tee output.log
python plotSolution.py
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.08567758149e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # 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: 968. # The end time of the simulation (in internal units).
dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_max: 10. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 12. # Time between statistics output
# Parameters governing the snapshots
Snapshots:
basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 48. # Time difference between outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
h_max: 60. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions
InitialConditions:
file_name: Disc-Patch.hdf5 # The file to read
# External potential parameters
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
x_disc: 400.
x_trunc: 300.
x_max: 350.
timestep_mult: 0.03
growth_time: 5.
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
# Tom Theuns (tom.theuns@durham.ac.uk)
# 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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 sys
import numpy as np
import math
import random
# Generates a disc-patch in hydrostatic equilibrium
#
# See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
#
#
# Disc parameters are: surface density -- sigma
# scale height -- b
# gas adiabatic index -- gamma
#
# Problem parameters are: Ratio height/width of the box -- z_factor
# Size of the patch -- side_length
# Parameters of the gas disc
surface_density = 10.
scale_height = 100.
gas_gamma = 5./3.
# Parameters of the problem
x_factor = 2
side_length = 400.
numPart = 1000
# File
fileName = "Disc-Patch.hdf5"
####################################################################
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.08567758149e18
PROTON_MASS_IN_CGS = 1.672621898e-24
BOLTZMANN_IN_CGS = 1.38064852e-16
YEAR_IN_CGS = 3.15569252e7
# choice of units
unit_length_in_cgs = (PARSEC_IN_CGS)
unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
unit_velocity_in_cgs = (1e5)
unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
print "UnitLength_in_cgs: %.5e"%unit_length_in_cgs
print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
print "UnitTime_in_cgs: %.5e"%unit_time_in_cgs
print ""
# Derived units
const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
unit_length_in_cgs**-3
const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
unit_time_in_cgs**2
print "--- Some constants [internal units] ---"
print "G_Newton: %.5e"%const_G
print "m_proton: %.5e"%const_mp
print "k_boltzmann: %.5e"%const_kb
print ""
# derived quantities
temp = math.pi * const_G * surface_density * scale_height * const_mp / \
const_kb
u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
v_disp = math.sqrt(2 * u_therm)
soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
t_dyn = math.sqrt(scale_height / (const_G * surface_density))
t_cross = scale_height / soundspeed
print "--- Properties of the gas [internal units] ---"
print "Gas temperature: %.5e"%temp
print "Gas thermal_energy: %.5e"%u_therm
print "Dynamical time: %.5e"%t_dyn
print "Sound crossing time: %.5e"%t_cross
print "Gas sound speed: %.5e"%soundspeed
print "Gas 3D vel_disp: %.5e"%v_disp
print ""
# Problem properties
boxSize_x = side_length
boxSize_x *= x_factor
volume = boxSize_x
M_tot = surface_density * math.tanh(boxSize_x / (2. * scale_height))
density = M_tot / volume
entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
print "--- Problem properties [internal units] ---"
print "Box: %.1f"%boxSize_x
print "Volume: %.5e"%volume
print "Total mass: %.5e"%M_tot
print "Density: %.5e"%density
print "Entropy: %.5e"%entropy
print ""
####################################################################
# Now create enough copies to fill the volume in x
pos = np.zeros((numPart, 3))
h = np.zeros(numPart) + 2. * boxSize_x / numPart
for i in range(numPart):
pos[i, 0] = (i + 0.5) * boxSize_x / numPart
# Compute further properties of ICs
mass = M_tot / numPart
print "--- Particle properties [internal units] ---"
print "Number part.: ", numPart
print "Part. mass: %.5e"%mass
print ""
# Create additional arrays
u = np.ones(numPart) * u_therm
mass = np.ones(numPart) * mass
vel = np.zeros((numPart, 3))
ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
####################################################################
# Create and write output file
#File
file = h5py.File(fileName, 'w')
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = [boxSize_x, 1., 1.]
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, 0, 0, 0, 0, 0]
grp.attrs["Dimension"] = 1
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 1
# write gas particles
grp0 = file.create_group("/PartType0")
ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=pos)
ds = grp0.create_dataset('Velocities', (numPart, 3), 'f')
ds = grp0.create_dataset('Masses', (numPart,), 'f', data=mass)
ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)