How to Impose nonhomogeneous Neumann Boundary Condition in the DG Formulation

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 $$$$x_{L} = x_{1} and let $$D_{k} = [x_{k},x_{k+1}]$$. Consider the weak formulation of the problem $$$$\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}}$$$$ 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 $$$$u_{h}^{k} = \sum_{j=1}^{p+1} \alpha_{j}^{k}\phi_{j}^{k}$$$$ With the contribution of all elements, we get the following ODE $$$$M \dot{\alpha}-r = b-l$$$$ 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 $$$$f^{*} = \{f(u)\}+\frac{C}{2}[\hspace{-0.6mm}[u]\hspace{-0.6mm}]$$$$ 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 $$$$\frac{u_{h}^{1}-u_{h}^{0}}{h} \approx \frac{\partial u}{\partial x} = g_{1}(t)$$$$ 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?
• 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. Jan 1 at 23:41

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.