Consider the following partial differential equation
\begin{align}
\frac{\partial u}{\partial t}+\frac{\partial f}{\partial x} &= g(x,t), \ \ x\in \Omega = [x_{L},x_{R}] \\
u(x,0) &= u_{0}(x) \\
\frac{\partial u}{\partial x}(x_{L},t) &= g_{1}(t)
\end{align}
where the flux function $f(u,x)$ is nonlinear. Let $\{x_{i}\}_{i=1}^{N+1}$ be a partition of $\Omega$ with $N$ elements such that
\begin{equation} x_{L} = x_{1}<x_{2}<\cdots<x_{N}<x_{N+1} = x_{R} \end{equation}
and let $D_{k} = [x_{k},x_{k+1}]$. Consider the weak formulation of the problem
\begin{equation}
\int_{D_{k}} \Bigg[ \frac{\partial u}{\partial t}\phi_{i}^{k} -f(u,x)\frac{d \phi_{i}^{k}}{d x}\Bigg] = \int_{D_{k}} g(x,t)\phi_{i}^{k}-\Big[f^{*}\phi_{i}^{k} \Big]_{x_{k}}^{x_{k+1}}
\end{equation}
where $f^{*}$ is called numerical flux. If we consider approximations with polynomials of degree $p$, the approximation $u_{h}^{k}$ of $u$ in the element $D_{k}$ is given by
\begin{equation}
u_{h}^{k} = \sum_{j=1}^{p+1} \alpha_{j}^{k}\phi_{j}^{k}
\end{equation}
With the contribution of all elements, we get the following ODE
\begin{equation}
M \dot{\alpha}-r = b-l
\end{equation}
where the entries of the vectors $b$ and $r$ in the element $D_{k}$ are given by
\begin{align}
[b_{k}]_{i} &= \int_{D_{k}} g(x,t)\phi_{i}^{k} dx \\
[r_{k}]_{i} &= \int_{D_{k}} f(u,x)\frac{d \phi_{i}^{k}}{d x} dx
\end{align}
The numerical flux is given by
\begin{equation}
f^{*} = \{f(u)\}+\frac{C}{2}[\hspace{-0.6mm}[u]\hspace{-0.6mm}]
\end{equation}
where $C = \max{ |f'(u)| }$, $\{\cdot\}$ is the average and $[\hspace{-0.6mm}[\cdot]\hspace{-0.6mm}]$ is the jump.
The numerical flux vector $l$ in the element $k$ has entries
\begin{align}
\Big[f^{*}\phi_{i}^{k} \Big]_{x_{k}}^{x_{k+1}} =& \Bigg[ \frac{1}{2}f(u_{h}^{k},x_{k+1})+ \frac{1}{2}f(u_{h}^{k+1},x_{k+1})+\frac{C}{2}u_{h}^{k}-\frac{C}{2}u_{h}^{k+1} \Bigg] \phi_{i}^{k} (x_{k+1}) \\&-\Bigg[ \frac{1}{2}f(u_{h}^{k-1},x_{k})+ \frac{1}{2}f(u_{h}^{k},x_{k})+\frac{C}{2}u_{h}^{k-1}-\frac{C}{2}u_{h}^{k} \Bigg] \phi_{i}^{k} (x_{k})
\end{align}
Case 1: $g_{1}(t) = 0$. In the first element, i.e, for $k=1$, I take $u_{h}^{0} = u_{h}^{1}$ and my implementation works fine.
Case 2: $g_{1}(t) \neq 0$. For $k=1$, I use the following approximation \begin{equation} \frac{u_{h}^{1}-u_{h}^{0}}{h} \approx \frac{\partial u}{\partial x} = g_{1}(t) \end{equation} and then I solve for $u_{h}^{0}$, that is $u_{h}^{0} = u_{h}^{1}-hg_{1}(t)$, where $h$ is the mesh size, but this doesn't work.
My questions are:
- Do you know books or papers that explain how to solve this PDE with the DG method?
- Do you know how should I take $\{f(u)\}$, $[\hspace{-0.6mm}[u]\hspace{-0.6mm}]$ in order to impose the Neumann boundary condition?
