Approximately well-balanced Discontinuous Galerkin methods using bases
  enriched with Physics-Informed Neural Networks
This work concerns the enrichment of Discontinuous Galerkin (DG) bases, so
that the resulting scheme provides a much better approximation of steady
solutions to hyperbolic systems of balance laws. The basis enrichment leverages
a prior - an approximation of the steady solution - which we propose to compute
using a Physics-Informed Neural Network (PINN). To that end, after presenting
the classical DG scheme, we show how to enrich its basis with a prior.
Convergence results and error estimates follow, in which we prove that the
basis with prior does not change the order of convergence, and that the error
constant is improved. To construct the prior, we elect to use parametric PINNs,
which we introduce, as well as the algorithms to construct a prior from PINNs.
We finally perform several validation experiments on four different hyperbolic
balance laws to highlight the properties of the scheme. Namely, we show that
the DG scheme with prior is much more accurate on steady solutions than the DG
scheme without prior, while retaining the same approximation quality on
unsteady solutions.
Emmanuel Franck, Victor Michel-Dansac, Laurent Navoret
Approximately well-balanced Discontinuous Galerkin methods using bases enriched with Physics-Informed Neural Networks

###### Abstract

This work concerns the enrichment of Discontinuous Galerkin (DG) bases, so that the resulting scheme provides a much better approximation of steady solutions to hyperbolic systems of balance laws. The basis enrichment leverages a prior - an approximation of the steady solution - which we propose to compute using a Physics-Informed Neural Network (PINN). To that end, after presenting the classical DG scheme, we show how to enrich its basis with a prior. Convergence results and error estimates follow, in which we prove that the basis with prior does not change the order of convergence, and that the error constant is improved. To construct the prior, we elect to use parametric PINNs, which we introduce, as well as the algorithms to construct a prior from PINNs. We finally perform several validation experiments on four different hyperbolic balance laws to highlight the properties of the scheme. Namely, we show that the DG scheme with prior is much more accurate on steady solutions than the DG scheme without prior, while retaining the same approximation quality on unsteady solutions.

## 1 Objective and model

In the last decades, much work has been devoted to proposing numerical methods for hyperbolic systems with source terms, which correctly capture stationary solutions of the system, as well as perturbations of flows around these steady states. If the perturbation is smaller than the scheme error, traditional numerical schemes are not able to provide a good approximation of the perturbed steady solution. To address such an issue, a first possibility is to refine the mesh in space. However, for small perturbations, this would greatly increase the computational overhead. To avoid this, schemes specifically dedicated to capturing stationary solutions have been introduced. They are called _well-balanced schemes_.

There are two families of well-balanced (WB) schemes: exactly and approximately WB schemes. Exactly WB schemes give an exact representation of the equilibria. Such schemes are usually developed for subclasses of steady solutions, especially for complex balance laws, or multidimensional problems. For instance, first- and second-order accurate exactly WB schemes have been developed for the shallow water equations [4, 31, 34, 35] or the Euler equations with gravity [28, 49]. High-order exactly well-balanced schemes include [22, 37, 21, 36, 8, 9, 24] with finite volume methods or related approaches, or [53, 11] with discontinuous Galerkin methods. The second family, approximately WB schemes, consist in ensuring a better approximation of the equilibria compared to traditional numerical schemes. This better approximation can be under the form of a better order of accuracy [18, 20, 23] or a better error constant [1, 16]. Both families of WB schemes may incur significant additional computational cost compared to traditional schemes, due to the extensive modifications necessary to ensure the WB property, especially for complex systems and equilibria.

In this work, we focus on providing an approximately WB scheme for the following parametric partial differential equation (PDE):

\[\begin{cases}\partial_{t}u+\partial_{x}F_{\mu_{1}}(u)=S_{\mu_{2}}(u),\\ u(t=0,x)=u_{0}(x),\end{cases} \tag{1.1}\]
with \(\mu_{1}\) and \(\mu_{2}\) the parameters of the PDE. We set \(\mu=\{\mu_{1},\mu_{2}\}\), and we assume that \(\mu\in\mathbb{P}\subset\mathbb{R}^{m}\). In (1.1), the unknown function is \(u\); \(F_{\mu_{1}}\) is called the physical flux function, while \(S_{\mu_{2}}\) is the source term. We assume that the equation is hyperbolic, that is to say that the Jacobian matrix of \(F_{\mu_{1}}\) is diagonalizable with real eigenvalues. Our goal will be to construct an approximately well-balanced approach for the general steady state \(\partial_{x}F_{\mu_{1}}(u)=S_{\mu_{2}}(u)\). The combination of learning and numerical methods (known as Scientific Machine Learning) has produced good results for hyperbolic PDEs. Examples include work on the design of limiters or shock detection [41, 7, 55], artificial viscosity [19, 42, 54, 10], or numerical schemes [5].

The approach proposed in this paper is also based on the hybridization of classical approaches and neural networks. We endeavor to improve the classical Discontinuous Galerkin (DG) method, which usually relies on a discontinuous approximation of the solution in a suitable polynomial basis. More information on the DG method can be found in [26, 39] for instance. A natural way of improving the traditional DG method to improve the accuracy on some family of solutions is to enrich the basis with a prior. This is for example the case of the Trefftz method [30, 6, 12, 27], or the non-polynomial bases studied in [56].

To perform this basis enrichment, we use a learning-based offline computation with a neural network to build a prior which approximates a parametrized family of equilibria. To that end, we use Physics-Informed Neural Networks (PINNs), see e.g. [40, 13], and parametric neural networks [48]. This prior is then introduced into the Discontinuous Galerkin basis, to increase the accuracy of the scheme around this family of equilibria. Note that the prior construction could be handled without the use of neural networks, but we will show that the neural network approach is more efficient. This framework could require significant offline calculation cost (depending on the problem), but will generate a very small additional cost in the online phase, i.e., when actually using the modified scheme. This method enhances the DG basis functions with a prior provided by a neural network. Similar techniques based on neural networks have already been successfully implemented for other applications. In [3], for elliptic problems, the authors use a network to provide a finite element basis that is dynamically enriched. In [46, 47], the authors show that random neural networks can be used as DG bases, and can be more accurate than classical ones for a sufficiently large number of basis function in each cell.

The paper is constructed as following. First, we assume that we know a prior (an approximation) of a family of equilibria, and we introduce the modification of the DG basis. Theoretical results show that this modification does not change the order of accuracy of the method, but decreases the error constant close to steady solutions. Then, we introduce the learning methods that will enable us to build our prior for a family of equilibria, and finally we perform numerical experiments, in one and two space dimensions, on several linear and nonlinear systems of balance laws. A conclusion ends this paper.

## 2 Modified Discontinuous Galerkin scheme

This section is devoted to the presentation of the modified DG scheme. We start by quickly introducing the classical DG scheme in Section2.1, and then move on to proposing the modification in Section2.2. Theoretical convergence results related to this modification will be presented in Section3. In this section and the following one, we write the scheme in the case of a scalar and one-dimensional PDE, but the method is easily extendable to systems and to higher dimensions.

### Classical Discontinuous Galerkin scheme

The goal of this section is to present the classical DG scheme in order to discretize the PDE (1.1). To that end, we discretize the space domain \(\Omega\subset\mathbb{R}^{d}\) in cells \(\Omega_{k}=(x_{k-1/2},x_{k+1/2})\) of size \(\Delta x_{k}\), and of centers \(x_{k}\).

The idea behind the classical DG scheme is to first compute the weak form of the considered PDE, and then to locally approximate the solution in each cell, by projecting it onto a finite-dimensional vector space \(V_{h}\). We consider a space \(V_{h}\) of dimension \(q+1\):

\[V_{h}=\operatorname{Span}\left(\phi_{k,0},\ldots,\phi_{k,q}\right).\]

Note that the space \(V_{h}\) can be different for each cell \(k\).

The first assumption of DG scheme is to approximate the solution \(u\) to the PDE, in each cell, with a value in \(V_{h}\):

\[\forall k,\quad u\big{|}_{\Omega_{k}}(t,x)\simeq u_{k}(t,x).\]
Since \(u_{k}\in V_{h}\), we can write

\[u_{k}(t,x)\coloneqq\sum_{j=0}^{q}u_{k,j}(t)\phi_{k,j}(x). \tag{2.1}\]

To obtain the DG scheme, we first write the weak form of the equation in each cell:

\[\int_{\Omega_{k}}\partial_{t}u(t,x)\phi(x)\,dx+\int_{\Omega_{k}}\partial_{x}F_{ \mu_{1}}(u(t,x))\phi(x)\,dx=\int_{\Omega_{k}}S_{\mu_{2}}(u(t,x))\phi(x)\,dx. \tag{2.2}\]

with \(\phi(x)\) a smooth test function. Performing an integration by parts, this form is equivalent to

\[\partial_{t}\left(\int_{\Omega_{k}}u\,\phi\right)-\int_{\Omega_{k}}F_{\mu_{1} }(u)\,\partial_{x}\phi+\left[F_{\mu_{1}}(u)\,\phi\right]^{x_{k+1/2}}_{x_{k-1/2 }}=\int_{\Omega_{k}}S_{\mu_{2}}(u)\,\phi. \tag{2.3}\]

We now plug the DG representation (2.1) in the weak form (2.3), using \(\phi_{k,i}\) as test function, for any \(i\in\{0,\ldots,q\}\).

1. We begin with the first term: \[\int_{\Omega_{k}}u\,\phi_{k,i}=\sum_{j=0}^{q}\left(\int_{\Omega_{k}}u_{k,j}(t )\phi_{k,j}(x)\phi_{k,i}(x)\,dx\right)=\sum_{j=0}^{q}u_{k,j}(t)\left(\int_{ \Omega_{k}}\phi_{k,j}(x)\phi_{k,i}(x)\,dx\right).\] To handle the integral in the expression above, we introduce the following quadrature formula, with weights \(w_{k,p}\) and points \(x_{k,p}\), valid for any smooth function \(\phi\): \[\int_{\Omega_{k}}\phi(x)\,dx\simeq\sum_{p=1}^{N_{q}}w_{k,p}\,\phi(x_{k,p}).\] We assume that the first and last quadrature points coincide with the cell boundaries, i.e. \(x_{k,1}=x_{k-1/2}\) and \(x_{k,N_{q}}=x_{k+1/2}\). In practice, we use the well-known Gauss-Lobatto quadrature rule, see e.g. [2] for more information. Equipped with this quadrature formula, we introduce \[M_{k,i,j}=\sum_{p=1}^{N_{q}}w_{k,p}\,\phi_{k,j}(x_{k,p})\,\phi_{k,i}(x_{k,p}) \simeq\int_{\Omega_{k}}\phi_{k,j}\phi_{k,i},\] so that the first term of (2.3) becomes \[\int_{\Omega_{k}}u(t,x)\phi_{k,i}(x)\,dx\simeq\sum_{j=0}^{q}M_{k,i,j}u_{k,j}(t).\]
2. Using the same techniques, the second term is approximated in the following way: \[\int_{\Omega_{k}}F_{\mu_{1}}(u)\,\partial_{x}\phi_{k,i}\simeq\sum_{p=1}^{N_{ q}}\left[w_{k,p}F_{\mu_{1}}\left(\sum_{j=0}^{q}u_{k,j}(t)\phi_{k,j}(x_{k,p}) \right)\partial_{x}\phi_{k,i}(x_{k,p})\right].\]
3. We note that the third term reduces to \[\left[F_{\mu_{1}}(u)\,\phi_{k,i}\right]^{x_{k+1/2}}_{x_{k-1/2}}=F_{\mu_{1}} \left(u_{k}(t,x_{k+\frac{1}{2}})\right)\phi_{k,i}(x_{k+\frac{1}{2}})-F_{\mu_{ 1}}\left(u_{k}(t,x_{k-\frac{1}{2}})\right)\phi_{k,i}(x_{k-\frac{1}{2}}),\] where the physical flux \(F_{\mu_{1}}\) has to be approximated at the cell boundaries. To that end, like the well-known finite volumes method, the DG method requires the introduction of a consistent numerical flux \[G_{\mu_{1}}(u_{L},u_{R})\quad\text{ such that }G_{\mu_{1}}(u,u)=F_{\mu_{1}}(u).\] This numerical flux is then used to approximate the interface flux, as follows \[F_{\mu_{1}}\left(u_{k}(t,x_{k+\frac{1}{2}})\right)\simeq G_{\mu_{1}}\left(u_{ k}(t,x_{k+\frac{1}{2}}),u_{k+1}(t,x_{k+\frac{1}{2}})\right).\]
4. Finally, for the last term, we use a straightforward application of the quadrature rule: \[\int_{\Omega_{k}}S_{\mu_{2}}(u)\,\phi\simeq\sum_{p=1}^{N_{q}}\left[w_{k,p}S_{\mu_{ 2}}\left(\sum_{j=0}^{q}u_{k,j}(t)\phi_{k,j}(x_{k,p})\right)\phi_{k,i}(x_{k,p}) \right].\]

Gathering all these terms, we show that, in each cell, the DG scheme can be written as an ordinary differential equation, where the interface flux term couples the cell \(\Omega_{k}\) with its neighbors:

\[\mathcal{M}_{k}\,\partial_{t}u_{k}(t)-\mathcal{F}_{\mu_{1}}(u_{k})+\mathcal{G }_{\mu_{1}}(u_{k-1},u_{k},u_{k+1})=\mathcal{S}_{\mu_{2}}(u_{k}).\]

Now that we have recalled the classical DG space discretization, we have all the tools we need to introduce a modification to this discretization that will enable us to provide an approximately WB scheme.

### Enrichment of the modal DG basis

There are many vector spaces able to represent the solution in each cell. For instance, nodal DG schemes [26] use Lagrange polynomials or other polynomials based on nodes chosen within each cell. Legendre polynomials or Taylor expansions around the cell centers lead to _modal DG schemes_. In this work, we focus on the Taylor basis, given on each cell \(\Omega_{k}\) by

\[V_{h}=\mathrm{Span}\left(\phi_{k,0},\phi_{k,1},\phi_{k,2},\ldots,\phi_{k,q} \right)=\mathrm{Span}\left(1,(x-x_{k}),\frac{1}{2}(x-x_{k})^{2},\ldots,\frac {1}{q!}(x-x_{k})^{q}\right). \tag{2.4}\]

In the remainder of this section, we assume that we have access to a prior on the equilibrium, denoted by \(u_{\theta}(x,\mu)\). Obtaining such a prior is discussed in Section4. For the moment, suffice it to say that \(u_{\theta}\) provides an approximation of the steady solution for \(x\in\Omega\) and for \(\mu\) in some parameter space \(\mathbb{P}\) to be defined.

Given the prior \(u_{\theta}\), we modify the local basis \(V_{h}\) to incorporate the prior: for that, we propose two possibilities.

* The _additive correction_\(V_{h}^{+}\) consists in replacing the first element of \(V_{h}\) by the prior: \[V_{h}^{+}=\mathrm{Span}\left(\phi_{k,0}^{+},\phi_{k,1}^{+},\phi_{k,2}^{+}, \ldots,\phi_{k,q}^{+}\right)=\mathrm{Span}\left(u_{\theta}(x,\mu),(x-x_{k}),\ldots,\frac{1}{q!}(x-x_{k})^{q}\right).\] (2.5)
* The _multiplicative correction_\(V_{h}^{*}\) consists in multiplying each element of \(V_{h}\) by the prior: \[V_{h}^{*}=\mathrm{Span}\left(\phi_{k,0}^{*},\phi_{k,1}^{*},\phi_{k,2}^{*}, \ldots,\phi_{k,q}^{*}\right)=\mathrm{Span}\left(u_{\theta}(x,\mu),(x-x_{k}) \,u_{\theta}(x,\mu),\ldots,\frac{1}{q!}(x-x_{k})^{q}u_{\theta}(x,\mu)\right).\] (2.6)

A first remark is that, if the prior is exactly equal to the steady solution, then it can be exactly represented by an element of \(V_{h}^{+}\) or \(V_{h}^{*}\) (namely, the first one) in each cell, which is not the case for the classical space \(V_{h}\). However, whether the prior is exact or not, the method will only be of interest if the projector onto the modified vector space is accurate (or even exact in the case of an exact prior). The second point to note is that, unlike conventional DG approaches, the bases are not polynomial. We must therefore ensure that this does not hinder the convergence of the DG method. In the next section, we follow Yuan and Shu's work [56] to study the convergence of the modified DG method, and provide error estimates.

## 3 Error estimates

In this section, we prove some convergence results on the modified DG scheme. We assume that our prior \(u_{\theta}\) is \(p\) times continuously differentiable, i.e., that it has differentiability class \(\mathcal{C}^{p}\), with \(p\geqslant q+1\). This hypothesis is compatible with the construction of the prior from Section4.

In [56], the authors study the convergence of the DG scheme for non-polynomial bases. They show that, if the non-polynomial basis can be represented in a specific way by a polynomial basis, then the convergence of the
local and global projection operators is not hampered. Using some stability results (given in [56] for the transport equation) together with these estimations, convergence can be recovered.

These theoretical results will be split in two parts. To begin with, in Section3.1, by prove that the bases proposed in Section2.2 fit into the hypotheses of [56], which ensures convergence. However, this study is insufficient to show that the better the prior, the more accurate the modified DG scheme. To that end, in Section3.2, we derive the projector estimates in the case of \(V_{h}^{*}\), in order to show the potential gains of the method.

### Convergence in non-polynomial DG bases

In [56], the authors prove the following lemma.

**Lemma 3.1**.: _Consider an approximation vector space \(V_{h}\) with local basis \((v_{k,0},\ldots,v_{k,q})\), which may depend on the cell \(\Omega_{k}\). If there exists constant real numbers \(a_{j\ell}\) and \(b_{j}\) independent of the size of the cell \(\Delta x_{k}\) such that, in each cell \(\Omega_{k}\),_

\[\forall j\in\{0,\ldots,q\},\quad\left|v_{k,j}(x)-\sum_{\ell=0}^{q}a_{j\ell}(x -x_{k})^{\ell}\right|\leq b_{j}(\Delta x_{k})^{q+1}, \tag{3.1}\]

_then for any function \(u\in H^{q+1}(\Omega_{k})\), there exists \(v_{h}\in V_{h}\) and a constant real number \(C\) independent of \(\Delta x_{k}\), such that_

\[\|v_{h}-u\|_{L^{\infty}(\Omega_{k})}\leq C\|u\|_{H^{q+1}(\Omega_{k})}(\Delta x _{k})^{q+\frac{1}{2}}.\]

Using this result, the authors show that the global projection error in the DG basis converges with an error in \((\Delta x_{k})^{q+1}\) in the Sobolev norm \(H^{q+1}\), and later prove the convergence of the whole scheme using a monotone flux for a scalar equation. In the remainder of this section, we prove that the two new bases proposed in Section2.2 satisfy the assumptions of Lemma3.1. Using these results together with the proofs of [56], we will obtain that both bases lead to a convergent scheme.

**Proposition 3.2**.: _If the prior \(u_{\theta}(x;\mu)\) has differentiability class \(\mathcal{C}^{q+1}(\mathbb{R})\) with respect to \(x\), then the approximation space \(V_{h}^{+}\) satisfies the assumption (3.1)._

Proof.: Since the prior is \(\mathcal{C}^{q+1}(\mathbb{R})\), we can write its Taylor series expansion around the cell center \(x_{k}\). Namely, there exists a constant \(c\in[x_{k-1/2},x_{k+1/2}]\) such that

\[u_{\theta}(x)=u_{\theta}(x_{k})+(x-x_{k})u_{\theta}^{\prime}(x_{k})+\cdots+ \frac{1}{q!}(x-x_{k})^{q}u^{(q)}(x_{k})+\frac{(x-x_{k})^{q+1}}{(q+1)!}u^{(q+1 )}(c). \tag{3.2}\]

With that expansion, we can write our basis \(V_{h}^{+}\) with respect to the classical modal basis \(V_{h}\) as follows:

\[\begin{pmatrix}u_{\theta}(x)\\ (x-x_{k})\\ \vdots\\ (x-x_{k})^{q}\end{pmatrix}=\underbrace{\begin{pmatrix}u_{\theta}(x_{k})&u_{ \theta}^{\prime}(x_{k})&\ldots&\frac{1}{q!}u_{\theta}^{(q)}(x_{k})\\ 0&1&\ldots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\ldots&1\end{pmatrix}}\begin{pmatrix}1\\ (x-x_{k})\\ \vdots\\ (x-x_{k})^{q}\end{pmatrix}+(x-x_{k})^{q+1}\underbrace{\begin{pmatrix}\frac{u^{ (q+1)}(c)}{(q+1)!}\\ 0\\ \vdots\\ 0\end{pmatrix}}_{b_{+}}.\]

We remark that the matrix \(A_{+}\) and the vector \(b_{+}\) are independent of \(\Delta x_{k}\). Hence, assumption (3.1) is verified, and Lemma3.1 can be applied. 

**Proposition 3.3**.: _If the prior \(u_{\theta}(x;\mu)\) has differentiability class \(\mathcal{C}^{q+1}(\mathbb{R})\) with respect to \(x\), then the approximation space \(V_{h}^{*}\) satisfies the assumption (3.1)._
Proof.: The proof follows the same lines as the proof of the previous proposition. Namely, (3.2) is still satisfied since the prior is \(\mathcal{C}^{q+1}(\mathbb{R})\). Then, the basis \(V_{h}^{*}\) is written with respect to the classical modal basis \(V_{h}\) as follows:

\[\begin{pmatrix}u_{\theta}(x)\\ (x-x_{k})\,u_{\theta}(x)\\ \vdots\\ (x-x_{k})^{q}\,u_{\theta}(x)\end{pmatrix}=\underbrace{\begin{pmatrix}u_{ \theta}(x_{k})&u_{\theta}^{\prime}(x_{k})&\dots&\frac{u^{(q)}(x_{k})}{q!}\\ 0&u_{\theta}(x_{k})&\dots&\frac{u_{\theta}^{(q-1)}(x_{k})}{(q-1)!}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&u_{\theta}(x_{k})\end{pmatrix}}_{\mathcal{A}_{*}}\begin{pmatrix}1 \\ (x-x_{k})\\ (x-x_{k})^{q}\end{pmatrix}+(x-x_{k})^{q+1}\underbrace{\begin{pmatrix}\frac{u_ {\theta}^{q+1}(c)}{q!}\\ \vdots\\ 1\end{pmatrix}}_{b_{*}}\]

Just like before, the matrix \(A_{*}\) and the vector \(b_{*}\) are independent of \(\Delta x_{k}\). Hence, assumption (3.1) is verified, and Lemma3.1 can be applied. 

These two propositions show that, if the prior is sufficiently smooth, we can apply the results of [56], which shows the convergence of the method. However, this approach does not give an estimation of the error with respect to the quality of the prior. Indeed, we expect the modified DG scheme to be more accurate when the prior is closer to the solution. Obtaining such an estimate is the objective of the following section.

### Estimate with prior dependency

The goal of this section is to refine the error estimates from Section3.1 for a specific modified basis. We consider the case of \(V_{h}^{*}\), since it is easier to write the projector onto the classical basis. This will enable us to quantify the gains that can be expected when using this new basis. The case of \(V_{h}^{+}\) is more complicated, since the projector is harder to write. Nevertheless, we will show in the numerical experiments from Section5 that both modified bases exhibit similar behavior.

Recall that the basis \(V_{h}^{*}\) is obtained by multiplying each element of \(V_{h}\) by the prior. Therefore, its basis functions are given by \(\phi_{k,j}^{*}=\phi_{k,j}u_{\theta}\) for each cell \(\Omega_{k}\) and for \(j\in\{0,\dots,q\}\).

**Lemma 3.4**.: _Assume that the prior \(u_{\theta}\) satisfies_

\[u_{\theta}(x;\mu)^{2}>m^{2}>0,\quad\forall x\in\Omega,\quad\forall\mu\in \mathbb{P}.\]

_For a given cell \(\Omega_{k}\), for any function \(u\in H^{q+1}(\Omega_{k})\), the \(L^{2}\) projector onto \(V_{h}^{*}\), denoted by \(P_{h}\) and such that \(P_{h}(u)\in V_{h}^{*}\), satisfies the inequality_

Proof.: The proof uses a strategy similar to [56]. We consider the cell \(\Omega_{k}\). For any smooth function \(f\) defined on \(\Omega_{k}\), we define the operator \(T\) by

\[T(f)=\left(\sum_{j=0}^{q}f^{(j)}(x_{k})\frac{1}{j!}(x-x_{k})^{j}\right)\]

and the operator \(T_{\theta}\) by

\[T_{\theta}(f)=\left(\sum_{j=0}^{q}\left(\frac{f}{u_{\theta}}\right)^{(j)}(x_{k };\mu)\frac{1}{j!}(x-x_{k})^{j}\right)u_{\theta}(x;\mu).\]

For simplicity, we no longer explicitly write the dependence in \(\mu\) in this proof. Let \(u\in H^{q+1}(\Omega_{k})\). Using \(T_{\theta}\), we write the following estimation:

\[\|u-P_{h}(u)\|_{L^{\infty}(\Omega_{k})}\leq\|u-T_{\theta}(u)\|_{L^{\infty}( \Omega_{k})}+\|T_{\theta}(u)-P_{h}(u)\|_{L^{\infty}(\Omega_{k})}\eqqcolon N _{1}+N_{2}. \tag{3.3}\]

To complete the proof, we need to estimate both terms \(N_{1}\) and \(N_{2}\).
We start with the estimation of \(N_{1}\). We obtain, according to the relationship between \(T\) and \(T_{\theta}\),

\[N_{1}=\|u-T_{\theta}(u)\|_{L^{\infty}(\Omega_{k})}=\left\|\frac{u}{u_{\theta}}u_{ \theta}-T\left(\frac{u}{u_{\theta}}\right)u_{\theta}\right\|_{L^{\infty}( \Omega_{k})}\leq\left\|\frac{u}{u_{\theta}}-T\left(\frac{u}{u_{\theta}}\right) \right\|_{L^{\infty}(\Omega_{k})}\left\|u_{\theta}\right\|_{L^{\infty}(\Omega_ {k})}. \tag{3.4}\]

We can now use an intermediate result from [56]: for all \(f\) smooth enough, the Taylor formula and the Cauchy-Schwartz inequality, followed by a direct computation, gives

\[\|f-T(f)\|_{L^{\infty}(\Omega_{k})} =\sup_{x\in\Omega_{k}}\left|\int_{x_{k}}^{x}f^{(q+1)}(\xi)\frac{( x-\xi)^{q}}{q!}d\xi\right|\] \[\leq\sup_{x\in\Omega_{k}}\left[\left(\int_{x_{k}}^{x}\left|f^{(q+ 1)}(\xi)\right|^{2}d\xi\right)^{\frac{1}{2}}\left(\int_{x_{k}}^{x}\left|\frac {(x-\xi)^{q}}{q!}\right|^{2}d\xi\right)^{\frac{1}{2}}\right],\] \[\lesssim|f|_{H^{q+1}(\Omega_{k})}\left(\Delta x_{k}\right)^{q+ \frac{1}{2}}. \tag{3.5}\]

Going back to \(N_{1}\) and plugging (3.5) into the estimate (3.4), we obtain

\[N_{1}\lesssim\left|\frac{u}{u_{\theta}}\right|_{H^{q+1}(\Omega_{k})}\left( \Delta x_{k}\right)^{q+\frac{1}{2}}\left\|u_{\theta}\right\|_{L^{\infty}( \Omega_{k})} \tag{3.6}\]

Now, we proceed with estimating \(N_{2}\), the second term of (3.3). The \(L^{2}\) projector \(P_{h}\) onto \(V_{h}^{*}\) is defined by

\[P_{h}(u)=\sum_{j=0}^{q}\frac{\alpha_{j}}{(\Delta x_{k})^{j}}\phi_{j}^{*},\]

with \(\alpha=(\alpha_{j})_{j\in\{1,\dots,q\}}=(M^{*})^{-1}b\), where

\[M_{j\ell}^{*}=\int_{\Omega_{k}}\frac{\phi_{j}^{*}(x)}{(\Delta x_{k})^{j}} \frac{\phi_{\ell}^{*}(x)}{(\Delta x_{k})^{\ell}}dx,\quad\text{and}\quad b_{j} =\int_{\Omega_{k}}u(x)\frac{\phi_{j}^{*}(x)}{(\Delta x_{k})^{j}}dx. \tag{3.7}\]

We are now ready to start estimating \(N_{2}\). Note that

\[N_{2}=\|T_{\theta}(u)-P_{h}(u)\|_{L^{\infty}(\Omega_{k})} =\sup_{x\in\Omega_{k}}\left|\sum_{j=0}^{q}\left(\left(\frac{u}{u_ {\theta}}\right)^{(j)}(x_{k})-\frac{\alpha_{j}}{(\Delta x_{k})^{j}}\right) \frac{1}{j!}(x-x_{k})^{j}u_{\theta}(x)\right|\] \[\leq\sup_{x\in\Omega_{k}}\left|\sum_{j=0}^{q}\left((\Delta x_{k}) ^{j}\left(\frac{u}{u_{\theta}}\right)^{(j)}(x_{k})-\alpha_{j}\right)\frac{1}{j! }\frac{(x-x_{k})^{j}}{(\Delta x_{k})^{j}}\right|\left\|u_{\theta}\right\|_{L^{ \infty}(\Omega_{k})}.\]

Using the Cauchy-Schwartz inequality on the sum, and bounding the resulting polynomial on the cell, we obtain the estimate

\[N_{2}\lesssim\left[\sum_{j=0}^{q}\left((\Delta x_{k})^{j}\left(\frac{u}{u_{ \theta}}\right)^{(j)}(x_{k})-\alpha_{j}\right)^{2}\right]^{\frac{1}{2}}\left\|u _{\theta}\right\|_{L^{\infty}(\Omega_{k})}=\left\|\delta-\alpha\right\|_{2} \left\|u_{\theta}\right\|_{L^{\infty}(\Omega_{k})} \tag{3.8}\]

where the vector \(\delta=(\delta_{j})_{j\in\{0,\dots,q\}}\) is defined by

\[\delta_{j}=(\Delta x_{k})^{j}\left(\frac{u}{u_{\theta}}\right)^{(j)}(x_{k}).\]

Recalling the definition \(\alpha=(M^{*})^{-1}b\), we obtain

\[\left\|\delta-\alpha\right\|_{2}=\left\|(M^{*})^{-1}(M^{*}\delta-b)\right\|_{2 }\leq\left\|(M^{*})^{-1}\right\|_{2}\left\|M^{*}\delta-b\right\|_{2}. \tag{3.9}\]
We first take care of the term in \(M^{*}\delta-b\). We have

\[\left\|M^{*}\delta-b\right\|_{2}^{2} =\sum_{j=0}^{q}\left[\sum_{\ell=1}^{q}M^{*}_{j\ell}\delta_{\ell}-b_ {j}\right]^{2}\] \[=\sum_{j=0}^{q}\left[\sum_{\ell=0}^{q}\left(\int_{\Omega_{k}} \frac{\phi_{j}^{*}(x)}{(\Delta x_{k})^{j}}\frac{\phi_{\ell}^{*}(x)}{(\Delta x_ {k})^{\ell}}dx\right)(\Delta x_{k})^{\ell}\left(\frac{u}{u_{\theta}}\right)^{( \ell)}(x_{k})-\int_{\Omega_{k}}u(x)\frac{\phi_{j}^{*}(x)}{(\Delta x_{k})^{j}} dx\right]^{2}=:\sum_{j=0}^{q}\Xi_{j}^{2}\]

We denote the summand by \(\Xi_{j}\), and we use the definition of the basis to obtain

\[\forall j\in\{0,\dots,q\},\quad\Xi_{j} :=\sum_{\ell=0}^{q}(\Delta x_{k})^{\ell}\left(\frac{u}{u_{\theta} }\right)^{(\ell)}(x_{k})\int_{\Omega_{k}}\frac{\phi_{j}^{*}(x)}{(\Delta x_{k })^{j}}\frac{\phi_{\ell}^{*}(x)}{(\Delta x_{k})^{\ell}}dx-\int_{\Omega_{k}}u(x )\frac{\phi_{j}^{*}(x)}{(\Delta x_{k})^{j}}dx\] \[=\sum_{\ell=0}^{q}(\Delta x_{k})^{\ell}\left(\frac{u}{u_{\theta} }\right)^{(\ell)}(x_{k})\int_{\Omega_{k}}\frac{\phi_{j}(x)}{(\Delta x_{k})^{j }}\frac{\phi_{\ell}(x)}{(\Delta x_{k})^{\ell}}u_{\theta}^{2}(x)dx-\int_{\Omega _{k}}\frac{u(x)}{u_{\theta}(x)}\frac{\phi_{j}(x)}{(\Delta x_{k})^{j}}u_{\theta }^{2}(x)dx\] \[=\int_{\Omega_{k}}\left(\sum_{\ell=0}^{q}\left(\frac{u}{u_{\theta }}\right)^{(\ell)}(x_{k})\phi_{\ell}(x)-\frac{u(x)}{u_{\theta}(x)}\right)\frac {\phi_{j}(x)}{(\Delta x_{k})^{j}}u_{\theta}^{2}(x)dx. \tag{3.10}\]

Using a Taylor expansion, we obtain, for all \(j\in\{1,\dots,q\}\),

\[\Xi_{j}= -\int_{\Omega_{k}}\left(\int_{x_{k}}^{x}\left(\frac{u}{u_{\theta} }\right)^{(q+1)}(\xi)\,\frac{(x-\xi)^{q}}{q!}d\xi\right)\frac{\phi_{j}(x)}{( \Delta x_{k})^{j}}u_{\theta}^{2}(x)dx,\]

from which we get the following upper bound

\[\forall j\in\{1,\dots,q\},\quad\left|\Xi_{j}\right|\leq\sup_{x\in\Omega_{k}} \left|\int_{x_{k}}^{x}\left(\frac{u}{u_{\theta}}\right)^{(q+1)}(\xi)\,\frac{(x -\xi)^{q}}{q!}d\xi\right|\left|\int_{\Omega_{k}}\frac{\phi_{j}(x)}{(\Delta x_{ k})^{j}}u_{\theta}^{2}(x)dx\right|.\]

Using the same ingredients as in the computation of (3.5) for the leftmost term and bounding the rightmost term by the \(L^{\infty}\) norm of the prior and by noting that the classical basis functions are bounded, we obtain the estimate

\[\forall j\in\{1,\dots,q\},\quad\left|\Xi_{j}\right|\lesssim\left|\frac{u}{u_ {\theta}}\right|_{H^{q+1}(\Omega_{k})}\left(\Delta x_{k}\right)^{q+\frac{1}{2}} \left(\Delta x_{k}\right)\left\|u_{\theta}^{2}\right\|_{L^{\infty}(\Omega_{k})}.\]

Going back to what we had set out to prove, we get

\[\left\|M^{*}\delta-b\right\|_{2}=\left(\sum_{j=0}^{q}\left|\Xi_{j}\right|^{2} \right)^{\frac{1}{2}}\lesssim\left|\frac{u}{u_{\theta}}\right|_{H^{q+1}(\Omega _{k})}\left(\Delta x_{k}\right)^{q+\frac{1}{2}}\left(\Delta x_{k}\right)\left\| u_{\theta}^{2}\right\|_{L^{\infty}(\Omega_{k})}. \tag{3.11}\]

Plugging (3.11) into (3.9) and then into (3.8), we get

\[N_{2}\lesssim\left\|(M^{*})^{-1}\right\|_{2}\left|\frac{u}{u_{\theta}}\right| _{H^{q+1}(\Omega_{k})}\left(\Delta x_{k}\right)^{q+\frac{1}{2}}\left(\Delta x_ {k}\right)\left\|u_{\theta}^{2}\right\|_{L^{\infty}(\Omega_{k})}\left\|u_{ \theta}\right\|_{L^{\infty}(\Omega_{k})}. \tag{3.12}\]

Finally, we note that, for any \(y\in\mathbb{R}^{q+1}\), given the expression (3.7) of \(M^{*}\),

\[(M^{*}y,M^{*}y)=\int_{\Omega_{k}}\left(\sum_{j=0}^{q}\frac{\phi_{j}^ {*}(x)}{(\Delta x_{k})^{j}}y_{j}\right)^{2}dx =\int_{\Omega_{k}}\left(\sum_{j=0}^{q}\frac{\phi_{j}(x)}{(\Delta x _{k})^{j}}y_{j}\right)^{2}u_{\theta}(x)^{2}dx\] \[\geqslant m^{2}\int_{\Omega_{k}}\left(\sum_{j=0}^{q}\frac{\phi_{j}( x)}{(\Delta x_{k})^{j}}y_{j}\right)^{2}dx=m^{2}\left(My,My\right),\]
where \(M\) is the mass matrix associated with the classical basis functions

\[M_{j\ell}=\int_{\Omega_{k}}\frac{\phi_{j}(x)}{(\Delta x_{k})^{j}}\frac{\phi_{\ell }(x)}{(\Delta x_{k})^{\ell}}dx=\frac{\Delta x_{k}}{1+j+\ell}=\Delta x_{k}H_{j \ell},\]

where \(H=(H_{j\ell})_{j\ell}\) is the Hilbert matrix. Then we deduce the following inequality

\[\|(M^{*})^{-1}\|_{2}\leqslant\frac{1}{m^{2}}\|M^{-1}\|_{2}=\frac{1}{m^{2}} \frac{1}{\Delta x_{k}}\|H^{-1}\|_{2}. \tag{3.13}\]

Combining (3.12) and (3.13), we obtain

\[N_{2}\lesssim\left|\frac{u}{u_{\theta}}\right|_{H^{q+1}(\Omega_{k})}(\Delta x _{k})^{q+\frac{1}{2}}\frac{\left\|u_{\theta}^{2}\right\|_{L^{\infty}(\Omega_{k })}}{m^{2}}\left\|u_{\theta}\right\|_{L^{\infty}(\Omega_{k})}. \tag{3.14}\]

We get, from (3.6) and (3.14), the expected result. 

The above proof relies on the smoothness of the prior. This may seem counter-intuitive in a hyperbolic context. However, since the prior will be obtained from a neural network in Section 4, this smoothness assumption becomes reasonable.

**Lemma 3.5**.: _We make the same assumptions as in the previous lemma, and still consider the vector space \(V_{h}^{*}\). For any function \(u\in H^{q+1}(\Omega)\),_

\[\left\|u-P_{h}(u)\right\|_{L^{2}(\Omega)}\lesssim\left|\frac{u}{u_{\theta}} \right|_{H^{q+1}(\Omega)}(\Delta x_{k})^{q+1}\left\|u_{\theta}\right\|_{L^{ \infty}(\Omega)}.\]

Proof.: We begin by stating the definition of the discrete \(L^{2}\) norm: by assuming that

\[\Omega=\bigcup_{k=1}^{N}\Omega_{k},\]

we obtain

\[\left\|u-P_{h}(u)\right\|_{L^{2}(\Omega)}^{2}\leqslant\sum_{k=1}^{N}\Delta x _{k}\left\|u-P_{h}(u)\right\|_{L^{\infty}(\Omega_{k})}^{2}.\]

Using the result from Lemma 3.4, we get

\[\left\|u-P_{h}(u)\right\|_{L^{2}(\Omega)}^{2} \lesssim\sum_{k=1}^{N}\Delta x_{k}\left(\left|\frac{u}{u_{\theta }}\right|_{H^{q+1}(\Omega_{k})}(\Delta x_{k})^{q+\frac{1}{2}}\,\left(1+\frac{ \left\|u_{\theta}^{2}\right\|_{L^{\infty}(\Omega_{k})}}{m^{2}}\right)\left\|u _{\theta}\right\|_{L^{\infty}(\Omega_{k})}\right)^{2}\] \[\lesssim\sum_{k=1}^{N}\left|\frac{u}{u_{\theta}}\right|_{H^{q+1} (\Omega_{k})}^{2}(\Delta x_{k})^{2q+2}\,\left(1+\frac{\left\|u_{\theta}^{2} \right\|_{L^{\infty}(\Omega_{k})}}{m^{2}}\right)^{2}\left\|u_{\theta}\right\|_ {L^{\infty}(\Omega_{k})}^{2}.\]

We assume that there exists \(\delta_{-}\), \(\delta_{+}\) and \(\Delta x\) such that, for all \(k\in\{1,\ldots,N\}\), \(\delta_{-}\Delta x\leqslant\Delta x_{k}\leqslant\delta_{+}\Delta x\). Then, since \(\left\|u_{\theta}\right\|_{L^{\infty}(\Omega_{k})}\leqslant\left\|u_{\theta} \right\|_{L^{\infty}(\Omega)}\), we obtain

\[\left\|u-P_{h}(u)\right\|_{L^{2}(\Omega)}^{2} \lesssim(\Delta x)^{2q+2}\sum_{k=1}^{N}\left|\frac{u}{u_{\theta}} \right|_{H^{q+1}(\Omega_{k})}^{2}\left(1+\frac{\left\|u_{\theta}^{2}\right\|_ {L^{\infty}(\Omega)}}{m^{2}}\right)^{2}\left\|u_{\theta}\right\|_{L^{\infty}( \Omega_{k})}^{2}\] \[\lesssim(\Delta x)^{2q+2}\left(1+\frac{\left\|u_{\theta}^{2} \right\|_{L^{\infty}(\Omega_{k})}}{m^{2}}\right)^{2}\left\|u_{\theta}\right\|_ {L^{\infty}(\Omega)}^{2}\sum_{k=1}^{N}\left|\frac{u}{u_{\theta}}\right|_{H^{q+ 1}(\Omega_{k})}^{2}.\]

The proof is concluded by recognizing the \(H^{q+1}(\Omega)\) seminorm.
This global error estimate Lemma 3.5 shows that projection error onto the basis \(V_{h}^{*}\) is bounded by

\[\left|\frac{u}{u_{\theta}}\right|_{H^{q+1}(\Omega)}.\]

This bound is equal to zero if the prior is exact, since it is nothing but the \((q+1)^{\text{th}}\) derivative of the constant function equal to one. This estimate also proves that the closer the prior is to the solution, the smaller the bound of the projection error. However, to obtain an even smaller bound, we need the prior and the solution to be close in the sense of the \(H^{q+1}(\Omega)\) seminorm. This means that the prior must be constructed in such a way that it also gives a good approximation of the derivatives of the solution.

As a summary, we have shown that the \(L^{2}\) projection error tends to zero when the prior tends to the solution. This result gives an idea of the expected gains in error ensured by using the modified basis \(V_{h}^{*}\). The final convergence error depends on this projection error, as has been shown in [56]. The proof to obtain the final convergence result is the same as in [56].

For the additive basis \(V_{h}^{+}\), such error estimates are harder to obtain, since the projection in the new basis is harder to write with respect to the traditional one. We expect an error bounded by a term in \(\|u-u_{\theta}\|_{H^{q+1}(\Omega)}\), which would enable us to draw similar conclusions as for the multiplicative basis \(V_{h}^{*}\). Namely, the error would also tend to zero when the prior tends to the solution, and the derivatives of the prior would need to be close to the derivatives of the solution. Proving this result is out of the scope of this paper, even though it should be ensured by the results of [56]. However, we will extensively study the behavior of the additive basis in Section 5.

## 4 Prior construction and algorithm

Equipped with the modified bases from Section 2.2 and with the theoretical results from Section 3, what is left to do is to propose a way to obtain a suitable prior \(u_{\theta}\).

Note that the approach described in Section 2 will be interesting if the prior \(u_{\theta}\) is a good approximation of the steady solution to (1.1) for a wide range of parameters. In addition, according to Section 3.2, the derivatives of the prior must also be good approximations of the derivatives of the steady solution.

This means that we wish to capture large families of solutions, i.e. we want to be able to calculate an approximation for several parameters. For example, assuming that (1.1) depends on 4 physical parameters leads to \(\mu\in\mathbb{R}^{4}\), and considering a problem in two space dimensions, leads to \(x\in\mathbb{R}^{2}\). Therefore, we are looking for a prior \(u_{\theta}(x;\mu)\), where \(u_{\theta}\in\mathcal{C}^{q}(\mathbb{R}^{2}\times\mathbb{R}^{4},\mathbb{R})\). Approaching such a function using polynomials defined on a mesh would be a very difficult task, especially if the space or parameter domains have a complex geometry. Neural networks have demonstrated their ability to approximate functions in fairly high dimensions, notably thanks to their intrinsic regularity. PINNs are a mesh-free approach to solving PDEs using neural networks. Their properties make them good candidates for approaching solutions to high-dimensional problems.

To build our prior, we propose to solve the parametric steady problem with PINNs. To that end, we now briefly introduce this method in Section 4.1, and we show how to compute and store the prior. Then, our algorithm is summarized in Section 4.2.

### Parametric PINNs

Note that the steady solutions to (1.1) are given by

\[\partial_{x}F_{\mu_{1}}(u)=S_{\mu_{2}}(u).\]

This is nothing but a parametric elliptic problem. Therefore, we introduce PINNs for the following generic boundary value problems (BVPs):

\[\begin{cases}\mathcal{D}(u,x;\mu)=f(x,\mu)&\text{ for }x\in\Omega,\\ u(t,x)=g(x,\mu),&\text{ for }x\in\partial\Omega,\end{cases} \tag{4.1}\]

where \(\mathcal{D}\) is a differential operator involving the solution \(u\) and its space derivatives, and with \(\mu\) some physical parameters. We recall that \(\mu\in\mathbb{P}\subset\mathbb{R}^{m}\). PINNs use the fact that classical fully-connected neural networks are
smooth functions of their inputs, as long as their activation functions are also smooth, to approximate the solution to (4.1). Contrary to traditional numerical schemes such as the DG method, where the degrees of freedom encode some explicit modal or nodal values of the solutions, the degrees of freedom of PINNs representation are the weights \(\theta\) of the neural network, and so do not explicitly represent the solution. Equipped with both of these remarks, the idea behind PINNs is to plug the network, which represents the solution to (4.1), into the equation. Then, the degrees of freedom (i.e. the weights \(\theta\) of the network) are found by minimizing a loss function. Since the neural network is differentiable, the derivatives can be exactly computed. In our case, the PINN is thus a smooth neural network that takes as input the space variable \(x\) and the parameter vector \(\mu\), which we denote by \(u_{\theta}(x;\mu)\).

Thanks to these definitions, solving the PDE can be rewritten as the following minimization problem:

\[\min_{\theta}\mathcal{J}(\theta),\quad\text{where }\mathcal{J}(\theta)= \mathcal{J}_{r}(\theta)+\mathcal{J}_{b}(\theta)+\mathcal{J}_{\text{data}}( \theta). \tag{4.2}\]

In (4.2), we have introduced three different terms: the residual loss function \(\mathcal{J}_{r}\), the boundary loss function \(\mathcal{J}_{b}\), and the data loss function \(\mathcal{J}_{\text{data}}\). For parameters \(\mu\in\mathbb{P}\), the residual loss function is defined by

\[\mathcal{J}_{r}(\theta)=\int_{\mathbb{P}}\int_{\Omega}\big{\|} \mathcal{D}(u_{\theta},x;\mu)-f(x;\mu)\big{\|}_{2}^{2}\,dxd\mu, \tag{4.3}\]

while the boundary loss function is given by

\[\mathcal{J}_{b}(\theta)=\int_{\mathbb{P}}\int_{\partial\Omega} \big{\|}u_{\theta}(x,\mu)-g(x,\mu)\big{\|}_{2}^{2}\,dxd\mu. \tag{4.4}\]

Finally, to define the data loss function, we assume that we know the exact solution to (4.1) at some points \(x_{i}\) and for some parameters \(\mu_{i}\), and we set

\[\mathcal{J}_{\text{data}}(\theta)=\sum_{i}\big{\|}u_{\theta}(x_{i},\mu_{i})-u (x_{i},\mu_{i})\big{\|}_{2}^{2}.\]

In practice, the integrals in (4.3) and (4.4) are approximated using a Monte-Carlo method. This method relies on sampling a certain number of so-called "collocation points" in order to approximate the integrals. Then, the minimization problem on \(\theta\) is solved using a gradient-type method, which corresponds to the learning phase.

The main advantage of PINNs is that they are mesh-free and less sensitive to dimension than classical methods. Indeed, neural networks easily deal with large input dimensions, and the Monte-Carlo method converges independently of the dimension. Consequently, PINNs are particularly well-suited to solving parametric PDEs such as (4.1). Thanks to that, we do not solve for a single equilibrium but rather for families of equilibria indexed by the parameters \(\mu\).

Traditional PINNs use this method to approximate both (4.3) and (4.4). However, for the boundary conditions, we elected to use another approach, which makes it possible to completely eliminate \(\mathcal{J}_{b}\) from the minimization algorithm. The idea is to define the approximate solution through a boundary operator \(\mathcal{B}\), which can for instance be a multiplication by a function which satisfies the boundary condition. We obtain

\[\widetilde{u}_{\theta}(x;\mu)=\mathcal{B}\big{(}u_{\theta},x;\mu\big{)},\]

with \(u_{\theta}\) the neural network and \(\mathcal{B}\) a simple operator such as \(\widetilde{u}_{\theta}\) exactly satisfies the boundary conditions. Using \(\widetilde{u}_{\theta}\), the residual loss becomes, instead of (4.3):

\[\mathcal{J}_{r}(\theta)=\int_{\mathbb{P}}\int_{\Omega}\big{\|} \mathcal{D}(\widetilde{u}_{\theta},x;\mu)-f(x;\mu)\big{\|}_{2}^{2}\,dxd\mu. \tag{4.5}\]

Examples of such functions \(\mathcal{B}\) are provided in Section5.

With this approach, we have presented one method for offline construction of our prior for a family of equilibria. Note that it is possible to further enhance this prior with data from previous simulations, thanks to the loss function \(\mathcal{J}_{\text{data}}\). Even though training PINNs may be harder than training traditional purely data-driven neural networks, they are much more efficient as priors. Indeed, the error estimates of Section3.2 show that the error
depends on the \(q\)-th derivative of the ratio between the prior and the solution. Therefore, to obtain a small error, it is important for the prior to provide a good approximation of not only the steady solution, but also of its derivatives. Since the PINN loss (4.5) inherently contains derivatives of \(u_{\theta}\), the resulting trained PINN will be more efficient in this respect. Note that a purely data-driven network could also be interesting if the data contains information on the derivatives.

### Algorithm

Now that we have discussed the strategy we use to obtain our prior, we give some details on the offline and online algorithms that we developed to construct the modified DG bases in practice. We start by describing the offline step, where the families of priors are computed. Then, we move on to an online algorithm, explaining how to construct the DG bases using the prior, and how to apply them to the actual DG time iterations.

```
1:space domain \(\Omega\), parameter set \(\mathbb{P}\), initial neural network \(u_{\theta_{0}}(x;\mu)\), \(\eta\) learning rate, \(N\) number of collocation points, \(n_{\text{epochs}}\) number of training epochs
2:trained neural network \(u_{\theta}(x;\mu)\)
3:initialize the weights: \(\theta=\theta_{0}\)
4:for\(n\leq n_{\text{epochs}}\)do
5: sample \(N\) values of \(x\) in \(\Omega\) and \(\mu\) in \(\mathbb{P}\)
6: compute the loss function \(\mathcal{J}(\theta)\)
7: update \(\theta\) using the gradient of \(\mathcal{J}(\theta)\): \(\theta=\theta-\eta\nabla_{\theta}J(\theta)\)
8:endfor
```

**Algorithm 1**Offline part: neural network training

In practice, we do not use a classical gradient descent to update the weights, but rather the Adam algorithm. Moreover, sampling is done through a uniform law on the space and parameter domains. It would also be possible use non-uniform sampling like in [52] for instance, but we elected to use uniform sampling for the sake of simplicity. Note that Algorithm 1 does not contain solution data in its inputs. Indeed, almost all numerical experiments from Section 5 do not require data on the solution. This avoids the cost of data production, which would otherwise require sampling the exact solution if it is known, or using a numerical scheme otherwise.

```
1:prior \(u_{\theta}\), degree \(N_{q}\) of the Gauss-Lobatto quadrature rule, initial data \(u_{0}\), space mesh \(\Omega_{h}\), parameters \(\mu\), \(n_{t}\) number of time steps
2:numerical solution \(u_{k}(t,x)\) on each cell \(\Omega_{k}\)
3:use the mesh \(\Omega_{h}\) to obtain all quadrature points \(x_{k,p}\) in each cell \(\Omega_{k}\)
4:evaluate the prior at each point \(x_{k,p}\): we obtain \(\tilde{u}_{k,p}\coloneqq u_{\theta}(x_{k,p};\mu)\)
5: reconstruct \(u_{k}(0,x)\) using \(\tilde{u}_{k,p}\)
6:for\(n\leq n_{t}\)do
7: construct the mass matrix \(\mathcal{M}\), the nonlinear flux \(\mathcal{F}\), the interface flux \(\mathcal{G}\) and the source term \(\mathcal{S}\) using \(\tilde{u}_{k,p}\) and the quadrature rule
8: update the solution \(u_{k}\) at the next time step, using \(u_{k}\) at the previous time step as well as the terms computed in the previous step
9:endfor
```

**Algorithm 2**Online part: using the neural network in the DG scheme

In this second step, the additional computational overhead associated to our method, compared to the classical DG scheme, comes from two distinct sources. The first one is a preprocessing phase, where we evaluate the prior on the quadrature points (step 2 of Algorithm 2). Even though such networks have been made to be quickly evaluated on GPUs, this evaluation step remains fast on CPUs. The second source of computational cost is associated to the quadrature rule. Indeed, in some cases, we will require a quadrature rule with a higher degree than the traditional DG scheme. The classical approach is to use \(N_{q}=q\) quadrature points for bases made of \(q\) polynomial functions, since the quadrature is exact for polynomials of degree \(q\). However, in our case, our basis
is non-polynomial. Hence, to have a good approximation of the integral of the prior, we may need to increase the degree of the quadrature. In most cases, this increase is slight; for a few test cases, especially to approximate functions with large derivatives, we will need to use fine quadrature rules.

## 5 Applications and numerical results

This section is dedicated to a validation of the approach on several parametric hyperbolic systems of balance laws: the linear advection equation in Section 5.1, the 1D shallow water equations in Section 5.2, the Euler-Poisson system in Section 5.3, and the 2D shallow water equations in Section 5.4. In the first two cases, there exist some exact well-balanced schemes in the literature. However, for the Euler-Poisson system and the 2D shallow water equations, exact (or even approximate) WB schemes are either not available or very complicated to implement. Code replicating the experiments of Section 5.1 is freely available [33] on GitHub1.

Footnote 1: [https://github.com/Victor-MichelDansac/DG-PINNs.git](https://github.com/Victor-MichelDansac/DG-PINNs.git)

In this section, we denote by \(K\) the number of cells. We test both bases \(V_{h}^{*}\) and \(V_{h}^{+}\) at first, showing that both display similar results. To cut down on the number of tables, we then only present the results for the additive basis \(V_{h}^{+}\).

Moreover, the time step \(\Delta t\) is given by

\[\Delta t=C_{\mathrm{CFL}}\,C_{\mathrm{RK}}\,\frac{\min_{k\in\{1,\ldots,K\}} \Delta x_{k}}{\lambda}, \tag{5.1}\]

where \(\lambda\) is the maximal wave speed of the system, \(C_{\mathrm{CFL}}\) is a CFL (Courant-Friedrichs-Lewy) coefficient, and \(C_{\mathrm{RK}}\) is the stability coefficient associated to the time discretization. All experiments are run using a strong stability-preserving Runge-Kutta (SSPRK) time discretization of the correct order. The time discretizations, with their associated stability coefficients \(C_{\mathrm{RK}}\), are collected in Table 1. To determine \(C_{\mathrm{CFL}}\), we run a study of the stability condition for the first experiment; this study is not repeated for the other experiments, since the new bases do not influence the stability condition.

### Linear advection

We first consider the case of a linear advection equation with a source term, on the space domain \(\Omega=(0,1)\). The equation is given as follows:

\[\left\{\begin{aligned} \partial_{t}u+\partial_{x}u& =s(u;\mu),\qquad\text{ for }x\in\Omega\\ u(t=0,x)&=u_{\mathrm{ini}}(x;\mu),\\ u(t,x=0)&=u_{0},\end{aligned}\right. \tag{5.2}\]

Here, the parameter vector \(\mu\) is made of three elements:

\[\mu=\begin{pmatrix}\alpha\\ \beta\\ u_{0}\end{pmatrix}\in\mathbb{P}\subset\mathbb{R}^{3},\quad\alpha\in\mathbb{R }_{+},\quad\beta\in\mathbb{R}_{+},\quad u_{0}\in\mathbb{R}_{+}^{*}.\]

The source term depends on \(\mu\) as follows:

\[s(u;\mu)=\alpha u+\beta u^{2},\]

\begin{table}
\begin{tabular}{c c c c c} \hline \hline number \(q\) of basis elements & 0 & 1 & 2 & 3 \\ \hline time discretization & explicit Euler & SSPRK2 [45] & SSPRK3(5) [45] & SSPRK4(10) [29] \\ stability coefficient \(C_{\mathrm{RK}}\) & 1 & 1 & 2.65 & 3 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Stability coefficients \(C_{\mathrm{RK}}\) of the high-order time discretizations used in the numerical experiments, with respect to the number \(q\) of basis elements.
and straightforward computations show that the associated steady solutions take the form

\[u_{\text{eq}}(x;\mu)=\frac{\alpha u_{0}}{(\alpha+\beta u_{0})e^{-\alpha x}-\beta u _{0}}. \tag{5.3}\]

To compute the time step, we take \(\lambda=1\) in (5.1), since the advection velocity in (5.2) is equal to \(1\). The first paragraph of this section shows how to choose \(C_{\text{CFL}}\) to complete the determination of the time step. Unless otherwise stated, we prescribe Dirichlet boundary conditions consisting in the steady solution.

To obtain a suitable prior \(u_{\theta}\), we train a PINN with parameters \(\theta\). To avoid cumbersome penalization of boundary conditions, we define \(\widetilde{u}_{\theta}\) using a boundary operator \(\mathcal{B}\), as follows:

\[\widetilde{u}_{\theta}(x;\mu)=\mathcal{B}(u_{\theta},x;\mu)=u_{0}+xu_{\theta} (x;\mu),\]

so that the boundary condition \(\widetilde{u}_{\theta}(0;\mu)=u_{0}\) is automatically satisfied by \(\widetilde{u}_{\theta}\). The parameter space \(\mathbb{P}\) is chosen such that the steady solution is well-defined, and we take

\[\mathbb{P}=[0.5,1]\times[0.5,1]\times[0.1,0.2]. \tag{5.4}\]

Thanks to the boundary operator \(\mathcal{B}\), the loss function only concerns the ODE residue, and we set

\[\mathcal{J}(\theta)=\left\|\partial_{x}\widetilde{u}_{\theta}-\alpha \widetilde{u}_{\theta}-\beta\widetilde{u}_{\theta}^{2}\right\|_{2}^{2}.\]

We use a neural network with \(5\) fully connected hidden layers, and around \(1200\) trainable parameters. Training takes about \(4\) minutes on a dual NVIDIA K80 GPU, until the loss is equal to about \(10^{-6}\). For this experiment, we increased the order of the quadrature compared to the baseline for the case with one basis function. Indeed, we take \(n_{Q}=\max(q+2,3)\), to ensure a sufficient precision when integrating the prior.

In this section, we compare four strategies: the basis \(V_{h}\) (2.4), the basis \(V_{h}^{*}\) with multiplicative prior (2.6), the basis \(V_{h}^{+}\) with additive prior (2.5), and the basis \(V_{h}^{\text{ex},+}\) which uses the exact steady solution (5.3) as a prior. First, we study the stability condition in Section 5.1.1. Then, we tackle the approximation of a steady solution without perturbation in Section 5.1.2 and with perturbation in Section 5.1.3. Finally, the approximation of an unsteady solution is computed in Section 5.1.4.

#### 5.1.1 Study of the stability condition

The very first experiment we run aims at making sure that the new bases do not alter the stability condition of the DG scheme. To that end, we slowly increase \(C_{\text{CFL}}\) until the time step \(\Delta t\) is too large for the scheme to be stable. For this experiment, the initial condition is made of the steady solution

\[u_{\text{ini}}(x;\mu)=u_{\text{eq}}(x;\mu), \tag{5.5}\]

and the final time is \(T=0.5\). Table 2 contains the optimal values of \(C_{\text{CFL}}\) (larger values leading to instabilities) obtained with the four bases and for \(q\in\{0,1,2,3\}\). We observe that the new bases do not change the stability condition, except for \(V_{h}^{\text{ex},+}\) with \(q=1\), which is slightly more stable. This study will not be repeated for other experiments, since it would yield similar results. In practice, we take \(C_{\text{CFL}}=0.1\) to ensure stability.

\begin{table}
\begin{tabular}{l c c c c} \hline \hline \(q\) & basis \(V_{h}\) & basis \(V_{h}^{*}\) & basis \(V_{h}^{+}\) & basis \(V_{h}^{\text{ex},+}\) \\ \hline \(0\) & \(1.250\) & \(1.250\) & \(1.250\) & \(1.250\) \\ \(1\) & \(0.399\) & \(0.399\) & \(0.399\) & \(0.416\) \\ \(2\) & \(0.209\) & \(0.209\) & \(0.209\) & \(0.209\) \\ \(3\) & \(0.185\) & \(0.185\) & \(0.185\) & \(0.185\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Maximal values of \(C_{\text{CFL}}\) obtained for the four bases and for a number of basis elements \(q\in\{0,1,2,3\}\).
#### 5.1.2 Steady solution

We now study the approximation of a steady solution, with and without perturbation. The goal of this section is to check whether the prior indeed makes it possible to decrease the error compared to the usual modal basis. For this experiment, the initial condition remains (5.5), and the final time is \(T=0.1\).

As a first step, the values of the parameters \(\mu\) are set to the midpoints of the intervals making up the parameter space (5.4). The \(L^{2}\) errors between the exact and approximate solutions are collected in Table 3.

In this case, we expect both \(V_{h}^{*}\) and \(V_{h}^{+}\) to show similar behavior. Moreover, we expect the basis \(V_{h}^{\mathrm{ex},+}\) to provide an exactly well-balanced scheme, up to machine precision. To that end, only for \(V_{h}^{\mathrm{ex},+}\), we take \(n_{Q}=\max(q+2,5)\), to ensure that the quadrature of the exact prior is also exact, up to machine precision.

We observe that the bases with and without prior allow a convergence of the correct order, i.e. of the same order as the number of basis elements. Moreover, we observe a consistent gain for all mesh resolutions, for a given size of the modal basis, which is lower the larger the size of the basis. Bases \(V_{h}^{*}\) and \(V_{h}^{+}\) seem to have comparable performance, with \(V_{h}^{*}\) being somewhat better for large values of \(q\), and \(V_{h}^{+}\) taking the lead for small values of \(q\). Finally, we observe that the basis \(V_{h}^{\mathrm{ex},+}\) is indeed able to provide a solution that is exact up to machine precision, thus validating the exact well-balanced property of the scheme using this basis.

As a second step, to refine this study, we now consider \(10^{3}\) parameters, randomly sampled from the parameter space (5.4). For \(q\in\{0,1,2,3\}\) and \(K=10\) discretization cells, we compute the minimum, average and maximum gains obtained with both bases \(V_{h}^{*}\) and \(V_{h}^{+}\). These values are reported in Table 4. We observe, on average, a significant gain in all cases, with larger gains obtained for smaller values of \(q\). Furthermore, the minimum gain is always greater than one. Like in the previous experiment, we observe that, even though both bases display similar behavior and very good results, \(V_{h}^{+}\) behaves better than \(V_{h}^{*}\) for small values of \(q\), and vise versa. Consequently, and to limit the number of tables in the remainder of this section, we perform all subsequent experiments with the basis \(V_{h}^{+}\).

#### 5.1.3 Perturbed steady solution

We now test the scheme on a perturbed steady solution, For this experiment, the initial condition is similar to (5.5), but with a perturbation. Indeed, we take

\[u_{\mathrm{ini}}(x;\mu)=\big{(}1+\varepsilon\sin(2\pi x)\big{)}u_{\mathrm{eq} }(x;\mu),\]

where \(\varepsilon\) is taken nonzero or zero, to control the strength of the perturbation. The final time is \(T=2\), and we study the impact of the perturbation by taking \(\varepsilon\in\{10^{-4},10^{-2},1\}\), and \(K=10\) discretization cells. The results are collected in Figure 1. We observe two different states: first, while the perturbation is being dissipated, the errors with the two bases are similar. Then, we note that the introduction of the prior has made it possible for the approximate solution to converge towards a final solution that is closer to the exact, unperturbed steady solution.

#### 5.1.4 Unsteady solution

Next, we seek to confirm that our proposed basis does not deteriorate the approximation of unsteady solutions. To that end, we consider an unsteady solution of the homogeneous problem, i.e. a solution to (5.2) with \(s(u;\mu)=0\). We take the following initial condition:

\[u_{0}(x)=0.1\left(1+\exp\left(-100(x-0.5)^{2}\right)\right),\]

so that \(u(t,x)=u_{0}(x-t)\). The final time is set to \(T=1\), and periodic boundary conditions are prescribed.

We compute the approximate solution with the two bases, for several values of \(q\). The results are collected in Table 5. We note that the basis with prior does not affect the approximate solution for \(q\geq 1\), while the results are slightly worse with the prior for \(q=0\). To improve the results here, one could introduce a space-time basis in a space-time discontinuous Galerkin method; this will be the object of future work.
\begin{table}

\end{table}
Table 3: Advection equation: errors, orders of accuracy, and gain obtained when approximating a steady solution for bases without prior (basis \(V_{h}\)), with a PINN prior (bases \(V_{h}^{*}\) and \(V_{h}^{+}\)), and with an exact prior (basis \(V_{h}^{\text{ex},+}\)).
### Shallow water equations

After studying a scalar linear advection equation in Section 5.1, we now turn to a nonlinear system of conservation laws. Namely, we tackle the shallow water equations

\[\begin{cases}\partial_{t}h+\partial_{x}Q=0,\\ \partial_{t}Q+\partial_{x}\left(\frac{Q^{2}}{h}+\frac{1}{2}g\frac{1}{h} \overline{\overline{\varepsilon}}\right)=-gh\partial_{x}Z(x;\alpha,\beta), \end{cases} \tag{5.6}\]

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline  & \multicolumn{3}{c}{gains in basis \(V_{h}^{*}\)} & \multicolumn{3}{c}{gains in basis \(V_{h}^{+}\)} \\ \cline{2-6} \(q\) & minimum & average & maximum & minimum & average & maximum \\ \hline
0 & 63.46 & 735.08 & 4571.89 & 63.46 & 735.08 & 4571.89 \\
1 & 32.22 & 149.38 & 450.74 & 26.01 & 190.08 & 830.20 \\
2 & 6.20 & 54.16 & 118.45 & 5.92 & 45.47 & 313.07 \\
3 & 1.55 & 19.54 & 108.10 & 1.56 & 13.69 & 184.17 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Advection equation: statistics of the gains obtained for the approximation of a steady solution in bases \(V_{h}^{*}\) and \(V_{h}^{+}\) with respect to basis \(V_{h}\).

Figure 1: Advection equation: errors, with respect to time, for the approximation of a perturbed steady solution for bases with and without prior.
where \(h>0\) is the water height, \(Q\) the water discharge, \(g=9.81\) the gravity constant, and where the parameterized topography function is

\[Z(x;\alpha,\beta)=\beta\omega\left(\alpha\left(x-\frac{1}{2}\right)\right). \tag{5.7}\]

In (5.7), the function \(\omega\in\{\omega_{g},\omega_{c}\}\) is either a Gaussian bump function

\[\omega_{g}(x)=\frac{1}{4}e^{-50x^{2}} \tag{5.8}\]

or a compactly supported bump function, with parameter \(\mathfrak{G}=0.15\):

\[\omega_{c}(x)=\begin{cases}\exp\left(1-\frac{1}{1-\left(\frac{x}{\mathfrak{G }}\right)^{2}}\right)&\text{ if }|x|<\mathfrak{G},\\ 0&\text{ otherwise.}\end{cases} \tag{5.9}\]

Unless otherwise mentioned, the final physical time is \(T=0.05\), and the space domain is \(\Omega=(0,1)\). For each experiment, Dirichlet boundary conditions corresponding to the steady solution are prescribed.

The steady solutions are given by cancelling the time derivatives in (5.6), and we get the following characterization:

\[Q_{\text{eq}}=\text{constant}=:Q_{0}\qquad\text{and}\qquad\left(1-\frac{Q_{0} ^{2}}{gh_{\text{eq}}(x;\mu)^{3}}\right)\partial_{x}h_{\text{eq}}(x;\mu)+ \partial_{x}Z(x;\alpha,\beta)=0. \tag{5.10}\]

To solve the nonlinear ODE on \(h\), we impose \(h=h_{0}\) at some point in space. Without loss of generality, we restrict the study to the case \(Q_{0}>0\). This leads us to a family of steady solutions with four parameters, and thus a parameter vector \(\mu\) made of four elements:

\[\mu=\begin{pmatrix}\alpha\\ \beta\\ h_{0}\\ Q_{0}\end{pmatrix}\in\mathbb{P}\subset\left(\mathbb{R}_{+}^{*}\right)^{4}.\]

\begin{table}

\end{table}
Table 5: Advection equation: errors, orders of accuracy, and gain obtained when approximating an unsteady solution for bases with and without prior.
Hence, we compute \(\Delta t\) in (5.1) by taking \(\lambda=\frac{Q_{0}}{h_{0}}+\sqrt{gh_{0}}\).

Depending on the values of these parameters, the Froude number

\[\text{Fr}=\sqrt{\frac{Q^{2}}{gh^{3}}}\]

controls the so-called flow regime for the steady solution. They can be in three distinct regimes: subcritical (\(\text{Fr}<1\) everywhere), supercritical (\(\text{Fr}>1\) everywhere) or transcritical (\(\text{Fr}=1\) somewhere in the domain). Each regime has its own parameter space for \(h_{0}\) and \(Q_{0}\), described later, but in all cases we take, unless otherwise stated,

\[0.5\leq\alpha\leq 1.5\quad;\qquad 0.5\leq\beta\leq 1.5. \tag{5.11}\]

To approximate the steady water height within this parameter space, we use a fully-connected PINN with about \(4000\) trainable parameters. Its result \(\widetilde{h}_{\theta}\) is modified through a boundary function \(\mathcal{B}\) that will be defined for each regime. The loss function is once again made only of the steady ODE, and we minimize

\[\mathcal{J}(\theta)=\left\|\left(1-\frac{Q_{0}^{2}}{g\widetilde{h}_{\theta}(x; \mu)^{3}}\right)\partial_{x}\widetilde{h}_{\theta}(x;\mu)+\partial_{x}Z(x; \alpha,\beta)\right\|.\]

Training takes about \(5\) minutes on a dual NVIDIA K80 GPU, and lasts until the loss is about \(10^{-4}\), depending on the regime.

#### 5.2.1 Subcritical flow

We start with a subcritical flow, where the parameter space for \(h_{0}\) and \(Q_{0}\) is:

\[2\leq h_{0}\leq 3\quad;\qquad 3\leq Q_{0}\leq 4. \tag{5.12}\]

To strongly enforce the boundary conditions, the prior \(\widetilde{h}_{\theta}\) is obtained as follows from the result \(h_{\theta}\) of the PINN:

\[\widetilde{h}_{\theta}(x;\mu)=\mathcal{B}(h_{\theta},x;\mu)=h_{0}+Z(x;\alpha, \beta)\,h_{\theta}(x;\mu). \tag{5.13}\]

To test the preservation of the steady solution, we set the initial water height to \(h_{\text{eq}}\).

A goal of this section is to better understand the differences between the two topography functions: the Gaussian bump (5.8) and the compactly supported bump (5.9). It is well-known that compactly supported functions exhibit large derivatives close to the support, see for instance [44]. As a consequence, to get a good approximation of these derivatives when computing integrals involving the PINN, we take \(n_{Q}=q+6\) when \(\omega=\omega_{c}\). Note that this choice is also motivated by the results in [44], where the authors had to take larger polynomial degrees to observe the correct orders of convergence. The Gaussian topography also suffers from the same drawback, but to a lesser extent, and we take \(n_{Q}=q+3\) when \(\omega=\omega_{g}\) when integrating the result of the PINN.

For the compactly supported topography, the results are reported in Table 6; for the Gaussian topography, the results are reported in Table 7.

As a conclusion of this first test case, we observe that using a Gaussian topography compared to a compactly supported topography leads to a more stable order of accuracy, but with lower gains, except for small values of \(K\) where the compactly supported topography is not well-approximated. The most important point is that the Gaussian topography requires a lower order quadrature to converge. These results are in line with [44]. As a consequence, we use the Gaussian topography in the remainder of this section.

Like in the previous section, we now consider \(10^{3}\) parameters in \(\mathbb{P}\), and we compute the minimum, average and maximum gains for \(q\in\{0,1,2\}\). To that end, we take \(K=20\) discretization cells. The results are reported in Table 8, where we observe that the average gains are substantial, whatever the value of \(q\), and that the minimum gain is always greater than \(1\).
#### 5.2.2 Supercritical flow

We now turn to a supercritical flow. In this case, the remaining parameters \(h_{0}\) and \(Q_{0}\) are taken such that:

\[0.5\leq h_{0}\leq 0.75\quad;\qquad 4\leq Q_{0}\leq 5. \tag{5.14}\]

The boundary conditions are enforced using the same expression (5.13) as in the subcritical case. We check the approximate preservation of the steady solution by taking the initial water height equal to the steady solution.

The results are displayed in Table 9, and we note that the gains are in line with the subcritical case, from Table 7.

Furthermore, in Table 10, we display some statistics on the gains obtained by using the prior, in the same configuration as for the subcritical regime. We draw similar conclusions to the subcritical case.

#### 5.2.3 Transcritical flow

The last steady experiment we study is the preservation of a transcritical steady solution. Such steady solutions are significantly harder to capture. Indeed, when \(\mathrm{Fr}=1\), the steady ODE (5.10) yields \(\partial_{x}Z=0\), and therefore the derivative of the steady water height is not defined using only (5.10). This is a well-known issue when approximating transcritical steady solutions, see for instance [14, 24].

\begin{table}

\end{table}
Table 6: Shallow water system, compactly supported topography (5.9): errors, orders of accuracy, and gain obtained when approximating a subcritical steady solution for bases with and without prior.
However, in our case, we provide a simple PINN by defining the water height such that the Froude number is equal to \(1\) at the top of the topography bump, i.e. at \(x=1/2\), where \(\partial_{x}Z=0\). When \(\mathrm{Fr}=1\), the water height becomes equal to \(h_{c}(\mu)=Q_{0}^{2/3}g^{-1/3}\), and we fix this value for the prior evaluated at \(x=1/2\). This eliminates \(h_{0}\) as a degree of freedom, and we choose \(2\leq Q_{0}\leq 3\). Moreover, we take \(0.75\leq\alpha\leq 1.25\) for this regime.

Then, to ensure a correct treatment of the boundary conditions and to obtain the correct value of \(\widetilde{h}_{\theta}\) at the

\begin{table}

\end{table}
Table 7: Shallow water system, Gaussian topography (5.8): errors, orders of accuracy, and gain obtained when approximating a subcritical steady solution for bases with and without prior.

\begin{table}

\end{table}
Table 8: Shallow water system, Gaussian topography (5.8): statistics of the gains obtained for the approximation of a subcritical steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).
top of the bump, we take

\[\widetilde{h}_{\theta}(x;\mu)=h_{R}(\mu)+\left(1-\tanh\left(15\left(x-\frac{1}{2} \right)\right)\right)\frac{h_{L}(\mu)-h_{R}(\mu)}{2}\ h_{\theta}(x;\mu).\]

In this expression, \(h_{L}(\mu)\) and \(h_{R}(\mu)\) are the left and right boundary conditions. Since we consider a smooth

\begin{table}

\end{table}
Table 10: Shallow water system, Gaussian topography (5.8): statistics of the gains obtained for the approximation of a supercritical steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).

\begin{table}

\end{table}
Table 9: Shallow water system, Gaussian topography (5.8): errors, orders of accuracy, and gain obtained when approximating a supercritical steady solution for bases with and without prior.
steady solution, relations (5.10) lead to

\[E(h,x;\mu)\coloneqq\frac{Q_{0}^{2}}{2h^{2}}+g\big{(}h+Z(x;\alpha,\beta)\big{)}= \text{constant}.\]

Since \(Z(0.5;\alpha,\beta)=\beta/4\), we obtain that \(h_{L}(\mu)>h_{R}(\mu)\) are the two solutions of the following equation, with unknown \(h\):

\[\frac{Q_{0}^{2}}{2h^{2}}+gh=E\bigg{(}h_{c}(\mu),\frac{1}{2};\mu\bigg{)}.\]

Table 11 contains the errors, order of convergence and the gains. We observe that the gains are lower than in the other two cases, but that was to be expected since the transcritical solution comes from a singular ODE, and it is harder for the PINN to approximate its solutions.

Finally, we report in Table 12 the minimum, average and maximum gains obtained by using the basis \(V_{h}^{+}\) instead of the basis \(V_{h}\). We draw the same conclusions as in the other two regimes, even though the gains are, on average, lower. This was to be expected, since the transcritical regime is harder to capture than the subcritical and supercritical ones, and therefore that the prior is of lower quality. Nevertheless, the gains remain substantial for all values of \(q\).

\begin{table}

\end{table}
Table 11: Shallow water system, Gaussian topography (5.8): errors, orders of accuracy, and gain obtained when approximating a transcritical steady solution for bases with and without prior.
#### 5.2.4 Perturbation of a steady flow

This last experiment related to the shallow water equations concerns a perturbed steady flow. We only perform this study on the subcritical flow, but the other regimes behave the same. We take \(\varepsilon\in\{5\cdot 10^{-k}\}_{k\in\{1,2,3\}}\) and \(20\) space cells, and set the initial water height to \(h(0,x;\mu)=\big{(}1+\varepsilon\sin(2\pi x)\big{)}h_{\text{eq}}(x;\mu)\).

The errors on \(h\) with respect to time are displayed in Figure 2, until the final physical time \(T=1\). Like in Section 5.1, with the prior, the error decreases to a much lower level than without the prior. This good behavior was expected since the prior makes it possible for the enhanced DG scheme to achieve higher accuracy on steady solutions.

### Euler-Poisson equations in spherical geometry

We now consider the Euler-Poisson equations in spherical geometry. This system is used in astrophysics, for instance, where it serves to model stars held together by gravitation, see e.g. [15, 17, 32]. They are given by

\[\begin{cases}\partial_{t}\rho+\partial_{r}Q=-\dfrac{2}{r}Q,\\ \partial_{t}Q+\partial_{r}\left(\dfrac{Q^{2}}{\rho}+p\right)=- \dfrac{2}{r}\dfrac{Q^{2}}{\rho}-\rho\partial_{r}\phi,\\ \partial_{t}E+\partial_{r}\left(\dfrac{Q}{\rho}(E+p)\right)=- \dfrac{2}{r}\dfrac{Q}{\rho}(E+p)-Q\partial_{r}\phi,\\ \dfrac{1}{r^{2}}\partial_{rr}(r^{2}\phi)=4\pi G\rho,\end{cases} \tag{5.15}\]

where \(G\) is a gravity constant, fixed to \(G=1\) in our applications, and where we take \(p\) as a function of \(\rho\), \(Q\) and \(E\) through a pressure law to be specified. In (5.15), \(\rho\) is the density, \(Q\) is the momentum, \(E\) is the energy, and \(\phi\) is the gravitational potential. Unless otherwise mentioned, the boundary conditions are of Dirichlet type, with the value of the steady solution prescribed of the boundaries.

The space domain is \(r\in(0,1)\). The apparent singularity at \(r=0\) is resolved by imposing suitable boundary conditions, namely \(\rho(0)=1\) and \(\partial_{r}\rho(0)\) given according to the pressure law. Indeed, the assumption that there is no gravity at \(r=0\) leads to \(\partial_{r}p(0)=0\), which makes it possible to determine \(\partial_{r}\rho(0)\). For more information on the boundary conditions and on the DG discretization of (5.15), the reader is referred to [57].

The steady solutions at rest are given by

\[\begin{cases}Q=0,\\ \partial_{r}p+\rho\partial_{r}\phi=0,\\ \partial_{rr}(r^{2}\phi)=4\pi r^{2}G\rho.\end{cases}\]

For the steady solutions, we shall distinguish two cases for the pressure law: a polytropic pressure law, and a temperature-dependent pressure law.

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline  & \multicolumn{2}{c}{minimum gain} & \multicolumn{2}{c}{average gain} & \multicolumn{2}{c}{maximum gain} \\ \cline{2-6} \cline{5-6} \(q\) & \(h\) & \(Q\) & \(h\) & \(Q\) & \(h\) & \(Q\) \\ \cline{2-6} \(0\) & \(35.82\) & \(26.19\) & \(254.53\) & \(177.02\) & \(928.03\) & \(668.73\) \\ \(1\) & \(5.51\) & \(4.73\) & \(30.83\) & \(38.69\) & \(134.83\) & \(142.11\) \\ \(2\) & \(4.55\) & \(6.16\) & \(16.49\) & \(24.29\) & \(96.95\) & \(109.94\) \\ \hline \hline \end{tabular}
\end{table}
Table 12: Shallow water system, Gaussian topography (5.8): statistics of the gains obtained for the approximation of a transcritical steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).
#### 5.3.1 Polytropic pressure law

In this case, we introduce two parameters \(\kappa\) and \(\gamma\), so the parameter vector \(\mu\) is composed of two elements:

\[\mu=\begin{pmatrix}\kappa\\ \gamma\end{pmatrix}\in\mathbb{P}\subset\mathbb{R}^{2},\quad\kappa\in\mathbb{R} _{+},\quad\gamma\in(1,+\infty).\]

Equipped with this parameter vector, we define the polytropic pressure law

\[p(\rho;\mu)=\kappa\rho^{\gamma},\]

and the steady solutions are then given as solutions to the following nonlinear second-order ordinary differential equation:

\[\frac{d}{dr}\left(r^{2}\kappa\gamma\rho^{\gamma-2}\frac{d\rho}{dr}\right)=4 \pi r^{2}G\rho.\]

Figure 2: Shallow water equations, compactly supported topography: errors, with respect to time, for the approximation of a perturbed subcritical steady solution for bases with and without prior.
In general, this ODE does not have analytic solutions. However, it turns out that, for specific values of \(\gamma\), there exists an analytic solution to this ODE. For instance, with \(\gamma=2\), we obtain

\[\rho(r)=\frac{\sin(\alpha r)}{\alpha r},\qquad\text{with}\quad\alpha=\sqrt{\frac {2\pi G}{\kappa}}.\]

Regarding the boundary conditions, the condition \(\partial_{r}p(0)=0\) leads to \(\partial_{r}\rho(0)=0\) for this pressure law. In this case, we take \(\lambda=1+\sqrt{\gamma}\) in (5.1) to compute the time step \(\Delta t\).

To obtain a prior \(\rho_{\theta}\), as usual, we train a PINN with about 1400 trainable parameters on 7 fully connected layers. The boundary conditions are taken into account by setting

\[\widetilde{\rho}_{\theta}(r;\mu)=1+r^{2}\rho_{\theta}(r;\mu),\]

where \(\rho_{\theta}\) is the result of the PINN. The PINN is trained on the parameter space

\[\mathbb{P}=[2,5]\times[1.5,3.5], \tag{5.16}\]

with only the physics-based loss function corresponding to the steady solution:

\[\mathcal{L}_{\theta}=\left\|\frac{d}{dr}\left(r^{2}\kappa\gamma\widetilde{ \rho}_{\theta}^{\gamma-2}\frac{d\widetilde{\rho_{\theta}}}{dr}\right)-4\pi r^ {2}G\widetilde{\rho}_{\theta}\right\|.\]

In addition, the prior for \(q\) is set to \(Q_{\theta}=1\) since we wish to approximate a constant momentum. Finally, the prior for \(E\) is set to \(E_{\theta}=p(\widetilde{\rho}_{\theta};\mu)/(\gamma-1)\). Training takes about 5 minutes on a dual NVIDIA K80 GPU, until the loss is equal to about \(5\cdot 10^{-5}\). In the DG discretization, the degree of the quadrature formula is the usual \(n_{Q}=q+2\): there is no need to further increase the order of the quadrature rule in this case.

We first collect, in Table 13, the results of the approximation in both bases (with and without prior), for \(\kappa=2\) and \(\gamma=2.5\), and until the final time \(T=0.01\). As usual, the observed gain is larger for smaller number of basis elements. We observe a slight superconvergence on the momentum \(Q\) when using the prior with \(q=0\). For these values of \(\kappa\) and \(\gamma\), gains on the density are not very large for \(q=2\), but this is compensated by larger gains on the energy.

To extend this study, we compute the statistics over the whole parameter space (5.16) by uniformly sampling \(10^{3}\) values and taking 10 cells in the mesh. The results are reported in Table 14. Just like before, the average gain is substantial, while the minimum rarely falls below 1. Moreover, note that the gains recorded in Table 13 correspond to a rather bad set of parameters compared to the average.

#### 5.3.2 Temperature-dependent pressure law

In this case, we take a given smooth temperature function \(T(r,\mu)\) parameterized by \(\mu\), where the parameter vector \(\mu\) is composed of two elements:

\[\mu=\begin{pmatrix}\kappa\\ \alpha\end{pmatrix}\in\mathbb{P}\subset\mathbb{R}^{2},\quad\kappa\in\mathbb{R} _{+},\quad\alpha\in\mathbb{R}_{+}.\]

This allows us to define the parameterized temperature function \(T(r;\alpha)=e^{-\alpha r}\), and so we get the following temperature-based pressure law:

\[p(\rho;\mu)=\kappa\rho T.\]

For this pressure law, the steady solutions are given by the following nonlinear second-order ODE:

\[\frac{d}{dr}\left(r^{2}\kappa\frac{T}{\rho}\frac{d\rho}{dr}\right)+\frac{d}{ dr}\left(r^{2}\kappa\frac{dT}{dr}\right)=4\pi r^{2}G\rho,\]

and the boundary condition \(\partial_{r}p(0)=0\) leads to \(\partial_{r}\rho(0)=\alpha\). For this pressure law, we also take \(\lambda=1+\sqrt{\gamma}\) in (5.1) to compute \(\Delta t\).
The prior \(\rho_{\theta}\) is obtained _via_ a PINN with the same characteristics as in the polytropic case, and whose result is still denoted by \(\rho_{\theta}\). To impose the boundary conditions, this time, we set

\[\widetilde{\rho}_{\theta}(r;\mu)=1+\alpha r+r^{2}\rho_{\theta}(r;\mu).\]

The parameter space is

\[\mathbb{P}=[2,5]\times[0.5,1.5], \tag{5.17}\]

\begin{table}

\end{table}
Table 13: Euler-Poisson system, polytropic pressure law: errors, orders of accuracy, and gain obtained when approximating a steady solution for bases with and without prior.

\begin{table}

\end{table}
Table 14: Euler-Poisson system, polytropic pressure law: statistics of the gains obtained for the approximation of a steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).
and the PINN is trained using only the physics-based loss function

\[\mathcal{L}_{\theta}=\left\|\frac{d}{dr}\left(r^{2}\kappa\frac{T}{\tilde{\rho}_{ \theta}}\frac{d\widetilde{\rho}_{\theta}}{dr}\right)+\frac{d}{dr}\left(r^{2} \kappa\frac{dT}{dr}\right)=4\pi r^{2}G\widetilde{\rho}_{\theta}\right\|.\]

Training takes about \(5\) minutes on a dual NVIDIA K80 GPU, until the loss is equal to about \(5\cdot 10^{-4}\). The priors \(Q_{\theta}\) and \(E_{\theta}\) are then defined in the same way as in the polytropic case. In this case, we also take \(n_{Q}=q+2\).

As is becoming usual, we first report, in Table 15, the results of the approximation in both bases (with and without prior). The final time is set to \(T=0.01\), and we take \(\kappa=3.5\) and \(\alpha=0.5\). As usual, using the prior provides significant gains, especially for low values of \(q\). Compared to the polytropic case, gains are consistently better for the large values of \(q\).

To understand gains on the whole parameter space (5.17), we uniformly sample \(10^{3}\) values of \(\kappa\) and \(\alpha\) and take a mesh made of \(10\) cells. We compute the minimum, average and maximum gains. These values are reported in Table 16. For this pressure law, the minimum gain is always larger than \(1\), and we obtain consistently large average gains, even for \(q=2\).

#### 5.3.3 Spherical blast wave

The goal of this last test case is to show that our prior does not negatively affect the capability of the scheme to capture discontinuous solutions. Let us emphasize that numerical viscosity is not an object of this study, and

\begin{table}

\end{table}
Table 15: Euler-Poisson system, temperature-based pressure law: errors, orders of accuracy, and gain obtained when approximating a steady solution for bases with and without prior.
therefore that we have not used any regularization procedure. Consequently, results will show some oscillations.

This experiment is nothing but a Riemann problem in spherical geometry, inspired by the experiments in [50]. As such, the initial condition is piecewise constant on the space domain \(r\in(0,0.4)\), as follows:

\[\rho(0,x)=\begin{cases}2&\text{if }r<0.2,\\ 1&\text{otherwise};\end{cases}\qquad Q(0,x)=0;\qquad p(0,x)=\begin{cases}2& \text{if }r<0.2,\\ 1&\text{otherwise}.\end{cases}\]

For this experiment, the pressure law is the standard ideal gas law

\[p=(\gamma-1)\left(E-\frac{1}{2}\frac{Q^{2}}{\rho}\right),\]

and we take the gas constant \(\gamma\) equal to \(1.4\). The experiment is run until the final time \(T=0.1\), and with Neumann boundary conditions. We take \(25\) discretization cells, and we use a basis made of \(3\) elements. Moreover, the source term is deactivated: we set \(\phi=0\), and we merely consider the Euler equations in spherical geometry, without gravity effects.

The results are depicted in Figure 3, where we compare the two bases (with and without prior, blue and orange lines respectively) to a reference solution (green line). We observe very good agreement with the reference solution, even though oscillations are present, as expected. We also note that the graphs for the solutions with and without prior are superimposed, which means that the quality of the approximation of this discontinuous solution has not been degraded by the introduction of the prior in the basis.

\begin{table}
\begin{tabular}{l r r r r r r r r} \hline \hline  & \multicolumn{4}{c}{minimum gain} & \multicolumn{4}{c}{average gain} & \multicolumn{4}{c}{maximum gain} \\ \cline{2-9} \cline{5-10} \(q\) & \(\rho\) & \(Q\) & \(E\) & \(\rho\) & \(Q\) & \(E\) & \(\rho\) & \(Q\) & \(E\) \\ \hline
0 & 13.30 & 1.05 & 16.24 & 151.96 & 1.88 & 150.63 & 600.13 & 2.91 & 473.83 \\
1 & 6.30 & 7.53 & 5.40 & 72.63 & 77.20 & 51.09 & 321.20 & 302.58 & 257.19 \\
2 & 3.35 & 3.45 & 2.20 & 18.96 & 22.58 & 13.56 & 55.47 & 63.45 & 47.83 \\ \hline \hline \end{tabular}
\end{table}
Table 16: Euler-Poisson system, temperature-based pressure law: statistics of the gains obtained for the approximation of a steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).

Figure 3: Euler equations, temperature-based pressure law: statistics for the approximation of a steady solution for bases with and without prior.
### Shallow water equations in two space dimensions

The last system considered in this series of experiments is the two-dimensional shallow water system. It is given by

\[\begin{cases}\partial_{t}h+\mathbf{\nabla}\cdot\mathbf{Q}=0,\\ \partial_{t}\mathbf{Q}+\mathbf{\nabla}\cdot\left(\frac{\mathbf{Q}\otimes\mathbf{Q}}{h}+\frac{1} {2}gh^{2}\text{Id}\right)=-gh\mathbf{\nabla}Z(\mathbf{x};\mu),\end{cases}\]

where \(g\) is the gravity constant, \(\text{Id}\) is the \(2\times 2\) identity matrix, \(h\) is the water height, \(\mathbf{Q}\) is the water discharge, and \(Z\) is the topography. For this system, \(\Delta t\) is computed by setting \(\lambda=2+\sqrt{\gamma}\) in (5.1).

The space variable \(\mathbf{x}=(x_{1},x_{2})\) belongs to the space domain \(\Omega=[-3,3]^{2}\), and we introduce three parameters:

\[\mu=\begin{pmatrix}\alpha\\ \Gamma\\ r_{0}\end{pmatrix}\in\mathbb{P}\subset\mathbb{R}^{3},\quad\alpha\in\mathbb{R} ^{*}_{+},\quad\Gamma\in\mathbb{R}^{*}_{+},\quad r_{0}\in\mathbb{R}^{*}_{+}.\]

This enables us to define the topography as the following Gaussian bump function, with \(r=\|\mathbf{x}\|\):

\[Z(\mathbf{x};\mu)=\Gamma\exp\left(\alpha(r_{0}^{2}-r^{2})\right),\]

see for instance [43] for a similar test case. On this topography, we consider the following steady solution:

\[\begin{cases}h_{\text{eq}}(\mathbf{x};\mu)=2-Z(\mathbf{x};\mu)-\frac{\Gamma}{8\alpha g }Z(\mathbf{x};\mu)^{4},\\ \mathbf{Q}_{\text{eq}}(\mathbf{x};\mu)=-\mathbf{x}^{\perp}\,h_{\text{eq}}(\mathbf{x};\mu)\,u_ {\text{eq}}(\mathbf{x};\mu),\end{cases} \tag{5.18}\]

with \(\mathbf{x}^{\perp}=(-x_{2},x_{1})\) and where

\[u_{\text{eq}}(\mathbf{x};\mu)=\alpha Z(\mathbf{x};\mu)^{2}.\]

To obtain a relevant prior, we approximate \(h_{\text{eq}}\) and \(u_{\text{eq}}\), using a different PINN for each of the two functions. The results of the PINN are denoted by \(h_{\theta}\) and \(u_{\theta}\), and we define the priors \(\widetilde{h}_{\theta}\) and \(\widetilde{u}_{\theta}\) as follows, to include the boundary conditions:

\[\widetilde{h}_{\theta}(\mathbf{x};\mu)=2-Z(\mathbf{x};\mu)\,h_{\theta}(\mathbf{x};\mu)^{ 2}\quad\text{and}\quad\widetilde{u}_{\theta}(\mathbf{x};\mu)=Z(\mathbf{x};\mu)\,u_{ \theta}(\mathbf{x};\mu).\]

Another possibility would be to strongly impose the divergence-free constraint, by learning a potential and taking the prior \(\mathbf{Q}_{\theta}(\mathbf{x};\mu)\) as the curl of this potential. However, we elected not to do so, since the current strategy was able to train faster. The parameter space is

\[\mathbb{P}=[0.25,0.75]\times[0.1,0.4]\times[0.5,1.25].\]

The loss function is made in equal parts of the now usual PDE loss, and of the minimization with respect to data. Data is regenerated at each epoch, and helps to avoid falling in a local minimum corresponding to a lake at rest, where \(\widetilde{h}_{\theta}+Z=\text{constant}\) and \(u_{\theta}=0\). Each PINN has about 2500 parameters, and training takes about 10 minutes on an NVIDIA V100 GPU, until the loss functions reaches about \(4\times 10^{-7}\). This prior is integrated with a quadrature formula of degree \(n_{Q}=q+3\): we needed to increase the usual quadrature degree by 1 to obtain the best possible approximation.

#### 5.4.1 Approximation of a steady solution

We take the steady solution (5.18) as the initial condition to test the approximate well-balanced property. The experiments are run until the final physical time \(T=0.01\). We prescribe Dirichlet boundary conditons consisting in the value of the steady solution.

First, we take the parameters as the center of the parameter cube \(\mathbb{P}\). The results are collected in Table 17, and we note that, as expected, the presence of the prior makes it possible to reach much lower errors, especially for the water height \(h\).

In addition, we provide some statistics over the whole parameter space \(\mathbb{P}\), computed on a mesh with \(25\times 25\) cells, in Table 18. We not that, on average, the gains are substantial. However, note that the minimum gains may be smaller than 1, which denotes a loss of precision due to the prior. This happens in around 0.75% of cases, so we obtain an improvement in an overwhelming majority of cases.
#### 5.4.2 Perturbed steady solution

We now compare the new basis to the classical one when the initial condition is a perturbed steady solution. To that end, the initial water height is set to

\[h(0,\mathbf{x};\mu)=h_{\text{eq}}(\mathbf{x};\mu)-0.02\exp\left(-2((x_{1}+2)^{2}+(x_{2}+ 2)^{2})\right),\]

thus creating a bump-shaped perturbation whose center is located at \((-2,-2)\). For simplicity, we still use the value of the steady solution as Dirichlet boundary conditions. Moreover, we set the parameters \(\mu\) to be the center of the cube \(\mathbb{P}\), and we take \(q=1\) with \(16^{2}\) discretization cells.

\begin{table}

\end{table}
Table 17: Shallow water equations in two space dimensions: errors, orders of accuracy, and gain obtained when approximating a steady solution for bases with and without prior.

\begin{table}

\end{table}
Table 18: Shallow water equations in two space dimensions: statistics of the gains obtained for the approximation of a steady solution in basis \(V_{h}^{+}\) with respect to basis \(V_{h}\).
The pointwise difference between \(h\) and \(h_{\text{eq}}\) is displayed in Figure 4. We observe that the prior-enriched basis \(V_{h}^{+}\) (right panels) is able to capture the perturbation much better than the classical basis \(V_{h}\) (left panels). Indeed, the background, underlying steady solution has been smeared by basis \(V_{h}\), while is is preserved with much greater resolution by basis \(V_{h}^{+}\).

## 6 Conclusion

In this work, we proposed a Discontinuous Galerkin scheme whose basis has been enriched by neural networks to ensure an approximate well-balance property for a generic PDE and a generic equilibrium. The offline phase of the algorithm consists in learning a family of equilibria using parametric PINNs. Then, during the online phase, the trained network is used to enrich the DG basis and to approximate the solution to the PDE.

The results show significant gains in accuracy compared with the conventional DG method, particularly for low-dimensional approximation spaces. To obtain the same accuracy, we can significantly reduce the number of cells and use larger time steps. The method has been validated on a wide range of PDEs and equilibria, showing that it is a general-purpose approach. Furthermore, it makes it possible to handle complicated equilibria, on complex geometries, which are rarely treated by conventional WB schemes, especially in two space dimensions. The cost of training the network is low, as is the cost of inference. The main additional cost of the method comes from the quadrature rule, whose order has to be increased to ensure a good approximation of the integral of the prior. In most cases, this increase in order is not very important, and the gain between our approach and the classical ones remains significant.

There are several possible ways of extending our approach. From an application point of view, we wish to deal with more difficult equilibria, such as equilibria for the magnetohydrodynamics in tokamaks. From a methodological point of view, we would like to improve the determination of the prior by replacing parametric PINNs with physics-informed neural operators [51, 25] in order to widen the family of equilibria that can be considered. The other approach is to extend the method with time-dependent priors, in order to increase the accuracy of the scheme around families of unsteady solutions. To that end, we wish to move on to space-time DG methods, see e.g. [38].
