diff --git a/examples/Noh_1D/makeIC.py b/examples/Noh_1D/makeIC.py index 3e060e0e7422972f71db8e38a89d0776a7ea6ae1..176f3517455db7a8b0994ac7d1e65fb9cb7419d4 100644 --- a/examples/Noh_1D/makeIC.py +++ b/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 diff --git a/examples/Noh_1D/plotSolution.py b/examples/Noh_1D/plotSolution.py index 01cec16a8e3696e2f31433ec756b9be281a7ebcf..f4916af6e6066d21f76c28b5acef41e1907a83fd 100644 --- a/examples/Noh_1D/plotSolution.py +++ b/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) diff --git a/examples/Noh_2D/makeIC.py b/examples/Noh_2D/makeIC.py index 9b8604b87768b69fd7c21793e9e5580625a90a82..f7239fa3cd188637df929f86451d20a9978bd1f5 100644 --- a/examples/Noh_2D/makeIC.py +++ b/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 diff --git a/examples/Noh_2D/plotSolution.py b/examples/Noh_2D/plotSolution.py index 75b53fc76fb17a934946549d25f5895d6d79714c..70ded94b9bc3b921de5a0ab62df1d632bdf13e1c 100644 --- a/examples/Noh_2D/plotSolution.py +++ b/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)