diff --git a/examples/current_sheet_2d/Config.sh b/examples/current_sheet_2d/Config.sh index cf1487b4e7440a8345f2484db4164f4930f37301..17bc564920f0c0b42432075a6b6dceee134014d9 100644 --- a/examples/current_sheet_2d/Config.sh +++ b/examples/current_sheet_2d/Config.sh @@ -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 diff --git a/examples/current_sheet_2d/check.py b/examples/current_sheet_2d/check.py index 838ba0a3ebea3d7ffcbf4e7e2dcd9bc6e7ded3c1..c3561330338ade8855708c42aefb618b0153ee38 100644 --- a/examples/current_sheet_2d/check.py +++ b/examples/current_sheet_2d/check.py @@ -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: diff --git a/examples/current_sheet_2d/param.txt b/examples/current_sheet_2d/param.txt index 514e76e63ee6a4c3cd590056c8f86cac810f28f1..ee2723f17d670087ed1db95e435d75826e2c0759 100644 --- a/examples/current_sheet_2d/param.txt +++ b/examples/current_sheet_2d/param.txt @@ -1,10 +1,10 @@ %% 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 diff --git a/examples/gresho_2d/Config.sh b/examples/gresho_2d/Config.sh index e8b941cc7d4c8c84d6e1aa16f21a89eeaba59229..0d39caf4d7b80968cf0102223e4ebc7d4a8dd7a9 100644 --- a/examples/gresho_2d/Config.sh +++ b/examples/gresho_2d/Config.sh @@ -1,13 +1,11 @@ #!/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) diff --git a/examples/gresho_2d/check.py b/examples/gresho_2d/check.py index 331f709c5cb9a6bd4f87f1f2db0d7ca1bcfc0130..b260ae467d6070c77e5e0e18acbab5f2765f4f31 100644 --- a/examples/gresho_2d/check.py +++ b/examples/gresho_2d/check.py @@ -1,4 +1,4 @@ -""" @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) diff --git a/examples/gresho_2d/create.py b/examples/gresho_2d/create.py index d5f9d5d35be95bb26c59cd307fd98df7efe16776..df68e21687dda831e478905ba6bd1ab3c618be17 100644 --- a/examples/gresho_2d/create.py +++ b/examples/gresho_2d/create.py @@ -1,4 +1,4 @@ -""" @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' diff --git a/examples/gresho_2d/param.txt b/examples/gresho_2d/param.txt index b850b472920b408513be909f41bfe585e37f98ce..13f6457de8a49b069ce06052896f989798994ce6 100644 --- a/examples/gresho_2d/param.txt +++ b/examples/gresho_2d/param.txt @@ -1,10 +1,10 @@ -%% ./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 diff --git a/examples/noh_2d/Config.sh b/examples/noh_2d/Config.sh index 77b11d30e7afa64d7fd6330e089e8bab99fc99d8..70222e470d79f7bf7311caba4d56b24dbdc8e7c0 100644 --- a/examples/noh_2d/Config.sh +++ b/examples/noh_2d/Config.sh @@ -1,6 +1,6 @@ #!/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) diff --git a/examples/noh_2d/check.py b/examples/noh_2d/check.py index e9d178d61d323d5049e2c25011c9f81d77be77c1..2fb8030bb4a3d62596583fca8e643b7cf977ee8c 100644 --- a/examples/noh_2d/check.py +++ b/examples/noh_2d/check.py @@ -1,4 +1,4 @@ -""" @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) diff --git a/examples/noh_2d/create.py b/examples/noh_2d/create.py index 7a5e1aedbffa6bd16b18e352fe635a2e91e67ee2..90906a09d6770d8143d5a0a872026cddd9d09853 100644 --- a/examples/noh_2d/create.py +++ b/examples/noh_2d/create.py @@ -1,4 +1,4 @@ -""" @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 diff --git a/examples/noh_2d/param.txt b/examples/noh_2d/param.txt index d78cf695b3ba286cf65f1bd0eea8bf4e09bd0230..699cbe669df5f17283b84bb697224c36b10e36d4 100644 --- a/examples/noh_2d/param.txt +++ b/examples/noh_2d/param.txt @@ -1,10 +1,10 @@ -%% 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 diff --git a/examples/yee_2d/Config.sh b/examples/yee_2d/Config.sh index 40824217a024da9bf60baca654706542289d3a0d..c68322c41ea2fa374d9501e18114c9eefa52c01e 100644 --- a/examples/yee_2d/Config.sh +++ b/examples/yee_2d/Config.sh @@ -1,6 +1,6 @@ #!/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) diff --git a/examples/yee_2d/check.py b/examples/yee_2d/check.py index 36508d0fc9520cd09e6a315492e6e0d4e14c10d8..3c76f915ffc74dae50c7584c5c94fa64e2b016db 100644 --- a/examples/yee_2d/check.py +++ b/examples/yee_2d/check.py @@ -1,16 +1,26 @@ -""" @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) diff --git a/examples/yee_2d/create.py b/examples/yee_2d/create.py index 37f1f65deee660d0198025240708e1fab7f05732..23f0532a44c85660ae43b06aa684de9e04f11d94 100644 --- a/examples/yee_2d/create.py +++ b/examples/yee_2d/create.py @@ -1,4 +1,4 @@ -""" @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) diff --git a/examples/yee_2d/param.txt b/examples/yee_2d/param.txt index 9c01a0dbf13bfcadccb4bf092f4112b3600913ba..3cfd812d156ea2de257814b95013e55496518415 100644 --- a/examples/yee_2d/param.txt +++ b/examples/yee_2d/param.txt @@ -1,10 +1,10 @@ -%% ./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