Commit 68a6f07c authored by Rainer Weinberger's avatar Rainer Weinberger

new example polytrope_1d_spherical and minor changes in existing ones

parent 0d96f291
......@@ -24,7 +24,7 @@ Redshifts = [1, 0]
status = 0
CompareAgainstReferenceRun = True ## comparison for small L50m32 box; deactivate this when comparing against self-created ICs
createFigures = True
createFigures = False
for i_file, z in enumerate(Redshifts):
""" try to read in snapshot """
......
......@@ -10,7 +10,7 @@ import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff
createFigures = True
createFigures = False
simulation_directory = str(sys.argv[1])
print("examples/noh_2d/check.py: checking simulation output in directory " + simulation_directory)
......@@ -100,7 +100,7 @@ while True:
L1_dens = np.average(abs_delta_dens, weights=CellVolume[i_compare] )
L1_max = DeltaMaxAllowed * time
print("examples/Noh_2d/check snapshot %d: DEBUG: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
print("examples/Noh_2d/check snapshot %d: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
if L1_dens > L1_max:
print("examples/Noh_2d/check: ERROR: snaphshot %d: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
......
......@@ -10,7 +10,7 @@ import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff
createFigures = True
createFigures = False
simulation_directory = str(sys.argv[1])
print("examples/noh_3d/check.py: checking simulation output in directory " + simulation_directory)
......
#!/bin/bash # this line only there to enable syntax highlighting in this file
## ./examples/polytrope_1d_spherical/Config.sh
## config file for 1d polytrope probelm in a spherically symmetric setup
#--------------------------------------- Basic operation mode of code
ONEDIMS # 1d simulation
ONEDIMS_SPHERICAL # 1d spherically symmetric simulation
GAMMA=2 # Adiabatic index of gas
#--------------------------------------- Mesh motion and regularization
VORONOI_STATIC_MESH # static mesh
#--------------------------------------- Gravity treatment; default: no gravity
SELFGRAVITY # gravitational intraction between simulation particles/cells
#--------------------------------------- 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
READ_MASS_AS_DENSITY_IN_INPUT # Reads the mass field in the IC as density
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
#--------------------------------------- output
OUTPUTACCELERATION # output gravitational acceleration
OUTPUT_PRESSURE_GRADIENT # output pressure gradient
""" @package ./examples/polytrope_1d_spherical/check.py
Code that checks results of 1d polytrope test problem
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 # plot results
createFigures = False
simulation_directory = str(sys.argv[1])
print("examples/polytrope_1d_spherical/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])
## maximum L1 error after one propagation; based on experience
DeltaMaxAllowed = 0.001 * FloatType(NumberOfCells/256.0)**-2
""" Initial conditions as reference """
i_snap = 0
directory = simulation_directory+"/output/"
filename = "snap_%03d.hdf5" % (i_snap)
data = h5py.File(directory+filename, "r")
Pos_ref = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)[:,0]
Density_ref = np.array(data["PartType0"]["Density"], dtype = FloatType)
Velocity_ref = np.array(data["PartType0"]["Velocities"], dtype = FloatType)[:,0]
Uthermal_ref = np.array(data["PartType0"]["InternalEnergy"], dtype = FloatType)
Accel_ref = np.array(data["PartType0"]["Acceleration"], dtype = FloatType)[:,0]
GradPress_ref = np.array(data["PartType0"]["PressureGradient"], dtype = FloatType)[:,0]
while True:
i_snap +=1
""" compare to snapshots """
filename = "snap_%03d.hdf5" % (i_snap)
try:
data = h5py.File(directory+filename, "r")
except:
break
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)[:,0]
Density = np.array(data["PartType0"]["Density"], dtype = FloatType)
Velocity = np.array(data["PartType0"]["Velocities"], dtype = FloatType)[:,0]
Uthermal = np.array(data["PartType0"]["InternalEnergy"], dtype = FloatType)
Accel = np.array(data["PartType0"]["Acceleration"], dtype = FloatType)[:,0]
GradPress = np.array(data["PartType0"]["PressureGradient"], dtype = FloatType)[:,0]
""" plots """
if createFigures:
fig, ax = plt.subplots(ncols=1, nrows=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, Density, 'r', label='final')
ax[0].plot(Pos_ref, Density_ref, 'b--', label='IC')
ax[0].set_ylabel(r"density")
ax[0].fill_between([0.1,0.7],[0.0,0.0],[1.01,1.01], color='k',alpha=0.2)
ax[1].plot(Pos, Density/Density_ref, 'r')
ax[1].plot(Pos_ref, [1.0]*Pos_ref.shape[0], 'b--')
ax[1].set_ylabel(r"rel. density")
ax[1].set_ylim([0.99,1.01])
ax[1].fill_between([0.1,0.7],[0.99,0.99],[1.01,1.01], color='k',alpha=0.2)
ax[2].plot(Pos, Velocity, 'r')
ax[2].plot(Pos_ref, [0.0]*Pos_ref.shape[0], 'b--')
ax[2].set_ylabel(r"velocity")
ax[2].set_ylim([-0.01,0.01])
ax[2].fill_between([0.1,0.7],[-0.01,-0.01],[0.01,0.01], color='k',alpha=0.2)
ax[3].plot(Pos[:-1], Accel[:-1] - GradPress[:-1] / Density[:-1], 'r')
ax[3].plot(Pos_ref[:-1], Accel_ref[:-1] - GradPress_ref[:-1] / Density_ref[:-1], 'b--')
ax[3].set_ylim([-0.5,0.1])#
ax[3].set_ylabel(r"net accel.")
ax[3].set_xlabel(r"radius")
ax[3].set_xlim([0.0,0.8])
ax[3].fill_between([0.1,0.7],[-0.5,-0.5],[0.1,0.1], color='k',alpha=0.2)
ax[0].legend(loc=1, frameon=False)
fig.align_ylabels(ax[:])
fig.savefig(simulation_directory+"profiles_%03d.pdf"%i_snap)
""" compare to ICs by interpolating to IC positions """
delta_dens = np.interp(Pos_ref, Pos, Density) - Density_ref
delta_vel = np.interp(Pos_ref, Pos, Velocity) - Velocity_ref
delta_uthermal = np.interp(Pos_ref, Pos, Uthermal) - Uthermal_ref
## L1 norm
i_check, = np.where(Pos_ref < 0.7)
L1_dens = np.average( np.abs(delta_dens[i_check]) )
L1_vel = np.average( np.abs(delta_vel[i_check]) )
L1_uthermal = np.average( np.abs(delta_uthermal[i_check]) )
""" printing results """
print("examples/polytrope_1d_spherical/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_uthermal)
print("\t tolerance: %g for %d cells" % (DeltaMaxAllowed, NumberOfCells) )
""" criteria for failing the test """
if L1_dens > DeltaMaxAllowed or L1_vel > DeltaMaxAllowed or L1_uthermal > DeltaMaxAllowed:
sys.exit(1)
""" normal exit """
sys.exit(0)
""" @package examples/polytrope_1d_spherical/create.py
Code that creates 1d polytrope test problem;
supposed to be as siple as possible and self-contained
created by Rainer Weinberger, last modified 04.03.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/polytrope_1d_spherical/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)
NumberOfCells = IntType(256)
## initial state
G = FloatType(1.0)
density_0 = FloatType(1.0)
PolytropicIndex = FloatType(1.0)
velocity_0 = FloatType(0.0)
gamma = FloatType( (1.0 + PolytropicIndex) / PolytropicIndex) ## note: this has to be consistent with the parameter settings for Arepo
pressure_0 = FloatType(1.)/gamma
""" set up grid: uniform radial 1d grid """
## spacing
dx = Boxsize / FloatType(NumberOfCells)
## position of first and last cell
pos_first, pos_last = FloatType(0.1) + 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 = FloatType(4.) * np.pi / FloatType(3.) * ( (Pos[:,0]+FloatType(0.5)*dx)**3 - (Pos[:,0]-FloatType(0.5)*dx)**3 )
""" set up hydrodynamical quantitites """
Rscale = FloatType(0.8) * Boxsize / np.pi
Radius = Pos[:,0]
theta = np.sinc(Radius/Rscale/np.pi) ## np.sinc(x) = np.sin(np.pi x)/(np.pi x) or 1 if x=0
theta[theta<1.e-12] = FloatType(1.e-12)
K = Rscale**2 * FloatType(4.0) * np.pi * G / ( (PolytropicIndex+FloatType(1.0) ) * density_0**(FloatType(1.0)/PolytropicIndex - FloatType(1.0) ) )
## explicit n=1 case
K = Rscale**2 * FloatType(2.0) * np.pi * G
## hydrodynamic quantities
Density = density_0 * theta**PolytropicIndex
Velocity = np.zeros([NumberOfCells, 3], dtype=FloatType)
Velocity[:,0] = velocity_0
Pressure = K * Density**gamma
Uthermal = Pressure / Density / (gamma - FloatType(1.0) )
## density in mass filed in input
Mass = Density
""" 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/polytrope_1d_spherical/param.txt
% parameter file for 1d polytrope test problem
InitCondFile ./run/examples/polytrope_1d_spherical/IC
ICFormat 3
OutputDir ./run/examples/polytrope_1d_spherical/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 0.125
UnitVelocity_in_cm_per_s 1.0
UnitLength_in_cm 1.0
UnitMass_in_g 1.0
GravityConstantInternal 1.0
ErrTolIntAccuracy 0.1
ErrTolTheta 0.1
ErrTolForceAcc 0.1
MaxSizeTimestep 0.125
MinSizeTimestep 1e-10
CourantFac 0.3
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 32
MaxNumNgbDeviation 2
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.01
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
CoreRadius 0.1
CoreMass 4.12455e-3
......@@ -12,8 +12,8 @@ import matplotlib.pyplot as plt ## needs to be active for plotting!
from Riemann import * ## Riemann-solver
createFigures = True
forceExitOnError=False ## exits immediately when tolerance is exceeded
createFigures = False
forceExitOnError = False ## exits immediately when tolerance is exceeded
""" check functiions """
def CheckL1Error(Pos, W, W_L, W_R, gamma, position_0, time):
......
......@@ -16,6 +16,7 @@ TESTS=""
TESTS+="wave_1d "
TESTS+="shocktube_1d "
TESTS+="interacting_blastwaves_1d "
TESTS+="polytrope_1d_spherical "
## 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