Skip to content
Snippets Groups Projects
Commit cdeb54c5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated ICs for the Gresho-Chan vortex test-case

parent 849ad16a
No related branches found
No related tags found
1 merge request!216Hydrodynamics in 2D + New test cases
......@@ -21,93 +21,83 @@
import h5py
import random
from numpy import *
import sys
# Generates a swift IC file for the Gresho Vortex in a periodic box
# Parameters
inputFile = "glass_256.hdf5"
periodic= 1 # 1 For periodic box
factor = 3
boxSize = [ 1.0 , 1.0, 1.0/factor ]
L = 120 # Number of particles along one axis
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
rho = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
fileName = "greshoVortex.hdf5"
vol = boxSize[0] * boxSize[1] * boxSize[2]
fileOutputName = "greshoVortex.hdf5"
#---------------------------------------------------
numPart = L**3 / factor
mass = boxSize[0]*boxSize[1]*boxSize[2] * rho / numPart
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype='L')
partId=0
for i in range(L):
for j in range(L):
for k in range(L/factor):
index = i*L*L/factor + j*L/factor + k
x = i * boxSize[0] / L + boxSize[0] / (2*L)
y = j * boxSize[0] / L + boxSize[0] / (2*L)
z = k * boxSize[0] / L + boxSize[0] / (2*L)
r2 = (x - boxSize[0] / 2)**2 + (y - boxSize[1] / 2)**2
r = sqrt(r2)
coords[index,0] = x
coords[index,1] = y
coords[index,2] = z
v_phi = 0.
if r < 0.2:
v_phi = 5.*r
elif r < 0.4:
v_phi = 2. - 5.*r
else:
v_phi = 0.
v[index,0] = -v_phi * (y - boxSize[0] / 2) / r
v[index,1] = v_phi * (x - boxSize[0] / 2) / r
v[index,2] = 0.
m[index] = mass
h[index] = eta * boxSize[0] / L
P = P0
if r < 0.2:
P = P + 5. + 12.5*r2
elif r < 0.4:
P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
else:
P = P + 3. + 4.*log(2.)
u[index] = P / ((gamma - 1.)*rho)
ids[index] = partId + 1
partId = partId + 1
fileInput = h5py.File(inputFile, "r")
header = fileInput["/Header"]
boxSize = header.attrs["BoxSize"][0]
coords = fileInput["/PartType0/Coordinates"][:,:]
h = fileInput["/PartType0/SmoothingLength"][:]
ids = fileInput["/PartType0/ParticleIDs"][:]
numPart = size(ids)
m = ones(shape(ids)) * boxSize**2 / numPart
u = zeros(shape(m))
v = zeros(shape(coords))
for i in range(numPart):
x = coords[i,0]
y = coords[i,1]
z = coords[i,2]
r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
r = sqrt(r2)
v_phi = 0.
if r < 0.2:
v_phi = 5.*r
elif r < 0.4:
v_phi = 2. - 5.*r
else:
v_phi = 0.
v[i,0] = -v_phi * (y - boxSize / 2) / r
v[i,1] = v_phi * (x - boxSize / 2) / r
v[i,2] = 0.
P = P0
if r < 0.2:
P = P + 5. + 12.5*r2
elif r < 0.4:
P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
else:
P = P + 3. + 4.*log(2.)
u[i] = P / ((gamma - 1.)*rho)
#File
file = h5py.File(fileName, 'w')
fileOutput = h5py.File(fileOutputName, 'w')
# Header
grp = file.create_group("/Header")
grp = fileOutput.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["NumFileOutputsPerSnapshot"] = 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]
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp = fileOutput.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp = fileOutput.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.
......@@ -115,20 +105,20 @@ 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 = fileOutput.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
ds = grp.create_dataset('Masses', (numPart, 1), 'f')
ds[()] = m.reshape((numPart,1))
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
ds[()] = h.reshape((numPart,1))
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds[()] = u.reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
ds[()] = ids[:]
ds[()] = ids.reshape((numPart,1))
file.close()
fileOutput.close()
###############################################################################
# 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 sys
import random
from numpy import *
# Generates a swift IC file containing a perturbed cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1.
L = int(sys.argv[1]) # Number of particles along one axis
rho = 1. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
pert = 0.1 # Perturbation scale (in units of the interparticle separation)
fileName = "perturbedPlane.hdf5"
#---------------------------------------------------
numPart = L**2
mass = boxSize**2 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
#Generate particles
coords = zeros((numPart, 3))
v = zeros((numPart, 3))
m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype='L')
for i in range(L):
for j in range(L):
index = i*L + j
x = i * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
y = j * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
z = 0
coords[index,0] = x
coords[index,1] = y
coords[index,2] = z
v[index,0] = 0.
v[index,1] = 0.
v[index,2] = 0.
m[index] = mass
h[index] = 1.23485 * boxSize / L
u[index] = internalEnergy
ids[index] = index
#--------------------------------------------------
#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, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total"] = numPart
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#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")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids + 1
file.close()
# 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: 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_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: perturbedPlane # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # 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.1 # 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: ./perturbedPlane.hdf5 # The file to read
......@@ -37,8 +37,8 @@
#define const_max_u_change 0.1f
/* Dimensionality of the problem */
#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_3D
#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment