Commit ad6dadd8 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated plots and text decribing the kernel functions.

parent 4194ac5d
......@@ -80,9 +80,9 @@ tests/testRiemannHLLC
tests/testMatrixInversion
theory/latex/swift.pdf
theory/kernel/kernels.pdf
theory/kernel/kernel_derivatives.pdf
theory/kernel/kernel_definitions.pdf
theory/SPH/kernel/kernels.pdf
theory/SPH/kernel/kernel_derivatives.pdf
theory/SPH/kernel/kernel_definitions.pdf
theory/paper_pasc/pasc_paper.pdf
m4/libtool.m4
......
......@@ -18,7 +18,7 @@
#
##############################################################################
# Computes the analytical solution of the 2D Sedov blast wave.
# Computes the analytical solution of the 3D Sedov blast wave.
# The script works for a given initial box and dumped energy and computes the solution at a later time t.
# Parameters
......
......@@ -218,7 +218,7 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
\begin{figure}
\includegraphics[width=\columnwidth]{kernel_derivatives.pdf}
\caption{The first and secon derivatives of the kernel functions
\caption{The first and second derivatives of the kernel functions
available in \swift for a mean inter-particle separation $\langle
x\rangle=1.5$ and a resolution $\eta=1.2348$. A Gaussian kernel
with the same smoothing length is shown for comparison.}
......@@ -229,14 +229,13 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
\section{Kernel Derivatives}
The derivatives of the kernel function have relatively simple
expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}.
expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}:
\begin{eqnarray*}
\vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
\vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|}, \\
\frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) +
\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right].
\end{eqnarray*}
Note that for all the kernels listed above, $f'(0) = 0$.
\end{document}
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import matplotlib
matplotlib.use("Agg")
from pylab import *
......@@ -22,9 +38,9 @@ params = {
'text.usetex': True,
'figure.figsize' : (4.15,4.15),
'figure.subplot.left' : 0.14,
'figure.subplot.right' : 0.99 ,
'figure.subplot.bottom' : 0.08 ,
'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.hspace' : 0. ,
'lines.markersize' : 6,
......@@ -139,6 +155,10 @@ arrow(H_WendlandC6, 0.12*maxY , 0., -0.12*maxY*0.9, fc='y', ec='y', length_inclu
plot([h, h], [0., 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')
# Show sigma
plot([h/2, h/2], [0., maxY], 'k:', linewidth=0.5)
text(h/2, maxY*0.05, "$\\sigma\\equiv h/2$", rotation=90, backgroundcolor='w', ha='center', va='bottom')
# Show <x>
plot([dx, dx], [0., maxY], 'k:', linewidth=0.5)
text(dx, maxY*0.35, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
......@@ -163,6 +183,9 @@ plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
# Show h
plot([h, h], [0., 1.], 'k:', linewidth=0.5)
# Show sigma
plot([h/2, h/2], [0., 1.], 'k:', linewidth=0.5)
# Show <x>
plot([dx, dx], [0., 1.], 'k:', linewidth=0.5)
......@@ -272,13 +295,16 @@ maxY = d_Gaussian(h/2, h)
# Show h
plot([h, h], [2*maxY, 0.1], 'k:', linewidth=0.5)
# Show sigma
plot([h/2, h/2], [2*maxY, 0.1], 'k:', linewidth=0.5)
# Show <x>
plot([dx, dx], [2*maxY, 0.1], 'k:', linewidth=0.5)
xlim(0., 2.5*h)
gca().xaxis.set_ticklabels([])
ylim(1.2*maxY, -0.1*maxY)
ylim(1.1*maxY, -0.1*maxY)
xlabel("$r$", labelpad=0)
ylabel("$\\partial W(r,h)/\\partial r$", labelpad=0.5)
legend(loc="lower right")
......@@ -288,12 +314,14 @@ legend(loc="lower right")
subplot(212)
maxY = d2_Gaussian(h,h)
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)
text(h, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
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, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
text(dx, -3.*maxY, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
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, d2W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$")
......@@ -304,7 +332,7 @@ plot(xx, d2W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$")
plot(xx, d2W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
xlim(0., 2.5*h)
ylim(-3.2*maxY, 1.4*maxY)
ylim(-3.2*maxY, 1.3*maxY)
xlabel("$r$", labelpad=0)
ylabel("$\\partial^2 W(r,h)/\\partial r^2$", labelpad=0.5)
......
#!/bin/bash
python kernels.py
pdflatex kernel_definitions.tex
pdflatex kernel_definitions.tex
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