I got stuck with Hestaven/Warburton's dG-FEM Matlab code. Starting with the file AdvecRHS1D.m, we see in line 11

du(:) = (u(vmapM)-u(vmapP)).*(a*nx(:)-(1-alpha)*abs(a*nx(:)))/2;

which calculates the jump $[u]$ between adjacent elements and multiplies it with some factor, giving

$$du=[u]\frac{an-(1-\alpha)|an|}{2}=[u]a\frac{n-(1-\alpha)}{2}, n\in \{-1,1\},\\ a>,\ 1\geq\alpha\geq 0$$

which looks very similar to the numerical flux described as

$$(au)^*=\{\{u\}\} + a\frac{1-\alpha}{2}[u]$$

if we leave out the average $\{\{u\}\}$. Yet these terms are obviously not the same.

Then, for the computation of the right hand side of the PDE, this happens:

rhsu = -a*rx.*(Dr*u) + LIFT*(Fscale.*(du));

which is supposed to calculate the integral

$$\left[l^k(x)(au^k_h) - (au)^* \right]_{x_l^k}^{x_r^k} = \oint_{\delta D^k}\hat{n}(au_h - (au)^*)l_i(r)$$

But I don't understand how this line of code is expressing this.

EDIT: I found a paper, see p.5, that gives some explanation on the calculation of rhsu. But still, it remains pretty unclear to me.


We can break down the code

rhsu = -a*rx.*(Dr*u) + LIFT*(Fscale.*(du));

into the two parts




The first part is simply taking a derivative in reference space by multiplying with the matrix Dr, then transforming into physical space by multiplying with the Jacobian rx and then the coefficient a. Next we need to do the face integral which is

$$ \oint_{\delta D^k}\hat{n}(au_h - (au)^*)l_i(r).$$

In this case, LIFT just returns the matrix $$M^{-1} \mathcal{E}$$ which was previously defined in the book, du was previously defined representing $$(au_h - (au)^*).$$

In short the DG semi-discretization is

$$ \frac{du}{dt} = -a M^{-1} S u - M^{-1} \mathcal{E} N (au^* - E(au)) $$

which is what rhsu gives with $M^{-1} S = rx*(Dr)$, $LIFT = M^{-1} \mathcal{E}$, $N (au^* - E(au)) = -du$.

| cite | improve this answer | |
  • 1
    $\begingroup$ Thanks! can you elaborate on what LIFT does? $\endgroup$ – dba Jul 8 '19 at 19:10
  • 1
    $\begingroup$ As I mentioned, LIFT just returns the matrix $M^{-1} \mathcal{E}$. $\mathcal{E}$ is the matrix that is equivalent to $l_i(r)$ and $M^{-1}$ is obtained because we invert $M \frac{du}{dt}$. $\endgroup$ – Vikram Jul 9 '19 at 9:15
  • $\begingroup$ Thanks! How does $\mathcal{E}$ refer to $l_i$? In the code it is defined as a $N_p \times 2$ matrix. With the first and the last entry as 1, 0 everywhere else. The Integral is defined as $(u_h - u^*)|_{r_{np}}e_{N_p} - (u_h - u^*)|_{r_1}e_{1}$. Meaning, that it must be vector valued with the flux values of the boundary in the first and last component. Is that correct? $\endgroup$ – dba Jul 9 '19 at 14:55
  • 1
    $\begingroup$ No, $\mathcal{E}$ has nothing to do with flux values. It is just a restriction operator with two rows. The first row is for interpolating to the left face and the second row for the right face. This is exactly equivalent to $l_i(r)$ with $r = x_r$ and $r=x_l$. For Gauss Lobatto points, each row has 0 everywhere else except at the first and last points, however this will not be the case for Gauss points. $\endgroup$ – Vikram Jul 9 '19 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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