Commit 012082c1 authored by Rainer Weinberger's avatar Rainer Weinberger

updated 1d examples; other minor comment corrections

parent 8afad291
......@@ -8,9 +8,16 @@ created by Rainer Weinberger, last modified: 04.03.2019
import sys # needed for exit codes
import numpy as np # scientific computing package
import h5py # hdf5 format
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
""" check functions """
def CheckL1Error(Pos, W, gamma, i_snap):
......@@ -28,9 +35,9 @@ def CheckL1Error(Pos, W, gamma, i_snap):
"""
try:
xx, rho, v, pres = np.loadtxt("./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap ).T
xx, rho, v, pres = np.loadtxt(simulation_directory+"/reference_%03d.txt" % i_snap ).T
except:
print("Check_L1_Error: could not find: ./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap)
print("Check_L1_Error: could not find: %s/reference_%03d.txt" % (simulation_directory,i_snap))
return 1
W_exact = np.array([rho, v, pres]).T
......@@ -89,22 +96,28 @@ def PlotSimulationData(Pos, W, gamma, i_snap, simulation_directory):
try:
## exact data, from 20,000 cell fixed grid simulation; downsampled by factor of 20
xx, rho, v, pres = np.loadtxt("./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap).T
xx, rho, v, pres = np.loadtxt(simulation_directory+"/reference_%03d.txt" % i_snap).T
except:
print("PlotSimulationData: could not find: ./examples/interacting_blastwaves_1d/reference_%03d.txt" % i_snap)
print("PlotSimulationData: could not find: %s/reference_%03d.txt" % (simulation_directory,i_snap))
return 1
## plot:
fig, ax = plt.subplots(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, W[:,0], ls="", marker='o')
ax[0].plot(xx, rho, color="k")
ax[1].plot(Pos, W[:,1], ls="", marker='o')
ax[1].plot(xx, v, color="k")
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), ls="", marker='o')
ax[2].plot(xx, pres / rho / (gamma - 1.0), color="k")
ax[3].plot(Pos, W[:,2], ls="", marker='o')
ax[3].plot(xx, pres, color="k")
ax[0].plot(Pos, W[:,0], 'r+', label="Arepo cells")
ax[1].plot(Pos, W[:,1], 'r+')
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), 'r+')
ax[3].plot(Pos, W[:,2], 'r+')
ax[0].plot(Pos, W[:,0], 'b', label="Arepo")
ax[1].plot(Pos, W[:,1], 'b')
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), 'b')
ax[3].plot(Pos, W[:,2], 'b')
ax[0].plot(xx, rho, color="k", lw=0.7, label="Exact solution")
ax[1].plot(xx, v, color="k", lw=0.7)
ax[2].plot(xx, pres / rho / (gamma - 1.0), color="k", lw=0.7)
ax[3].plot(xx, pres, color="k", lw=0.7)
for i_plot in np.arange(4):
ax[i_plot].set_ylim([ min[i_plot], max[i_plot] ])
......@@ -118,13 +131,15 @@ def PlotSimulationData(Pos, W, gamma, i_snap, simulation_directory):
ax[3].set_ylabel(r"pressure")
fig.align_ylabels(ax[:])
fig.savefig(simulation_directory+"/output/figure_%03d.pdf" % (i_file))
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
fig.savefig(simulation_directory+"/plots/figure_%03d.pdf" % (i_file))
plt.close(fig)
return 0
simulation_directory = str(sys.argv[1])
print("examples/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
print("wave_1d: checking simulation output in directory " + simulation_directory)
##parameters
Dtype = np.float64 # double precision: np.float64, for single use np.float32
......@@ -155,16 +170,16 @@ internalEnergy = np.array(data["PartType0"]["InternalEnergy"], dtype = Dtype)
W = np.array([density, vel[:,0], (gamma-1.0)*internalEnergy*density], dtype=Dtype).T ## shape: (n,3)
""" plot data if you want """
if createFigures:
if makeplots:
ReturnFlag += PlotSimulationData(position[:,0], W, gamma, i_file, simulation_directory)
""" perform checks """
ReturnFlag += CheckL1Error(position[:,0], W, gamma, i_file)
if ReturnFlag == 0:
print("examples/interacting_blastwaves_1d/check.py: success!")
print("interacting_blastwaves_1d: success!")
else:
print("examples/interacting_blastwaves_1d/check.py: failed!")
print("interacting_blastwaves_1d: failed!")
sys.exit(ReturnFlag)
......@@ -7,11 +7,11 @@ CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
%% initial conditions
InitCondFile IC
InitCondFile ./IC
ICFormat 3
%% output options
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
......@@ -8,15 +8,22 @@ created by Rainer Weinberger, last modified 12.03.2019
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
plt.rcParams['text.usetex'] = True
createFigures = False
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
CreateReferenceSolution = False
simulation_directory = str(sys.argv[1])
print("examples/mhd_shocktube_1d/check.py: checking simulation output in directory " + simulation_directory)
print("mhd_shocktube_1d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
......@@ -59,10 +66,10 @@ while True:
if CreateReferenceSolution:
checkData = verificationData[::40,:]
np.savetxt('examples/mhd_shocktube_1d/data_%03d.txt'%i_file,checkData)
np.savetxt(simulation_directory+'/data_%03d.txt'%i_file,checkData)
status = 0
else:
checkData = np.loadtxt('examples/mhd_shocktube_1d/data_%03d.txt'%i_file)
checkData = np.loadtxt(simulation_directory+'/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])
......@@ -75,26 +82,39 @@ while True:
delta_u = u - u_ref
delta_B = absB - absB_ref
if createFigures:
if makeplots:
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].plot(x[:,0], rho, 'b', zorder=3, label="Arepo")
ax[0].plot(x[:,0], rho, 'r+', zorder=2, label="Arepo cells")
ax[0].plot(checkData[:,0], checkData[:,1], c='k', zorder=1, lw=0.7, label="Exact solution")
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].plot(x[:,0], vel[:,0], 'b', zorder=3)
ax[1].plot(x[:,0], vel[:,0], 'r+', zorder=2)
ax[1].plot(checkData[:,0], checkData[:,2], c='k', zorder=1, lw=0.7)
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].plot(x[:,0], u, 'b', zorder=3)
ax[2].plot(x[:,0], u, 'r+', zorder=2)
ax[2].plot(checkData[:,0], checkData[:,3], c='k', zorder=1, lw=0.7)
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].plot(x[:,0], absB, 'b', zorder=3)
ax[3].plot(x[:,0], absB, 'r+', zorder=2)
ax[3].plot(checkData[:,0], checkData[:,4], c='k', zorder=1, lw=0.7)
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)
ax[0].legend( loc='upper right', frameon=False, fontsize=8 )
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
fig.savefig(simulation_directory+'/plots/figure_%03d.pdf'%i_file)
......@@ -123,4 +143,4 @@ while True:
i_file += 1
exit(status)
\ No newline at end of file
exit(status)
%% examples/mhd_shocktube_1d.txt
% parameter file for 1d mhd shocktube problem
InitCondFile IC
ICFormat 3
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
......@@ -8,9 +8,16 @@ created by Rainer Weinberger, last modified 04.03.2019
import sys # needed for exit codes
import numpy as np # scientific computing package
import h5py # hdf5 format
import matplotlib.pyplot as plt # plot results
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/polytrope_1d_spherical/check.py: checking simulation output in directory " + simulation_directory)
......@@ -62,29 +69,31 @@ while True:
""" plots """
if createFigures:
if makeplots:
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].plot(Pos, Density, 'b', label='evolved profile')
ax[0].plot(Pos, Density, 'r+', label='Arepo cells')
ax[0].plot(Pos_ref, Density_ref, 'k--', label='initial profile', lw=0.7)
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[0].set_ylim( 0., 1. )
ax[1].plot(Pos, Density/Density_ref, 'r')
ax[1].plot(Pos_ref, [1.0]*Pos_ref.shape[0], 'b--')
ax[1].plot(Pos, Density/Density_ref, 'b')
ax[1].plot(Pos_ref, [1.0]*Pos_ref.shape[0], 'k--', lw=0.7)
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].plot(Pos, Velocity, 'b')
ax[2].plot(Pos_ref, [0.0]*Pos_ref.shape[0], 'k--', lw=0.7)
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].plot(Pos[:-1], Accel[:-1] - GradPress[:-1] / Density[:-1], 'b')
ax[3].plot(Pos_ref[:-1], Accel_ref[:-1] - GradPress_ref[:-1] / Density_ref[:-1], 'k--', lw=0.7)
ax[3].set_ylim([-0.5,0.1])#
ax[3].set_ylabel(r"net accel.")
......@@ -96,7 +105,9 @@ while True:
fig.align_ylabels(ax[:])
fig.savefig(simulation_directory+"profiles_%03d.pdf"%i_snap)
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
fig.savefig(simulation_directory+"plots/profiles_%03d.pdf"%i_snap)
""" compare to ICs by interpolating to IC positions """
......
%% ./examples/polytrope_1d_spherical/param.txt
% parameter file for 1d polytrope test problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
......@@ -8,12 +8,20 @@ created by Rainer Weinberger, last modified: 19.02.2019
import sys # needed for exit codes
import numpy as np # scientific computing package
import h5py # hdf5 format
import os # file specific calls
import matplotlib.pyplot as plt ## needs to be active for plotting!
plt.rcParams['text.usetex'] = True
from Riemann import * ## Riemann-solver
createFigures = False
forceExitOnError = False ## exits immediately when tolerance is exceeded
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
forceExitOnError=False ## exits immediately when tolerance is exceeded
""" check functiions """
def CheckL1Error(Pos, W, W_L, W_R, gamma, position_0, time):
......@@ -188,44 +196,55 @@ def PlotSimulationData(Pos, W, W_L, W_R, gamma, position_0, time, simulation_dir
fig, ax = plt.subplots(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)
nskip = 0
i_show = np.arange(nskip, len(density)-nskip)
ax[0].plot(Pos, W[:,0], ls="", marker='o')
ax[0].plot(xx, W_exact[:,0], color="k")
ax[1].plot(Pos, W[:,1], ls="", marker='o')
ax[1].plot(xx, W_exact[:,1], color="k")
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), ls="", marker='o')
ax[2].plot(xx, W_exact[:,2]/W_exact[:,0]/(gamma - 1.0), color="k")
ax[3].plot(Pos, W[:,2], ls="", marker='o')
ax[3].plot(xx, W_exact[:,2], color="k")
ax[0].plot(xx, W_exact[:,0], color="k", lw=0.7, label="Exact solution")
ax[1].plot(xx, W_exact[:,1], color="k", lw=0.7)
ax[2].plot(xx, W_exact[:,2]/W_exact[:,0]/(gamma - 1.0), color="k", lw=0.7)
ax[3].plot(xx, W_exact[:,2], color="k", lw=0.7)
ax[0].plot(Pos, W[:,0], '+r', mec='r', mfc="None", label="Arepo cells")
ax[1].plot(Pos, W[:,1], '+r', mec='r', mfc="None")
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), '+r', mec='r', mfc="None")
ax[3].plot(Pos, W[:,2], '+r', mec='r', mfc="None")
ax[0].plot(Pos, W[:,0], 'b', label="Arepo")
ax[1].plot(Pos, W[:,1], 'b')
ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), 'b')
ax[3].plot(Pos, W[:,2], 'b')
ax[3].set_xlim([0,20])
for i_plot in np.arange(4):
ax[i_plot].set_ylim([ min[i_plot], max[i_plot] ])
ax[0].legend( loc='upper right', frameon=False, fontsize=8 )
## set labels
ax[3].set_xlabel(r"pos")
ax[0].set_ylabel(r"density")
ax[1].set_ylabel(r"velocity")
ax[2].set_ylabel(r"spec. int. energy")
ax[3].set_ylabel(r"pressure")
fig.align_ylabels(ax[:])
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
print(simulation_directory+"/output/figure_%03d.pdf" % (i_file) )
fig.savefig(simulation_directory+"/output/figure_%03d.pdf" % (i_file))
print(simulation_directory+"/plots/figure_%03d.pdf" % (i_file) )
fig.savefig(simulation_directory+"/plots/figure_%03d.pdf" % (i_file))
plt.close(fig)
return 0
simulation_directory = str(sys.argv[1])
print("examples/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
print("wave_1d: checking simulation output in directory " + simulation_directory)
Dtype = np.float64 # double precision: np.float64, for single use np.float32
## open initial conditiions to get parameters
directory = simulation_directory+""
directory = simulation_directory+"/"
filename = "IC.hdf5"
try:
data = h5py.File(directory+filename, "r")
except:
print( "Could not find file {0}".format(directory+filename) )
sys.exit(1)
IC_position = np.array(data["PartType0"]["Coordinates"], dtype = np.float64)
IC_mass = np.array(data["PartType0"]["Masses"], dtype = np.float64)
......@@ -270,7 +289,7 @@ while True:
""" plot data if you want """
if createFigures:
if makeplots:
ReturnFlag += PlotSimulationData(position[:,0], W, W_L, W_R, gamma, position_0, time, simulation_directory)
if ReturnFlag > 0 and forceExitOnError:
print('ERROR: something went wrong in plot')
......
......@@ -7,11 +7,11 @@ CpuTimeBetRestartFile 9000
TimeLimitCPU 90000
%% initial conditions
InitCondFile IC
InitCondFile ./IC
ICFormat 3
%% output options
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
......@@ -18,6 +18,7 @@ FORCE_EQUAL_TIMESTEPS # variable but global timestep
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)
......
......@@ -8,13 +8,21 @@ created by Rainer Weinberger, last modified 19.2.2019 -- comments welcome
import sys # system specific calls
import numpy as np # scientific computing package
import h5py # hdf5 format
import os # file specific calls
simulation_directory = str(sys.argv[1])
print("examples/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
print("wave_1d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32 # integer type
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
""" open initial conditiions to get parameters """
try:
data = h5py.File(simulation_directory + "/IC.hdf5", "r")
......@@ -52,10 +60,12 @@ while True:
""" get simulation data """
## simulation data
Pos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType)
Pos = np.array(data["PartType0"]["CenterOfMass"], 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)
Uthermal = np.array(data["PartType0"]["InternalEnergy"], dtype = FloatType)
Volume = Mass / Density
""" calculate analytic solution at new cell positions """
Density_ref = np.full(Pos.shape[0], density_0, dtype=FloatType)
......@@ -69,18 +79,18 @@ while True:
""" compare data """
## density
abs_delta_dens = np.abs(Density - Density_ref) / Density_ref
L1_dens = np.average(abs_delta_dens)
L1_dens = np.average(abs_delta_dens, weights=Volume)
## velocity, here, use absolute error (velocity should be zero!)
abs_delta_vel = np.abs(Velocity - Velocity_ref)
L1_vel = np.average(abs_delta_vel)
## velocity, here, use absolute error (velocity should be zero! check only x-vel, taking all components dilutes the norm!)
abs_delta_vel = np.abs(Velocity - Velocity_ref)[:,0]
L1_vel = np.average(abs_delta_vel, weights=Volume)
## internal energy
abs_delta_utherm = np.abs(Uthermal - Uthermal_ref) / Uthermal_ref
L1_utherm = np.average(abs_delta_utherm)
L1_utherm = np.average(abs_delta_utherm, weights=Volume)
""" printing results """
print("examples/wave_1d/check.py: L1 error of " + filename +":")
print("wave_1d: L1 error of " + filename +":")
print("\t density: %g" % L1_dens)
print("\t velocity: %g" % L1_vel)
print("\t specific internal energy: %g" % L1_utherm)
......@@ -89,6 +99,40 @@ while True:
""" criteria for failing the test """
if L1_dens > DeltaMaxAllowed or L1_vel > DeltaMaxAllowed or L1_utherm > DeltaMaxAllowed:
sys.exit(1)
if makeplots and i_file > 0:
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
# only import matplotlib if needed
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
f = plt.figure( figsize=(3.5,3.5) )
ax = plt.axes( [0.19, 0.12, 0.75, 0.75] )
refpos = np.linspace( 0, Boxsize, 257 )
ax.plot( refpos, density_0 * (1. + delta * np.sin(2.0 * np.pi * refpos / Boxsize)), 'k', lw=0.7, label="Analytical solution" )
ax.plot( Pos[:,0], Density, 'o-r', mec='r', mfc="None", label="Arepo" )
ax.set_xlim( 0, Boxsize )
ax.set_xlabel( "x" )
ax.set_ylabel( "Density deviation" )
ax.legend( loc='upper right', frameon=False, fontsize=8 )
ax.set_title( "$\mathrm{wave\_1d:}\ \mathrm{N}=%d,\ \mathrm{L1}=%4.1e$" % (NumberOfCells,L1_dens), loc='right', size=8 )
plt.ticklabel_format( axis='y', style='sci', scilimits=(0,0) )
f.savefig( simulation_directory+"plots/density.pdf" )
f = plt.figure( figsize=(3.5,3.5) )
ax = plt.axes( [0.19, 0.12, 0.75, 0.75] )
refpos = np.linspace( 0, Boxsize, 257 )
ax.plot( refpos, np.ones(257)*velocity_0, 'k', lw=0.7, label="Analytical solution" )
ax.plot( Pos[:,0], Velocity[:,0], 'o-r', mec='r', mfc="None", label="Arepo" )
ax.set_xlim( 0, Boxsize )
ax.set_xlabel( "x" )
ax.set_ylabel( "Velocity" )
ax.legend( loc='lower center', frameon=False, fontsize=8 )
ax.set_title( "$\mathrm{wave\_1d:}\ \mathrm{N}=%d,\ \mathrm{L1}=%4.1e$" % (NumberOfCells,L1_vel), loc='right', size=8 )
plt.ticklabel_format( axis='y', style='sci', scilimits=(0,0) )
f.savefig( simulation_directory+"plots/velocity.pdf" )
i_file += 1
""" normal exit """
......
......@@ -22,7 +22,10 @@ IntType = np.int32 # integer type
## computational domain
Boxsize = FloatType(1.0)
NumberOfCells = IntType(32)
if len(sys.argv) > 3:
NumberOfCells = IntType(sys.argv[3])
else:
NumberOfCells = IntType(32)
## initial state
density_0 = FloatType(1.0)
......
%% ./examples/wave_1d/param.txt
% parameter file for 1d wave propagation problem
InitCondFile IC
InitCondFile ./IC
ICFormat 3
OutputDir output
OutputDir ./output/
SnapshotFileBase snap
SnapFormat 3
NumFilesPerSnapshot 1
......
......@@ -238,7 +238,7 @@ static int subfind_so_potegy_loctree_treebuild(int npart, int start)
/* seems like we're dealing with particles at identical (or extremely close)
* locations. Shift subnode index to allow tree construction. Note: Multipole moments
* of tree are still correct, but one should MAX_TREE_LEVEL large enough to have
* DomainLen/2^MAX_TREE_LEEL < gravitational softening length
* DomainLen/2^MAX_TREE_LEVEL < gravitational softening length
*/
for(int j = 0; j < 8; j++)
{
......@@ -345,7 +345,9 @@ static int subfind_so_potegy_loctree_treebuild(int npart, int start)
/*! \brief Walk the tree and update node data recursively.
*
* This routine computes the multipole moments for a given internal node and
* all its subnodes using a recursive computation.
* all its subnodes using a recursive computation. Note that this switches
* the information stored in LocNodes[no].u from suns to d!
*
*
* \param[in] no Node index.
* \param[in] sib Sibling index.
......
......@@ -46,7 +46,7 @@ do
mkdir -p ${RUNDIR}
## copy Config and parameter file to run directory
cp ${DIR}/* ${RUNDIR}
cp -r ${DIR}/* ${RUNDIR}
## create ICs in run directory
echo ${DIR}
......
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