\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Eighth Mississippi State - UAB Conference on Differential Equations and
Computational Simulations.
{\em Electronic Journal of Differential Equations},
Conf. 19 (2010),  pp. 65--73.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{65}
\title[\hfilneg EJDE-2010/Conf/19/\hfil A hybrid scheme]
{A hybrid semi-primitive shock capturing scheme for conservation laws}

\author[R. K. Dubey\hfil EJDE/Conf/19 \hfilneg]
{Ritesh Kumar Dubey} 

\address{Ritesh Kumar Dubey \newline
  Institut de   Math\'ematiques de Toulouse,
  Universit\'e Paul  Sabatier, 
  118 route de Narbonne,
  F-31062 Toulouse Cedex 9, Fance}
\email{riteshkd@gmail.com}

\thanks{Published September 25, 2010.}
\subjclass[2000]{35L65, 65M06, 65M12}
\keywords{Hyperbolic conservation laws; primitive scheme;
  upwind difference; \hfill\break\indent hybrid schemes}

\begin{abstract}
  A hybrid semi-primitive shock capturing scheme is
  presented for hyperbolic problems. Upwind based construction
  is done using explicit information on the wave propagation direction
  associated with the problem.
  This scheme captures the shock waves at right location but
  shows unphysical sonic expansion shock. This phenomena
  of unphysical expansion shock in the presence of expensive sonic point
  is not surprising and it is common for the wave based upwind schemes.
  A hybrid scheme approach using an iteration criteria based on sonic
  entropy fix is proposed to avoid such expansion shocks.
  Numerical results for  scalar test problems are presented which
  show that proposed scheme captures the shock accurately.
\end{abstract}

\maketitle
\numberwithin{equation}{section}

\section{Introduction}

In this work we consider the hyperbolic differential equations in the
following primitive (non-conservative) form,
\begin{equation}
  \frac{\partial u}{\partial t} +   \alpha(u) \frac{\partial (u)}{\partial x}
  =0, \quad (x,t)\in \mathbb{R}\times \mathbb{R}^{+},\label{eq1a}
\end{equation}
together with the appropriate boundary and initial conditions. In
\eqref{eq1a}, the unknown variable $u \in \mathbb{R}$ can be a
physically conserved variable and in such situation the
characteristic speed $\alpha(u)=\frac{\partial f}{\partial u}$
where, $f=f(u)$ be the physical flux and \eqref{eq1a} can be
written in conservative form as
\begin{equation}
  \frac{\partial u}{\partial t} +   \frac{\partial f(u)}{\partial x}
  =0.\label{eq1b}
\end{equation}
It is well known that even if the initial condition is smooth
enough, various kind of discontinuity e.g, shock wave, contact
discontinuity and expansion fan may arise in the solution of such
hyperbolic problems. Looking at literature to solve hyperbolic
problems numerically one can experience an imbalance in between
the number of available conservative schemes and primitive
(non-conservative) schemes. The main reason for this imbalance
seems to be the possible occurrence of shock waves in the solution
of \eqref{eq1b} and it is shown in \cite{hou} that conservative
form is mandatory for capturing the shock at right location while
the primitive form captures the shock at wrong position and
primitive schemes often converge to wrong week solution. In fact,
most of the numerical schemes developed for hyperbolic problems
are conservative. A state of art detail on various conservative
numerical schemes for solving \eqref{eq1b} can be found in
\cite{laney,leveque1,leveque2002,Toro1,wesseling}. On the other
hand primitive or semi-primitive schemes work well for the smooth
solution, discontinuities such as contact, shear waves and mild
shocks \cite{Toro1} but these methods are simply discarded.
Recently few attempts are  made to construct primitive schemes.
Primitive schemes based on upwinding are given in
\cite{Toro2,karni} and centered schemes are proposed in
\cite{siviglia}.


In this work our main aim is to construct a primitive scheme which
can capture the shock waves at right location. In order to do this
we look for a primitive scheme which changes into conservative
form around discontinuities so that it can capture shock
accurately as suggested in \cite{hou}. The scheme is designed
based on upwinding by using the explicit information on
characteristic speed $\alpha= f'(u)$ of primitive form
\eqref{eq1a}.   Thus constructed semi-primitive scheme captures
the shock wave correctly but in the presence of expensive sonic
point where wave speed changes its sign, produces entropy
violating solution. In order to overcome this, we use the idea of
hybrid schemes and use entropy satisfying  Lax-Friedrichs in the
presence of sonic point. To decide upon iterations a criteria is
given using an entropy fix approach.

\section{Primitive scheme formulation \label{sec2}}

We define semi-primitive (semi-conservative) schemes as a scheme which remains
primitive (non-conservative) for the region of smooth solution and
locally corrected into conservative form around discontinuities. For
the simplicity of presentation, consider the wave equation
\begin{equation}
  \frac{\partial u}{\partial t} +   a\frac{\partial u}{\partial x}
  =0,\quad  0\neq a\in\mathbb{R}.\label{eq3}
\end{equation}
We discretize the domain with fixed space and time step sizes
$\Delta x$, $\Delta t$ respectively. Integrating \eqref{eq3} over
the rectangle
$[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[t^{n},t^{n+1}]$ and
introducing the definitions of the spatial and temporal cell
averages
\begin{equation}
U_{i}^{n} = \frac{1}{\Delta
  x}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}u(x,t^{n})dx,\quad
U_{i \pm \frac{1}{2}}  =  \frac{1}{\Delta
  t}\int_{t^{n}}^{t^{n+1}}u(x_{i\pm\frac{1}{2}},t)dt,\label{eq4}
 \end{equation}
one can obtain a primitive difference scheme of the form
\begin{equation}
U_{i}^{n+1}=U_{i}^{n}- a\lambda\left(
U_{i+\frac{1}{2}}-U_{i-\frac{1}{2}}\right) \label{eq5}
\end{equation}
%
where  $\lambda = \frac{\Delta t}{\Delta x}$,  $U_{i\pm\frac{1}{2}}$
are intermediate states defined on the cell interfaces $x_{i \pm
  \frac{1}{2}}$.
Note that for constant wave speed $a$, the difference scheme
\eqref{eq5} can be written in conservative form,
\begin{equation}
U_{i}^{n+1}=U_{i}^{n}- \lambda\left(
H_{i+\frac{1}{2}}-H_{i-\frac{1}{2}}\right) \label{eq6}
\end{equation}
where $H_{i \pm \frac{1}{2}}  = aU^{n}_{i \pm\frac{1}{2}}$ is the
numerical flux function. It is well-known that the first order
upwind approximation for \eqref{eq3} can be obtained by using the
following intermediate states $U_{i\pm\frac{1}{2}}$ in \eqref{eq5}
\begin{equation}
 U_{i+\frac{1}{2}}=  \begin{cases}
 U_{i},  & a\geq0,\\
 U_{i+1},  & a\leq0.
 \end{cases} \label{eq7}
\end{equation}
Consider the primitive form \eqref{eq1b}  of hyperbolic problem
\eqref{eq1a}
\begin{equation}
  \frac{\partial u}{\partial t} +   \alpha(u)\frac{\partial u}{\partial x}
  =0, \label{eq8}
\end{equation}
On comparison with \eqref{eq3} and utilizing the explicit
information of the characteristic speed $\alpha(u)$ we can
construct a primitive upwind scheme so that it respects physical
hyperbolicity property associated with \eqref{eq1a}. Analogous to
first order accurate primitive difference scheme for linear
equation \eqref{eq3}, we propose a first order primitive upwind
scheme for non-linear equation \eqref{eq8} which can be written as
\begin{equation}
U_{i}^{n+1}=U_{i}^{n}-\lambda \bar{s_{i}}\big\{
{U}_{i+\frac{1}{2}}-{U}_{i-\frac{1}{2}}\big\} \label{eq9}
\end{equation}
%
where ${U}_{i\pm\frac{1}{2}}$ are intermediate states on the cell
interface at $x_{i\pm\frac{1}{2}}$ and defined similar to
\eqref{eq7} as follows,
\[
 U_{i+\frac{1}{2}}=  \begin{cases}
  U_{i},  & a_{i+\frac{1}{2}}\geq0,\\
  U_{i+1},  & a_{i+\frac{1}{2}}\leq0,
  \end{cases},\quad
 a_{i+\frac{1}{2}}=
\begin{cases}
 \frac{F_{i+1}-F_{i}}{U_{i+1} -U_{i}},  & U_{i+1} \neq U_{i},\\
\alpha(U_{i}),  & U_{i+1} = U_{i}.
 \end{cases}
 \]
 where  ${F}_{i} = f({U}_{i})$.


Numerical wave speed  $\bar{s}_{i}$ in \eqref{eq9} is the local
linearized approximation for $\alpha(u)$ in the $i_{th}$ cell. We
define the numerical wave speed using intermediate states in such
a way that the resulting scheme changes into conservative form in
the vicinity of shock as follows,
\begin{equation}
\bar{s}_{i}=  \begin{cases}
\frac{{F}_{i+\frac{1}{2}} - {F}_{i-\frac{1}{2}}}{{U}_{i+\frac{1}{2}} -
  {U}_{i-\frac{1}{2}}},& \text{if }{U}_{i+\frac{1}{2}} \neq
     {U}_{i-\frac{1}{2}},\\
f'_{i}=f'({U}_{i-\frac{1}{2}}), & \text{if }
 {U}_{i+\frac{1}{2}} =
{U}_{i-\frac{1}{2}},
\end{cases} \label{eq10}
\end{equation}

Note that, in the vicinity of shock discontinuity, proposed scheme
gives conservative approximation. For this, let shock
discontinuity be locally in the $i_{th}$ cell then
${U}_{i+\frac{1}{2}} \neq {U}_{i-\frac{1}{2}}$ and \eqref{eq9},
\eqref{eq10} yield a conservative approximation
\begin{equation}
U^{n+1}_{i} = U^{n}_{i} - \lambda(F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}).
\end{equation}
This shows that proposed semi-primitive scheme is
locally corrected by a conservative scheme, hence can capture shocks
accurately \cite{hou}. Presented numerical results also validate this
feature.

\section{Hybrid Scheme}

We have observed that the proposed semi-primitive scheme captures the
shock discontinuity correctly but in presence of expansive sonic point where
wave speed changes its $sign$, produces entropy violating solution
which shows the less dissipative nature of semi-primitive scheme
near expensive sonic point (see Fig. \ref{Fig3} Left). This phenomena is
not new for the wave propagation based upwind schemes and it is well
known in the literature that one need to use an entropy fix
in the solution region where
wave speed changes its sign \cite{hyman,kermani} . In order to apply proposed method for the
problem with
expensive sonic point we use hybrid schemes approach. The idea is to
use entropy satisfying Lax-Friedrichs centered scheme in the region of
sonic point. We propose an iteration approach based of the entropy
fix formula for residual distribution scheme given in \cite{sermeus}
as follows,
\begin{equation}
  U^{n+1}_{i}= \begin{cases}
        U_{i} -\lambda \bar{\alpha}_{i}\big(
        {U}_{i+\frac{1}{2}}-{U}_{i-\frac{1}{2}}\big),  &
        |\bar{s}_{i}| \geq\delta_{i}, \\
        \frac{(U_{i+1} +U_{i-1})}{2} -\frac{\lambda}{2}(f(U_{i+1}) -
        f(U_{i-1})),  & |\bar{s}_{i}| <\delta_{i}.
    \end{cases} \label{eq11}
\end{equation}
where
\begin{equation}
\delta_{i} = \gamma\max(0, (\tilde{a}_{i+\frac{1}{2}}
-\tilde{a}_{i-\frac{1}{2}})),\quad
0\leq \gamma\leq 1,\quad
\tilde{a}_{i} = f'(U_{i}).\label{eq12}
\end{equation}

\noindent{\bf Remarks}
(1) One can also use other entropy fix formulas
 for numerical computation \cite{hyman,kermani}.

(2) For non-linear problem \eqref{eq1b}, the proposed scheme is
  different from the classical conservative upwind scheme though for
  the linear problem \eqref{eq3} it can be written as classical upwind
  scheme. In fact, in smooth solution region the proposed scheme
  \eqref{eq9} may give primitive approximation to \eqref{eq1b}.

\section{Numerical Results \label{numerical}}

\noindent{\bf Remarks} (1) In all our numerical results, we denote
the results obtained using proposed semi-primitive scheme (without
entropy iteration) by {\it SP} scheme and with entropy iteration
that is hybrid scheme by {\it HSP} scheme. For all numerical
simulation of {\it HSP}, a fixed value $\gamma=\frac{1}{2}$ is
taken in \eqref{eq12}.

(2)The proposed {\it SP} scheme works very well for most of the
  case but fails to capture centered expansion fan. Hence for safer side
  {\it HSP} scheme should always be used. For some tests results are given
  only for {\it SP/HSP} as similar results are obtained with proposed {\it HSP/SP}
  scheme.

(3) One can easily show that proposed scheme is stable for linear
  problem \eqref{eq3} under the CFL condition $|a\lambda| \leq 1$ this
  is why we use CFL like stability condition $|\lambda
  \tilde{\alpha}|\leq 1$ for non-linear problems.

\subsection{Convex Flux function: Inviscid Burgers Equation}

We consider the inviscid Burgers equation
  \begin{equation}
 \frac{\partial u}{\partial t} + \frac{\partial }{\partial
   x}\Bigl(\frac{u^2}{2}\Bigr) =0, \quad  t>0, \label{n3}
\end{equation}
with periodic boundary conditions. It is well known that the solutions
of inviscid Burgers equation may contain shocks \cite{laney}. We
present three different test cases to show the accuracy and shock
capturing at correct location using the proposed scheme.

\subsubsection{Case 1: Pre-/Post-Shock:}

In this example, we take smooth sinusoidal initial condition,
\begin{equation}
u(x, 0)= \sin(x),\quad x\in[0,\,2\pi]. \label{n4}
\end{equation}
In this test case the final solution beginning from smooth initial conditions
contains a shock discontinuity after a breaking time $t=1.0$.

In Fig. \ref{Fig1}, approximate solutions obtained by proposed
semi-primitive scheme  {\it SP } are presented at pre-shock time
$t=0.25,\;0.75$ when the solution remains smooth and at post shock
time $t=1.5$ when the shock is completely developed. Graph shows that
{\it SP}  gives good approximation for the smooth solution and
captures the shock at right location.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1} %SinCis0p6Nis40.eps
\end{center}
\caption{Computed solution by semi-primitive (SP) scheme for various
    pre and post shock time with  $CFL=0.6$ }
 \label{Fig1}
\end{figure}


\subsubsection{Case 2: Moving Shock:}
Here we consider the Burgers equation \eqref{n3}  with the
 initial condition
\begin{equation}
  u(x,0) =   \begin{cases}
    2,  & \text{for } x \leq 0, \\
    -1, & \text{for } x > 0.
  \end{cases}
\end{equation}
%
This test case has no sonic point in its solution. The initial
shock discontinuity at $x=0$ propagates without losing its initial
shape with time. In Figure \ref{Fig2} we compare the solution obtained
with {\it SP} at different time level which shows that the scheme is
able to capture the shock at right location within $2$ grid
point. Similar results are obtained with {\it HSP}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2} % shock.eps
\end{center}
\caption{ Shock capturing feature of proposed
   SP scheme for $C=0.5,\;\Delta x =0.1$ at various time level }
 \label{Fig2}
\end{figure}


\subsubsection{Case 3: Sonic expansion fan in presence of sonic point:}

Consider the Burgers equation \eqref{n3} with the initial
condition
\begin{equation}
  u(x,0) =   \begin{cases}
    -1, &\text{for } x < 0, \\
    2,  &\text{for } x \geq 0.
  \end{cases}
\end{equation}
The exact solution is an expansion fan emanating from the sonic point
at $x=0.$ In Figure \ref{Fig3}, numerical solutions are compared with
exact reference solution at time $t=0.4$ on a uniform grid with $200$ points.
The solution computed by proposed {\it SP} scheme shows an expansion shock
while the entropy fix based {\it HSP} scheme removes the expansion
shock very satisfactorily.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3a} % expansionshock.eps
\includegraphics[width=0.7\textwidth]{fig3b} % expansionfan.eps
\end{center}
\caption{Sonic expansion: semi-primitive scheme
    (A) and hybrid semi-primitive scheme (B) for data $C=0.8,  \Delta x
    =0.02$ at time $t=0.4$}
 \label{Fig3}
\end{figure}

\subsection{Convex-Concave flux: Buckley Leverett Equation}

A more demanding test example for the scalar one-dimensional
problem is the Buckley-Leverett equation. This equation models the
two phase flows that arise oil-recovery problems and physically
represents a mixture of oil and water through the porous medium.
In conservative form \eqref{eq1b} the flux function for this
problem is given by a convex-concave (S-shaped) function
\begin{equation}
f(u) = \frac{u^2}{u^2 + \nu(1-u)^2}. \label{buckley}
\end{equation}
Here $\nu$ is viscosity ratio and $u$ represents the saturation of water and
lies between $0$ and $1$.

\subsubsection{One moving shock \cite{wesseling}}

Consider the flux \eqref{buckley} with $\nu =\frac{1}{2}$ and
initial condition
\begin{equation}
u(x,0)= \begin{cases}
  1, & x<0,\\
  0,  & x>0.
\end{cases}
\end{equation}
The solution involves one single moving shock followed by a rarefaction wave.
Numerical  results using proposed hybrid semi-primitive scheme ({\it HSP})
are presented  in Fig. \ref{Fig4}. {\it HSP} sharply captures the moving
shock at right location and rarefaction wave accurately.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig4} % HSC1shock.eps
\end{center}
\caption{Solution profile for data $C=0.4$,  $\Delta x =0.02$ at
time $t=0.5,\,1.0,\, 1.5$}
\label{Fig4}
\end{figure}

\subsubsection{Two moving shock \cite{flaherty}}

Consider the flux \eqref{buckley} with $\nu =\frac{1}{4}$ and
subject to initial condition
\begin{equation}
u(x,\,0)= \begin{cases}
  1, & -0.5\leq x\leq 0,\\
  0,  & elsewhere.
\end{cases}
\end{equation}
The solution involves two moving shocks, each followed by an
rarefaction wave. Numerical result is given in Fig. \ref{Fig5} which
shows that the {\it HSP} sharply captures both the fast and slow shocks. The
rarefaction waves are also accurately approximated by {\it HSP}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig5} % HSC2shock.eps
\end{center}
\caption{Solution profile using $C=0.4$,
  $\Delta x =0.01$ at time $t=0.5$}
\label{Fig5}
\end{figure}

\subsection{Flux Sine problem: Moving and Stationary Shock}

We consider third scalar Riemann problem with
non-convex sinusoidal flux function which is introduced by Leveque
in \cite{leveque2002}. It is given by
\begin{equation}
\frac{\partial u}{\partial t} +\frac{\partial (\sin(u))}{\partial
x} =0 \label{fluxsine}
\end{equation}
with the initial condition
\begin{equation}
u(x,0)= \begin{cases}
  \frac{\pi}{4}, & x<0,\\
  \frac{15\pi}{4},  & x>0.
\end{cases}
\end{equation}
The solution consists one stationary shock at the origin, one
moving shock followed by an expansion wave and one moving
expansion wave. In Fig \ref{Fig6} solution obtained by {\it HSP}
is shown for time $t=1.0$. Figure \ref{Fig6} shows that {\it HSP}
capture the stationary shock as well moving shocks nicely with a
nice resolution for the two expansion waves.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig6} % leveque.eps
\end{center}
 \caption{Solution profile using $C=0.6$,
     $\Delta x =0.02$ at time $t=1.0$}
\label{Fig6}
\end{figure}

\subsection*{Conclusion and future work}

In this work, we define the semi-primitive scheme by considering
the primitive form of conservation law.  Primitive scheme is
constructed in such a way that it becomes conservative in the vicinity
of discontinuity. Proposed scheme captures the shock with high
accuracy but yields rarefaction shocks in the presence of sonic
point. A hybrid scheme approach is proposed based on entropy fix
iteration. Presented numerical examples show that hybrid scheme captures shock
at right location and resolves the expansion fan. Extension for
the hyperbolic system and high order accurate hybrid semi-primitive
scheme is under investigation.

\begin{thebibliography}{00}

\bibitem{flaherty} Xin, J.; Flaherty, J. E.;
 Viscous stabilization of discontinuous
 Galerkin solutions of hyperbolic conservation laws, App. Num. Math. 56 (2006), 444-458.

\bibitem{hou} Hou, T.Y.; LeFloch, P.G.;
 Why non-conservative schemes converge to wrong solutions:
error analysis, Mathematics of Computation, 62(1994), 497-530.

\bibitem{hyman} Harten, A.; Hyman, J. M.;
 Self-Adjusting Grid Methods   for One-Dimensional Hyperbolic
Conservation Laws, J. of  Comput. Physics, 50(1983), 235-269.

\bibitem{karni} Karni, S.;
 Multicomponent flow calculations using a
  consistent primitive algorithm. J. Comp. Phys. 112( 1994), 31-43.

\bibitem{kermani} Kermani, M. J.; Plett, E. G.;
Modified Entropy  Correction Formula for the Roe Scheme, AIAA 2001- 0083.

\bibitem{laney} Laney, C. B.;
 Computational Gas-dynamics (Ist edn). Cambridge
  University press: Cambridge, 1998.

\bibitem{leveque1} Leveque, R. J.;
 Numerical methods for conservation
  Laws, Lectures in Mathematics, ETH Zurich, Birkhauser, 1999.

\bibitem{leveque2002} Leveque, R. J.;
  Finite volume methods for hyperbolic
  problems, Cambridge University Press, Cambridge, 2002.

\bibitem{sermeus}  Sermeus, K.; Deconinck, H.;
An entropy fix for  multi dimensional
upwind residual distribution schemes, Computers and
Fluids, 34(2005), 617-640.

\bibitem{siviglia} Toro, E. F.; Silviglia, A.;
 PRICE: primitive centered schemes for hyperbolic systems,
  Int. J. Num. Meth. Fluids, 42(2003), 1263-1291.

\bibitem{Toro1} Toro, E. F.;
  Riemann solvers and numerical methods for
  fluid dynamics, A practical introduction 2nd edition, Springer,
  1999.
\bibitem{Toro2} Toro, E. F.;
 Primitive, conservative and adaptive schemes
  for hyperbolic conservation laws. In Numerical Methods for Wave
  Propagation. Toro E. F., Clarke J. F. (eds). Kluwer Academic Publishers:
  Dordrecht, 1998; 323-385.

\bibitem{wesseling} Wesseling, P.;
 Principles of computational fluid dynamics,   Springer 2001.


\end{thebibliography}

\end{document}