Commit 2bb6eb54 authored by Rainer Weinberger's avatar Rainer Weinberger

added example mhd_shocktube_1d

parent ca675a67
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/mhd_shocktube_1d/Config.sh
## config file for 1d mhd shocktube probelm
#--------------------------------------- Basic operation mode of code
ONEDIMS # 1d simulation
REFLECTIVE_X=2 # open boundaries
#--------------------------------------- Magnetohydrodynamics
MHD # Master switch for magnetohydrodynamics
MHD_POWELL # Powell div(B) cleaning scheme for magnetohydrodynamics
MHD_POWELL_LIMIT_TIMESTEP # Timestep constraint due to Powell cleaning scheme
#--------------------------------------- Riemann solver
RIEMANN_HLLD # HLLD approximate Riemann solver (required to use for MHD)
#--------------------------------------- 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/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
""" @package ./examples/mhd_shocktube_1d/check.py
Code that checks results of 1d mhd shocktube problem
created by Rainer Weinberger, last modified 12.03.2019
"""
""" load libraries """
import sys ## load sys; needed for exit codes
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt
createFigures = False
CreateReferenceSolution = False
simulation_directory = str(sys.argv[1])
print("examples/mhd_shocktube_1d/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) ## 2d sim
""" loop over all output files """
status = 0
i_file = 1
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 """
x = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
vel = np.array(data["PartType0"]["Velocities"], dtype = FloatType)
rho = np.array(data["PartType0"]["Density"], dtype = FloatType)
u = np.array(data['PartType0']['InternalEnergy'], dtype = FloatType)
B = np.array(data["PartType0"]["MagneticField"], dtype = FloatType)
absB = np.sqrt(B[:,0]*B[:,0] + B[:,1]*B[:,1] + B[:,2]*B[:,2])
alphaB = (B[:,1]/B[:,2])
verificationData = np.array([x[:,0], rho, vel[:,0], u, absB]).T
if CreateReferenceSolution:
checkData = verificationData[::40,:]
np.savetxt('examples/mhd_shocktube_1d/data_%03d.txt'%i_file,checkData)
status = 0
else:
checkData = np.loadtxt('examples/mhd_shocktube_1d/data_%03d.txt'%i_file)
rho_ref = np.interp(x[:,0], checkData[:,0], checkData[:,1])
vel_ref = np.interp(x[:,0], checkData[:,0], checkData[:,2])
u_ref = np.interp(x[:,0], checkData[:,0], checkData[:,3])
absB_ref = np.interp(x[:,0], checkData[:,0], checkData[:,4])
delta_rho = rho - rho_ref
delta_vel = vel[:,0] - vel_ref
delta_u = u - u_ref
delta_B = absB - absB_ref
if createFigures:
fig, ax = plt.subplots(4, sharex=True, figsize=np.array([6.9,6.0]))
fig.subplots_adjust(left = 0.09, bottom = 0.09,right = 0.98, top = 0.98)
ax[0].plot(x[:,0], rho, marker='o', markersize=2, zorder=2)
ax[0].plot(checkData[:,0], checkData[:,1], c='red', zorder=1)
ax[0].set_ylabel(r'density')
ax[1].plot(x[:,0], vel[:,0], marker='o', markersize=2, zorder=2)
ax[1].plot(checkData[:,0], checkData[:,2], c='red', zorder=1)
ax[1].set_ylabel(r'vel')
ax[2].plot(x[:,0], u, marker='o', markersize=2, zorder=2)
ax[2].plot(checkData[:,0], checkData[:,3], c='red', zorder=1)
ax[2].set_ylabel(r'u$_{th}$')
ax[3].plot(x[:,0], absB, marker='o', markersize=2, zorder=2)
ax[3].plot(checkData[:,0], checkData[:,4], c='red', zorder=1)
ax[3].set_ylabel(r'abs(B)')
ax[3].set_xlabel(r'position')
ax[3].set_xlim([-0.1,2.6])
fig.align_ylabels(ax[:])
fig.savefig(simulation_directory+'/snap_%03d.pdf'%i_file)
res_scaling = 200. / np.float64(len(rho))
tolerance_rho = 0.02 * res_scaling
tolerance_vel = 0.04 * res_scaling
tolerance_u = 0.09 * res_scaling
tolerance_B = 0.05 * res_scaling
if np.std(delta_rho) > tolerance_rho:
status += 1
if np.std(delta_vel) > tolerance_vel:
status += 1
if np.std(delta_u) > tolerance_u:
status += 1
if np.std(delta_B) > tolerance_B:
status +=1
print('standard deviations of absolute error and tolerance (density, velocity, int. energy, magnetic field:')
print(np.std(delta_rho), tolerance_rho)
print(np.std(delta_vel), tolerance_vel)
print(np.std(delta_u), tolerance_u)
print(np.std(delta_B), tolerance_B)
i_file += 1
exit(status)
\ No newline at end of file
""" @package ./examples/mhd_shocktube_1d/create.py
Code that creates 1d mhd shocktube initial conditions
created by Rainer Weinberger, last modified 12.03.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 write initial conditions in hdf5 format
simulation_directory = str(sys.argv[1])
print("examples/mhd_shocktube_1d/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(2.5) # quadratic box
NumberOfCells = IntType(400)
alpha = np.pi
## parameters
density_L = FloatType(1.0)
velocity_L = FloatType(0.0)
pressure_L = FloatType(1.0)
b_L = np.array([1.0, 1.0, 0.0], dtype=FloatType)
density_R = FloatType(0.2)
velocity_R = FloatType(0.0)
pressure_R = FloatType(0.2)
b_R = np.array([1.0, np.cos(alpha), np.sin(alpha)], dtype=FloatType)
gamma = FloatType(5.0/3.0)
gamma_minus_one = FloatType(gamma - 1.0)
""" set up 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 magnetohydrodynamical quantitites """
## left state
Mass = np.full(NumberOfCells, density_L*dx, dtype=FloatType)
Velocity = np.zeros([NumberOfCells,3], dtype=FloatType)
Velocity[:,0] = velocity_L
Uthermal = np.full(NumberOfCells, (pressure_L/density_L/gamma_minus_one), dtype=FloatType)
Bfield = np.zeros([NumberOfCells, 3], dtype=FloatType)
for dim in np.arange(3):
Bfield[:,dim] = b_L[dim]
## right state
i_right, = np.where(Pos[:,0] > 1.0)
Mass[i_right] = density_R*dx
Velocity[i_right,0] = velocity_R
Uthermal[i_right] = pressure_R/density_R/gamma_minus_one
for dim in np.arange(3):
Bfield[:,dim] = b_R[dim]
""" 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 = Mass)
part0.create_dataset("Velocities", data = Velocity)
part0.create_dataset("InternalEnergy", data = Uthermal)
part0.create_dataset("MagneticField", data = Bfield)
## close file
IC.close()
This diff is collapsed.
This diff is collapsed.
%% examples/mhd_shocktube_1d.txt
% parameter file for 1d mhd shocktube problem
InitCondFile ./run/examples/mhd_shocktube_1d/IC
ICFormat 3
OutputDir ./run/examples/mhd_shocktube_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 2.5
PeriodicBoundariesOn 1
ComovingIntegrationOn 0
MaxMemSize 2500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 0.005
TimeBegin 0.0
TimeMax 0.4
TimeBetSnapshot 0.2
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.1
MinSizeTimestep 1e-10
CourantFac 0.3
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
DesNumNgb 64
MultipleDomains 2
TopNodeFactor 4
ActivePartFracForNewDomainDecomp 0.5
MaxNumNgbDeviation 2
TypeOfTimestepCriterion 1
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
......@@ -16,6 +16,7 @@ TESTS=""
TESTS+="wave_1d "
TESTS+="shocktube_1d "
TESTS+="interacting_blastwaves_1d "
TESTS+="mhd_shocktube_1d "
TESTS+="polytrope_1d_spherical "
## available 2d examples
......
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