Commit 1a8f00d5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Improved visual quality of kernel function plots

parent 3e463812
...@@ -26,31 +26,53 @@ from scipy.optimize import fsolve ...@@ -26,31 +26,53 @@ from scipy.optimize import fsolve
from matplotlib.font_manager import FontProperties from matplotlib.font_manager import FontProperties
import numpy import numpy
params = { params = {'axes.labelsize': 9,
'axes.labelsize': 10, 'axes.titlesize': 10,
'axes.titlesize': 8, 'font.size': 12,
'font.size': 10, 'legend.fontsize': 12,
'legend.fontsize': 9, 'xtick.labelsize': 9,
'xtick.labelsize': 10, 'ytick.labelsize': 9,
'ytick.labelsize': 10, 'text.usetex': True,
'xtick.major.pad': 2.5, 'figure.figsize' : (3.15,3.15),
'ytick.major.pad': 2.5, 'figure.subplot.left' : 0.17,
'text.usetex': True, 'figure.subplot.right' : 0.99 ,
'figure.figsize' : (4.15,4.15), 'figure.subplot.bottom' : 0.08 ,
'figure.subplot.left' : 0.14, 'figure.subplot.top' : 0.99 ,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.06,
'figure.subplot.top' : 0.99,
'figure.subplot.wspace' : 0. , 'figure.subplot.wspace' : 0. ,
'figure.subplot.hspace' : 0. , 'figure.subplot.hspace' : 0. ,
'lines.markersize' : 6, 'lines.markersize' : 6,
'lines.linewidth' : 1.5, 'lines.linewidth' : 3.,
'text.latex.unicode': True 'text.latex.unicode': True
} }
rcParams.update(params) rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']}) rc('font',**{'family':'sans-serif','sans-serif':['Times']})
# params = {
# 'axes.labelsize': 10,
# 'axes.titlesize': 8,
# 'font.size': 10,
# 'legend.fontsize': 9,
# 'xtick.labelsize': 10,
# 'ytick.labelsize': 10,
# 'xtick.major.pad': 2.5,
# 'ytick.major.pad': 2.5,
# 'text.usetex': True,
# 'figure.figsize' : (4.15,4.15),
# 'figure.subplot.left' : 0.14,
# 'figure.subplot.right' : 0.99,
# 'figure.subplot.bottom' : 0.06,
# 'figure.subplot.top' : 0.99,
# 'figure.subplot.wspace' : 0. ,
# 'figure.subplot.hspace' : 0. ,
# 'lines.markersize' : 6,
# 'lines.linewidth' : 1.5,
# 'text.latex.unicode': True
# }
# rcParams.update(params)
# rc('font',**{'family':'sans-serif','sans-serif':['Times']})
#Parameters #Parameters
eta = 1.2348422195325 # Resolution (Gives 48 neighbours for a cubic spline kernel) eta = 1.2348422195325 # Resolution (Gives 48 neighbours for a cubic spline kernel)
dx = 1.5#4 #2.7162 # Mean inter-particle separation dx = 1.5#4 #2.7162 # Mean inter-particle separation
...@@ -135,13 +157,13 @@ xx = linspace(0., 5*h, 1000) ...@@ -135,13 +157,13 @@ xx = linspace(0., 5*h, 1000)
maxY = 1.2*Gaussian(0, h) maxY = 1.2*Gaussian(0, h)
# Plot the kernels # Plot the kernels
plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm %14s\\quad H=\\infty}$"%("Gaussian~~~~~~")) plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm %14s~ H=\\infty}$"%("Gaussian~~~~~~"))
plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm %14s\\quad H=%4.3f}$"%("Cubic~spline~~", H_cubic)) plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm %14s~ H=%4.3f}$"%("Cubic~spline~~", H_cubic), lw=1.5)
plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm %14s\\quad H=%4.3f}$"%("Quartic~spline", H_quartic)) plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm %14s~ H=%4.3f}$"%("Quartic~spline", H_quartic), lw=1.5)
plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm %14s\\quad H=%4.3f}$"%("Quintic~spline", H_quintic)) plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm %14s~ H=%4.3f}$"%("Quintic~spline", H_quintic), lw=1.5)
plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C2~", H_WendlandC2)) plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm %14s~ H=%4.3f}$"%("Wendland~C2~", H_WendlandC2), lw=1.5)
plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C4~", H_WendlandC4)) plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm %14s~ H=%4.3f}$"%("Wendland~C4~", H_WendlandC4), lw=1.5)
plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C6~", H_WendlandC6)) plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm %14s~ H=%4.3f}$"%("Wendland~C6~", H_WendlandC6), lw=1.5)
# Indicate the position of H # Indicate the position of H
arrow(H_cubic, 0.12*maxY , 0., -0.12*maxY*0.9, fc='b', ec='b', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3) arrow(H_cubic, 0.12*maxY , 0., -0.12*maxY*0.9, fc='b', ec='b', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
...@@ -152,33 +174,33 @@ arrow(H_WendlandC4, 0.12*maxY , 0., -0.12*maxY*0.9, fc='m', ec='m', length_inclu ...@@ -152,33 +174,33 @@ arrow(H_WendlandC4, 0.12*maxY , 0., -0.12*maxY*0.9, fc='m', ec='m', length_inclu
arrow(H_WendlandC6, 0.12*maxY , 0., -0.12*maxY*0.9, fc='y', ec='y', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3) arrow(H_WendlandC6, 0.12*maxY , 0., -0.12*maxY*0.9, fc='y', ec='y', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
# Show h # Show h
plot([h, h], [0., maxY], 'k:', linewidth=0.5) plot([h, h], [0., 0.63*maxY], 'k:', linewidth=0.5)
text(h, maxY*0.35, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom') text(h, maxY*0.65, "$h\\equiv\\eta\\langle x\\rangle$", rotation=90, ha='center', va='bottom', fontsize=9)
# Show sigma # Show sigma
plot([h/2, h/2], [0., maxY], 'k:', linewidth=0.5) plot([h/2, h/2], [0., 0.63*maxY], 'k:', linewidth=0.5)
text(h/2, maxY*0.05, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom') text(h/2, maxY*0.65, "$\\sigma\\equiv h/2$", rotation=90, ha='center', va='bottom', fontsize=9)
# Show <x> # Show <x>
plot([dx, dx], [0., maxY], 'k:', linewidth=0.5) plot([dx, dx], [0., 0.63*maxY], 'k:', linewidth=0.5)
text(dx, maxY*0.35, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom') text(dx, maxY*0.65, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, ha='center', va='bottom', fontsize=9)
xlim(0., 2.5*h) xlim(0., 2.5*h)
ylim(0., maxY) ylim(0., maxY)
gca().xaxis.set_ticklabels([]) gca().xaxis.set_ticklabels([])
ylabel("$W(r,h)$", labelpad=1.5) ylabel("$W(r,h)$", labelpad=1.5)
legend(loc="upper right", handlelength=1.2, handletextpad=0.2) legend(loc="upper right", frameon=False, handletextpad=0.1, handlelength=1.2, fontsize=8)
# Same but now in log space # Same but now in log space
subplot(212, yscale="log") subplot(212, yscale="log")
plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$") plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$", lw=1.5)
plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$") plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$", lw=1.5)
plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$") plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$", lw=1.5)
plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$") plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$", lw=1.5)
plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$") plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$", lw=1.5)
plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$") plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$", lw=1.5)
plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$") plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$", lw=1.5)
# Show h # Show h
plot([h, h], [0., 1.], 'k:', linewidth=0.5) plot([h, h], [0., 1.], 'k:', linewidth=0.5)
...@@ -191,21 +213,22 @@ plot([dx, dx], [0., 1.], 'k:', linewidth=0.5) ...@@ -191,21 +213,22 @@ plot([dx, dx], [0., 1.], 'k:', linewidth=0.5)
# Show plot properties # Show plot properties
text(h/5., 1e-3, "$\\langle x \\rangle = %3.1f$"%(dx), va="top", backgroundcolor='w') text(h/5., 2e-3, "$\\langle x \\rangle = %3.1f$"%(dx), va="top", backgroundcolor='w', fontsize=9)
text(h/5.+0.06, 3e-4, "$\\eta = %5.4f$"%(eta), va="top", backgroundcolor='w') text(h/5.+0.06, 4e-4, "$\\eta = %5.4f$"%(eta), va="top", backgroundcolor='w', fontsize=9)
text(h/5.+0.06, 8e-5, "$h = \\eta\\langle x \\rangle = %5.4f$"%(h), va="top", backgroundcolor='w', fontsize=9)
# Show number of neighbours # Show number of neighbours
text(1.9*h, 2e-1/2.9**0, "$N_{\\rm ngb}=\\infty$", fontsize=10) text(1.75*h, 2e-1/2.9**0, "$N_{\\rm ngb}=\\infty$", fontsize=10)
text(1.9*h, 2e-1/2.9**1, "$N_{\\rm ngb}=%3.1f$"%(N_H_cubic), color='b', fontsize=9) text(1.75*h, 2e-1/2.9**1, "$N_{\\rm ngb}=%3.1f$"%(N_H_cubic), color='b', fontsize=9)
text(1.9*h, 2e-1/2.9**2, "$N_{\\rm ngb}=%3.1f$"%(N_H_quartic), color='c', fontsize=9) text(1.75*h, 2e-1/2.9**2, "$N_{\\rm ngb}=%3.1f$"%(N_H_quartic), color='c', fontsize=9)
text(1.9*h, 2e-1/2.9**3, "$N_{\\rm ngb}=%3.1f$"%(N_H_quintic), color='g', fontsize=9) text(1.75*h, 2e-1/2.9**3, "$N_{\\rm ngb}=%3.1f$"%(N_H_quintic), color='g', fontsize=9)
text(1.9*h, 2e-1/2.9**4, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC2), color='r', fontsize=9) text(1.75*h, 2e-1/2.9**4, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC2), color='r', fontsize=9)
text(1.9*h, 2e-1/2.9**5, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC4), color='m', fontsize=9) text(1.75*h, 2e-1/2.9**5, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC4), color='m', fontsize=9)
text(1.9*h, 2e-1/2.9**6, "$N_{\\rm ngb}=%3.0f$"%(N_H_WendlandC6), color='y', fontsize=9) text(1.75*h, 2e-1/2.9**6, "$N_{\\rm ngb}=%3.0f$"%(N_H_WendlandC6), color='y', fontsize=9)
xlim(0., 2.5*h) xlim(0., 2.5*h)
ylim(1e-5, 0.7) ylim(1e-5, 0.7)
xlabel("$r$", labelpad=0) xlabel("$r$", labelpad=-1)
ylabel("$W(r,h)$", labelpad=0.5) ylabel("$W(r,h)$", labelpad=0.5)
savefig("kernels.pdf") savefig("kernels.pdf")
...@@ -282,13 +305,13 @@ figure() ...@@ -282,13 +305,13 @@ figure()
subplot(211) subplot(211)
plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7) plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7)
plot(xx, d_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$") plot(xx, d_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$", lw=1.5)
plot(xx, dW_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$") plot(xx, dW_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$", lw=1.5)
plot(xx, dW_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$") plot(xx, dW_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$", lw=1.5)
plot(xx, dW_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$") plot(xx, dW_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$", lw=1.5)
plot(xx, dW_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$") plot(xx, dW_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$", lw=1.5)
plot(xx, dW_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$") plot(xx, dW_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$", lw=1.5)
plot(xx, dW_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$") plot(xx, dW_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$", lw=1.5)
maxY = d_Gaussian(h/2, h) maxY = d_Gaussian(h/2, h)
...@@ -307,7 +330,7 @@ gca().xaxis.set_ticklabels([]) ...@@ -307,7 +330,7 @@ gca().xaxis.set_ticklabels([])
ylim(1.1*maxY, -0.1*maxY) ylim(1.1*maxY, -0.1*maxY)
xlabel("$r$", labelpad=0) xlabel("$r$", labelpad=0)
ylabel("$\\partial W(r,h)/\\partial r$", labelpad=0.5) ylabel("$\\partial W(r,h)/\\partial r$", labelpad=0.5)
legend(loc="lower right") legend(loc="lower right", frameon=False, handletextpad=0.1, handlelength=1.2, fontsize=8)
...@@ -317,23 +340,23 @@ maxY = d2_Gaussian(h,h) ...@@ -317,23 +340,23 @@ maxY = d2_Gaussian(h,h)
plot([h/2, h/2], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5) plot([h/2, h/2], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
plot([h, h], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5) plot([h, h], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
plot([dx, dx], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5) plot([dx, dx], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
text(h/2, -3.*maxY, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom') text(h/2, -2.9*maxY, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom', fontsize=9)
text(h, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom') text(h, -2.9*maxY, "$h\\equiv\\eta\\langle x\\rangle$", rotation=90, backgroundcolor='w', ha='center', va='bottom', fontsize=9)
text(dx, -3.*maxY, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom') text(dx, -2.9*maxY, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom', fontsize=9)
plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7) plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7)
plot(xx, d2_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$") plot(xx, d2_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$")
plot(xx, d2W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$") plot(xx, d2W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$", lw=1.5)
plot(xx, d2W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$") plot(xx, d2W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$", lw=1.5)
plot(xx, d2W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$") plot(xx, d2W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$", lw=1.5)
plot(xx, d2W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$") plot(xx, d2W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$", lw=1.5)
plot(xx, d2W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$") plot(xx, d2W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$", lw=1.5)
plot(xx, d2W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$") plot(xx, d2W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$", lw=1.5)
xlim(0., 2.5*h) xlim(0., 2.5*h)
ylim(-3.2*maxY, 1.3*maxY) ylim(-3.2*maxY, 1.3*maxY)
xlabel("$r$", labelpad=0) xlabel("$r$", labelpad=-1)
ylabel("$\\partial^2 W(r,h)/\\partial r^2$", labelpad=0.5) ylabel("$\\partial^2 W(r,h)/\\partial r^2$", labelpad=0.5)
......
Supports Markdown
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