Commit 6763b415 authored by Rainer Weinberger's avatar Rainer Weinberger
Browse files

added example gresho_2d

parent e5aeefb1
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/gresho_2d/Config.sh
## config file for 2d Gresho vortex probelm
#--------------------------------------- Basic operation mode of code
TWODIMS # 2d simulation
#--------------------------------------- MPI/Threading Hybrid
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/gresho_2d/check.py
Code that checks results of 2d Gresho vortex problem
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 read snapshots
simulation_directory = str(sys.argv[1])
print("examples/gresho_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
density_0 = 1.0
velocity_0 = 3.0
gamma = 5.0/3.0
gamma_minus_one = gamma - 1.0
## maximum L1 error after one propagation; empirically based
DeltaMaxAllowed = 0.03 * (FloatType(CellsPerDimension) / 40.0)**-1.4
""" 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)
Velocity[:,0] -= velocity_0
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 )
i_select, = np.where(Radius != 0.0)
RotationVelocity = yPosFromCenter * Velocity[:,0] - xPosFromCenter * Velocity[:,1]
RotationVelocity[i_select] /= Radius[i_select]
RotationVelocity[Radius[:] == 0.0] = 0.0
""" calculate analytic solution """
## density
Density_ref = np.full(Density.shape, density_0, dtype=FloatType)
## different zones
i1, = np.where( Radius < 0.2 )
i2, = np.where( (Radius >= 0.2) & (Radius < 0.4) )
i3, = np.where( Radius >= 0.4 )
## velocity
RotationVelocity_ref = np.zeros(Velocity.shape[0], dtype=FloatType)
RotationVelocity_ref[i1] = 5.0 * Radius[i1]
RotationVelocity_ref[i2] = 2.0 - 5.0 * Radius[i2]
RotationVelocity_ref[i3] = 0.0
Velocity_ref = np.zeros(Velocity.shape, dtype=FloatType )
i_select, = np.where(Radius != 0.0)
Velocity_ref[i_select,0] = RotationVelocity_ref[i_select] * (Pos[i_select,1] - 0.5 * Boxsize) / Radius[i_select]
Velocity_ref[i_select,1] = -RotationVelocity_ref[i_select] * (Pos[i_select,0] - 0.5 * Boxsize) / Radius[i_select]
## specific internal energy
Pressure_ref = np.zeros(Uthermal.shape, dtype=FloatType)
Pressure_ref[i1] = 5.0 + 12.5 * Radius[i1] * Radius[i1]
Pressure_ref[i2] = 9.0 + 12.5 * Radius[i2] * Radius[i2] - 20.0 * Radius[i2] + 4.0 * np.log(Radius[i2] / 0.2)
Pressure_ref[i3] = 3.0 + 4.0 * np.log(2.0)
Uthermal_ref = Pressure_ref / density_0 / gamma_minus_one
""" compare data """
## density
abs_delta_dens = np.abs(Density - Density_ref)
L1_dens = np.average(abs_delta_dens, weights = Volume)
## velocity
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/Gresho_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
print("normal exit")
""" normal exit """
sys.exit(0)
""" @package ./examples/gresho_2d/create.py
Code that creates 2d Gresho vortex initial conditions
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/gresho_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(1.0)
CellsPerDimension = IntType(40)
## parameters
density_0 = 1.0
velocity_0 = 3.0 ## bulk velocity
gamma = 5.0/3.0
gamma_minus_one = gamma - 1.0
""" 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 """
## mass insetad of density
Density = np.full(NumberOfCells, density_0, dtype=FloatType)
## different zones
i1, = np.where( Radius < 0.2 )
i2, = np.where( (Radius >= 0.2) & (Radius < 0.4) )
i3, = np.where( Radius >= 0.4 )
## velocity
RotationVelocity = np.zeros( NumberOfCells, dtype=FloatType )
RotationVelocity[i1] = 5.0 * Radius[i1]
RotationVelocity[i2] = 2.0 - 5.0 * Radius[i2]
RotationVelocity[i3] = 0.0
Vel = np.zeros([NumberOfCells, 3], dtype=FloatType )
i_all_but_central = np.arange(1,NumberOfCells)
Vel[i_all_but_central,0] = RotationVelocity[i_all_but_central] * yPosFromCenter[i_all_but_central] / Radius[i_all_but_central]
Vel[i_all_but_central,1] = -RotationVelocity[i_all_but_central] * xPosFromCenter[i_all_but_central] / Radius[i_all_but_central]
Vel[:,0] += velocity_0
## specific internal energy
Pressure = np.zeros( NumberOfCells, dtype=FloatType)
Pressure[i1] = 5.0 + 12.5 * Radius[i1] * Radius[i1]
Pressure[i2] = 9.0 + 12.5 * Radius[i2] * Radius[i2] - 20.0 * Radius[i2] + 4.0 * np.log(Radius[i2] / 0.2)
Pressure[i3] = 3.0 + 4.0 * np.log(2.0)
Uthermal = Pressure / density_0 / gamma_minus_one
""" 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/gresho_2d/param.txt
% parameter file for 2d Gresho vortex problem
InitCondFile ./run/examples/gresho_2d/IC
ICFormat 3
OutputDir ./run/examples/gresho_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 1.0
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
MaxMemSize 2500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 0.5
TimeBegin 0.0
TimeMax 3.0
TimeBetSnapshot 3.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
MaxNumNgbDeviation 2
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.5
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
%% ./examples/yee_2d/param.txt %% ./examples/yee_2d/param.txt
% parameter file for 2d Yee vortex problem % parameter file for 2d Yee vortex problem
InitCondFile ./run/examples/yee_2d/IC InitCondFile ./run/examples/yee_2d/IC
ICFormat 3 ICFormat 3
OutputDir ./run/examples/yee_2d/output/ OutputDir ./run/examples/yee_2d/output/
SnapshotFileBase snap SnapshotFileBase snap
SnapFormat 3 SnapFormat 3
NumFilesPerSnapshot 1 NumFilesPerSnapshot 1
NumFilesWrittenInParallel 1 NumFilesWrittenInParallel 1
ResubmitOn 0 ResubmitOn 0
ResubmitCommand my-scriptfile ResubmitCommand my-scriptfile
OutputListFilename ol OutputListFilename ol
OutputListOn 0 OutputListOn 0
CoolingOn 0 CoolingOn 0
StarformationOn 0 StarformationOn 0
Omega0 0.0 Omega0 0.0
OmegaBaryon 0.0 OmegaBaryon 0.0
OmegaLambda 0.0 OmegaLambda 0.0
HubbleParam 1.0 HubbleParam 1.0
BoxSize 10.0 BoxSize 10.0
PeriodicBoundariesOn 1 PeriodicBoundariesOn 1
ComovingIntegrationOn 0 ComovingIntegrationOn 0
MaxMemSize 2500 MaxMemSize 2500
TimeOfFirstSnapshot 0.0
TimeOfFirstSnapshot 0.0 CpuTimeBetRestartFile 9000
CpuTimeBetRestartFile 9000 TimeLimitCPU 90000
TimeLimitCPU 90000
TimeBetStatistics 0.5
TimeBetStatistics 0.5 TimeBegin 0.0
TimeBegin 0.0 TimeMax 10.0
TimeMax 10.0 TimeBetSnapshot 10.0
TimeBetSnapshot 10.0
UnitVelocity_in_cm_per_s 1.0
UnitVelocity_in_cm_per_s 1.0 UnitLength_in_cm 1.0
UnitLength_in_cm 1.0 UnitMass_in_g 1.0
UnitMass_in_g 1.0 GravityConstantInternal 0.0
GravityConstantInternal 0.0
ErrTolIntAccuracy 0.1
ErrTolIntAccuracy 0.1 ErrTolTheta 0.1
ErrTolTheta 0.1 ErrTolForceAcc 0.1
ErrTolForceAcc 0.1
MaxSizeTimestep 1.0
MaxSizeTimestep 1.0 MinSizeTimestep 1e-5
MinSizeTimestep 1e-5 CourantFac 0.3
CourantFac 0.3
LimitUBelowThisDensity 0.0
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0 LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 64 DesNumNgb 64
......
...@@ -17,6 +17,7 @@ TESTS+="wave_1d " ...@@ -17,6 +17,7 @@ TESTS+="wave_1d "
TESTS+="shocktube_1d " TESTS+="shocktube_1d "
## available 2d examples ## available 2d examples
TESTS+="gresho_2d "
TESTS+="yee_2d " TESTS+="yee_2d "
## loop over all tests ## loop over all tests
......
Supports Markdown
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