Rôle du flux numérique dans la DG-FEM

J'apprends la théorie derrière les méthodes DG-FEM en utilisant le livre Hesthaven / Warburton et je suis un peu confus quant au rôle du «flux numérique». Je m'excuse s'il s'agit d'une question fondamentale, mais j'ai cherché et je n'ai pas trouvé de réponse satisfaisante.

Considérons l'équation d'onde scalaire linéaire:

$$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0$$ où le flux linéaire est donné comme$f(u) = au$.

Comme introduit dans le livre de Hesthaven, pour chaque élément , on se retrouve avec N équations, une pour chaque fonction de base, imposant que le résidu s'évanouit faiblement:$k$$N$

$$R_h(x,t) = \frac{\partial u_h}{\partial t} + \frac{\partial au_h}{\partial x}$$

$$\int_{D^k} R_h(x,t) \psi_n(x) \, dx = 0$$

Bien. Nous passons donc par l'intégration par parties une fois pour arriver à la «forme faible» (1) et nous intégrons par parties deux fois pour obtenir la «forme forte» (2). J'adopterai la forme intégrale de surface de Hesthaven mais trop généralisée en 1D:

(1)

$$\int_{D^k} \left( \frac{\partial u_h^k}{\partial t} \psi_n-au_h^k\frac{d \psi_n}{d x} \right)\, dx = - \int_{\partial D^k} \hat{n}\cdot(au_h)^*\psi_n \,dx \qquad 1 \leq n\leq N$$

(2)

$$\int_{D^k} R_h \psi_n\, dx = \int_{\partial D^k} \hat{n}\cdot \left( au_h^k-(au_h)^* \right)\psi_n \,dx \qquad 1 \leq n\leq N$$

Pourquoi choisissons-nous un flux numérique? Pourquoi n'utilisons-nous pas la valeur d' à la frontière de (1) au lieu d'utiliser un flux? Oui, il est vrai que la valeur de cette quantité peut être définie par multiplication entre les éléments, mais chaque équation ne dépasse que 1 élément D k , alors pourquoi est-ce important?$au_h^k$$D^k$

De plus, le terme limite de la deuxième intégration par parties donne clairement une quantité différente la deuxième fois en (2), ce qui n'a aucun sens pour moi. Nous faisons la même opération! Pourquoi les deux termes de la frontière ne s'annuleraient-ils pas, rendant (2) inutile? Comment avons-nous introduit de nouvelles informations?$au_h^k$

De toute évidence, il me manque quelque chose de crucial pour la méthode, et je voudrais y remédier. J'ai fait une analyse réelle et fonctionnelle, donc s'il y a une réponse plus théorique concernant la formulation, j'aimerais savoir!

Le flux numérique est choisi pour garantir que les informations du problème voyagent dans le sens des courbes caractéristiques de l'équation (remontée). Comme mentionné dans les commentaires, le flux numérique est nécessaire pour coupler les sous-problèmes définis sur chaque élément.

Une façon d'obtenir une intuition pour le rôle du flux numérique est de considérer l'exemple simple suivant.

Considérons l'équation d'advection scalaire (où pour simplifier $a=1$ )

$$\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0 \qquad\text{on $\Omega$},$$ où le domaine est donné par$\Omega = [0,1]$ . Comme il s'agit d'une équation hyperbolique et que les informations se propagent de gauche à droite, nous devons appliquer une condition aux limites à$x = 0$ (mais pas à$x = 1$ ). Pour le concret, supposons que nous appliquons la condition de Dirichlet$u(0,t) = g_D$ pour certains$g_D$ donnés.

Supposons maintenant on discrétise cette équation en utilisant la méthode DG, et nous utilisons deux éléments, $D_1 = [0,1/2]$ et $D_2 = [1/2,1]$ . On pourrait tout aussi bien discrétiser l'ensemble suivant de deux PDE couplés,

$$\begin{align*} \text{(PDE 1):}&& \quad v_t + v_x &= 0 \quad\text{on $D_1$},\\ \text{(PDE 2):}&& \quad w_t + w_x &= 0 \quad\text{on $D_2$}, \end{align*}$$ où nous allons coupler ces équations pour les rendre équivalentes à l'équation d'origine.

$D_1$$D_2$

$D_1$$v(0,t) = g_D$$D_2$$w(1/2,t) = v(1/2,t)$

$\psi$$D_k$

$$\begin{align*} \int_{\partial D_1} \hat{n} \cdot v \psi \, dx &= \left[ v\psi \right]_0^{1/2}\\ \int_{\partial D_2} \hat{n} \cdot w \psi \, dx &= \left[ w\psi\right]_{1/2}^1 \\ \end{align*}$$ In order to "weakly" enforce the boundary conditions, we replace $v$ and $w$ with the prescribed values at those points where boundary conditions are specified (i.e. the left endpoints of $D_1$ and $D_2$). This means we replace $v(0,t)$ by $g_D$ and $w(1/2,t)$ by $v(1/2,t)$ in the boundary integrals.

In other words, we define $u_h^* = g_D$ at $x = 0$, and $u_h^* = v(1/2,t)$ at $x=1/2$, and we recover exactly the standard upwind flux that is used in the DG method.

Looking at things this way, we can consider the numerical flux functions as weakly enforcing the boundary conditions on each element that are required to couple the equations in such a way that respects the characteristic structure of the equations.

For equations more complicated than constant-coefficient advection, information may not propagate always in the same direction, and so the numerical flux must be determined by solving (or approximating the solution to) a Riemann problem at the interface. This is discussed for linear problems in Section 2.4 of Hesthaven's book.

Very loosely speaking there are two things most discretization techniques need in order to converge to the actual solution of your PDE as you increase their approximation quality, regardless if you're using DG or not:

1. Consistency (If a function $u$ satisfies the PDE, then it also satisfies your weak formulation)
2. Stability (small changes in data result in small changes in answer)

The first steps of a DG derivation where you integrate by parts on each mesh element preserves (1) because you are starting with the PDE and only applying legal operations from there.

This does not give you (2) though. You can see this yourself by trying to assemble the matrix of a partially formulated DG weak form and looking at its eigenvalues - for time dependent problem we want them all on the left half plane, but without a proper numerical flux they will be everywhere. This leads to a solution that explodes exponentially in time even if the physical problem does not.

Thus you need to add terms to your formulation so that (2) is satisfied, but without hurting (1). This is tricky to do, but not impossible. You can replace function values with cell wall averages without hurting consistency, and you can always add the jump in at cell walls without hurting consistency (because for a solution $u$ with suitable smoothness property, the jump is just 0!)

The trick is to take combinations of jumps and averages and combine them in a way that your scheme is still consistent but also stable. After that a convergence theorem usually reveals itself.

This is the basics, but you can also often bring in additional physics into the numerical flux so that it doesn't simply satisfy these mathematical requirements but also plays nicely with conservation principles.

When you choose the test function equal to the trial function in the DG method you are creating an optimization problem. That is, you have a Galerkin rather than a Petrov-Galerkin method. You are seeking the time derivatives of the trial function amplitudes that will minimise the element residual in the L2 norm, and you make this miminimisation on the assumption of a given flux function at inflow.