3
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
4
$\begingroup$

We can break down the code 

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

into the two parts

-a*rx.*(Dr*u)

and

LIFT*(Fscale.*(du));

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$.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Thanks! can you elaborate on what LIFT does? $\endgroup$
    – dba
    Commented Jul 8, 2019 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
    Commented Jul 9, 2019 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
    Commented Jul 9, 2019 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
    Commented Jul 9, 2019 at 16:10

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.