Commit c94143dd authored by Rainer Weinberger's avatar Rainer Weinberger
Browse files

example wave_1d

parent b1f3f617
#!/bin/bash # this line only there to enable syntax highlighting in this file
## ./examples/wave_1d/Config.sh
## config file for 1d wave propagation probelm
#--------------------------------------- Basic operation mode of code
ONEDIMS # 1d simulation
#--------------------------------------- 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/wave_1d/check.py
Code that checks results of 1d wave propagation problem
created by Rainer Weinberger, last modified 19.2.2019 -- comments welcome
"""
""" 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/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32 # integer type
""" 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 = np.int32(data["Header"].attrs["NumPart_Total"][0])
""" maximum L1 error after one propagation; empirically motivated """
DeltaMaxAllowed = 5e-5 * FloatType(NumberOfCells)**-2
""" initial state -- copied from create.py """
density_0 = FloatType(1.0)
velocity_0 = FloatType(0.0)
pressure_0 = FloatType(3.0) / FloatType(5.0)
gamma = FloatType(5.0) / FloatType(3.0)
gamma_minus_one = gamma - FloatType(1.0)
delta = FloatType(1e-6) # relative density perturbation
uthermal_0 = pressure_0 / density_0 / gamma_minus_one
"""
loop over all output files; need to be at times when analytic
solution equals the initial conditions
"""
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 """
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Density = np.array(data["PartType0"]["Density"], dtype = FloatType)
Velocity = np.array(data["PartType0"]["Velocities"], dtype = FloatType)
Uthermal = np.array(data["PartType0"]["InternalEnergy"], dtype = FloatType)
""" calculate analytic solution at new cell positions """
Density_ref = np.full(Pos.shape[0], density_0, dtype=FloatType)
Velocity_ref = np.zeros(Pos.shape, dtype=FloatType)
Uthermal_ref = np.full(Pos.shape[0], uthermal_0, dtype=FloatType)
## perturbations
Density_ref *= FloatType(1.0) + delta * np.sin( FloatType(2.0) * FloatType(np.pi) * Pos[:,0] / Boxsize )
Velocity_ref[:,0] = velocity_0
Uthermal_ref *= (Density / density_0)**gamma_minus_one
""" compare data """
## density
abs_delta_dens = np.abs(Density - Density_ref) / Density_ref
L1_dens = np.average(abs_delta_dens)
## velocity, here, use absolute error (velocity should be zero!)
abs_delta_vel = np.abs(Velocity - Velocity_ref)
L1_vel = np.average(abs_delta_vel)
## internal energy
abs_delta_utherm = np.abs(Uthermal - Uthermal_ref) / Uthermal_ref
L1_utherm = np.average(abs_delta_utherm)
""" printing results """
print("examples/wave_1d/check.py: L1 error of " + filename +":")
print("\t density: %g" % L1_dens)
print("\t velocity: %g" % L1_vel)
print("\t specific internal energy: %g" % L1_utherm)
print("\t tolerance: %g for %d cells" % (DeltaMaxAllowed, NumberOfCells) )
""" criteria for failing the test """
if L1_dens > DeltaMaxAllowed or L1_vel > DeltaMaxAllowed or L1_utherm > DeltaMaxAllowed:
sys.exit(1)
i_file += 1
""" normal exit """
sys.exit(0)
""" @package examples/wave_1d/create.py
Code that creates 1d wave test problem;
supposed to be as siple as possible and self-contained
created by Rainer Weinberger, last modified 19.2.2019 -- comments welcome
"""
""" 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/wave_1d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
## ensure calculations happen with predefined precision
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32 # integer type
## computational domain
Boxsize = FloatType(1.0)
NumberOfCells = IntType(32)
## initial state
density_0 = FloatType(1.0)
velocity_0 = FloatType(0.0)
pressure_0 = FloatType(3.0) / FloatType(5.0)
gamma = FloatType(5.0) / FloatType(3.0)
gamma_minus_one = gamma - FloatType(1.0)
delta = FloatType(1e-6) # relative density perturbation
uthermal_0 = pressure_0 / density_0 / gamma_minus_one
""" set up grid: uniform 1d grid """
## spacing
dx = Boxsize / FloatType(NumberOfCells)
## position of first and last cell
pos_first, pos_last = FloatType(0.5) * dx, Boxsize - FloatType(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 """
## set up unperturbed system; density, velocity and specific internal energy
Density = np.full(Pos.shape[0], density_0, dtype=FloatType)
Velocity = np.zeros(Pos.shape, dtype=FloatType)
Uthermal = np.full(Pos.shape[0], uthermal_0, dtype=FloatType)
## perturbations
Density *= FloatType(1.0) + delta * np.sin( FloatType(2.0) * FloatType(np.pi) * Pos[:,0] / Boxsize )
Velocity[:,0] = velocity_0
Uthermal *= (Density / density_0)**gamma_minus_one
## 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/wave_1d/param.txt
% parameter file for 1d wave propagation problem
InitCondFile ./run/examples/wave_1d/IC
ICFormat 3
OutputDir ./run/examples/wave_1d/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 1.0
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
MaxMemSize 2500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 1.0
TimeBegin 0.0
TimeMax 1.0
TimeBetSnapshot 1.0
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
CellShapingSpeed 0.5
CellMaxAngleFactor 2.25
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 64
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.01
MaxNumNgbDeviation 2
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 0
SofteningTypeOfPartType2 0
SofteningTypeOfPartType3 0
SofteningTypeOfPartType4 0
SofteningTypeOfPartType5 0
InitGasTemp 0.0
MinGasTemp 0.0
MinEgySpec 0.0
MinimumDensityOnStartUp 0.0
#!/bin/bash # this line only there to enable syntax highlighting in this file
## shell script for code testing
#################################
## perform predefined examples ##
#################################
## Number of cores to compile and run the problem on
## use NUMBER_OF_TASKS=1 for 1d test problems!
NUMBER_OF_TASKS=1
## choose your tests
TESTS=""
## available 1d test cases
TESTS+="wave_1d "
## loop over all tests
for TEST in $TESTS
do
DIR="./examples/"${TEST}"/"
RUNDIR="./run/examples/"${TEST}"/"
## clean up
rm -rf ./run
## create run directory
mkdir ./run
mkdir ./run/examples
mkdir ${RUNDIR}
## create ICs in run directory
echo ${DIR}
python ${DIR}/create.py ${RUNDIR}
((return_value=$?)) ## get return value
if [ $return_value != 0 ] ## check return value
then echo "ERROR: test.sh:\t" $DIR "\t python create.py failed!"
exit $return_value
fi
## copy Config and parameter file to run directory
cp ${DIR}/Config.sh ${RUNDIR}
cp ${DIR}/param.txt ${RUNDIR}
## compile Arepo
make -j ${NUMBER_OF_TASKS} CONFIG=${RUNDIR}/Config.sh BUILD_DIR=${RUNDIR}/build EXEC=${RUNDIR}/Arepo
((return_value=$?)) ## get return value
if [ $return_value != 0 ] ## check return value
then echo "ERROR: test.sh:\t" $DIR "\t make failed!"
exit $return_value
fi
## execute test simulation
mpiexec -np ${NUMBER_OF_TASKS} ${RUNDIR}/Arepo ${RUNDIR}/param.txt
((return_value=$?)) ## get return value
if [ $return_value != 0 ] ## check return value
then echo "ERROR: test.sh:\t" $DIR "\t execution failed!"
exit $return_value
fi
## check result in example directory, this also creates some check plots
python ${DIR}/check.py ${RUNDIR}
((return_value=$?)) ## get return value
if [ $return_value != 0 ] ## check return value
then echo "ERROR: test.sh: test failed!"
exit $return_value
else echo "test.sh:\t" $DIR "\t test passed!"
echo "cleaning up..."
fi
## clean up
rm -rf ./run
done
echo "tests"
echo ${TESTS}
echo "passed!"
exit ${return_value}
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