Better looking output of solution script for Noh problems.

examples/Noh_1D/makeIC.py

@@ -20,7 +20,7 @@
 import h5py
 from numpy import *
 
-# Generates a swift IC file for the Noh problem test in a periodic cubic box
+# Generates a swift IC file for the 1D Noh problem test in a periodic box
 
 # Parameters
 numPart = 1001

examples/Noh_1D/plotSolution.py

@@ -160,6 +160,7 @@ ylim(-0.05, 0.2)
 subplot(236, frameon=False)
 
 text(-0.49, 0.9, "Noh problem with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
 plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
 text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
 text(-0.49, 0.4, scheme, fontsize=10)

examples/Noh_2D/makeIC.py

@@ -21,7 +21,7 @@
 import h5py
 from numpy import *
 
-# Generates a swift IC file for the 2D Sod Shock in a periodic box
+# Generates a swift IC file for the 2D Noh problem in a periodic box
 
 # Parameters
 gamma = 5./3.          # Gas adiabatic index

examples/Noh_2D/plotSolution.py

@@ -29,6 +29,7 @@ v0 = 1
 import matplotlib
 matplotlib.use("Agg")
 from pylab import *
+from scipy import stats
 import h5py
 
 # Plot parameters
@@ -56,7 +57,6 @@ rc('font',**{'family':'sans-serif','sans-serif':['Times']})
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
 sim = h5py.File("noh_%03d.hdf5"%snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
@@ -79,15 +79,30 @@ rho = sim["/PartType0/Density"][:]
 r = np.sqrt((x-1)**2 + (y-1)**2)
 v = -np.sqrt(vx**2 + vy**2)
 
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+# Analytic solution
 N = 1000  # Number of points
 x_min = -1
 x_max = 1
 
-#x += x_min
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
 x_s = np.arange(0, 2., 2./N) - 1.
 rho_s = np.ones(N) * rho0
 P_s = np.ones(N) * rho0
@@ -118,8 +133,9 @@ figure()
 
 # Velocity profile --------------------------------
 subplot(231)
-plot(r, v, '.', color='r', ms=4.0)
+plot(r, v, '.', color='r', ms=0.5, alpha=0.2)
 plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
 ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
 xlim(0, 0.5)
@@ -127,8 +143,9 @@ ylim(-1.2, 0.4)
 
 # Density profile --------------------------------
 subplot(232)
-plot(r, rho, '.', color='r', ms=4.0)
+plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
 plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
 ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
 xlim(0, 0.5)
@@ -136,8 +153,9 @@ ylim(0.95, 19)
 
 # Pressure profile --------------------------------
 subplot(233)
-plot(r, P, '.', color='r', ms=4.0)
+plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
 plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
 ylabel("${\\rm{Pressure}}~P$", labelpad=0)
 xlim(0, 0.5)
@@ -145,8 +163,9 @@ ylim(-0.5, 11)
 
 # Internal energy profile -------------------------
 subplot(234)
-plot(r, u, '.', color='r', ms=4.0)
+plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
 plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
 ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
 xlim(0, 0.5)
@@ -154,8 +173,9 @@ ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
 subplot(235)
-plot(r, S, '.', color='r', ms=4.0)
+plot(r, S, '.', color='r', ms=0.5, alpha=0.2) 
 plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
 ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
 xlim(0, 0.5)
@@ -165,6 +185,7 @@ ylim(-0.05, 0.2)
 subplot(236, frameon=False)
 
 text(-0.49, 0.9, "Noh problem with  $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
 plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
 text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
 text(-0.49, 0.4, scheme, fontsize=10)