Commit f9dfcf8c authored by Rainer Weinberger's avatar Rainer Weinberger

updated 2d examples

parent 012082c1
......@@ -33,4 +33,4 @@ HAVE_HDF5 # needed when HDF5 I/O support is desir
DEBUG # enables core-dumps, should this be standard?
#--------------------------------------- Output options
OUTPUT_BFIELD_GRADIENT # gradient estimates of B-field
\ No newline at end of file
OUTPUT_BFIELD_GRADIENT # gradient estimates of B-field
......@@ -8,12 +8,19 @@ created by Rainer Weinberger, last modified 13.03.2019
import sys ## system calls
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt
import os # file specific calls
import matplotlib.pyplot as plt # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
createFigures = True
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
simulation_directory = str(sys.argv[1])
print("examples/current_sheet_2d/check.py: checking simulation output in directory " + simulation_directory)
print("current_sheet_2d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
......@@ -68,10 +75,44 @@ while True:
Time = np.array(Time)
MagneticEnergy = np.array(MagneticEnergy)
if createFigures:
if makeplots:
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
fig, ax = plt.subplots(1)
ax.plot(Time, MagneticEnergy/MagneticEnergy[0])
fig.savefig(simulation_directory+'/MagneticEnergy.pdf')
fig.savefig(simulation_directory+'/plots/MagneticEnergy.pdf')
for i_file in [0,2,4,8]:
filename = "snap_%03d.hdf5" % (i_file)
try:
data = h5py.File(directory+filename, "r")
except:
break
VoronoiPos = np.array(data["PartType0"]["Coordinates"], dtype = np.float64)
MagneticField = np.array(data["PartType0"]["MagneticField"], dtype = np.float64)
Time = data["Header"].attrs["Time"]
Boxsize = data["Header"].attrs["BoxSize"]
NumberOfCells = np.int32(data["Header"].attrs["NumPart_Total"][0])
Nplot = 256
from scipy import spatial # needed for KDTree that we use for nearest neighbour search
Edges1d = np.linspace(0., Boxsize, Nplot+1, endpoint=True, dtype=np.float64)
Grid1d = 0.5 * (Edges1d[1:] + Edges1d[:-1])
xx, yy = np.meshgrid(Grid1d, Grid1d)
Grid2D = np.array( [xx.reshape(Nplot**2), yy.reshape(Nplot**2)] ).T
dist, cells = spatial.KDTree( VoronoiPos[:,:2] ).query( Grid2D, k=1 )
fig = plt.figure( figsize=(3.5,3.5), dpi=300 )
ax = plt.axes( [0,0,1,1] )
ax.streamplot( Grid1d, Grid1d, MagneticField[cells,0].reshape((Nplot,Nplot)), MagneticField[cells,1].reshape((Nplot,Nplot)), density=2., linewidth=0.5, arrowsize=1e-6, integration_direction='both', minlength=0.3 )
ax.set_xlim( 0, Boxsize )
ax.set_ylim( 0, Boxsize )
plt.text( 0.02, 0.95, "$N=%d^2,\ t=%3.1f$" % (np.int32(np.sqrt(NumberOfCells)), Time), color='k', transform=fig.transFigure, bbox={'boxstyle':'round', 'facecolor':'w'} )
fig.savefig( simulation_directory+"plots/streamlines_%03d.pdf" % (i_file), dpi=300 )
plt.close(fig)
## if everything is ok
if np.min(MagneticEnergy/MagneticEnergy[0]) > 0.92:
......
%% examples/current_sheet_2d.txt
% parameter file for 2d current sheet problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/gresho_2d/Config.sh
## 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
......@@ -22,6 +20,7 @@ TREE_BASED_TIMESTEPS # non-local timestep criterion (take 's
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_CENTER_OF_MASS # output centers of cells
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
......
""" @package ./examples/gresho_2d/check.py
""" @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
......@@ -8,9 +8,19 @@ created by Rainer Weinberger, last modified 21.02.2019 -- comments welcome
import sys # system specific calls
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import os # file specific calls
import matplotlib.pyplot as plt ## needs to be active for plotting!
plt.rcParams['text.usetex'] = True
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
simulation_directory = str(sys.argv[1])
print("examples/gresho_2d/check.py: checking simulation output in directory " + simulation_directory)
print("Gresho_2d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
......@@ -46,7 +56,9 @@ while True:
""" get simulation data """
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Time = FloatType(data["Header"].attrs["Time"])
Pos = np.array(data["PartType0"]["CenterOfMass"], dtype = FloatType)
VoronoiPos = 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)
......@@ -102,12 +114,67 @@ while True:
L1_utherm = np.average(abs_delta_utherm, weights = Volume)
""" printing results """
print("examples/Gresho_2d/check.py: L1 error of " + filename +":")
print("Gresho_2d: 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) )
if makeplots:
## plot:
fig = plt.figure( figsize=np.array([7.0,3.5]), dpi=300 )
Nplot = 256
from scipy import spatial # needed for KDTree that we use for nearest neighbour search
Edges1d = np.linspace(0., Boxsize, Nplot+1, endpoint=True, dtype=FloatType)
Grid1d = 0.5 * (Edges1d[1:] + Edges1d[:-1])
xx, yy = np.meshgrid(Grid1d, Grid1d)
Grid2D = np.array( [xx.reshape(Nplot**2), yy.reshape(Nplot**2)] ).T
vor = spatial.Voronoi( VoronoiPos[:,:2] )
dist, cells = spatial.KDTree( VoronoiPos[:,:2] ).query( Grid2D, k=1 )
ax = plt.axes( [0.05,0.15,0.35,0.7] )
pc = ax.pcolormesh( Edges1d, Edges1d, RotationVelocity[cells].reshape((Nplot,Nplot)), rasterized=True )
cax = plt.axes( [0.42,0.15,0.01,0.7] )
plt.colorbar( pc, cax=cax )
ax.set_title( 'Rotation velocity' )
spatial.voronoi_plot_2d( vor, ax=ax, line_colors='w', show_points=False, show_vertices=False, line_width=0.5 )
ax.set_xlim( 0, Boxsize )
ax.set_ylim( 0, Boxsize )
error = RotationVelocity[cells]-RotationVelocity_ref[cells]
ax = plt.axes( [0.53,0.15,0.35,0.7] )
pc = ax.pcolormesh( Edges1d, Edges1d, error.reshape((Nplot,Nplot)), rasterized=True )
cax = plt.axes( [0.90,0.15,0.01,0.7] )
plt.colorbar( pc, cax=cax )
ax.set_title( 'Rotation velocity error' )
spatial.voronoi_plot_2d( vor, ax=ax, line_colors='w', show_points=False, show_vertices=False, line_width=0.5 )
ax.set_xlim( 0, Boxsize )
ax.set_ylim( 0, Boxsize )
plt.text( 0.5, 0.92, "$t=%4.1f,\ L_1=%5.1e$" % (Time,L1_dens), ha='center', size=12, transform=fig.transFigure )
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
print(simulation_directory+"/plots/figure_%03d.pdf" % (i_file) )
fig.savefig(simulation_directory+"/plots/figure_%03d.pdf" % (i_file), dpi=300)
plt.close(fig)
fig = plt.figure( figsize=np.array([3.5,3.5]), dpi=300 )
ax = plt.axes( [0.19, 0.12, 0.75, 0.75] )
ax.plot( [0.,0.2,0.4], [0.,1.,0.], 'k--', label='Analytic solution' )
ax.plot( Radius, RotationVelocity, 'o', label='Arepo cells', mec='b', mfc="None", mew=0.3 )
ax.set_ylabel( "$v_\phi$" )
ax.set_xlabel( "$r$" )
ax.set_ylim( 0, 1 )
ax.legend( loc='upper right', frameon=False, fontsize=8 )
ax.set_title( "$\mathrm{Gresho\_2d:}\ \mathrm{N}=%d^2,\ \mathrm{L1}=%4.1e$" % (np.int32(CellsPerDimension),L1_vel), loc='right', size=8 )
plt.ticklabel_format( axis='y', style='sci', scilimits=(0,0) )
fig.savefig( simulation_directory+"plots/velocity_%03d.pdf" % (i_file) )
plt.close(fig)
""" criteria for failing the test """
if L1_dens > DeltaMaxAllowed or L1_vel > DeltaMaxAllowed or L1_utherm > DeltaMaxAllowed:
sys.exit(-1)
......
""" @package ./examples/gresho_2d/create.py
""" @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
......@@ -10,7 +10,7 @@ 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)
print("examples/Gresho_2d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
......
%% ./examples/gresho_2d/param.txt
%% ./examples/Gresho_2d/param.txt
% parameter file for 2d Gresho vortex problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/noh_2d/Config.sh
## examples/Noh_2d/Config.sh
## config file for 2d Noh probelm
#--------------------------------------- Basic operation mode of code
......@@ -19,6 +19,7 @@ TREE_BASED_TIMESTEPS # non-local timestep criterion (take 's
#---------------------------------------- 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)
......
""" @package ./examples/noh_2d/check.py
""" @package ./examples/Noh_2d/check.py
Code that checks results of 2d Noh problem
created by Rainer Weinberger, last modified 23.02.2019
......@@ -8,12 +8,19 @@ created by Rainer Weinberger, last modified 23.02.2019
import sys # system specific calls
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff
import os # file specific calls
import matplotlib.pyplot as plt # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
createFigures = False
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
simulation_directory = str(sys.argv[1])
print("examples/noh_2d/check.py: checking simulation output in directory " + simulation_directory)
print("Noh_2d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
......@@ -50,9 +57,11 @@ while True:
break
""" get simulation data """
time = FloatType(data["Header"].attrs["Time"])
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Pos = np.array(data["PartType0"]["CenterOfMass"], dtype = FloatType)
VoronoiPos = 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)
......@@ -74,25 +83,65 @@ while True:
Density_ref[i_postshock] = 16.0 ## for 2d
#### plot density profile
if createFigures:
fig, ax = plt.subplots(3, sharex=True, figsize=np.array([6.9,6.0]) )
fig.subplots_adjust(left = 0.13, bottom = 0.09,right = 0.98, top = 0.98)
if makeplots:
fig = plt.figure( figsize=np.array([7.5,6.0]), dpi=300 )
i_sorted = np.argsort(Radius)
ax[0].scatter(Radius, Density, rasterized=True, s=0.2)
ax[1].scatter(Radius, vRad, rasterized=True, s=0.2)
ax[2].scatter(Radius, Uthermal, rasterized=True, s=0.2)
if r_shock == 0:
r_shock = 0.1 * Boxsize
i_sorted = np.argsort(Radius)
ax[0].plot(Radius[i_sorted], Density_ref[i_sorted], label="analytic solution")
ax[0].set_xlim(0,0.8)
ax[0].set_ylim(0,22)
ax[2].set_xlabel(r"radius")
ax[0].set_ylabel(r"density")
ax[1].set_ylabel(r"radial velocity")
ax[2].set_ylabel(r"spec. internal energy")
Nplot = 512
from scipy import spatial # needed for KDTree that we use for nearest neighbour search and Voronoi mesh
Edges1d = np.linspace(0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.4*r_shock, Nplot+1, endpoint=True, dtype=FloatType)
Grid1d = 0.5 * (Edges1d[1:] + Edges1d[:-1])
xx, yy = np.meshgrid(Grid1d, Grid1d)
Grid2D = np.array( [xx.reshape(Nplot**2), yy.reshape(Nplot**2)] ).T
dist, cells = spatial.KDTree( VoronoiPos[:,:2] ).query( Grid2D, k=1 )
ax = plt.axes( [0.08,0.70,0.52,0.25] )
ax.plot(Radius[i_sorted], Density_ref[i_sorted], 'k', label="Analytic solution")
ax.plot(Radius, Density, 'b.', rasterized=True, ms=0.2)
ax.set_ylabel(r"Density")
ax.set_ylim( 0, 20 )
ax = plt.axes( [0.65,0.70,0.20,0.25] )
pc = ax.pcolormesh( Edges1d, Edges1d, Density[cells].reshape((Nplot,Nplot)), rasterized=True, cmap=plt.get_cmap('viridis'), vmin=0, vmax=20. )
cax = plt.axes( [0.88,0.70,0.02,0.25] )
plt.colorbar( pc, cax=cax )
ax.set_xlim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
ax.set_ylim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
ax = plt.axes( [0.08,0.38,0.52,0.25] )
ax.plot(Radius, vRad, 'b.', rasterized=True, ms=0.2)
ax.set_ylabel(r"Radial velocity")
ax.set_ylim( -1.1, 1.1 )
ax = plt.axes( [0.65,0.38,0.20,0.25] )
pc = ax.pcolormesh( Edges1d, Edges1d, vRad[cells].reshape((Nplot,Nplot)), rasterized=True, cmap=plt.get_cmap('plasma'), vmin=-1.1, vmax=1.1 )
cax = plt.axes( [0.88,0.38,0.02,0.25] )
plt.colorbar( pc, cax=cax )
ax.set_xlim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
ax.set_ylim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
fig.align_ylabels(ax[:])
fig.savefig(directory+"/figure_%03d.pdf"%i_file)
ax = plt.axes( [0.08,0.07,0.52,0.25] )
ax.plot(Radius, Uthermal, 'b.', rasterized=True, ms=0.2)
ax.set_xlabel(r"Radius")
ax.set_ylabel(r"Spec. internal energy")
ax.set_ylim( -0.05, 0.8 )
ax = plt.axes( [0.65,0.07,0.20,0.25] )
pc = ax.pcolormesh( Edges1d, Edges1d, Uthermal[cells].reshape((Nplot,Nplot)), rasterized=True, cmap=plt.get_cmap('magma'), vmin=-0.05, vmax=0.8 )
cax = plt.axes( [0.88,0.07,0.02,0.25] )
plt.colorbar( pc, cax=cax )
ax.set_xlim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
ax.set_ylim( 0.5*Boxsize-1.3*r_shock, 0.5*Boxsize+1.3*r_shock )
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
print(simulation_directory+"/plots/figure_%03d.pdf" % (i_file) )
fig.savefig(simulation_directory+"/plots/figure_%03d.pdf" % (i_file), dpi=300)
plt.close(fig)
#### check against analytic solution
i_compare, = np.where(Radius < 0.8) ## only check inner region; boundary has spourious effects in this testcase
......@@ -100,13 +149,14 @@ while True:
L1_dens = np.average(abs_delta_dens, weights=CellVolume[i_compare] )
L1_max = DeltaMaxAllowed * time
print("examples/Noh_2d/check snapshot %d: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
print("Noh_2d: snapshot %d: DEBUG: 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) )
if L1_dens > L1_max and time > 0.:
print("Noh_2d: ERROR: snaphshot %d: L1_dens = %g, DeltaMaxAllowed = %g!" % (i_file, L1_dens, L1_max) )
sys.exit(1)
i_file += 1
## if everything is ok
sys.exit(0)
""" @package ./examples/noh_2d/create.py
""" @package ./examples/Noh_2d/create.py
Code that creates 2d Noh test problem initial conditions
created by Rainer Weinberger, last modified 24.02.2019
......@@ -10,7 +10,7 @@ import numpy as np # scientific computing package
import h5py # hdf5 format
simulation_directory = str(sys.argv[1])
print("examples/noh_2d/create.py: creating ICs in directory " + simulation_directory)
print("examples/Noh_2d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
......@@ -95,4 +95,4 @@ part0.create_dataset("Velocities", data = Velocity)
part0.create_dataset("InternalEnergy", data = Uthermal)
## close file
IC.close()
IC.close()
\ No newline at end of file
%% examples/noh_2d/param.txt
%% examples/Noh_2d/param.txt
% parameter file for 2d Noh problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
#!/bin/bash # this line only there to enable syntax highlighting in this file
## examples/yee_2d/Config.sh
## examples/Yee_2d/Config.sh
## config file for 2d Yee vortex probelm
......@@ -21,6 +21,7 @@ TREE_BASED_TIMESTEPS # non-local timestep criterion (take 's
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_CENTER_OF_MASS # output centers of cells
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
......
""" @package ./examples/yee_2d/check.py
""" @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
created by Rainer Weinberger, last modified 07.01.2018 -- 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
import os # file specific calls
import matplotlib.pyplot as plt ## needs to be active for plotting!
plt.rcParams['text.usetex'] = True
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
simulation_directory = str(sys.argv[1])
print("examples/yee_2d/check.py: checking simulation output in directory " + simulation_directory)
print("Yee_2d: 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
......@@ -45,7 +55,9 @@ while True:
""" get simulation data """
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Time = FloatType(data["Header"].attrs["Time"])
Pos = np.array(data["PartType0"]["CenterOfMass"], dtype = FloatType)
VoronoiPos = 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)
......@@ -84,12 +96,54 @@ while True:
L1_utherm = np.average(abs_delta_utherm, weights=Volume)
""" printing results """
print("examples/Yee_2d/check.py: L1 error of " + filename +":")
print("Yee_2d: 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) )
if makeplots:
## plot:
fig = plt.figure( figsize=np.array([7.0,3.5]), dpi=300 )
Nplot = 256
from scipy import spatial # needed for KDTree that we use for nearest neighbour search
Edges1d = np.linspace(0., Boxsize, Nplot+1, endpoint=True, dtype=FloatType)
Grid1d = 0.5 * (Edges1d[1:] + Edges1d[:-1])
xx, yy = np.meshgrid(Grid1d, Grid1d)
Grid2D = np.array( [xx.reshape(Nplot**2), yy.reshape(Nplot**2)] ).T
vor = spatial.Voronoi( VoronoiPos[:,:2] )
dist, cells = spatial.KDTree( VoronoiPos[:,:2] ).query( Grid2D, k=1 )
ax = plt.axes( [0.05,0.15,0.35,0.7] )
pc = ax.pcolormesh( Edges1d, Edges1d, Density[cells].reshape((Nplot,Nplot)), rasterized=True )
cax = plt.axes( [0.42,0.15,0.01,0.7] )
plt.colorbar( pc, cax=cax )
ax.set_title( 'Density' )
spatial.voronoi_plot_2d( vor, ax=ax, line_colors='w', show_points=False, show_vertices=False, line_width=0.5 )
ax.set_xlim( 0, Boxsize )
ax.set_ylim( 0, Boxsize )
error = (Density[cells]-Density_ref[cells])/Density_ref[cells]
ax = plt.axes( [0.53,0.15,0.35,0.7] )
pc = ax.pcolormesh( Edges1d, Edges1d, error.reshape((Nplot,Nplot)), rasterized=True )
cax = plt.axes( [0.90,0.15,0.01,0.7] )
plt.colorbar( pc, cax=cax )
ax.set_title( 'Relative density error' )
spatial.voronoi_plot_2d( vor, ax=ax, line_colors='w', show_points=False, show_vertices=False, line_width=0.5 )
ax.set_xlim( 0, Boxsize )
ax.set_ylim( 0, Boxsize )
plt.text( 0.5, 0.92, "$t=%4.1f,\ L_1=%5.1e$" % (Time,L1_dens), ha='center', size=12, transform=fig.transFigure )
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
print(simulation_directory+"/plots/figure_%03d.pdf" % (i_file) )
fig.savefig(simulation_directory+"/plots/figure_%03d.pdf" % (i_file), dpi=300)
plt.close(fig)
""" criteria for failing the test """
if L1_dens > DeltaMaxAllowed or L1_vel > DeltaMaxAllowed or L1_utherm > DeltaMaxAllowed:
sys.exit(1)
......
""" @package ./examples/yee_2d/create.py
""" @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
......@@ -11,7 +11,7 @@ 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)
print("examples/Yee_2d/create.py: creating ICs in directory " + simulation_directory)
""" initial condition parameters """
FilePath = simulation_directory + '/IC.hdf5'
......@@ -20,7 +20,10 @@ FloatType = np.float64 # double precision: np.float64, for single use np.float3
IntType = np.int32
Boxsize = FloatType(10.0)
CellsPerDimension = IntType(50)
if len(sys.argv) > 3:
CellsPerDimension = IntType(sys.argv[3])
else:
CellsPerDimension = IntType(50)
#### parameters
Tinf = FloatType(1.0)
......
%% ./examples/yee_2d/param.txt
%% ./examples/Yee_2d/param.txt
% parameter file for 2d Yee vortex problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
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