Trefftz discontinuous Galerkin basis functions for a class of Friedrichs systems coming from linear transport

Abstract

The Trefftz discontinuous Galerkin (TDG) method provides natural well-balanced (WB)and asymptotic preserving (AP) discretization, since exact solutions are used locally in the basis functions. However, one difficult point may be the construction of such solutions, which is a necessary first step in order to apply the TDG scheme. This work deals with the construction of solutions to Friedrichs systems with relaxation with application to the spherical harmonics approximation of the transport equation (the so-called P_N models). Various exponential and polynomial solutions are constructed. Two numerical tests on the P₃ model illustrate the good accuracy of the TDG method. They show that the exponential solutions lead to accurate schemes to capture boundary layers on a coarse mesh and that a combination of exponential and polynomial solutions is efficient in a regime with vanishing absorption coefficient.

Appendices

Appendix A: Time dependent solutions of equation (1)

We give some possible ways to get time dependent solutions to the P_N model. These solutions can also be used as basis functions for the TDG method when considering space-time meshes. Many other time-dependent solutions can be constructed starting from [5, 10].

• A first possibility is to consider one-dimensional solution under the form as follows:

  $$ \mathbf{v}(t,x) = \mathbf{q}(t,x) e^{\lambda x}, $$

  where $\mathbf {q}(t,x) \in \mathbb {R}^{m}$ is polynomial vector in x and t. A concrete example is given in [33] [Proposition 4.2], for the case of the P₁ model. Then, one can use the rotational invariance of the P_N model and gets the following solutions:

  $$ \mathbf{v}(t,\mathbf{x}) = U_{\theta}\mathbf{q}(t,x \cos \theta + y \sin \theta) e^{\lambda (x \cos \theta + y \sin \theta)}. $$

  (29)

  Another possibility is to search directly for two-dimensional solutions under the form as follows:

  $$ \mathbf{v}(t,\mathbf{x}) = \mathbf{p}(t,\mathbf{x}) e^{\lambda (x \cos \theta + y \sin \theta)}, $$

  (30)

  where $\mathbf {p}(t,\mathbf {x}) \in \mathbb {R}^{m}$ is polynomial vector in x, y ant t. Note that the functions obtained with (30) may differ from the functions (29). A complete example is given in [32], [Section 5.2.2] for the P₁ model.

• Another possibility is to consider time dependent solutions under the form as follows:

  $$ \mathbf{v}(t,\mathbf{x}) = \mathbf{g}(\mathbf{x})e^{\alpha t}, $$

  (31)

  with $\alpha \in \mathbb {R}$. One can inject this solution in the P_N model (1). One gets after removing the exponentials as follows:

  $$ \left( A_{1} \partial_{x} + A_{2} \partial_{y} + (R+ \varepsilon\alpha I_{m})\right) \mathbf{g}(\mathbf{x}) = \mathbf{0}, $$

  where I_m is the identity matrix of $\mathbb {R}^{m \times m}$. The function g(x) is very similar to the stationary solutions already calculated in Section 3.1. The matrix R is just replaced by the matrix $\widetilde {R}:=R+\varepsilon \alpha I_{m}$.

  For example, if one takes α such that σ_a + εα > 0, then g(x) is one of the exponential solutions (10). In particular, if α > 0, the functions (31) naturally degenerate toward nontrivial time-dependent solutions when $\sigma _{a} \rightarrow 0$. This is one asset of the functions (31).

Appendix B: A variational Trefftz discontinuous Galerkin method

The objective of this Appendix is to provide the basic material for the construction of the TDG discussed in this work. The notations are taken from [33]. We write v = (v₁,...,v_m)^T where ^T denotes the transpose. The domain Ω is partitioned: let $\mathcal {T}_{S,h}$ be a partition of Ω_S into polyhedral subdomains Ω_S,k. For a stationary problem, we set $\mathcal {T}_{h} = \mathcal {T}_{S,h}$. For time-dependent problems, we consider $N \in \mathbb {N}^{*}$ time steps [t_i,t_i+ 1], 0 = t₀ < t₁ < ... < t_N = T. Then, $\mathcal {T}_{h}$ is the partition of Ω made of the polyhedral subdomains ${\Omega }_{k^{n}}={\Omega }_{S,k} \times [t_{n},t_{n+1}]$; here, the index k is related to $\mathcal {T}_{S,h}$ and n to the time step. In both cases, we use the general notation Ω_k when designing a cell of $\mathcal {T}_{h}$. With this notations, the time variable is treated as the space variables in the domain Ω. One must therefore be careful that $\mathcal {T}_{h}$ denotes a spatial mesh for stationary model and a space-time mesh for time-dependent model.

The matrix R is constant in each cell. The broken Sobolev space is as follows:

$$ H^{1}(\mathcal{T}_{h}) := \{\mathbf{v} \in L^{2}({\Omega}), \ \mathbf{v}_{|{\Omega}_{k}} \in H^{1}({\Omega}_{k}) \ \forall {\Omega}_{k} \in \mathcal{T}_{h} \}. $$

In the following, we assume $\mathbf {u} \in H^{1}(\mathcal {T}_{h})$. For convenience, we rewrite the model under the form $L = {\sum }_{i} A_{i}\partial _{i} + R$ and Lu = 0. We consider the adjoint operator $L^{*} = - {\sum }_{i} A_{i}\partial _{i} + R$. Multiplying the equation by $\mathbf {v} \in H^{1}(\mathcal {T}_{h})$ and integrating gives on Ω as follows:

$$ \sum\limits_{k} {\int}_{{\Omega}_{k}} \mathbf{v}_{k}^{T} L \mathbf{u}_{k} = 0, $$

(32)

where $\mathbf {v}_{k} = \mathbf {v}_{|{\Omega }_{k}}$, $\mathbf {u}_{k} = \mathbf {u}_{|{\Omega }_{k}}$. Integrating by parts one gets the following:

$$ \sum\limits_{k} {\int}_{{\Omega}_{k}} \left( L^{*} \mathbf{v}_{k}\right)^{T} \mathbf{u}_{k} + \sum\limits_{k} {\int}_{\partial{\Omega}_{k}}\mathbf{v}_{k}^{T}M_{k} \mathbf{u}_{k} = 0, $$

where ∂Ω_k is the contour of the element Ω_k with an outward unit normal $\mathbf {n}_{k} = (n_{t},n_{x_{1}},...,n_{x_{d}})^{T}$, $M = A_{0} n_{t} + {\sum }_{i} A_{i} n_{i}$ and M_k = M(n_k(x)). Denoting Σ_kj = ∂Ω_k ∩ ∂Ω_j (n_k + n_j = 0 on Σ_kj), one can writes the following:

$$ \begin{array}{@{}rcl@{}} \sum\limits_{k} {\int}_{{\Omega}_{k}} \left( L^{*} \mathbf{v}_{k}\right)^{T} \mathbf{u}_{k} &+& \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} \left( \mathbf{v}^{T} M \mathbf{u}\right)_{k} + \left( \mathbf{v}^{T} M \mathbf{u}\right)_{j} \\ &+& \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{+} \mathbf{u}_{k} = - \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{-} \mathbf{g} . \end{array} $$

For u satisfying the equation, the normal flux is continuous in a weak sense as follows:

$$ M_{k} \mathbf{u}_{k} = - M_{j} \mathbf{u}_{j} = f_{\text{kj}}(\mathbf{u}_{k},\mathbf{u}_{j}), \quad \text{on } {\Sigma}_{\text{kj}} $$

(33)

where f_kj is a numerical flux on Σ_kj defined below. One has the following:

$$ \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (\mathbf{v}^{T} M \mathbf{u})_{k} + (\mathbf{v}^{T} M \mathbf{u})_{j} = \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (\mathbf{v}_{k} - \mathbf{v}_{j})^{T} f_{\text{kj}}(\mathbf{u}_{k},\mathbf{u}_{j}) . $$

Because M is symmetric one can decompose M under the form M = M⁺ + M⁻ where M⁺ is a non-negative matrix and M⁻ is a non-positive matrix. In the following, we will consider the upwind flux $f_{\text {kj}}(\mathbf {u}_{k},\mathbf {u}_{j}) = M_{\text {kj}}^{+}\mathbf {u}_{k} + M_{\text {kj}}^{-}\mathbf {u}_{j},$ where $M_{\text {kj}} = M_{k|{\Sigma }_{\text {kj}}}$. Finally, one has the following:

$$ \begin{array}{@{}rcl@{}} \sum\limits_{k} {\int}_{{\Omega}_{k}} \left( L^{*} \mathbf{v}_{k}\right)^{T} \mathbf{u}_{k} &+& \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (\mathbf{v}_{k} - \mathbf{v}_{j})^{T} (M_{\text{kj}}^{+} \mathbf{u}_{k} + M_{\text{kj}}^{-} \mathbf{u}_{j}) \end{array} $$

(34)

$$ \begin{array}{@{}rcl@{}} &+& \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{+} \mathbf{u}_{k} = -\sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{-} \mathbf{g} . \end{array} $$

(35)

We define the bilinear form $a_{\text {DG}}:H^{1}(\mathcal {T}_{h}) \times H^{1}(\mathcal {T}_{h}) \rightarrow \mathbb {R}$ and the linear form $l:H^{1}(\mathcal {T}_{h}) \rightarrow \mathbb {R}$ as follows:

$$ \begin{array}{@{}rcl@{}} a_{\text{DG}}(\mathbf{u},\mathbf{v}) &=& \sum\limits_{k} {\int}_{{\Omega}_{k}} (L^{*} \mathbf{v}_{k})^{T} \mathbf{u}_{k} + \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (\mathbf{v}_{k} - \mathbf{v}_{j})^{T} (M_{\text{kj}}^{+} \mathbf{u}_{k} + M_{\text{kj}}^{-} \mathbf{u}_{j})\\ &&+ \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{+} \mathbf{u}_{k} , \quad \mathbf{u},\mathbf{v} \in H^{1}(\mathcal{T}_{h}),\\ l(\mathbf{v}) &=& -\sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{-} \mathbf{g} , \quad \mathbf{v} \in H^{1}(\mathcal{T}_{h}). \end{array} $$

(36)

One can rewrite (34) as $a_{\text {DG}}(\mathbf {u},\mathbf {v}) = l(\mathbf {v}), \ \forall \mathbf {v} \in H^{1}(\mathcal {T}_{h})$. The standard discontinuous Galerkin method is based on the use of polynomial basis functions [14, 31]. Define $\mathbb {P}_{q}^{d}$, the space of polynomials of d variables, of total degree at most q and the broken polynomial space as follows:

$$ \mathbb{P}_{q}^{d}(\mathcal{T}_{h}) := \{\mathbf{v} \in L^{2}({\Omega}), \mathbf{v}_{|{\Omega}_{k}} \in \mathbb{P}_{q}^{d} \ \forall {\Omega}_{k} \in \mathcal{T}_{h} \} \subset H^{1}(\mathcal{T}_{h}). $$

Definition B.1

Assume $P_{m}(\mathcal {T}_{h})$ is a finite subspace of $H^{1}(\mathcal {T}_{h})$, for example, $P_{m}(\mathcal {T}_{h})=\mathbb {P}_{q}^{d}(\mathcal {T}_{h})$. The standard upwind discontinuous Galerkin method for Friedrichs systems is formulated as follows:

$$ \left\{\begin{array}{ll} \text{find } \mathbf{u}_{h} \in P_{m}(\mathcal{T}_{h}) \text{ such that }\\ a_{\text{DG}}(\mathbf{u}_{h},\mathbf{v}_{h}) =l(\mathbf{v}_{h}), \quad \forall \mathbf{v}_{h} \in P_{m}(\mathcal{T}_{h}). \end{array}\right. $$

(37)

For Trefftz methods, one takes basis functions which are exact solutions to the equations in each cell as follows:

$$ V(\mathcal{T}_{h}) = \{\mathbf{v} \in H^{1}(\mathcal{T}_{h}), L \mathbf{v}_{k} = \mathbf{0} \quad \forall {\Omega}_{k} \in \mathcal{T}_{h}\} \subset H^{1}(\mathcal{T}_{h}). $$

The space $V(\mathcal {T}_{h})$ is a genuine subspace of $H^{1}(\mathcal {T}_{h})$ except in the case L = 0 which is of no interest. Starting from the bilinear form a_DG from (36), one notices that the volume term can be written for all functions in $V(\mathcal {T}_{h})$ as follows:

$$ {\int}_{{\Omega}_{k}} \left( L^{*} \mathbf{v}_{k}\right)^{T} \mathbf{u}_{k} = 2 {\int}_{{\Omega}_{k}} \mathbf{v}_{k}^{T} R\mathbf{u}_{k} , \quad \forall \mathbf{u},\mathbf{v} \in V(\mathcal{T}_{h}). $$

(38)

One can therefore define a bilinear form $a_{T}:V(\mathcal {T}_{h}) \times V(\mathcal {T}_{h}) \rightarrow \mathbb {R}$ as follows:

$$ \begin{array}{@{}rcl@{}} a_{T}(\mathbf{u},\mathbf{v}) &=& \sum\limits_{k} 2 {\int}_{{\Omega}_{k}} \mathbf{v}_{k}^{T} R\mathbf{u}_{k} + \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (\mathbf{v}_{k} - \mathbf{v}_{j})^{T} (M_{\text{kj}}^{+} \mathbf{u}_{k} + M_{\text{kj}}^{-} \mathbf{u}_{j})\\ &&+ \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{+} \mathbf{u}_{k} , \quad \mathbf{u},\mathbf{v} \in V(\mathcal{T}_{h}). \end{array} $$

(39)

Thanks to an integration by part, one gets an equivalent formulation of the bilinear form a_T(⋅,⋅) as follows:

$$ a_{T}(\mathbf{u},\mathbf{v}) = - \sum\limits_{k} \sum\limits_{j<k} {\int}_{{\Sigma}_{\text{kj}}} (M_{\text{kj}}^{-}\mathbf{v}_{k} + M_{kj}^{+} \mathbf{v}_{j})^{T} (\mathbf{u}_{k} - \mathbf{u}_{j}) -\!\! \sum\limits_{k} {\int}_{{\Sigma}_{\text{kk}}} \mathbf{v}_{k}^{T} M_{k}^{-} \mathbf{u}_{k} ,\!\!\!\!\!\! \quad \mathbf{u},\mathbf{v} \in\! V(\mathcal{T}_{h}). $$

(40)

The main advantage to the formulation (40) compare to (39) is that there is no volume term; thus, it can save some computational time and may be easier to code. Note also that in the formulation (40), the relaxation term R completely disappears. This is not a problem for the TDG method because the information about R will be contained in the basis functions. For the bilinear form $l:V(\mathcal {T}_{h}) \rightarrow \mathbb {R}$, it is the same as in (36), that is, $l(\mathbf {v}) = -{\sum }_{k} {\int \limits }_{{\Sigma }_{kk}} \mathbf {v}_{k}^{T} M_{k}^{-} \mathbf {g}$ for all $\mathbf {v} \in V(\mathcal {T}_{h})$.

Definition B.2

Assume $V_{m}(\mathcal {T}_{h})$ is a finite subspace of $V(\mathcal {T}_{h})$. The TDG method is as follows:

$$ \left\{\begin{array}{ll} \text{find } \mathbf{u}_{h} \in V_{m}(\mathcal{T}_{h}) \text{ such that }\\ a_{T}(\mathbf{u}_{h},\mathbf{v}_{h}) = l(\mathbf{v}_{h}), \quad \forall \mathbf{v}_{h} \in V_{m}(\mathcal{T}_{h}). \end{array}\right. $$

(41)