Commit 6db14388 authored by Rainer Weinberger's avatar Rainer Weinberger
Browse files

added example noh_3d

parent e8b51260
......@@ -50,7 +50,6 @@ while True:
break
""" get simulation data """
time = FloatType(data["Header"].attrs["Time"])
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
......@@ -74,7 +73,6 @@ while True:
i_postshock, = np.where(Radius < r_shock)
Density_ref[i_postshock] = 16.0 ## for 2d
#### plot density profile
if createFigures:
fig, ax = plt.subplots(3, sharex=True, figsize=np.array([6.9,6.0]) )
......@@ -109,7 +107,6 @@ while True:
sys.exit(1)
i_file += 1
## if everything is ok
sys.exit(0)
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/noh_3d/Config.sh
## config file for 3d Noh probelm
#--------------------------------------- Basic operation mode of code
REFLECTIVE_X=2 # in-/outflow boundary conditions in x direction
REFLECTIVE_Y=2 # in-/outflow boundary conditions in y direction
REFLECTIVE_Z=2 # in-/outflow boundary conditions in z direction
#--------------------------------------- Mesh motion and regularization
REGULARIZE_MESH_CM_DRIFT # Mesh regularization; Move mesh generating point towards center of mass to make cells rounder.
REGULARIZE_MESH_CM_DRIFT_USE_SOUNDSPEED # Limit mesh regularization speed by local sound speed
REGULARIZE_MESH_FACE_ANGLE # Use maximum face angle as roundness criterion in mesh regularization
#--------------------------------------- Time integration options
TREE_BASED_TIMESTEPS # non-local timestep criterion (take 'signal speed' into account)
#---------------------------------------- Single/Double Precision
DOUBLEPRECISION=1 # Mode of double precision: not defined: single; 1: full double precision 2: mixed, 3: mixed, fewer single precisions; unless short of memory, use 1.
INPUT_IN_DOUBLEPRECISION # initial conditions are in double precision
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired; should this be standard?
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps, should this be standard?
""" @package ./examples/noh_3d/check.py
Code that checks results of 3d Noh problem
created by Rainer Weinberger, last modified 24.02.2019
"""
""" load libraries """
import sys ## load sys; needed for exit codes
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff
createFigures = True
simulation_directory = str(sys.argv[1])
print("examples/noh_3d/check.py: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
## open initial conditiions to get parameters
try:
data = h5py.File(simulation_directory + "/IC.hdf5", "r")
except:
print("could not open initial conditions!")
exit(1)
Boxsize = FloatType(data["Header"].attrs["BoxSize"])
NumberOfCells = IntType(data["Header"].attrs["NumPart_Total"][0])
CellsPerDimension = NumberOfCells**(1./3.) ## 3d sim
## parameters for initial state
density_0 = 1.0
velocity_radial_0 = -1.0 ## radial inflow velocity
pressure_0 = 1.0e-4
gamma = 5./3. ## note: this has to be consistent with the parameter settings for Arepo!
utherm_0 = pressure_0 / ( gamma - 1.0 ) / density_0
## maximum L1 error after one propagation; empirically based
DeltaMaxAllowed = 0.25
""" loop over all output files """
i_file = 0
while True:
""" try to read in snapshot """
directory = simulation_directory+"/output/"
filename = "snap_%03d.hdf5" % (i_file)
try:
data = h5py.File(directory+filename, "r")
except:
break
""" get simulation data """
time = FloatType(data["Header"].attrs["Time"])
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Density = np.array(data["PartType0"]["Density"], dtype = FloatType)
Mass = np.array(data["PartType0"]["Masses"], dtype = FloatType)
Velocity = np.array(data["PartType0"]["Velocities"], dtype = FloatType)
Uthermal = np.array(data["PartType0"]["InternalEnergy"], dtype = FloatType)
CellVolume = Mass / Density
xPosFromCenter = (Pos[:,0] - 0.5 * Boxsize)
yPosFromCenter = (Pos[:,1] - 0.5 * Boxsize)
zPosFromCenter = (Pos[:,2] - 0.5 * Boxsize)
Radius = np.sqrt( xPosFromCenter**2 + yPosFromCenter**2 + zPosFromCenter**2 )
vRad = Velocity[:,0] * xPosFromCenter / Radius \
+ Velocity[:,1] * yPosFromCenter / Radius \
+ Velocity[:,2] * zPosFromCenter / Radius
r_shock = 1./3. * time
## exact solution, 2d
Radius[Radius == 0] = 1e-10
Density_ref = density_0 * (1.0 + time/Radius) * (1.0 + time/Radius) ## in 3d
i_postshock, = np.where(Radius < r_shock)
Density_ref[i_postshock] = 64.0 ## for 3d
#### plot profiles
if createFigures:
fig, ax = plt.subplots(3, sharex=True, figsize=np.array([6.9,6.0]) )
fig.subplots_adjust(left = 0.13, bottom = 0.09,right = 0.98, top = 0.98)
ax[0].scatter(Radius, Density, rasterized=True, s=0.2)
ax[1].scatter(Radius, vRad, rasterized=True, s=0.2)
ax[2].scatter(Radius, Uthermal, rasterized=True, s=0.2)
i_sorted = np.argsort(Radius)
ax[0].plot(Radius[i_sorted], Density_ref[i_sorted], label="analytic solution")
ax[0].set_xlim(0,0.8)
ax[0].set_ylim(0,70)
ax[2].set_xlabel(r"radius")
ax[0].set_ylabel(r"density")
ax[1].set_ylabel(r"radial velocity")
ax[2].set_ylabel(r"spec. internal energy")
fig.align_ylabels(ax[:])
fig.savefig(directory+"/figure_%03d.pdf"%i_file)
#### check against analytic solution
i_compare, = np.where(Radius < 0.8) ## only check inner region; boundary has spourious effects in this testcase
abs_delta_dens = np.abs(Density[i_compare] - Density_ref[i_compare]) / Density_ref[i_compare]
L1_dens = np.average(abs_delta_dens, weights=CellVolume[i_compare] )
L1_max = DeltaMaxAllowed * time
print("examples/Noh_3d/check snapshot %d: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
if L1_dens > L1_max:
sys.exit(1)
i_file += 1
## if everything is ok
sys.exit(0)
""" @package ./examples/noh_3d/create.py
Code that creates 3d Noh test problem initial conditions
created by Rainer Weinberger, last modified 24.02.2019
"""
""" load libraries """
import sys # system specific calls
import numpy as np # scientific computing package
import h5py # hdf5 format
simulation_directory = str(sys.argv[1])
print("examples/noh_3d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
Boxsize = FloatType(6.0)
CellsPerDimension = IntType(30)
NumberOfCells = CellsPerDimension * CellsPerDimension * CellsPerDimension
## initial state
density_0 = 1.0
velocity_radial_0 = -1.0 ## radial inflow velocity
pressure_0 = 1.0e-4
gamma = 5./3. ## note: this has to be consistent with the parameter settings for Arepo!
utherm_0 = pressure_0 / ( gamma - 1.0 ) / density_0
""" set up grid: cartesian 3d grid """
## spacing
dx = Boxsize / FloatType(CellsPerDimension)
## position of first and last cell
pos_first, pos_last = 0.5 * dx, Boxsize - 0.5 * dx
## set up grid
Grid1d = np.linspace(pos_first, pos_last, CellsPerDimension, dtype=FloatType)
xx, yy, zz = np.meshgrid(Grid1d, Grid1d, Grid1d)
Pos = np.zeros([NumberOfCells, 3], dtype=FloatType)
Pos[:,0] = xx.reshape(NumberOfCells)
Pos[:,1] = yy.reshape(NumberOfCells)
Pos[:,2] = zz.reshape(NumberOfCells)
## calculate distance from center
xPosFromCenter = (Pos[:,0] - 0.5 * Boxsize)
yPosFromCenter = (Pos[:,1] - 0.5 * Boxsize)
zPosFromCenter = (Pos[:,2] - 0.5 * Boxsize)
Radius = np.sqrt( xPosFromCenter**2 + yPosFromCenter**2 + zPosFromCenter**2 )
""" set up hydrodynamical quantitites """
## mass insetad of density
Mass = np.full(NumberOfCells, density_0*dx*dx*dx, dtype=FloatType)
## velocity
Velocity = np.zeros([NumberOfCells,3], dtype=FloatType)
Velocity[:,0] = velocity_radial_0 * xPosFromCenter / Radius
Velocity[:,1] = velocity_radial_0 * yPosFromCenter / Radius
Velocity[:,2] = velocity_radial_0 * zPosFromCenter / Radius
## specific internal energy
Uthermal = np.full(NumberOfCells, utherm_0, dtype=FloatType)
""" write *.hdf5 file; minimum number of fields required by Arepo """
IC = h5py.File(simulation_directory+'/IC.hdf5', 'w')
## create hdf5 groups
header = IC.create_group("Header")
part0 = IC.create_group("PartType0")
## header entries
NumPart = np.array([NumberOfCells, 0, 0, 0, 0, 0], dtype=IntType)
header.attrs.create("NumPart_ThisFile", NumPart)
header.attrs.create("NumPart_Total", NumPart)
header.attrs.create("NumPart_Total_HighWord", np.zeros(6, dtype=IntType) )
header.attrs.create("MassTable", np.zeros(6, dtype=IntType) )
header.attrs.create("Time", 0.0)
header.attrs.create("Redshift", 0.0)
header.attrs.create("BoxSize", Boxsize)
header.attrs.create("NumFilesPerSnapshot", 1)
header.attrs.create("Omega0", 0.0)
header.attrs.create("OmegaB", 0.0)
header.attrs.create("OmegaLambda", 0.0)
header.attrs.create("HubbleParam", 1.0)
header.attrs.create("Flag_Sfr", 0)
header.attrs.create("Flag_Cooling", 0)
header.attrs.create("Flag_StellarAge", 0)
header.attrs.create("Flag_Metals", 0)
header.attrs.create("Flag_Feedback", 0)
if Pos.dtype == np.float64:
header.attrs.create("Flag_DoublePrecision", 1)
else:
header.attrs.create("Flag_DoublePrecision", 0)
## copy datasets
part0.create_dataset("ParticleIDs", data=np.arange(1, NumberOfCells+1) )
part0.create_dataset("Coordinates", data=Pos)
part0.create_dataset("Masses", data=Mass)
part0.create_dataset("Velocities", data=Velocity)
part0.create_dataset("InternalEnergy", data=Uthermal)
## close file
IC.close()
%% examples/noh_3d/param.txt
% parameter file for 3d Noh problem
InitCondFile ./run/examples/noh_3d/IC
ICFormat 3
OutputDir ./run/examples/noh_3d/output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
NumFilesWrittenInParallel 1
ResubmitOn 0
ResubmitCommand my-scriptfile
OutputListFilename ol
OutputListOn 0
CoolingOn 0
StarformationOn 0
Omega0 0.0
OmegaBaryon 0.0
OmegaLambda 0.0
HubbleParam 1.0
BoxSize 6.0
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
MaxMemSize 2500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 0.005
TimeBegin 0.0
TimeMax 2.0
TimeBetSnapshot 0.5
UnitVelocity_in_cm_per_s 1.0
UnitLength_in_cm 1.0
UnitMass_in_g 1.0
GravityConstantInternal 0.0
ErrTolIntAccuracy 0.1
ErrTolTheta 0.1
ErrTolForceAcc 0.1
MaxSizeTimestep 0.5
MinSizeTimestep 1e-5
CourantFac 0.3
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 64
MaxNumNgbDeviation 2
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.1
TypeOfTimestepCriterion 0
TypeOfOpeningCriterion 1
GasSoftFactor 0.01
SofteningComovingType0 0.1
SofteningComovingType1 0.1
SofteningComovingType2 0.1
SofteningComovingType3 0.1
SofteningComovingType4 0.1
SofteningComovingType5 0.1
SofteningMaxPhysType0 0.1
SofteningMaxPhysType1 0.1
SofteningMaxPhysType2 0.1
SofteningMaxPhysType3 0.1
SofteningMaxPhysType4 0.1
SofteningMaxPhysType5 0.1
SofteningTypeOfPartType0 0
SofteningTypeOfPartType1 1
SofteningTypeOfPartType2 1
SofteningTypeOfPartType3 1
SofteningTypeOfPartType4 1
SofteningTypeOfPartType5 1
InitGasTemp 0.0
MinGasTemp 0.0
MinEgySpec 0.0
MinimumDensityOnStartUp 0.0
CellShapingSpeed 0.5
CellMaxAngleFactor 2.25
......@@ -21,6 +21,9 @@ TESTS+="gresho_2d "
TESTS+="noh_2d "
TESTS+="yee_2d "
## available 3d examples
TESTS+="noh_3d "
## loop over all tests
for TEST in $TESTS
do
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment