Commit 0d96f291 authored by Rainer Weinberger's avatar Rainer Weinberger

added example interacting_blastwaves_1d

parent 24e430d1
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/interacting_blastwaves_1d/Config.sh
## config file for 1d interacting blastwaves probelm
#--------------------------------------- Basic operation mode of code
ONEDIMS # 1d simulation
REFLECTIVE_X=1 # X Boundary; 1: Reflective, 2: Inflow/Outflow; not active: periodic
GAMMA=1.4 # Adiabatic index of gas; 5/3 if not set
#--------------------------------------- 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
FORCE_EQUAL_TIMESTEPS # variable but global timestep
#---------------------------------------- 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_IN_DOUBLEPRECISION # snapshot files will be written in double precision
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
""" @package ./examples/interacting_blastwaves_1d/check.py
Code that checks results of 1d interacting blastwaves
created by Rainer Weinberger, last modified: 04.03.2019
"""
""" load libraries """
import sys # needed for exit codes
import numpy as np # scientific computing package
import h5py # hdf5 format
import matplotlib.pyplot as plt ## needs to be active for plotting!
createFigures = False
""" check functions """
def CheckL1Error(Pos, W, gamma, i_snap):
"""
Compare the hydrodynamical quantities of the simulation with the high
resolution simulation solution, calculate the L1 error and check whether
avarge L1 error is acceptable
\param[in] Pos: shape (n,) array with 1d positions
\param[in] W: shape (n,3) array with density, velocity and pressure of each cell
\param[in] gamma: adiabatic index of gas
\param[in] i_snap: index of snapshot
\return status: zero if everything is within tolerances, otherwise 1
"""
try:
xx, rho, v, pres = np.loadtxt("./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap ).T
except:
print("Check_L1_Error: could not find: ./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap)
return 1
W_exact = np.array([rho, v, pres]).T
## match indices
dx = xx[1] - xx[0]
index = np.array( (Pos - 0.5*dx) / dx, dtype=np.int32 )
residual_x = Pos - xx[index]
W_predict = np.zeros([len(index), 3], dtype=np.float64)
for k in np.arange(3):
W_predict[:,k] = (1.0- residual_x/dx) * W_exact[index,k] + residual_x/dx * W_exact[index+1,k]
norm = np.abs(W_predict)
norm[:,1] = 1 ## use absolute error in velocity component
## calculate L1 norm
delta = np.abs(W_predict - W) / norm
L1 = np.average(delta, axis=0)
## tolarance value; found empirically, fist order convergence!
val_max = 0.05 * 400.0 / np.float(Pos.shape[0])
L1MaxAllowed = np.array([val_max, 4.0*val_max, val_max])
if np.any(L1 > L1MaxAllowed):
print("CheckL1Error: ERROR: L1 error too large: %g %g %g; tolerance %g %g %g"% (L1[0], L1[1], L1[2], L1MaxAllowed[0], L1MaxAllowed[1], L1MaxAllowed[2]) )
return 1
else:
print("CheckL1Error: L1 error fine: %g %g %g; tolerance %g %g %g"% (L1[0], L1[1], L1[2], L1MaxAllowed[0], L1MaxAllowed[1], L1MaxAllowed[2]) )
return 0
def PlotSimulationData(Pos, W, gamma, i_snap, simulation_directory):
"""
Plot density, velocity, specific internal energy and pressure of
Simulation and high resolution solution result on top
\param[in] Pos: shape (n,) array with 1d positions
\param[in] W: shape (n,3) array with density, velocity and pressure of each cell
\param[in] gamma: adiabatic index of gas
\param[in] i_snap: index of snapshot
\param[in] simulation_directory: path to simulation
\return status: zero if everything went fine, 1: could not find exact solution data
"""
## set limits for plot
min = np.zeros(4)
max = np.zeros(4)
min[0] = np.min( W[:,0] ) - 0.3
max[0] = np.max( W[:,0] ) + 1.0
min[1] = np.min( W[:,1] ) - 0.3
max[1] = np.max( W[:,1] ) + 1.0
min[2] = np.min( W[:,2] / W[:,0] / (gamma - 1.0) ) - 30.0
max[2] = np.max( W[:,2] / W[:,0] / (gamma - 1.0) ) + 30.0
min[3] = np.min( W[:,2] ) - 30.0
max[3] = np.max( W[:,2] ) + 30.0
try:
## exact data, from 20,000 cell fixed grid simulation; downsampled by factor of 20
xx, rho, v, pres = np.loadtxt("./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap).T
except:
print("PlotSimulationData: could not find: ./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap)
return 1
## plot:
fig, ax = plt.subplots(4, 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].plot(Pos, W[:,0], ls="", marker='o')
ax[0].plot(xx, rho, color="k")
ax[1].plot(Pos, W[:,1], ls="", marker='o')
ax[1].plot(xx, v, color="k")
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), ls="", marker='o')
ax[2].plot(xx, pres / rho / (gamma - 1.0), color="k")
ax[3].plot(Pos, W[:,2], ls="", marker='o')
ax[3].plot(xx, pres, color="k")
for i_plot in np.arange(4):
ax[i_plot].set_ylim([ min[i_plot], max[i_plot] ])
ax[i_plot].set_xlim([0.5,0.9])
## set labels
ax[3].set_xlabel(r"pos")
ax[0].set_ylabel(r"density")
ax[1].set_ylabel(r"velocity")
ax[2].set_ylabel(r"spec. int. energy")
ax[3].set_ylabel(r"pressure")
fig.align_ylabels(ax[:])
fig.savefig(simulation_directory+"/output/figure_%03d.pdf" % (i_file))
plt.close(fig)
return 0
simulation_directory = str(sys.argv[1])
print("examples/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
##parameters
Dtype = np.float64 # double precision: np.float64, for single use np.float32
gamma = 1.4 ## needs to be identical to Config.sh =(5./3.) if not specified there!
""" analyze snap_001.hdf5 """
i_file = 1
ReturnFlag = 0
## read in data
directory = simulation_directory+"/output/"
filename = "snap_%03d.hdf5" % (i_file)
print("try to open "+directory+filename)
## open hdf5 file
try:
data = h5py.File(directory+filename, "r")
except:
pass
print("analyzing "+ filename)
## get data from snapshot
time = np.float( data["Header"].attrs["Time"] )
position = np.array(data["PartType0"]["Coordinates"], dtype = Dtype)
density = np.array(data["PartType0"]["Density"], dtype = Dtype)
vel = np.array(data["PartType0"]["Velocities"], dtype = Dtype)
internalEnergy = np.array(data["PartType0"]["InternalEnergy"], dtype = Dtype)
## convert to more useful data structure
W = np.array([density, vel[:,0], (gamma-1.0)*internalEnergy*density], dtype=Dtype).T ## shape: (n,3)
""" plot data if you want """
if createFigures:
ReturnFlag += PlotSimulationData(position[:,0], W, gamma, i_file, simulation_directory)
""" perform checks """
ReturnFlag += CheckL1Error(position[:,0], W, gamma, i_file)
if ReturnFlag == 0:
print("examples/interacting_blastwaves_1d/check.py: success!")
else:
print("examples/interacting_blastwaves_1d/check.py: failed!")
sys.exit(ReturnFlag)
""" @package examples/interacting_blastwaves_1d/create.py
Code that creates 1d interacting blastwave initial conditions
supposed to be as siple as possible and self-contained
created by Rainer Weinberger, 04.05.2019
"""
""" load libraries """
import sys # system calls
import numpy as np # scientific computing package
import h5py # hdf5 format
simulation_directory = str(sys.argv[1])
print("examples/interacting_blast_waves_1d/create.py: creating ICs in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
Boxsize = FloatType(1.0)
NumberOfCells = IntType(400)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
## Riemann problems
position_0 = 0.1 # initial position of left discontinuity
position_1 = 0.9 # initial position of right discontinuity
density_0 = 1.0
density_1 = 1.0
density_2 = 1.0
velocity_0 = 0.0
velocity_1 = 0.0
velocity_2 = 0.0
pressure_0 = 1000.0
pressure_1 = 0.1
pressure_2 = 100.0
gamma = 1.4 ## note: this has to be consistent with the parameter settings for Arepo!
uthermal_0 = pressure_0 / ( gamma - 1.0 ) / density_0
uthermal_1 = pressure_1 / ( gamma - 1.0 ) / density_1
uthermal_2 = pressure_2 / ( gamma - 1.0 ) / density_2
""" set up grid: uniform 1d grid """
## spacing
dx = Boxsize / FloatType(NumberOfCells)
## position of first and last cell
pos_first, pos_last = 0.5 * dx, Boxsize - 0.5 * dx
## set up grid
Pos = np.zeros([NumberOfCells, 3], dtype=FloatType)
Pos[:,0] = np.linspace(pos_first, pos_last, NumberOfCells, dtype=FloatType)
Volume = np.full(NumberOfCells, dx, dtype=FloatType)
""" set up hydrodynamical quantitites """
## fill with central state
Density = np.full(Pos.shape[0], density_1, dtype=FloatType)
Velocity = np.zeros(Pos.shape, dtype=FloatType)
Velocity[:,0] = velocity_1
Uthermal = np.full(Pos.shape[0], uthermal_1, dtype=FloatType)
i_left, = np.where( Pos[:,0] < position_0)
i_right, = np.where( Pos[:,0] > position_1)
## left and right states
Density[i_left] = density_0
Velocity[i_left,0] = velocity_0
Uthermal[i_left] = uthermal_0
Density[i_right] = density_2
Velocity[i_right,0] = velocity_2
Uthermal[i_right] = uthermal_2
## mass instead of density needed for input
Mass = Density * Volume
""" write *.hdf5 file; minimum number of fields required by Arepo """
IC = h5py.File(FilePath, 'w') # open/create file
## create hdf5 groups
header = IC.create_group("Header") # create header group
part0 = IC.create_group("PartType0") # create particle group for gas cells
## write 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)
## write cell data
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()
""" normal exit """
sys.exit(0)
%% examples/interacting_blastwaves_1d/param.txt
% parameter file for 1d interacting blast wave problem
%% system options
MaxMemSize 2500
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
%% initial conditions
InitCondFile ./run/examples/interacting_blastwaves_1d/IC
ICFormat 3
%% output options
OutputDir ./run/examples/interacting_blastwaves_1d/output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
NumFilesWrittenInParallel 1
%% resubmit opitions
ResubmitOn 0
ResubmitCommand my-scriptfile
OutputListFilename ol
OutputListOn 0
%% simulation mode
CoolingOn 0
StarformationOn 0
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
%% Cosmological parameters
Omega0 0.0
OmegaBaryon 0.0
OmegaLambda 0.0
HubbleParam 1.0
%% Simulation parameters
BoxSize 1.0
TimeOfFirstSnapshot 0.0
TimeBetStatistics 0.019
TimeBegin 0.0
TimeMax 0.038
TimeBetSnapshot 0.038
%% Units
UnitVelocity_in_cm_per_s 1.0
UnitLength_in_cm 1.0
UnitMass_in_g 1.0
GravityConstantInternal 0.0
%% settings for gravity
ErrTolIntAccuracy 0.1
ErrTolTheta 0.1
ErrTolForceAcc 0.1
%% timestepping
MaxSizeTimestep 0.019
MinSizeTimestep 1e-10
%% moving mesh
CellShapingSpeed 0.5
CellMaxAngleFactor 2.25
TypeOfTimestepCriterion 0
TypeOfOpeningCriterion 1
%% hydrodynamics
CourantFac 0.3
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 64
MaxNumNgbDeviation 2
InitGasTemp 0.0
MinGasTemp 0.0
MinEgySpec 0.0
MinimumDensityOnStartUp 0.0
%% domain decomposition
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.01
%% gravitational softening
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 0
SofteningTypeOfPartType2 0
SofteningTypeOfPartType3 0
SofteningTypeOfPartType4 0
SofteningTypeOfPartType5 0
This diff is collapsed.
......@@ -15,6 +15,7 @@ TESTS=""
## available 1d examples
TESTS+="wave_1d "
TESTS+="shocktube_1d "
TESTS+="interacting_blastwaves_1d "
## available 2d examples
TESTS+="gresho_2d "
......
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