Commit 0444abba authored by Rainer Weinberger's avatar Rainer Weinberger

Updated examples to the published version of the paper

parent f278c3f0
......@@ -351,7 +351,7 @@ build: $(EXEC)
$(EXEC): $(OBJS)
$(LINKER) $(OPTIMIZE) $(OBJS) $(LIBS) -o $(EXEC)
lib$(LIBRARY).a: $(filter-out $(BUILD_DIR)/main.o,$(OBJS))
lib$(LIBRARY).a: $(filter-out $(BUILD_DIR)/main/main.o,$(OBJS))
$(AR) -rcs lib$(LIBRARY).a $(OBJS)
clean:
......
......@@ -44,9 +44,9 @@ if os.path.exists( "{0}/../../../music".format(sys.path[0]) ):
call(['cp','-r',"{0}/../../../music".format(sys.path[0]),simulation_directory+'/music'])
else:
print('CREATE: did not fine dir {0}/../../music, trying to clone from bitbucket.'.format(sys.path[0]))
status = call(['hg', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music'])
status = call(['git', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music'])
if status != 0:
print('CREATE: ERROR: hg clone failed!')
print('CREATE: ERROR: git clone failed!')
sys.exit(status)
cwd = os.getcwd()
os.chdir(simulation_directory+'/music/')
......
""" @package ./examples/current_sheet_2d/check.py
Code that checks results of 2d current sheet problem
created by Rainer Weinberger, last modified 13.03.2019
created by Rainer Weinberger, last modified 03.07.2020
"""
""" load libraries """
......@@ -34,16 +34,6 @@ Boxsize = FloatType(data["Header"].attrs["BoxSize"])
NumberOfCells = np.int32(data["Header"].attrs["NumPart_Total"][0])
CellsPerDimension = np.sqrt(NumberOfCells) ## 2d sim
## parameters for initial state
density_0 = 1.0
velocity_radial_0 = -1.0 ## radial inflow velocity
pressure_0 = 1.0e-4
gamma = 5./3. ## note: this has to be consistent with the parameter settings for Arepo!
utherm_0 = pressure_0 / ( gamma - 1.0 ) / density_0
## maximum L1 error after one propagation; based on tests
DeltaMaxAllowed = 0.1 * (FloatType(CellsPerDimension) / 150.0)**-1
Time = []
MagneticEnergy = []
......@@ -61,13 +51,16 @@ while True:
""" get simulation data """
B = np.array(data["PartType0"]["MagneticField"], dtype = FloatType)
Mass = np.array(data["PartType0"]["Masses"], dtype = FloatType)
Density = np.array(data["PartType0"]["Density"], dtype = FloatType)
Volume = Mass/Density
B2 = B[:,0]*B[:,0] + \
B[:,1]*B[:,1] + \
B[:,2]*B[:,2]
Time.append( FloatType(data["Header"].attrs["Time"]) )
MagneticEnergy.append( np.sum(B2) / 8.0 / np.pi )
MagneticEnergy.append( np.sum(B2 * Volume) / 8.0 / np.pi )
i_file += 1
......@@ -75,6 +68,9 @@ while True:
Time = np.array(Time)
MagneticEnergy = np.array(MagneticEnergy)
data= np.array([Time, MagneticEnergy]).T
np.savetxt(simulation_directory+'/magnetic_energy_vs_time.txt', data)
if makeplots:
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
......@@ -115,7 +111,7 @@ if makeplots:
plt.close(fig)
## if everything is ok
if np.min(MagneticEnergy/MagneticEnergy[0]) > 0.92:
if np.min(MagneticEnergy/MagneticEnergy[0]) > 0.7:
sys.exit(0)
else:
sys.exit(1)
\ No newline at end of file
sys.exit(1)
......@@ -3,7 +3,7 @@ Code that creates 2d current sheet initial conditions
literature reference: Gardiner and Stone (2005), JCoPh.205..509G
created by Rainer Weinberger, last modified 13.03.2019 -- comments welcome
created by Rainer Weinberger, last modified 03.07.2020 -- comments welcome
"""
""" load libraries """
......@@ -31,7 +31,7 @@ vpert = 0.1
pressure_0 = 0.1
gamma = FloatType(5.0/3.0)
gamma_minus_one = FloatType(gamma - 1.0)
b0 = FloatType(1.0)
b0 = FloatType(np.sqrt(4.0 * np.pi))
""" set up grid """
......
......@@ -124,7 +124,7 @@ def PlotSimulationData(Pos, W, gamma, i_snap, simulation_directory):
ax[i_plot].set_xlim([0.5,0.9])
## set labels
ax[3].set_xlabel(r"pos")
ax[3].set_xlabel(r"position")
ax[0].set_ylabel(r"density")
ax[1].set_ylabel(r"velocity")
ax[2].set_ylabel(r"spec. int. energy")
......@@ -139,7 +139,7 @@ def PlotSimulationData(Pos, W, gamma, i_snap, simulation_directory):
return 0
simulation_directory = str(sys.argv[1])
print("wave_1d: checking simulation output in directory " + simulation_directory)
print("interacting_blastwaves_1d: checking simulation output in directory " + simulation_directory)
##parameters
Dtype = np.float64 # double precision: np.float64, for single use np.float32
......
-1.0623633232 1.0000000000 0.0000000000 0.0000000000 0.0000000000 1.0000000000 1.0000000000 0.0000000000 1.0000000000
-0.7082423488 1.0000000000 0.0000000000 0.0000000000 0.0000000000 1.0000000000 1.0000000000 0.0000000000 1.0000000000
-0.6715409130 0.9640808586 0.0642650283 -0.0308081080 0.0000000000 1.0000000000 0.9468580270 0.0000000000 0.9408544290
-0.6353895336 0.9294519019 0.1275314856 -0.0626154426 0.0000000000 1.0000000000 0.8947854689 0.0000000000 0.8852070567
-0.5997767072 0.8960667875 0.1898128326 -0.0955160345 0.0000000000 1.0000000000 0.8436668070 0.0000000000 0.8328509799
-0.5646916186 0.8638808379 0.2511220299 -0.1296218774 0.0000000000 1.0000000000 0.7933768598 0.0000000000 0.7835915331
-0.5301232705 0.8328509799 0.3114715056 -0.1650681126 0.0000000000 1.0000000000 0.7437772171 0.0000000000 0.7372455645
-0.4960603918 0.8029356877 0.3708731151 -0.2020203061 0.0000000000 1.0000000000 0.6947112374 0.0000000000 0.6936407547
-0.4624913156 0.7740949272 0.4293380879 -0.2406849513 0.0000000000 1.0000000000 0.6459968167 0.0000000000 0.6526149762
-0.4294038152 0.7462901020 0.4868769612 -0.2813251383 0.0000000000 1.0000000000 0.5974155685 0.0000000000 0.6140156908
-0.3967848765 0.7194840023 0.5434994912 -0.3242848793 0.0000000000 1.0000000000 0.5486959577 0.0000000000 0.5776993822
-0.3646203767 0.6936407547 0.5992145375 -0.3700287563 0.0000000000 1.0000000000 0.4994856721 0.0000000000 0.5435310224
-0.2405925815 0.6936407547 0.5992145375 -0.3700287563 0.0000000000 1.0000000000 0.4994856721 0.0000000000 0.5435310224
-0.2405917815 0.6936407547 0.5992145375 -1.5507859988 0.1486044609 1.0000000000 -0.4839091543 0.1237653709 0.5435310224
-0.1513521286 0.6936407547 0.5992145375 -1.5507859988 0.1486044609 1.0000000000 -0.4839091543 0.1237653709 0.5435310224
-0.1513513286 0.7682366135 0.5042900674 -1.4363673835 0.1193405749 1.0000000000 -0.4063220123 0.1039215607 0.6446935393
0.2017156269 0.7682366135 0.5042900674 -1.4363673835 0.1193405749 1.0000000000 -0.4063220123 0.1039215607 0.6446935393
0.2017164268 0.3644655094 0.5042900671 -1.4363673865 0.1193405756 1.0000000000 -0.4063220141 0.1039215612 0.6446935394
0.5511846187 0.3644655094 0.5042900671 -1.4363673865 0.1193405756 1.0000000000 -0.4063220141 0.1039215612 0.6446935394
0.5511854187 0.1729621047 -0.4630389023 -0.4864779563 -0.1236046479 1.0000000000 -0.7087891172 0.1812810260 0.1569995009
0.7765834587 0.1729621047 -0.4630389023 -0.4864779563 -0.1236046479 1.0000000000 -0.7087891172 0.1812810260 0.1569995009
0.7765842587 0.1729621047 -0.4630389023 -0.4492235434 0.0640352632 1.0000000000 -0.7242827480 0.1032440019 0.1569995009
1.0453886576 0.1729621047 -0.4630389023 -0.4492235434 0.0640352632 1.0000000000 -0.7242827480 0.1032440019 0.1569995009
1.0720849906 0.1754926228 -0.4181953008 -0.3991472028 0.0568970539 1.0000000000 -0.7512180877 0.1070835415 0.1608464367
1.0989639531 0.1780601637 -0.3730326565 -0.3504444547 0.0499546456 1.0000000000 -0.7779870345 0.1108993623 0.1647876334
1.1260269931 0.1806652688 -0.3275489601 -0.3030108818 0.0431931537 1.0000000000 -0.8046279407 0.1146969314 0.1688254006
1.1532755626 0.1833084879 -0.2817421790 -0.2567549978 0.0365995373 1.0000000000 -0.8311747912 0.1184810932 0.1729621047
1.1807111175 0.1859903786 -0.2356102576 -0.2115961184 0.0301622952 1.0000000000 -0.8576579447 0.1222561752 0.1772001698
1.2083351183 0.1887115066 -0.1891511172 -0.1674626627 0.0238712237 1.0000000000 -0.8841047230 0.1260260720 0.1815420796
1.2361490305 0.1914724460 -0.1423626564 -0.1242907820 0.0177172213 1.0000000000 -0.9105398866 0.1297943132 0.1859903786
1.2641543249 0.1942737793 -0.0952427514 -0.0820232442 0.0116921299 1.0000000000 -0.9369860197 0.1335641180 0.1905476736
1.2923524779 0.1971160974 -0.0477892555 -0.0406085154 0.0057886035 1.0000000000 -0.9634638472 0.1373384408 0.1952166354
1.3207457719 0.2000000000 0.0000000000 0.0000000000 0.0000000000 1.0000000000 -0.9899924966 0.1411200081 0.2000000000
1.9811184578 0.2000000000 0.0000000000 0.0000000000 0.0000000000 1.0000000000 -0.9899924966 0.1411200081 0.2000000000
""" @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
created by Rainer Weinberger, last modified 03.07.2020
"""
""" load libraries """
......@@ -19,13 +19,13 @@ if len(sys.argv) > 2:
makeplots = True
else:
makeplots = False
CreateReferenceSolution = False
simulation_directory = str(sys.argv[1])
print("mhd_shocktube_1d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
## open initial conditiions to get parameters
try:
......@@ -34,7 +34,7 @@ except:
print("could not open initial conditions!")
exit(-1)
Boxsize = FloatType(data["Header"].attrs["BoxSize"])
NumberOfCells = np.int32(data["Header"].attrs["NumPart_Total"][0])
NumberOfCells = IntType(data["Header"].attrs["NumPart_Total"][0])
CellsPerDimension = np.sqrt(NumberOfCells) ## 2d sim
......@@ -58,89 +58,80 @@ while True:
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])
absB = np.sqrt(B[:,0]*B[:,0] + B[:,1]*B[:,1] + B[:,2]*B[:,2]) / np.sqrt(4. * np.pi)
verificationData = np.array([x[:,0], rho, vel[:,0], u, absB]).T
if CreateReferenceSolution:
checkData = verificationData[::40,:]
np.savetxt(simulation_directory+'/data_%03d.txt'%i_file,checkData)
status = 0
else:
checkData = np.loadtxt(simulation_directory+'/data_%03d.txt'%i_file)
checkData = np.loadtxt(simulation_directory+'/alpha_3.dat')
checkData[:,0] += 1.0
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[:,8]/checkData[:,1]*3./2.)
absB_ref = np.interp(x[:,0], checkData[:,0], np.sqrt(checkData[:,5]**2+checkData[:,6]**2+checkData[:,7]**2) )
delta_rho = rho - rho_ref
delta_vel = vel[:,0] - vel_ref
delta_u = u - u_ref
delta_B = absB - absB_ref
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)
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])
ax[0].plot(x[:,0], rho, 'b', zorder=3, label="Arepo")
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], 'b', zorder=3)
ax[1].plot(checkData[:,0], checkData[:,2], c='k', zorder=1, lw=0.7)
ax[1].set_ylabel(r'vel')
delta_rho = rho - rho_ref
delta_vel = vel[:,0] - vel_ref
delta_u = u - u_ref
delta_B = absB - absB_ref
ax[2].plot(x[:,0], u, 'b', zorder=3)
ax[2].plot(checkData[:,0], checkData[:,8]/checkData[:,1]*3./2., c='k', zorder=1, lw=0.7)
ax[2].set_ylabel(r'u$_{th}$')
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, '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], '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, '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, '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'$\left|\mathbf{B}\right|$')
ax[3].set_xlabel(r'position')
ax[3].set_xlim([-0.1,2.6])
fig.align_ylabels(ax[:])
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)
ax[3].plot(x[:,0], absB, 'b', zorder=3)
ax[3].plot(checkData[:,0], np.sqrt(checkData[:,5]**2+checkData[:,6]**2+checkData[:,7]**2), c='k', zorder=1, lw=0.7)
ax[3].set_ylabel(r'$\left|\mathbf{B}\right|$')
ax[3].set_xlabel(r'position')
ax[3].set_xlim([-0.1,2.6])
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
fig.align_ylabels(ax[:])
ax[0].legend( loc='upper right', frameon=False, fontsize=8 )
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
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
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)
fig.savefig(simulation_directory+'/plots/figure_%03d.pdf'%i_file)
res_scaling = 300. / np.float64(len(rho))
tolerance_rho = 0.04 * res_scaling
tolerance_vel = 0.15 * res_scaling
tolerance_u = 0.2 * res_scaling
tolerance_B = 0.03 * res_scaling
i_compare = np.where((x<0.6) | (x>0.9))[0]
if np.std(delta_rho[i_compare]) > tolerance_rho:
status += 1
if np.std(delta_vel[i_compare]) > tolerance_vel:
status += 1
if np.std(delta_u[i_compare]) > tolerance_u:
status += 1
if np.std(delta_B[i_compare]) > tolerance_B:
status +=1
print('standard deviations of absolute error and tolerance (density, velocity, int. energy, magnetic field):')
print(np.std(delta_rho[i_compare]), tolerance_rho)
print(np.std(delta_vel[i_compare]), tolerance_vel)
print(np.std(delta_u[i_compare]), tolerance_u)
print(np.std(delta_B[i_compare]), tolerance_B)
i_file += 1
exit(status)
......@@ -2,7 +2,7 @@
Code that creates 1d mhd shocktube initial conditions
created by Rainer Weinberger, last modified 12.03.2019 -- comments welcome
created by Rainer Weinberger, last modified 03.07.2020 -- comments welcome
"""
""" load libraries """
......@@ -19,20 +19,20 @@ 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)
NumberOfCells = IntType(300)
alpha = np.pi
alpha = 3.0
## 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)
b_L = np.array([1.0, 1.0, 0.0], dtype=FloatType) * np.sqrt(4.0 * np.pi)
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)
b_R = np.array([1.0, np.cos(alpha), np.sin(alpha)], dtype=FloatType) * np.sqrt(4.0 * np.pi)
gamma = FloatType(5.0/3.0)
gamma_minus_one = FloatType(gamma - 1.0)
......@@ -66,7 +66,7 @@ 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]
Bfield[i_right,dim] = b_R[dim]
""" write *.hdf5 file; minimum number of fields required by Arepo """
......
This diff is collapsed.
This diff is collapsed.
......@@ -37,7 +37,7 @@ TimeLimitCPU 90000
TimeBetStatistics 0.005
TimeBegin 0.0
TimeMax 0.4
TimeBetSnapshot 0.2
TimeBetSnapshot 0.4
UnitVelocity_in_cm_per_s 1.0
UnitLength_in_cm 1.0
......
......@@ -102,6 +102,7 @@ while True:
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_xlim( 0.0, 0.9 )
ax.set_ylim( 0, 20 )
ax = plt.axes( [0.65,0.70,0.20,0.25] )
......@@ -114,6 +115,7 @@ while True:
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_xlim( 0.0, 0.9 )
ax.set_ylim( -1.1, 1.1 )
ax = plt.axes( [0.65,0.38,0.20,0.25] )
......@@ -127,6 +129,7 @@ while True:
ax.plot(Radius, Uthermal, 'b.', rasterized=True, ms=0.2)
ax.set_xlabel(r"Radius")
ax.set_ylabel(r"Spec. internal energy")
ax.set_xlim( 0.0, 0.9 )
ax.set_ylim( -0.05, 0.8 )
ax = plt.axes( [0.65,0.07,0.20,0.25] )
......
......@@ -105,6 +105,7 @@ while True:
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_xlim( 0.0, 0.9 )
ax.set_ylim( 0, 70 )
ax = plt.axes( [0.65,0.70,0.20,0.25] )
......@@ -117,6 +118,7 @@ while True:
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_xlim( 0.0, 0.9 )
ax.set_ylim( -1.1, 1.1 )
ax = plt.axes( [0.65,0.38,0.20,0.25] )
......@@ -130,6 +132,7 @@ while True:
ax.plot(Radius, Uthermal, 'b.', rasterized=True, ms=0.2)
ax.set_xlabel(r"Radius")
ax.set_ylabel(r"Spec. internal energy")
ax.set_xlim( 0.0, 0.9 )
ax.set_ylim( -0.05, 0.8 )
ax = plt.axes( [0.65,0.07,0.20,0.25] )
......
......@@ -93,10 +93,10 @@ def CheckTotalVariation(Pos, W, W_L, W_R, gamma, position_0, time):
MaxRatioAllowed = 1.01
if np.any( TotalVariationSim / TotalVariationExact > MaxRatioAllowed):
print("CheckTotalVariation: ERROR: TotalVariation Sim/Exact: %g %g %g, tolerance: %g" % (TotalVariationSim[0] / TotalVariationExact[0], TotalVariationSim[1] / TotalVariationExact[1], TotalVariationSim[2] / TotalVariationExact[2], MaxRatioAllowed) )
return 1
return 1, TotalVariationSim / TotalVariationExact
else:
print("CheckTotalVariation: TotalVariation Sim/Exact fine: %g %g %g, tolerance: %g" % (TotalVariationSim[0] / TotalVariationExact[0], TotalVariationSim[1] / TotalVariationExact[1], TotalVariationSim[2] / TotalVariationExact[2], MaxRatioAllowed) )
return 0
return 0, TotalVariationSim / TotalVariationExact
def CheckWidthOfDiscontinuities(Pos, W, W_L, W_R, gamma, position_0, time):
......@@ -116,6 +116,7 @@ def CheckWidthOfDiscontinuities(Pos, W, W_L, W_R, gamma, position_0, time):
"""
xx, W_exact, PosOfCharacteristics = RiemannProblem(Pos, position_0, W_L, W_R, gamma, time)
jump_per_cell_list = np.zeros(5)
for i, pos_char in enumerate(PosOfCharacteristics):
ReturnFlag = 0
## PosOfCharacteristics will give different values in index 0 and 1 (3 and 4)
......@@ -139,6 +140,7 @@ def CheckWidthOfDiscontinuities(Pos, W, W_L, W_R, gamma, position_0, time):
i_low = np.full(3, i_sorted[0], dtype=np.int32)
i_high = np.full(3, i_sorted[0], dtype=np.int32)
max_jump_per_cell = np.zeros(3, dtype=np.float64)
## check by how many cells 5th to 95th percentile of jump are sampled
for j in np.arange(3):
......@@ -146,21 +148,30 @@ def CheckWidthOfDiscontinuities(Pos, W, W_L, W_R, gamma, position_0, time):
i_low[j] -= 1
while W[i_low[j],j] > percentile_05[j] and W[i_low[j],j] < percentile_95[j]:
i_high[j] += 1
# max fractional change between two cells
jump_cells = np.arange(np.min((i_low[j], i_high[j])), np.max((i_low[j], i_high[j]))+2)
dW = W[jump_cells, j]-W[jump_cells-1, j]
max_jump_per_cell[j] = np.max( dW / jump[j])
## sufficient for exact Riemann solver
MaxNumerOfCells = 4
if(i == 2):
print("CheckWidthOfDiscontinuities: density jump at contact discontinuity resolved by %d cells (5th to 95th precentile), tolerance: %d" \
% (i_high[0]-i_low[0], MaxNumerOfCells) )
print("Contact density jump %g percent per cell"%(max_jump_per_cell[0]*100))
else:
print("CheckWidthOfDiscontinuities: density, velocity and pressure jump at shock resolved by %d, %d and %d cells (5th to 95th precentile), tolerance: %d" \
% (i_high[0]-i_low[0], i_high[1]-i_low[1], i_high[2]-i_low[2], MaxNumerOfCells) )
print("Shock jump %g %g %g percent per cell"%(max_jump_per_cell[0]*100, max_jump_per_cell[1]*100, max_jump_per_cell[2]*100))
if np.any(i_high-i_low > MaxNumerOfCells):
print("CheckWidthOfDiscontinuities: ERROR: discontinuity too wide!")
ReturnFlag += 1
return ReturnFlag
jump_per_cell_list[i] = max_jump_per_cell[0] #only record density jump
return ReturnFlag, jump_per_cell_list
def PlotSimulationData(Pos, W, W_L, W_R, gamma, position_0, time, simulation_directory):
"""
......@@ -235,7 +246,7 @@ def PlotSimulationData(Pos, W, W_L, W_R, gamma, position_0, time, simulation_dir
return 0
simulation_directory = str(sys.argv[1])
print("wave_1d: checking simulation output in directory " + simulation_directory)
print("shocktube_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
......@@ -264,6 +275,11 @@ position_0 = 0.5 * (IC_position[i_0[0]-1,0]+IC_position[i_0[0],0]) ## discontinu
""" loop over all output files """
i_file = -1
ReturnFlag = 0
time_arr = np.empty(shape=(0,1), dtype=np.float64)
jump_arr = np.empty(shape=(0,5), dtype=np.float64) # 5 characteristics
tv_arr = np.empty(shape=(0,3), dtype=np.float64)
while True:
i_file += 1
......@@ -301,16 +317,27 @@ while True:
print('ERROR: exceeding tolerance!')
sys.exit(ReturnFlag)
ReturnFlag += CheckTotalVariation(position[:,0], W, W_L, W_R, gamma, position_0, time)
flag, tv = CheckTotalVariation(position[:,0], W, W_L, W_R, gamma, position_0, time)
ReturnFlag += flag
if ReturnFlag > 0 and forceExitOnError:
print('ERROR: exceeding tolerance!')
sys.exit(ReturnFlag)
ReturnFlag += CheckWidthOfDiscontinuities(position[:,0], W, W_L, W_R, gamma, position_0, time)
flag, jump_per_cell_list = CheckWidthOfDiscontinuities(position[:,0], W, W_L, W_R, gamma, position_0, time)
ReturnFlag += flag
if ReturnFlag > 0 and forceExitOnError:
print('ERROR: exceeding tolerance!')
sys.exit(ReturnFlag)
time_arr = np.concatenate((time_arr, np.array([time], ndmin=2).T), axis=0)
jump_arr = np.concatenate((jump_arr, np.array(jump_per_cell_list, ndmin=2)), axis=0)
tv_arr = np.concatenate((tv_arr, np.array(tv, ndmin=2)), axis=0)
jump_arr = np.array(jump_arr, dtype=np.float64)
tv_arr = np.array(tv_arr, dtype=np.float64)
data = np.array([time_arr[:,0], jump_arr[:,2], jump_arr[:,3], tv_arr[:,0], tv_arr[:,1], tv_arr[:,2]])
np.savetxt(simulation_directory+'/jumps_total_variation.txt', data.T)
if ReturnFlag == 0:
print("check.py: success!")
else:
......
......@@ -85,7 +85,7 @@
* \param[in] st_R Right hand side state.
* \param[out] st_face State at face.
*
* \return Central pressure.
* \return 0.
*/
double godunov_flux_3d(struct state *st_L, struct state *st_R, struct state_face *st_face)
{
......
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