diff --git a/advection/advection-higherorder.tex b/advection/advection-higherorder.tex index 1d76aa7..3c55246 100644 --- a/advection/advection-higherorder.tex +++ b/advection/advection-higherorder.tex @@ -4,8 +4,8 @@ \section{Advection and the finite-volume method} \label{ch:adv:sndorder} -In these notes, we will typically use a {\em finite-volume} discretization. Here we -explore this method for the +In these notes, we will typically use a {\em finite-volume} discretization. Here we +explore this method for the advection equation. First we rewrite the advection equation in {\em conservation form}: \begin{equation} @@ -13,7 +13,7 @@ \section{Advection and the finite-volume method} \label{eq:advect-cons} \end{equation} where $f(a) = ua$ is the flux of the quantity $a$. In conservation form, -the time derivative of a quantity is related to the divergence of +the time derivative of a quantity is related to the divergence of its flux. % figure created with figures/advection/fv-ghost.py @@ -36,9 +36,9 @@ \section{Advection and the finite-volume method} from ${x_{i-\myhalf}}$ to ${x_{i+\myhalf}}$, normalizing by the zone width, $\Delta x$: \begin{align} -\frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} a_t \, dx &= +\frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} a_t \, dx &= - \frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \frac{\partial}{\partial x} f(a) \, dx \\ -\frac{\partial}{\partial t} a_i &= +\frac{\partial}{\partial t} a_i &= - \frac{1}{\Delta x} \left \{ \left [f(a) \right ]_{i+\myhalf} - \left [f(a) \right ]_{i-\myhalf} \right \} \label{adv:eq:fvadv} \end{align} This is an evolution equation for the zone-average of $a$, and shows @@ -209,7 +209,7 @@ \section{Second-order predictor-corrector scheme} {\label{fig:fvadvect} Second-order finite volume advection showing the result of advecting a tophat profile through five periods with both unlimited and limited slopes. This calculation used 64 zones and -$\cfl=0.7$. \\ +$\cfl=0.7$. \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/advection.py}{advection.py}}} \end{figure} @@ -235,9 +235,9 @@ \subsection{Limiting} \frac{a_i - a_{i-1}}{\Delta x}, \frac{a_{i+1} - a_i}{\Delta x} \right ) \end{equation} instead of Eq.~\ref{eq:slopecentered}. -with +with \begin{equation} -\mathtt{minmod}(a,b) = \left \{ +\mathtt{minmod}(a,b) = \left \{ \begin{array}{ll} a & \mathit{if~} |a| < |b| \mathrm{~and~} a\cdot b > 0 \\ b & \mathit{if~} |b| < |a| \mathrm{~and~} a\cdot b > 0 \\ @@ -277,7 +277,7 @@ \subsection{Limiting} \includegraphics[width=0.75\linewidth]{rea-nolimit-start_003} \\ \includegraphics[width=0.75\linewidth]{rea-nolimit-start_004} \\ \includegraphics[width=0.75\linewidth]{rea-nolimit-start_005} \\ -\includegraphics[width=0.75\linewidth]{rea-nolimit-start_006} +\includegraphics[width=0.75\linewidth]{rea-nolimit-start_006} %\includegraphics[width=0.55\linewidth]{rea-nolimit-start_007} \\ %\includegraphics[width=0.55\linewidth]{rea-nolimit-start_008} \\ \caption[The effect of no limiting on initially discontinuous data]{\label{fig:limitingex}Initially discontinuous data evolved for several steps with @@ -310,7 +310,7 @@ \subsection{Limiting} \end{equation} Then the limited difference is \begin{equation} -\left . \frac{\partial a}{\partial x} \right |_i = +\left . \frac{\partial a}{\partial x} \right |_i = \left \{ \begin{array}{ll} \min \left [ \frac{| a_{i+1} - a_{i-1} |}{2 \Delta x}, @@ -324,7 +324,7 @@ \subsection{Limiting} % Note that a slightly different form of this limiter is presented in \cite{leveque:2002}, where all quantities are in a {\tt minmod}, which appears to limit a bit less. -This is second-order accurate for smooth flows. +This is second-order accurate for smooth flows. The main goal of a limiter is to reduce the slope near extrema. Figure~\ref{fig:limit} shows a finite-volume grid with the original @@ -372,18 +372,18 @@ \subsection{Limiting} \subsection{Reconstruct-evolve-average} Another way to think about these methods is as a reconstruction, -evolve, and average (R-E-A) process (see Figure~\ref{fig:rea}). +evolve, and average (R-E-A) process (see Figure~\ref{fig:rea}). We can write the conservative update as: \begin{align} -a_i^{n+1} &= a_i^n + \frac{\Delta t}{\Delta x} +a_i^{n+1} &= a_i^n + \frac{\Delta t}{\Delta x} (u a^{n+\myhalf}_{i-\myhalf} - u a^{n+\myhalf}_{i+\myhalf} ) \\ - &= a_i^n + \cfl (a^{n+\myhalf}_{i-\myhalf} - a^{n+\myhalf}_{i+\myhalf} ) + &= a_i^n + \cfl (a^{n+\myhalf}_{i-\myhalf} - a^{n+\myhalf}_{i+\myhalf} ) \end{align} If we take $u > 0$, then the Riemann problem will always choose the left state, so we can write this as: \begin{equation} -a_i^{n+1} = a_i^n + +a_i^{n+1} = a_i^n + \cfl \biggl [\underbrace{\left (a_{i-1}^n + \frac{1}{2} (1-\cfl) \Delta a_{i-1} \right )}_{a_{i-\myhalf,L}} - \underbrace{\left (a_{i}^n + \frac{1}{2} (1-\cfl) \Delta a_{i} \right )}_{a_{i+\myhalf,L}} \biggr ] \label{eq:rea_orig} @@ -394,20 +394,20 @@ \subsection{Reconstruct-evolve-average} \begin{equation} a_i(x) = a_i + \frac{\Delta a_i}{\Delta x} (x - x_i) \end{equation} -Consider zone $i$. +Consider zone $i$. If we are advecting with a CFL number $\cfl$, then that means that the fraction $\cfl$ of the zone immediately to the left of the $i-\myhalf$ interface will advect into zone $i$ over the timestep. And only the fraction $1-\cfl$ in zone $i$ -immediately to the right of the interface will stay in that zone. This -is indicated by the shaded regions in Figure~\ref{fig:rea}. +immediately to the right of the interface will stay in that zone. This +is indicated by the shaded regions in Figure~\ref{fig:rea}. The average of the quantity $a$ from zone $i-1$ that will advect into -zone $i$ is +zone $i$ is \begin{eqnarray} -\mathcal{I}_< &=& \frac{1}{\cfl \Delta x} +\mathcal{I}_< &=& \frac{1}{\cfl \Delta x} \int_{x_{i-\myhalf} - \cfl\Delta x}^{x_{i-\myhalf}} a_{i-1}(x) dx \\ % - &=& \frac{1}{\cfl \Delta x} + &=& \frac{1}{\cfl \Delta x} \int_{x_{i-\myhalf} - \cfl\Delta x}^{x_{i-\myhalf}} \left [ a_{i-1} + \frac{\Delta a_{i-1}}{\Delta x} (x - x_{i-1} ) \right ] dx \\ &=& a_{i-1} + \frac{1}{2} \Delta a_{i-1} (1-\cfl) @@ -416,16 +416,16 @@ \subsection{Reconstruct-evolve-average} And the average of the quantity $a$ in zone $i$ that will remain in zone $i$ is \begin{eqnarray} -\mathcal{I}_> &=& \frac{1}{(1-\cfl) \Delta x} +\mathcal{I}_> &=& \frac{1}{(1-\cfl) \Delta x} \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl) \Delta x} a_{i}(x) dx \\ % - &=& \frac{1}{(1-\cfl) \Delta x} - \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl)\Delta x} + &=& \frac{1}{(1-\cfl) \Delta x} + \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl)\Delta x} \left [ a_{i} + \frac{\Delta a_{i}}{\Delta x} (x - x_{i} ) \right ] dx \\ &=& a_{i} - \frac{1}{2} \Delta a_{i} \cfl \end{eqnarray} -The final part of the R-E-A procedure is to average the over the +The final part of the R-E-A procedure is to average the over the advected profiles in the new cell. The weighted average of the amount brought in from the left of the interface and that that remains in the cell is @@ -435,10 +435,10 @@ \subsection{Reconstruct-evolve-average} + (1-\cfl) \left [ a^n_i - \frac{1}{2} \Delta a_i \cfl \right ] \\ &= a^n_i + \cfl \left [a^n_{i-1} + \frac{1}{2}(1 - \cfl) \Delta a_{i-1} \right ] - \cfl \left [ a^n_i + \frac{1}{2} (1-\cfl) \Delta a_i \right ] -\end{align} +\end{align} This is identical to Eq.~\ref{eq:rea_orig}. This demonstrates that the R-E-A procedure is equivalent to our reconstruction, prediction of the -interface states, solving the Riemann problem, and doing the +interface states, solving the Riemann problem, and doing the conservative flux update. % this figure sequence is created by figures/advection/rea.py @@ -490,7 +490,7 @@ \subsection{Reconstruct-evolve-average} \begin{figure}[ht] \centering \includegraphics[width=0.8\linewidth]{fv-gaussian-limiters} \\[1em] -\includegraphics[width=0.8\linewidth]{fv-tophat-limiters} +\includegraphics[width=0.8\linewidth]{fv-tophat-limiters} \caption[Effect of different limiters on evolution] {\label{fig:limiter_panel} The effect of different limiters on the evolution of a Gaussian initial profile (top) and a tophat initial @@ -566,7 +566,7 @@ \section{Multi-dimensional advection} \caption[A 2-d grid with zone-centered indexes]{\label{fig:2dgrid} A simple 2-d grid with the zone-centered indexes. The $\times$s mark the left and right interface states at the upper edge of the $i,j$ zone in each - coordinate direction. For a finite-volume discretization, $a_{i,j}$ + coordinate direction. For a finite-volume discretization, $a_{i,j}$ represents the average of $a(x,y)$ over the zone.} \end{figure} @@ -581,21 +581,21 @@ \section{Multi-dimensional advection} define the average of $a$ in a zone by integrating it over the volume: \begin{equation} -a_{i,j} = \frac{1}{\Delta x \Delta y} - \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} +a_{i,j} = \frac{1}{\Delta x \Delta y} + \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} a(x,y,t) \, dx \, dy \end{equation} Integrating Eq.~\ref{eq:advect2d-cons} over $x$ and $y$, we have: \begin{align} -\frac{1}{\Delta x \Delta y} - \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} - \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} a_t \, dx \, dy = +\frac{1}{\Delta x \Delta y} + \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} + \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} a_t \, dx \, dy = &- \frac{1}{\Delta x \Delta y} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} (u a)_x \, dx \, dy \nonumber \\ &- \frac{1}{\Delta x \Delta y} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} - (v a)_y \, dx \, dy + (v a)_y \, dx \, dy \end{align} or using the divergence theorem, \begin{align} @@ -610,7 +610,7 @@ \section{Multi-dimensional advection} Now we integrate over time---the left side of our expression becomes just the difference between the new and old state. \begin{align} - a_{i,j}^{n+1} - a_{i,j}^n = + a_{i,j}^{n+1} - a_{i,j}^n = &- \frac{1}{\Delta x\Delta y} \int_{t^n}^{t^{n+1}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} \left \{ (u a)_{i+\myhalf,j} - (u a)_{i-\myhalf,j} \right \} dy dt \nonumber \\ &- \frac{1}{\Delta x\Delta y} \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} @@ -628,13 +628,13 @@ \section{Multi-dimensional advection} \item y-face \begin{equation} \langle (va)_{i,j+\myhalf}\rangle_{(t)} = \frac{1}{\Delta x \Delta t} - \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx dt + \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx dt \end{equation} \end{itemize} where $\langle . \rangle_{(t)}$ denotes the time-average over the face. For a second-order accurate method in time, we replace the average in time with the flux at the midpoint in time and the average over the face -with the flux at the center of the face: +with the flux at the center of the face: \begin{equation} \langle (ua)_{i+\myhalf,j} \rangle_{(t)} \approx (ua)_{i+\myhalf,j}^{n+\myhalf} \end{equation} @@ -647,7 +647,7 @@ \section{Multi-dimensional advection} For the advection problem, since $u$ and $v$ are constant, we need to find the interface states of $a$ alone. -There are two methods for computing these interface states, +There are two methods for computing these interface states, $a_{i\pm\myhalf,j}^{n+\myhalf}$ on $x$-interfaces and $a_{i,j\pm\myhalf}^{n+\myhalf}$ on $y$-interfaces: dimensionally split and unsplit. Dimensionally split methods are easier to code, since each dimension is operated on independent of the @@ -691,17 +691,17 @@ \subsection{Dimensionally split} case (Eq.~\ref{eq:statel} and \ref{eq:stater}). For example, the $a_{i+\myhalf,j,L}^{n+\myhalf}$ state is: \begin{eqnarray} -a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + - \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + - \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + +a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + + \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + + \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + \mathcal{O}(\Delta x^2) + \mathcal{O}(\Delta t^2) \nonumber \\ - &=& a_{i,j}^n + - \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + - \frac{\Delta t}{2} \left ( + &=& a_{i,j}^n + + \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + + \frac{\Delta t}{2} \left ( - u \left .\frac{\partial a}{\partial x} \right |_{i,j} \right ) + \ldots \nonumber \\ - &=& a_{i,j}^n + - \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) + &=& a_{i,j}^n + + \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) \left .\frac{\partial a}{\partial x} \right |_{i,j} + \ldots \label{eq:statels} \end{eqnarray} @@ -736,18 +736,18 @@ \subsection{Unsplit multi-dimensional advection} The construction of the $a_{i+\myhalf,j,L}^{n+\myhalf}$ interface state appears as \begin{eqnarray} -a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + - \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + - \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + +a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + + \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + + \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + \mathcal{O}(\Delta x^2) + \mathcal{O}(\Delta t^2) \nonumber \\ - &=& a_{i,j}^n + - \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + - \frac{\Delta t}{2} \left ( - - u \left .\frac{\partial a}{\partial x} \right |_{i,j} + &=& a_{i,j}^n + + \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + + \frac{\Delta t}{2} \left ( + - u \left .\frac{\partial a}{\partial x} \right |_{i,j} - v \left .\frac{\partial a}{\partial y} \right |_{i,j} \right ) + \ldots \nonumber \\ - &=& \underbrace{a_{i,j}^n + - \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) + &=& \underbrace{a_{i,j}^n + + \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) \left .\frac{\partial a}{\partial x} \right |_{i,j}}_{\hat{a}_{i+\myhalf,j,L}^{n+\myhalf}} \underbrace{- \frac{\Delta t}{2} v \left .\frac{\partial a}{\partial y} \right |_{i,j}}_{\mathrm{``transverse~flux~difference''}} + \ldots \label{eq:statelu} @@ -756,7 +756,7 @@ \subsection{Unsplit multi-dimensional advection} explicitly appearance of the ``transverse flux difference'' in the unsplit interface state. We rewrite this as: \begin{equation} -a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} +a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} - \frac{\Delta t}{2} v \left .\frac{\partial a}{\partial y} \right |_{i,j} \end{equation} Here, the $\hat{a}_{i+\myhalf,j,L}^{n+\myhalf}$ term is just the normal @@ -768,7 +768,7 @@ \subsection{Unsplit multi-dimensional advection} compute these one-dimensional $\hat{a}$'s at the left and right every interface in both coordinate directions. Note that these states are still one-dimensional, since we have not used any information from the -transverse direction in their computation. +transverse direction in their computation. \item Solve a Riemann problem at each of these interfaces: \begin{eqnarray} @@ -778,13 +778,13 @@ \subsection{Unsplit multi-dimensional advection} \hat{a}_{i,j+\myhalf,R}^{n+\myhalf}) \\ \end{eqnarray} These states are given the `$^T$' superscript since they are the states -that are used in computing the transverse flux difference. +that are used in computing the transverse flux difference. -\item Correct the +\item Correct the previously computed normal interface states (the $\hat{a}$'s) with the transverse flux difference: \begin{equation} -a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} +a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} - \frac{\Delta t}{2} v \frac{a^T_{i,j+\myhalf} - a^T_{i,j-\myhalf}}{\Delta y} \end{equation} Notice that the @@ -793,11 +793,11 @@ \subsection{Unsplit multi-dimensional advection} right state, it would be those that are transverse, but to the right of the interface: \begin{equation} -a_{i+\myhalf,j,R}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,R}^{n+\myhalf} +a_{i+\myhalf,j,R}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,R}^{n+\myhalf} - \frac{\Delta t}{2} v \frac{a^T_{i+1,j+\myhalf} - a^T_{i+1,j-\myhalf}}{\Delta y} \end{equation} \end{itemize} -A similar procedure happens at the $y$-interfaces. +A similar procedure happens at the $y$-interfaces. Figure~\ref{fig:unsplitstates} illustrates the steps involved in the construction of the $a_{i+\myhalf,j,L}^{n+\myhalf}$ state. \MarginPar{more panels illustrating the sequence?} @@ -821,7 +821,7 @@ \subsection{Unsplit multi-dimensional advection} Once all of the full states (normal prediction $+$ transverse flux difference) are computed to the left and right of all the interfaces -($x$ and $y$), we solve another Riemann problem to find the final +($x$ and $y$), we solve another Riemann problem to find the final state on each interface. \begin{equation} a_{i+\myhalf,j}^{n+\myhalf} = \mathcal{R}(a_{i+\myhalf,j,L}^{n+\myhalf}, @@ -829,7 +829,7 @@ \subsection{Unsplit multi-dimensional advection} \end{equation} The final conservative update is then done via Eq.~\ref{eq:update2du}. -See \cite{colella:1990} for more details on this unsplit method, +See \cite{colella:1990} for more details on this unsplit method, and \cite{saltzman:1994} for details of the extension to 3-d. Figures~\ref{fig:adv:gaussian-2d} and \ref{fig:adv:tophat-2d} show the @@ -916,12 +916,12 @@ \subsection{Method-of-lines in multi-dimensions} \item x-face: \begin{equation} \langle (ua)_{i+\myhalf,j} \rangle = \frac{1}{\Delta y} - \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} (ua)_{i+\myhalf,j}\, dy + \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} (ua)_{i+\myhalf,j}\, dy \end{equation} \item y-face \begin{equation} \langle (va)_{i,j+\myhalf} \rangle = \frac{1}{\Delta x} - \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx + \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx \end{equation} \end{itemize} These are similar to the expressions above, except there is no @@ -979,7 +979,7 @@ \subsection{Method-of-lines in multi-dimensions} \begin{equation} a^n_{i,j} = A^n e^{Ii\theta}e^{Ij\phi} \end{equation} -Next, to simplify things, we'll assume that $u = v$, and $\Delta x = \Delta y$, +Next, to simplify things, we'll assume that $u = v$, and $\Delta x = \Delta y$, then $\cfl = u \Delta t/ \Delta x = v \Delta t / \Delta y$. Inserting this into our difference equation gives: \begin{equation} @@ -1004,7 +1004,7 @@ \subsection{Method-of-lines in multi-dimensions} method in \cite{mccorquodalecolella} allows for $\Delta t \sum_{i=1}^d (|\Ub \cdot \mathbf{e}_d|)/\Delta x_d \lesssim 1.4$. Total variation diminishing Runge-Kutta methods are popular for this application~\cite{gottliebshu:1996}. -% run as: +% run as: % ./pyro.py advection_rk tophat inputs.tophat % ./plot.py -W 7.7 -H 6 -o tophat_rk4_cfl08.pdf advection_rk tophat_0040 % @@ -1019,7 +1019,7 @@ \subsection{Method-of-lines in multi-dimensions} the method-of-lines approach with 4th-order Runge-Kutta time integration. We use $u = v = 1$ for one period on a $32\times 32$ grid. On the left is $\cfl = 0.8$---this is clearly unstable. On the right is $\cfl = 0.4$, - which is stable. + which is stable. This was run as {\tt ./pyro.py advection\_rk tophat inputs.tophat}} \end{figure} @@ -1289,11 +1289,11 @@ \subsubsection{Implementation issues with WENO schemes} fourth order Runge-Kutta, and it is the error from this integrator that is dominating. -To confirm this, figure~\ref{fig:weno-converge-gaussian} uses the eigth order +To confirm this, figure~\ref{fig:weno-converge-sine} uses the eigth order Dormand-Price Runge-Kutta method with adaptive step size control (using the -\texttt{scipy.integrate.ode} routine). Here we see high order convergence for -all $r$, and in each case the convergence rate is $2 r - 2$ for \textbf{no -reason I can understand}. +\texttt{scipy.integrate.ode} routine). In this plot a sine wave is advected five +time around a periodic domain. Here we see high order convergence for +all $r$, and in each case the convergence rate is $2 r - 1$ as expected. \begin{figure}[t] \centering @@ -1308,9 +1308,9 @@ \subsubsection{Implementation issues with WENO schemes} \begin{figure}[t] \centering % figure generated by hydro_examples/advection/weno.py -\includegraphics[width=0.8\linewidth]{weno-converge-gaussian} +\includegraphics[width=0.8\linewidth]{weno-converge-sine} \caption[Very high order WENO convergence rates for linear advection] -{\label{fig:weno-converge-gaussian} WENO solutions for advecting a Gaussian one periods, using four different orders. An eighth order Dormand-Price Runge-Kutta method is used for time evolution. This minimizes the time integrator error and we see convergence at order $2 r - 2$ for all schemes, although the time integrator error eventually shows in the highest order case. \\ +{\label{fig:weno-converge-sine} WENO solutions for advecting a sine wave five periods, using three different orders. An eighth order Dormand-Price Runge-Kutta method is used for time evolution. This minimizes the time integrator error and we see convergence at order $2 r - 1$ for all schemes, although the highest order case shows limits (either from round off error, or the tolerance of the time integrator). \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/weno.py}{weno.py}}} \end{figure} % @@ -1341,18 +1341,18 @@ \section{Going further} \end{equation} The solution procedure is largely the same as described above. We write: \begin{align} -\phi_{i+\myhalf,j,L}^{n+\myhalf} &= \phi_{i,j}^n + +\phi_{i+\myhalf,j,L}^{n+\myhalf} &= \phi_{i,j}^n + \frac{\Delta x}{2} \frac{\partial \phi}{\partial x} + \frac{\Delta t}{2} \frac{\partial \phi}{\partial t} + \ldots \nonumber \\ % - &= \phi_{i,j}^n + + &= \phi_{i,j}^n + \frac{\Delta x}{2} \frac{\partial \phi}{\partial x} + \frac{\Delta t}{2} \left [ -(\phi u)_x -(\phi v)_y \right ]_{i,j} \nonumber\\ % - &= \underbrace{\phi_{i,j}^n + + &= \underbrace{\phi_{i,j}^n + \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u_{i,j} \right ) \frac{\partial \phi}{\partial x}}_{\hat{\phi}_{i+\myhalf,j,L}^{n+\myhalf}} - - \frac{\Delta t}{2} \left [\phi u_x \right ]_{i,j} + - \frac{\Delta t}{2} \left [\phi u_x \right ]_{i,j} - \frac{\Delta t}{2} \left [ (\phi v)_y \right ]_{i,j} \end{align} and upwinding is used to resolve the Riemann problem for both the @@ -1361,7 +1361,7 @@ \section{Going further} according to the velocity field in much the fashion shown here, and is described in \S~\ref{sec:lm:density}. - + For compressible hydrodynamics, we often have density-weighted quantities that we advect. This extension is described in \S~\ref{sec:euler:further}. @@ -1370,7 +1370,7 @@ \section{Going further} \section{\pyro\ experimentation} -To gain some experiences with these ideas, we can use the advection +To gain some experiences with these ideas, we can use the advection solver in \pyro\ (see Appendix~\ref{app:pyro} to get started). The \pyro\ advection solver implements the second-order unsplit advection algorithm described in the previous sections. To run @@ -1402,4 +1402,3 @@ \section{\pyro\ experimentation} and compare the results to the unsplit version. Pick a suitable norm and measure the convergence of the methods.} \end{exercise} - diff --git a/higher-order/weno-converge-sine-rk4.pdf b/higher-order/weno-converge-sine-rk4.pdf new file mode 100644 index 0000000..2811cef Binary files /dev/null and b/higher-order/weno-converge-sine-rk4.pdf differ diff --git a/higher-order/weno-converge-sine.pdf b/higher-order/weno-converge-sine.pdf new file mode 100644 index 0000000..d25f8e4 Binary files /dev/null and b/higher-order/weno-converge-sine.pdf differ