How to compute numerical fluxes in the local discontinuous galerkin method for poisson equation 1D

Some days ago I began to study the local discontinuous galerkin (LDG) method, this is my first time working with a discontinuous method, so I decided to solve the poisson equation in 1D to learn the new ideas, this is what I have done so far.

$$\mathbf{Weak \ Formulation}$$

Consider the following ordinary differential equation:

\begin{align} -\left( k u' \right)' & = f \ \ \text{in} \ [a,b] \tag{+}\label{ec} \\ u(a) & = u_{a} \\ u(b) & = u_{b} \end{align}

where $$k(x) > 0 \ \ \forall x \in [a,b]$$. Let $$\{x_{i}\}_{i=1}^{N}$$ be a partition of $$[a,b]$$ such that

$$a = x_{1}

Now we define $$I_{k} = [x_{k},x_{k+1}]$$ for $$k=1,...N-1$$. The LDG method starts introducing the auxiliar variable
$$q = ku'$$, and rewriting (\ref{ec}) as a system of first order differential equations as follows:

\begin{align} k^{-1}q & = u' \tag{1} \label{axq} \\ -q' & = f \tag{2} \label{ecldg} \end{align}

now we multiply (\ref{axq}) and (\ref{ecldg}) by test functions $$r$$ and $$v$$ (respectively) then integrating over an element $$I_{k}$$ and using integration by parts we get

\begin{align} \int_{I_{k}} k^{-1}qr- ur\big|_{x_{k}}^{x_{k+1}}+\int_{I_{k}} ur' & = 0 \\ \int_{I_{k}} qv'-qv\big|_{x_{k}}^{x_{k+1}} - \int_{I_{k}}fv & = 0 \end{align}

adding an stabilization term $$s$$ to the last equation we have

\begin{align} \int_{I_{k}} k^{-1}qr- \widehat{u}r\big|_{x_{k}}^{x_{k+1}}+\int_{I_{k}} ur' & = 0 \\ \int_{I_{k}} qv'-\widehat{q}v\big|_{x_{k}}^{x_{k+1}}+vs\big|_{x_{k}}^{x_{k+1}} - \int_{I_{k}}fv & = 0 \end{align}

$$\mathbf{Spatial \ Discretization}$$

let $$u^{k} = \sum_{j=1}^{p+1} U_{j}^{k}\phi_{j}^{k}$$ and $$q^{k} = \sum_{j=1}^{p+1} Q_{j}^{k}\phi_{j}^{k}$$ denote the approximation of $$u$$ and $$q$$ (respectively) in the element $$I_{k}$$, where $$\{ \phi_{j}^{k} \}_{j=1}^{p+1}$$ is the basis in this element and $$p$$ is the maximum degree of the polynomials in the basis, it follows that:

\begin{align} \int_{I_{k}} k^{-1}\left( \sum_{j=1}^{p+1} Q_{j}^{k}\phi_{j}^{k} \right)r- \widehat{u}r\big|_{x_{k}}^{x_{k+1}}+\int_{I_{k}} \left( \sum_{j=1}^{p+1} U_{j}^{k}\phi_{j}^{k} \right)r' & = 0 \\ \int_{I_{k}} \left( \sum_{j=1}^{p+1} Q_{j}^{k}\phi_{j}^{k} \right)v'-\widehat{q}v\big|_{x_{k}}^{x_{k+1}}+vs\big|_{x_{k}}^{x_{k+1}} - \int_{I_{k}}fv & = 0 \end{align}

simplifying and taking $$r = \phi_{i}^{k}$$ and $$v = \phi_{i}^{k}$$ for $$i=1,...,p+1$$ we obtain

\begin{align} \sum_{j=1}^{p+1} Q_{j}^{k} \int_{I_{k}} k^{-1} \phi_{j}^{k} \phi_{i}^{k}- \widehat{u}\phi_{i}^{k}\big|_{x_{k}}^{x_{k+1}}+ \sum_{j=1}^{p+1} U_{j}^{k} \int_{I_{k}} \phi_{j}^{k}\partial_{x}\phi_{i}^{k} & = 0 \tag{*} \label{la1} \\ \sum_{j=1}^{p+1} Q_{j}^{k} \int_{I_{k}} \phi_{j}^{k} \partial_{x}\phi_{i}^{k}-\widehat{q}\phi_{i}^{k}\big|_{x_{k}}^{x_{k+1}}+\phi_{i}^{k}s\big|_{x_{k}}^{x_{k+1}} - \int_{I_{k}}f\phi_{i}^{k} & = 0 \tag{**}\label{la2} \end{align}

$$\mathbf{Numerical \ Fluxes}$$

$$\widehat{u}$$ and $$\widehat{q}$$ in equations (\ref{la1}) and (\ref{la2}) are called numerical fluxes, and for $$\alpha_{k} \in [0,1]$$ are defined as follows

\begin{align} \widehat{u}(x_{k}) & = (1-\alpha_{k})u^{k-1}(x_{k}^{-})+\alpha_{k} u^{k}(x_{k}^{+}) \tag{3} \label{fx1} \\ \widehat{q}(x_{k}) & = \alpha_{k}q^{k-1}(x_{k}^{-})+(1-\alpha_{k}) q^{k}(x_{k}^{+}) \tag{4} \label{fx2} \end{align}

for $$k=2,...,N-1$$. For the boundary points the definition is

\begin{align} \widehat{u}(x_{1}) & = u_{a} \hspace{1.4cm} \ \widehat{u}(x_{N}) = u_{b} \\ \widehat{q}(x_{1}) & = q(x_{1}^{+}) \hspace{1cm} \widehat{q}(x_{N}) = q(x_{N}^{-}) \end{align}

The stabilization term $$s$$ in (\ref{la2}) is given by

$$$$s(x_{k}) = \frac{\eta_{k}}{l_{k}} \left( u^{k}(x_{k}^{+})-u^{k-1}(x_{k}^{-}) \right) \tag{5} \label{st}$$$$

for $$k=2,...,N-1$$. Where $$\eta_{k}>0$$ and $$l_{k} = \max\{ h_{k} , h_{k-1} \}$$ , $$h_{k} = x_{k+1}-x_{k}$$

$$\mathbf{Compute \ Blocks}$$

Let $$\widehat{I} = [-1,1]$$ and consider the mapping $$x:\widehat{I}\rightarrow I_{k}$$ such that

\begin{align} x(\widehat{x}) & = \frac{x_{k}}{2}(1-\widehat{x})+\frac{x_{k+1}}{2}(1+\widehat{x}) \\ x(\widehat{x}) & = \frac{x_{k}+x_{k+1}}{2}+\left(\frac{x_{k+1}-x_{k}}{2}\right)\widehat{x} \end{align}

I don't have problems computing the integrals in (\ref{la1}) and (\ref{la2}), these are given by

\begin{align} \int_{I_{k}} k^{-1} \phi_{i}^{k}\phi_{j}^{k} & = \int_{\widehat{I}} k^{-1}(x(\widehat{x}))\phi_{i}^{k}(x(\widehat{x}))\phi_{j}^{k}(x(\widehat{x})) \Big|\frac{dx}{d\widehat{x}}\Big| \\ & = \int_{\widehat{I}} k^{-1}(x(\widehat{x}))\widehat{\phi}_{i}(\widehat{x})\widehat{\phi}_{j}(\widehat{x}) \frac{h_{k}}{2} \end{align}

\begin{align} \int_{I_{k}} \frac{d\phi_{i}^{k}}{dx} \phi_{j}^{k} & = \int_{\widehat{I}} \frac{2}{h_{k}} \frac{d \widehat{\phi}_{i}}{d\widehat{x}} \widehat{\phi}_{j} \frac{h_{k}}{2} \\ & = \int_{\widehat{I}} \frac{d \widehat{\phi}_{i}}{d\widehat{x}} \widehat{\phi}_{j} \end{align}

and

\begin{align} \int_{I_{k}} f\phi_{i}^{k} = \int_{\widehat{I}} f(x(\widehat{x}))\widehat{\phi}_{i} \frac{h_{k}}{2} \end{align}

But I do have problems when I try to compute $$\widehat{u}(x_{k})$$ (\ref{fx1}), $$\widehat{q}(x_{k})$$ (\ref{fx2}) and $$s(x_{k})$$ (\ref{st})
in (\ref{la1}) and (\ref{la2}) because the expressions of the form $$u^{k}(x^{+}) = \lim_{\epsilon\rightarrow 0^{+}} u^{k}(x_{k}+\epsilon)$$ and $$u^{k}(x^{-}) = \lim_{\epsilon\rightarrow 0^{+}} u^{k}(x_{k}-\epsilon)$$ are new for me, I don't know how to deal with them.

$$\mathbf{ Questions: }$$

1. Could you explain how to compute the stabilization term and the numerical fluxes in the equations (\ref{la1}) and (\ref{la2})?

2. Do you know some references where I can study topics related to question 1.? ....thanks!

• Check lacan.upc.edu/dg2017 . There are lecture notes and MATLAB codes which may be helpful. – Abdullah Ali Sivas Jan 10 at 2:01
• It's not clear to me what the hats on u and q mean; why do you introduce them (just before discretization paragraph)? – Gabriele Jan 15 at 18:54