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:

  1. Do you know books or papers that explain how to solve this PDE with the DG method?
  2. Do you know how should I take $\{f(u)\}$, $[\hspace{-0.6mm}[u]\hspace{-0.6mm}]$ in order to impose the Neumann boundary condition?
  • $\begingroup$ Both cases do not make sense and are inconsistent. Keep in mind that your PDE is based on a polynomial approximation. Taking the difference of two points does not give you the polynomial gradient inside the element. $\endgroup$
    – ConvexHull
    Commented Jan 1, 2022 at 23:41

1 Answer 1


You cannot specify just a boundary condition on $\partial_x u$ at $x_L$.

Remember that in DG your solution is composed of piecewise discontinuous polynomials; at every interface you can have effectively any jump in $u$, which is separate from the slope of u to the left or right of the interface.

Unless $f(u,x)$ has any dependence on $\partial_x u$ (in which cases you need to use elliptic/higher derivative methods for DG such as LDG or interior penalty methods), it turns out that specifying $\partial_x u$ as a boundary condition makes no difference to the solution!

For example, suppose that $x_L = 0$, and consider the two following states for $u$: $$ u_1 = \begin{cases} x-1 & x < 0\\ x+1 & x > 0 \end{cases}\\ u_2 = \begin{cases} -x-1 & x < 0\\ x+1 & x > 0 \end{cases} $$ Using your definition of numerical flux, then we get in both cases $$ f^* = \frac{f(1,0) + f(-1,0)}{2} + C $$

Instead of specifying what $\partial_x u$ is, instead you can specify the amount of flux flowing into or out of your domain: $$ f(u,x) \cdot n^-|_{x_L} = g_2(t) $$ Then when it comes to calculating the numerical flux at the boundary, you just plug in the BC directly.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.