Skip to content
Snippets Groups Projects
Verified Commit b1f02225 authored by Florian Atteneder's avatar Florian Atteneder
Browse files

start notes on upwinding flux for weak boundary conditions

parent 7f66c8ad
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 428 additions and 0 deletions
tex/upwind-flux.tex 0 → 100644
+ 428
− 0
View file @ b1f02225
\documentclass[
10pt,
parskip=half
]{scrartcl}
\usepackage[margin=2cm]{geometry}
% basic packages
\usepackage[english]{babel} % language package
\usepackage{amsmath, amssymb, amstext} % math package
\usepackage[T1]{fontenc} % including important text signs and changing font of document
\usepackage[applemac]{inputenc} % input encoding package
% diverse packages for figures, tables, equations
\usepackage{graphicx} % include graphics
\usepackage{datetime} % date and time
\usepackage{physics} % useful math shortcuts
\usepackage{tensor} % enables multiple sub- and superscripts
\usepackage{enumitem} % change item lables for enumerates
\usepackage{xcolor, xspace} % advanced color and whitespace settings
\usepackage{csquotes} % for quotes, required by bilatex
\usepackage[backend=biber,style=ieee,dashed=false,citestyle=numeric-comp]{biblatex}
\usepackage{hyperref} % crossreferencing
\usepackage{empheq}
\usepackage[most]{tcolorbox}
\usepackage{rotating}
\usepackage{bbm}
%%% styled boxed equation environment
% use as
% \begin{empheq}[box=\mybox]{align}
% a^2 + b^2 = c^2
% \end{empheq}
\newtcbox{\mybox}[1][]{%
nobeforeafter, math upper, tcbox raise base,
enhanced, colframe=blue!30!black,
colback=gray!10, boxrule=1pt,
#1}
%%% bibliography
\addbibresource{../zotero.bib}
%%% hyperlink styling
\hypersetup{
draft=false,
colorlinks,
linkcolor={blue!50!black},
citecolor={green!50!black},
urlcolor={purple!80!black}
}
%%% document data
\title{Upwinding numerical flux for weak boundary conditions of nonlinear conservation laws}
\author{Florian Atteneder}
%%% Custom commands
\newcommand{\lie}{\mathcal{L}}
\newcommand{\coloredText}[2]{{\color{#1}{#2}}\xspace}
\newcommand{\todo}[1]{\coloredText{green!70!black}{[TODO: #1]}}
\newcommand{\wip}[1]{\coloredText{red!50!black}{[WIP: #1]}}
\newcommand{\question}[1]{\coloredText{purple}{[Q: #1]}}
\newcommand{\detail}[1]{\coloredText{cyan!70!black}{#1}}
\newcommand{\pdvt}{\partial_{t}}
\newcommand{\pdvx}{\partial_{x}}
\newcommand{\pdvr}{\partial_{r}}
\newcommand{\nonr}{\nonumber}
\newcommand{\sqrtg}{\sqrt{-g}}
\newcommand{\sqrtgamma}{\sqrt{\gamma}}
\newcommand{\sqrtsigma}{\sqrt{\sigma}}
\newcommand{\rmmin}{\textrm{min}}
\newcommand{\constoprim}{\textsc{cons2prim}}
\newcommand{\etal}{\textit{et al.}}
% use as \eg{}, \ie{}
\newcommand{\eg}{e.\,g.\@}
\newcommand{\ie}{i.\,e.\@}
\newcommand{\vs}{vs.\@}
% notation
\newcommand{\lmfp}{{l_\textrm{mfp}}}
\newcommand{\EE}{{\mathcal{E}}}
\newcommand{\FF}{{\mathcal{F}}}
\newcommand{\PP}{{\mathcal{P}}}
\newcommand{\NN}{{\mathcal{N}}}
\newcommand{\OO}{{\mathcal{O}}}
\newcommand{\MM}{{\mathcal{M}}}
\newcommand{\CC}{{\mathcal{C}}}
\newcommand{\SSS}{{\mathcal{S}}}
\newcommand{\eq}{\textrm{eq}}
\newcommand{\ideal}{\textrm{ideal}}
% \newcommand{\zz}[1]{{\overset{0}{#1}}}
% \newcommand{\oo}[1]{{\overset{1}{#1}}}
\newcommand{\zz}[1]{{{#1}^{(0)}}}
\newcommand{\oo}[1]{{{#1}^{(1)}}}
% state variables
\newcommand{\bD}{\bar{D}}
\newcommand{\bS}{\bar{S}}
\newcommand{\btau}{\bar{\tau}}
\newcommand{\hGa}{{\hat{\Gamma}}}
\newcommand{\bGa}{{\bar{\Gamma}}}
\newcommand{\Kn}{{\text{Kn}}}
\begin{document}
\maketitle
\tableofcontents
\section{Motivation}
\cite[TODO section nr]{koprivaImplementingSpectralMethods2009}
provides a construction of an upwinding numerical flux for
in order to impose boundary conditions on in- and outflow
boundaries.
The example is based on the linear wave equation and this
limits the obtained result to linear systems of equations.
Naively applying the result to nonlinear systems of conservation
laws yields a numerical flux function which is not consistent
with the initial flux of the PDE.
In this notes we extend the discussion from \cite{koprivaImplementingSpectralMethods2009}
to nonlinear systems of conservation laws and derive a
numerical flux function that is consistent with the PDE flux.
\section{Basics}
\todo{
Notation: Clarify that we only use vector arrow notation for the
numerical discretization.
}
Consider a (possibly nonlinear) hyperbolic system of conservation laws,
\begin{align}
\pdvt q + \partial_i f^i(q) = 0 \, ,
% \pdvt q + A^i \partial_i q = 0 \, ,
\label{eq:conservation-law}
\end{align}
We assume appropriate initial conditions and boundary conditions such that the problem
is well-posed.
When discussing boundary conditions it is useful to introduce
to rewrite the above wrt the unit normal vector $n^i$ to the
boundary under consideration.
This allows to distinguish flux components that are flowing normal and
tangential to a boundary.
In particular, the equation takes the form
\begin{align}
0 &= \pdvt q + n_j n^i \partial_i f^j(q)
+ (\tensor{\delta}{_j_i} - n_j n^i) \partial_i f^j(q)
\\
&= \pdvt q + \partial_n f^n(q) + \partial_m f^m(q)
\\
&= \pdvt q + A^n \partial_n q + A^m \partial_m q \, ,
\end{align}
where we introduce the normal and tangential derivatives
$\partial_n f^n = n_j n^i \partial_i f^j, \partial_m f^m = (\tensor{\delta}{_j^i} - n_j n^i)\partial_i f^j$, respectively,
as well as the associated normal and tangential Jacobians
$A^n = n_j n^i \partial_i f^j/\partial q$ and
$A^m = (\tensor{\delta}{_j^i} - n_j n^i) \partial_i f^j/\partial q$,
respectively.
$A^n$ is also refered to as the principle symbol.
\todo{Cite David's well-posedness notes.}
The present notes ignore the tangential part and any potential problems associated
with characteristic modes at the boundary \todo{was characteristic modes the right
term for those speeds that vanish?}.
Let $\psi$ denote a smooth test function.
The weak formulation of \eqref{eq:conservation-law} is given by
\begin{align}
\braket{\pdvt q}{\psi}_D - \braket{f^i(q)}{\partial_i \psi}_D + \eval{\qty[n_i f^{i,\ast} \psi]}_{\partial D} = 0 \, ,
\end{align}
where
\begin{align}
\braket{a}{b}_D = \int_D a(x) b(x)~\dd^3 x
\end{align}
is an inner product.
The weak formulation introduced the numerical flux functions $f^{i,\ast}$
that is associated with the physical flux $f^i$.
They are used to impose boundary conditions between next-neighbor subdomains
in a boundary conforming triangulation of the computational domain.
For practical purposes the choice of $f^{i,\ast}$ is crucial
for the stability properties of the resulting scheme.
To be more specific, consider a boundary that separates two subdomains with unit normal $n^-_i$
that points from the inner domain to the outer domain.
The state variable evaluated at the boundary in the inner and outer domain are refered
to as $q^-$ and $q^+$. The inwards pointing unit normal is $n^+_i = - n^-_i$.
The numerical flux function should then obey \textit{at least} the following two properties
\cite{hesthavenNodalDiscontinuousGalerkin2008}:
\begin{enumerate}
\item Conservation
\begin{align}
(n^-_i f^{i,\ast})(q^-, q^+, n^-_i) =
-(n^+_i f^{i,\ast})(q^+, q^-, n^+_i) \, .
\end{align}
\item Consistency \begin{align}
(n^-_i f^{i,\ast})(q, q, n^-_i) =
(n^-_i f^i)(q) \, .
\end{align}
\end{enumerate}
In the following we restrict ourselves to one spatial dimension, hence, the unit normal
is $n^-_i = \pm \delta_{i,x}$ and we can drop the $i$ indices from the normal vectors and fluxes.
However, all results straightforwardly extend to higher dimensions.
\section{Linear case}
In this section we superficially repeat the discussion from \cite{koprivaImplementingSpectralMethods2009}.
That is, we work with the general form \eqref{eq:conservation-law} with
$A^n = const.$, which holds for a linear system of conservation laws.
Note that in this case $n f(q) = A^n q$.
The idea of the upwinding numerical flux construction is to decompose the flow of the
conservation law into in- and outgoing waves wrt the boundary's unit normal
and then discard the ingoing part for an outflow boundary condition,
or discard the outgoing part for an inflow boundary condition.
This analysis makes use of the characteristic decomposition of the PDE.
To this end, one computes the eigenvalues and eigenvectors of the princple symbol $A$,
such that
\begin{align}
(T^n)^{-1} A^n T^n = \Lambda^n \, ,
\label{eq:diagonal}
\end{align}
where $T^n$ is now an orthogonal matrix that contains the right eigenvectors of $A^n$ as the columns,
and $\Lambda^n$ is a diagonal matrix of eigenvalues.
We now ignore the superscript $n$ for a moment.
Because $A = const.$, we also have $T = const.$ and we can introduce the characteristic
variables $v = T^{-1}$. Multiplying \eqref{eq:conservation-law} by $T^{-1}$ from the left
and inserting an identity $T T^{-1}$ we obtain the evolution equation for the characteristic
variables,
\begin{align}
0 &= T^{-1} \pdvt q + T^{-1} A \pdvx q
\\
&=
T^{-1} \pdvt q + T^{-1} A T T^{-1} \pdvx q
\\
&=
\pdvt T^{-1} q + T^{-1} A T \pdvx T^{-1} q
\\
&=
\pdvt v + \Lambda \pdvx v \, .
\end{align}
This is now a system of advection equations for $v$, which is due to $\Lambda$ being
constant and of diagonal shape. The analytical solution to this equation can be written
done immediately.
This allows one now to identify the eigenvalues as the corresponding speeds with which the
components of $v$ are advected by the conservation law.
The construction of the upwinding flux now starts from \eqref{eq:diagonal}.
For a hyperbolic conservation law all eigenvalues are real and so we can group the eigenvalues
wrt their sign. That is, we decompose the principle symbol as
\begin{align}
T \Lambda T^{-1}
= T ( \Lambda^{(+)} + \Lambda^{(0)} + \Lambda^{(-)} ) T^{-1}
= T \Lambda^{(+)} T^{-1} + T \Lambda^{(0)} T^{-1} + T \Lambda^{(-)} T^{-1}
= A^{(+)} + A^{(0)} + A^{(-)} = A \, .
\label{eq:decompose_principle_part}
\end{align}
Because the flux is linear, the flux itself can be decomposed as
\begin{align}
n^- f(q) = n^- f^{(+)}(q) + n^- f^{(0)}(q) + n^- f^{(-)}(q) = A^{(+)} q + A^{(0)} q + A^{(-)} q \, ,
\end{align}
where $\Lambda^{(+)}, \Lambda^{(-)}$ and $\Lambda^{(0)}$ contain now the speeds with which
waves are advected along the direction $n$, against it or not even advected at all, respectively.
Based on this decomposition one can now design a numerical flux that only allows
information to flow into the domain when considering an inflow boundary
condition, and similar for an outflow one.
That is, consider a subdomain $[x_l, x_r]$ where $x_r$ is an inflow boundary
with unit normal $n^- = 1$, \eg{} the flow satisfies $n^- f(q(x_r)) < 0$.
The numerical flux is now chosen such that only characteristic fields
with negative eigenvalues, $\Lambda^{n^-(-)}_{ii} < 0$,
are taken into account to impose the desired boundary condition weakly.
These waves are left moving, because $n^-$ is pointing to the right.
Let $q_{BC}$ be the desired boundary value.
The upwinding numerical flux function is then given by
\begin{align}
n^- f^\ast(q_r, q_{BC}, n^-)
&= n^- f^{(-)}(q_r) + n^- f^{(-)}(q_{BC})
\\
&= n^- f^{(-)}(q_r) + n^+ f^{(+)}(q_{BC})
\\
&= A^{n^-(-)} q_r + A^{n^+(+)} q_{BC} \, .
\end{align}
Here we used the fact the the flow $n^- f^{(-)} = A^{n^-(-)}$ along $n^- = 1$,
which corresponds to the left moving modes,
is equal to the flow $n^+ f^{(+)} = A^{n^+(+)}$ along $n^+ = -1$,
which also corresponds to the left moving modes.
This latter part of the rewriting above appears to be unnecessarily complicated,
however, doing so allows one to generalize the numerical flux function to arbitrary directions.
Another reason why this is often done is that one can now identify the
first term as the data that is coming from the interior of the considered subdomain,
and from which the normal vector of the boundary $n^-$ is outwards pointing.
Similarly, the second term corresponds to data coming from the exterior of the considered subdomain,
and the normal vector $n^+$ to the boundary as seen from the exterior side is now
also outwards pointing wrt the exterior.
The above numerical flux formula can be written in a more general form
using the interior and exterior states $q^-$ and $q^+$ as
\begin{align}
n^- f^\ast(q^-, q^+, n^-)
&= \frac{1}{2} \left( 1 - \frac{n^- f}{| n^- f |} \right) (n^- f^{(-)}(q^-) - n^- f^{(+)}(q^+))
\nonumber\\
&\quad + \frac{1}{2} \left( 1 + \frac{n^- f}{| n^- f |} \right) (n^- f^{(+)}(q^-) - n^- f^{(-)}(q^+)) \, .
\label{eq:num_flux_linear}
\end{align}
The prefactors $\frac{1}{2} \left( 1 \mp \frac{n^- f}{| n^- f |} \right)$ automatically
select the appropriate upwinding term based on whether the flow is incoming, $n^- f < 0$,
or outgoing, $n^- f > 0$.
In practice, these prefactors are omitted and one instead usese explicit if statements
to distinguish between in- and outflow, in part also to avoid potential divison by zero
problems when $n^- f = 0$.
\todo{Give an example. Maybe just the wave equation from \cite{koprivaImplementingSpectralMethods2009}.}
\todo{Give a counter example that illustrates that this construction fails
the consistency property when $f$ is a nonlinear function.}
\section{Nonlinear case}
The previous example illustrated that the numerical flux \eqref{eq:num_flux_linear} is not
suitable for imposing boundary conditions for nonlinear conservation laws,
because the consistencty property is violated.
The reason this fails is that $f$ is nonlinear, and so we no longer can
decompose it in analogy to how we decomposed the principle part \eqref{eq:decompose_principle_part}.
As we demonstrate below, \eqref{eq:num_flux_linear} can be generalized for nonlinear
fluxes by realizing that the upwinding decomposition can only be applied to suitably small
deviations from the average state variable between the boundary states.
To make the above claim precise, we first rewrite the nonlinear conservation law
in terms of an averaged state $q_0$ and a deviation $\delta q$ around it,
\begin{align}
q = q_0 + \delta q \, .
\label{eq:decomposition_state}
\end{align}
To be specific, consider the three subdomains
$[x^{k-1}_l, x^{k-1}_r], [x^{k}_l, x^{k}_r], [x^{k+1}_l, x^{k+1}_r]$,
where $x^{k-1}_r = x^k_l$ and $x^k_r = x^{k+1}_l$.
We assume that the global solution $q = q(x)$ is piecewise continuous and piecewise smooth
over the computational domain, \eg{} it is given by direct sum of local solutions,
\begin{align}
q(x) = \bigoplus_{k} q^k(x) \, ,
\label{eq:global-solution}
\end{align}
where each $q^k = q^k(x)$ is smooth.
Let $q^k_{l,r} = q(x^{k}_{l,r})$ and similarly for the other domains.
Note that in general $q^k_r \neq q^{k+1}_r$ and so on.
The averaged state and its deviation can now be defined in the $k$th subdomain by
\begin{align}
q^k_0(x) &= q^k(x) + \left( \frac{q^{k-1}_r + q^k_l}{2} - q^k_l \right) \alpha^k_l(x)
+ \left( \frac{q^{k+1}_l + q^k_r}{2} - q^k_r \right) \alpha^k_r(x) \, ,
\label{eq:local_sol_average}
\\
\delta q^k(x) &= \left( q^k_l - \frac{q^{k-1}_r + q^k_l}{2} \right) \alpha^k_l(x)
+ \left( q^k_r - \frac{q^{k+1}_l + q^k_r}{2} \right) \alpha^k_r(x) \, ,
\label{eq:local_sol_delta}
\end{align}
and similarly in other subdomains.
The $\alpha^k_l, \alpha^k_r$ shall satisfy
\begin{align}
\alpha^k_l(x_l) &= 1 \, ,
&
\alpha^k_l(x_r) &= 0 \, ,
\\
\alpha^k_r(x_l) &= 0 \, ,
&
\alpha^k_r(x_r) &= 1 \, .
\end{align}
In practice, for a nodal DG method in which the local solution is $q^k \in \mathbb{P}_N$,
\ie{} it as an $N$th order polynomial, the $\alpha^k_{l,r}$ can be taken to be $N$th order
Lagrange polynomials that are centered on the nodes $x^k_{l,r}$.
With this choice, we then also have $q^k_0, \delta q^k \in \mathbb{P}_N$.
Similarly to \eqref{eq:global-solution}, one can construct the global solutions $q_0$
and $\delta q$ from the local ones \eqref{eq:local_sol_average} and \eqref{eq:local_sol_delta},
\begin{align}
q_0(x) &= \bigoplus_{k} q_0^k(x) \, ,
&
\delta q(x) &= \bigoplus_{k} \delta q^k(x) \, .
\end{align}
It is apparent that $q_0$ is now also continuous across subdomains, but still only piecewise smooth.
Furthermore, the jumps in $\delta q$ across the subdomain boundaries like $x^k_r = x^{k+1}_l$
are given by $q^k_r - q^{k+1}_l$.
Using \eqref{eq:decomposition_state} we can rewrite \eqref{eq:conservation-law} as
\begin{align}
\pdvt q_0 + \pdvt \delta q + \pdvx f(q_0 + \delta q) = 0 \, .
\end{align}
We assume now that the jumps $\delta q$ between subdomains are sufficiently small such
that we can Taylor expand the flux around $q = q_0$, \ie{},
\begin{align}
f(q_0 + \delta q) - f(q_0) = f'(q_0) \delta q + \OO(\delta q^2) \, .
\end{align}
It is clear that in the linear case this expansion terminates after the first term on the RHS.
\todo{point out that we see here that we can only really control upwinding
of the deviations $\delta q$, because this is now a linear perturbation around $f(q_0)$.}
With this we arrive at
\begin{align}
\pdvt q_0 + \pdvx \delta f(q_0) + \pdvt \delta q + \pdvx f'(q_0) \delta q = \OO(\delta q^2) \, .
\end{align}
The corresponding weak formulation is given by
\begin{align}
\braket{q_0}{\psi}_D - \braket{f(q_0)}{\pdvx \psi}_D + \eval{\qty[n^- f^\ast(q_0) \psi]}_{\partial D}
+ \braket{\delta q}{\psi}_D - \braket{f'(q_0)\delta q}{\pdvx \psi}_D + \eval{\qty[n^- (f'(q_0) \delta q)^\ast \psi]}_{\partial D} \approx 0 \, .
\end{align}
On first sight it appears that one now needs to design two numerical flux functions,
$n^- f^\ast$ and $n^ (f(q_0)\delta q)^\ast$ in order to establish the next-neighbor coupling
and stability.
However, due to $q_0$ being continuous even across subdomain boundaries,
it is justified to use $n^- f^\ast = n^ f$, because the knowledge of the interior state $q_0^-$
of a boundary specifies the exterior state $q_0^+$, \ie{}, $q^-_0 = q^+_0$.
Consequently, the conservation and consistency property are each automatically satisfied.
It remains to discuss the design of $n^-(f'(q_0)\delta q)^\ast$ in which we now incorporate
the upwinding property.
\printbibliography
\end{document}
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