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

added example yee_2d

parent af0e98cb
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/yee_2d/Config.sh
## config file for 2d Yee vortex probelm
#--------------------------------------- Basic operation mode of code
TWODIMS # 2d simulation
GAMMA=1.4 # Adiabatic index of gas; 5/3 if not set
READ_MASS_AS_DENSITY_IN_INPUT # Reads the mass field in the IC as density
#--------------------------------------- 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_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/yee_2d/check.py
Code that checks results of 2d Yee vortex problem
created by Rainer Weinberger, last modified 21.02.2019 -- comments welcome
"""
#### load libraries
import sys ## load sys; needed for exit codes
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
simulation_directory = str(sys.argv[1])
print("examples/yee_2d/check.py: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
## 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])
CellsPerDimension = np.sqrt(NumberOfCells)
#### parameters
Tinf = FloatType(1.0)
beta = FloatType(5.0)
gamma = FloatType(1.4) ## note: this has to be consistent with the parameter settings for Arepo!
## maximum L1 error, empirically based
DeltaMaxAllowed = 0.005 * (FloatType(CellsPerDimension) / 50.0)**-2
""" 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 """
## 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)
Volume = Mass / Density
xPosFromCenter = (Pos[:,0] - 0.5 * Boxsize)
yPosFromCenter = (Pos[:,1] - 0.5 * Boxsize)
Radius = np.sqrt( xPosFromCenter**2 + yPosFromCenter**2 )
RotationVelocity = xPosFromCenter * Velocity[:,1] - yPosFromCenter * Velocity[:,0]
i_select = np.where(Radius > 0)[0]
RotationVelocity[i_select] /= Radius[i_select]
""" calculate analytic solution """
Temperature_ref = Tinf - ( ( gamma - 1.0 ) * beta * beta / ( 8.0 * np.pi * np.pi * gamma ) * np.exp( 1.0 - Radius*Radius ) )
exponent = 1.0 / ( gamma - 1.0 )
Density_ref = Temperature_ref**exponent
Vel_ref = np.zeros([NumberOfCells, 3], dtype=FloatType)
Vel_ref[:,0] = -0.5 * yPosFromCenter * beta / np.pi * np.exp( 0.5 * (1.0 - Radius * Radius) )
Vel_ref[:,1] = 0.5 * xPosFromCenter * beta / np.pi * np.exp( 0.5 * (1.0 - Radius * Radius) )
RotationVelocity_ref = xPosFromCenter * Vel_ref[:,1] - yPosFromCenter * Vel_ref[:,0]
RotationVelocity_ref /= Radius
Uthermal_ref = Temperature_ref / (gamma - 1.0)
""" compare data """
## density
abs_delta_dens = np.abs(Density - Density_ref)
L1_dens = np.average(abs_delta_dens, weights=Volume)
## velocity, here, use absolute error (velocity can be close to zero!)
abs_delta_vel = np.abs(RotationVelocity - RotationVelocity_ref)
L1_vel = np.average(abs_delta_vel, weights=Volume)
## internal energy
abs_delta_utherm = np.abs(Uthermal - Uthermal_ref)
L1_utherm = np.average(abs_delta_utherm, weights=Volume)
""" printing results """
print("examples/Yee_2d/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 per dimension" % (DeltaMaxAllowed, CellsPerDimension) )
""" 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/yee_2d/create.py
Code that creates 2d Yee vortex initial conditions
parameters are identical to Pakmor et al. (2016), MNRAS 455, 1134
created by Rainer Weinberger, last modified 21.02.2019 -- comments welcome
"""
#### load libraries
import sys # system specific calls
import numpy as np ## load numpy
import h5py ## load h5py; needed to write initial conditions in hdf5 format
simulation_directory = str(sys.argv[1])
print("examples/yee_2d/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(10.0)
CellsPerDimension = IntType(50)
#### parameters
Tinf = FloatType(1.0)
beta = FloatType(5.0)
gamma = FloatType(1.4) ## note: this has to be consistent with the parameter settings for Arepo!
""" set up grid: equidistant polar 2d grid; similar to Pakmor et al (2016) """
d_ring = Boxsize / CellsPerDimension
## place central cell
Radius = [0.0]
xPosFromCenter = [0.0]
yPosFromCenter = [0.0]
i_ring = 0
while d_ring * FloatType(i_ring) <= 0.71 * Boxsize:
i_ring += 1
## place i_ring cells at this distance
phi = np.random.uniform(0,2.0*np.pi) ## random starting angle
#print "i_ring", i_ring
for i_cell in np.arange( IntType(2.0*np.pi*i_ring) ):
radius_this_cell = d_ring * FloatType(i_ring)
#print radius_this_cell
xcoord = radius_this_cell * np.sin(phi)
ycoord = radius_this_cell * np.cos(phi)
## only include the cell if coordinates are within the box
if xcoord >= -0.5*Boxsize and xcoord < 0.5*Boxsize and ycoord >= -0.5*Boxsize and ycoord < 0.5*Boxsize:
Radius.append( radius_this_cell )
xPosFromCenter.append( xcoord )
yPosFromCenter.append( ycoord )
phi += (2.0*np.pi*i_ring) / IntType(2.0*np.pi*i_ring) / FloatType(i_ring)
## convert to numpy arrays now that number of cells is known
Radius = np.array(Radius, dtype=FloatType)
xPosFromCenter = np.array(xPosFromCenter, dtype=FloatType)
yPosFromCenter = np.array(yPosFromCenter, dtype=FloatType)
NumberOfCells = Radius.shape[0]
## set up structure for positions (in code coordinates, i.e. from 0 to Boxsize)
Pos = np.zeros([NumberOfCells, 3])
Pos[:,0] = xPosFromCenter + 0.5 * Boxsize
Pos[:,1] = yPosFromCenter + 0.5 * Boxsize
""" set up hydrodynamical quantitites """
Temperature = Tinf - ( ( gamma - 1.0 ) * beta * beta / ( 8.0 * np.pi * np.pi * gamma ) * np.exp( 1.0 - Radius*Radius ) )
exponent = 1.0 / ( gamma - 1.0 )
Density = Temperature**exponent
Vel = np.zeros([NumberOfCells, 3], dtype=FloatType)
Vel[:,0] = -0.5 * yPosFromCenter * beta / np.pi * np.exp( 0.5 * (1.0 - Radius * Radius) )
Vel[:,1] = 0.5 * xPosFromCenter * beta / np.pi * np.exp( 0.5 * (1.0 - Radius * Radius) )
Uthermal = Temperature / (gamma - 1.0)
""" 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 = Density)
part0.create_dataset("Velocities", data = Vel)
part0.create_dataset("InternalEnergy", data = Uthermal)
## close file
IC.close()
sys.exit(0)
%% ./examples/yee_2d/param.txt
% parameter file for 2d Yee vortex problem
InitCondFile ./run/examples/yee_2d/IC
ICFormat 3
OutputDir ./run/examples/yee_2d/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 10.0
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
MaxMemSize 2500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 0.5
TimeBegin 0.0
TimeMax 10.0
TimeBetSnapshot 10.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 1.0
MinSizeTimestep 1e-5
CourantFac 0.3
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
CellShapingSpeed 0.5
CellMaxAngleFactor 2.25
......@@ -9,12 +9,16 @@
## use NUMBER_OF_TASKS=1 for 1d test problems!
NUMBER_OF_TASKS=1
## choose your tests
## choose your examples
TESTS=""
## available 1d test cases
## available 1d examples
TESTS+="wave_1d "
TESTS+="shocktube_1d "
## available 2d examples
TESTS+="yee_2d "
## 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