Peran fluks numerik dalam DG-FEM

Saya belajar teori di balik metode DG-FEM menggunakan buku Hesthaven / Warburton dan saya agak bingung tentang peran 'fluks numerik.' Saya minta maaf jika ini adalah pertanyaan mendasar, tetapi saya telah mencari dan tidak menemukan jawaban yang memuaskan untuk itu.

Pertimbangkan persamaan gelombang skalar linier:

$$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0$$ mana fluks linier diberikan sebagai$f(u) = au$ .

Seperti yang diperkenalkan dalam buku Hesthaven, untuk setiap elemen , kita berakhir dengan persamaan N , satu untuk setiap fungsi basis, memastikan bahwa residu menghilang dengan lemah:$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$$

Baik. Jadi kita melalui integrasi dengan bagian sekali untuk sampai pada 'bentuk lemah' (1) dan berintegrasi dengan bagian dua kali untuk mendapatkan 'bentuk kuat' (2). Saya akan mengadopsi bentuk integral permukaan Hesthaven yang terlalu umum tetapi mudah digeneralisasi dalam 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$$

Mengapa kita memilih fluks numerik? Mengapa kita tidak menggunakan nilai pada batas di (1) daripada menggunakan fluks? Ya, memang benar bahwa nilai kuantitas ini dapat kalikan didefinisikan di elemen, tetapi masing-masing persamaan hanya lebih 1 elemen D $au_h^k$$D^k$ , jadi mengapa hal ini?

Lebih lanjut, istilah batas integrasi kedua dengan bagian jelas menghasilkan kuantitas yang berbeda kedua kalinya dalam (2), yang tidak masuk akal bagi saya. Kami melakukan operasi yang sama! Mengapa dua syarat batas tidak dibatalkan begitu saja, membuat (2) tidak berguna? Bagaimana kami memperkenalkan informasi baru?$au_h^k$

Jelas saya kehilangan sesuatu yang penting untuk metode ini, dan saya ingin memperbaikinya. Saya telah melakukan beberapa analisis nyata dan fungsional, jadi jika ada jawaban yang lebih berbasis teori mengenai formulasi, saya ingin tahu!

Fluks numerik dipilih untuk memastikan bahwa informasi dalam masalah bergerak ke arah kurva karakteristik persamaan (berliku). Seperti disebutkan dalam komentar, fluks numerik diperlukan untuk memasangkan subproblem yang didefinisikan pada setiap elemen.

Salah satu cara untuk mendapatkan intuisi untuk peran fluks numerik adalah dengan mempertimbangkan contoh sederhana berikut ini.

$a=1$

$$\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0 \qquad\text{on $\Omega$},$$ $\Omega = [0,1]$$x = 0$$x = 1$$u(0,t) = g_D$$g_D$ .

$D_1 = [0,1/2]$$D_2 = [1/2,1]$. We could equally-well be discretizing the following set of two coupled PDEs,

$$\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*}$$ where we will couple these equations to make them equivalent to the original equation.

To make the above equations well-posed, we need to enforce boundary conditions. As before, each equation is hyperbolic, and information is traveling from left to right. Therefore, we need to enforce a boundary condition for (PDE 1) on the left endpoint of $D_1$, and a boundary condition for (PDE 2) on the left endpoint of $D_2$.

The boundary condition on the left endpoint of $D_1$ must be chosen to be $v(0,t) = g_D$ in order to remain consistent with the original problem. We also look for smooth solutions, so the boundary condition on the left endpoint of $D_2$ must be chosen to enforce continuity. This condition reads $w(1/2,t) = v(1/2,t)$.

The DG method in this case chooses the numerical fluxes precisely to enforce the above boundary conditions. If we multiply by a test function $\psi$ and integrate by parts over each element $D_k$, we obtain boundary terms of the form

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