Commit 7dd21e62 authored by Rainer Weinberger's avatar Rainer Weinberger

Merge branch 'master' into 'master'

Master

See merge request vrs/arepo!1
parents 0444abba 39d4e583
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/alfven_wave_1d/Config.sh
## config file for 1d Alfven wave probelm
#--------------------------------------- Basic operation mode of code
ONEDIMS # 1d simulation
#--------------------------------------- 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_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_CENTER_OF_MASS # output centers of cells
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
This diff is collapsed.
""" @package examples/alfven_wave_1d/create.py
Code that creates 1d Alfven wave test problem;
supposed to be as simple and self-contained as possible
created by Alessandro Stenghel and Federico Marinacci,
last modified 13.7.2020 -- 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/alfven_wave_1d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/ics.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)
if len(sys.argv) > 3:
NumberOfCells = IntType(sys.argv[3])
else:
NumberOfCells = IntType(32)
## initial state
density_0 = FloatType(1.0)
velocity_0 = FloatType(0.0)
pressure_0 = FloatType(1.0)
gamma = FloatType(5.0) / FloatType(3.0)
gamma_minus_one = gamma - FloatType(1.0)
delta = FloatType(1e-6) # relative velocity perturbation
uthermal_0 = pressure_0 / density_0 / gamma_minus_one
bfield_0= FloatType(1.0)
k_z = 2*np.pi
omega = bfield_0*k_z/np.sqrt(density_0)
""" 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 magnetohydrodynamical 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)
Bfield = np.zeros(Pos.shape, dtype=FloatType)
## perturbations
Velocity[:,0] = velocity_0
Velocity[:,1] = delta*np.sin(k_z*Pos[:,0])
Velocity[:,2] = delta*np.cos(k_z*Pos[:,0])
Uthermal *= (Density / density_0)**gamma_minus_one
Bfield[:,0] = bfield_0
Bfield[:,1] = -k_z*bfield_0/omega*Velocity[:,1]
Bfield[:,2] = -k_z*bfield_0/omega*Velocity[:,2]
## 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)
part0.create_dataset("MagneticField", data=Bfield*np.sqrt(4*np.pi)) #conversion Lorentz System->Gaussian System
## close file
IC.close()
""" normal exit """
sys.exit(0)
%% examples/alfven_wave_1d.txt
% parameter file for 1d Alfven wave problem
InitCondFile ics
ICFormat 3
OutputDir 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 1500
TimeOfFirstSnapshot 0.0
CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
TimeBetStatistics 0.005
TimeBegin 0.0
TimeMax 2.0
TimeBetSnapshot 0.25
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
......@@ -171,8 +171,6 @@ int init(void)
{
P[i].Pos[1] = 0.0;
P[i].Pos[2] = 0.0;
P[i].Vel[1] = 0.0;
P[i].Vel[2] = 0.0;
}
#endif /* #ifdef ONEDIMS */
......@@ -180,7 +178,6 @@ int init(void)
for(i = 0; i < NumPart; i++)
{
P[i].Pos[2] = 0;
P[i].Vel[2] = 0;
}
#endif /* #ifdef TWODIMS */
......
......@@ -19,6 +19,7 @@ TESTS+="shocktube_1d "
TESTS+="interacting_blastwaves_1d "
TESTS+="mhd_shocktube_1d "
TESTS+="polytrope_1d_spherical "
TESTS+="alfven_wave_1d "
## 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