# A seamless, extended DG approach for hyperbolic-parabolic problems on
unbounded domains

Federico Vismara111Present address: CASA, Technische Universiteit Eindhoven,
The Netherlands, Tommaso Benacchio, Luca Bonaventura

###### Abstract

We propose and analyze a seamless extended Discontinuous Galerkin (DG)
discretization of hyperbolic-parabolic equations on semi-infinite domains. The
semi-infinite half line is split into a finite subdomain where the model uses
a standard polynomial basis, and a semi-unbounded subdomain where scaled
Laguerre functions are employed as basis and test functions. Numerical fluxes
enable the coupling at the interface between the two subdomains in the same
way as standard single domain DG interelement fluxes. A novel linear analysis
on the extended DG model yields stability constraints on the finite subdomain
grid size that get tighter for increasing values of the Péclet number. Errors
due to the use of different sets of basis functions on different portions of
the domain are negligible, as highlighted in numerical experiments with the
linear advection-diffusion and viscous Burgers’ equations. With an added
damping term on the semi-infinite subdomain, the extended framework is able to
efficiently simulate absorbing boundary conditions without additional
conditions at the interface. A few modes in the semi-infinite subdomain are
found to suffice to deal with outgoing single wave and wave train signals,
thus providing an appealing model for fluid flow simulations in unbounded
regions.

MOX – Modelling and Scientific Computing,
Dipartimento di Matematica, Politecnico di Milano
Via Bonardi 9, 20133 Milano, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Keywords: Laguerre functions; hyperbolic equations; parabolic equations;
discontinuous Galerkin methods; open boundary conditions.

AMS Subject Classification: 65M60, 65M70, 65Z99, 76M10, 76M22.

##  1 Introduction

The correct modelling of evolution problems over arbitrarily large regions has
a wide range of applications in computational physics and poses several still
unsolved challenges. An especially relevant application area is atmospheric
modelling, where the region of interest to forecasts - typically, the
troposphere and lower stratosphere - should not feature spurious reflections
of upwardly propagating waves generated by the computational model lid, see,
e.g., [9, 18, 17]. At the same time, as computational resources enable raising
the lid, an accurate description of upper atmosphere phenomena is of paramount
practical interest and the goal of fully integrated space weather models is
increasingly being discussed, see, e.g., [1, 16].

Approximations over arbitrarily large regions usually rely on the creation of
an artificial boundary separating the region of interest from an external
region. Analytical approaches, e.g., [12, 11, 15], attempt to impose
conditions at the artificial boundary in order to let outgoing perturbations
propagate without spurious reflections. However, these conditions can be
difficult to determine and computationally expensive, and may require _ad hoc_
information on the outgoing signal.

An alternative approach is provided by numerical techniques based on absorbing
(or sponge) layers. These are buffer regions placed next to the artificial
boundary where perturbations leaving the computational domain are damped to a
prescribed external solution by an artificial dissipation term. The choice of
the parameters to be employed in these regions, however, is non-trivial and
the corresponding buffers may be quite large and entail substantial
computational costs. A complete review of the proposed approaches for open
boundary conditions is beyond the scope of this paper. Comprehensive reviews
can be found, for example, in [2, 5, 6, 8, 13, 19].

In [6, 7], a numerical approach to open boundary conditions was proposed,
based on the use of scaled Laguerre functions [21, 22, 23, 25, 27] for
spectral approximations on the semi-infinite line, coupled to finite volume or
finite element discretizations on a finite domain. Different approaches were
used on either side of the artificial boundary, and only the hyperbolic case
was considered. In addition, the stability analysis in [7] only concerned the
discretization on the semi-infinite subdomain.

This paper expands and completes the previously proposed approach by
presenting two major improvements. First, the method is extended to
hyperbolic-parabolic conservation laws, thus showing that all relevant terms
of standard computational fluid dynamics models are amenable to a
discretization based on scaled Laguerre functions. Second, the hyperbolic-
parabolic model equations are discretized in a completely seamless way, using
a discontinuous Galerkin (DG) finite element formulation that relies on scaled
Laguerre functions as both basis and test functions on the semi-infinite
portion, along with Gauss-Laguerre-Radau quadrature rules for numerical
integration.

The resulting extended DG approach and its numerical stability on the entire
half line are analyzed in the case of a linear advection-diffusion equation,
by considering several options for the polynomial basis and quadrature rules.
Scaled Laguerre functions and Gauss-Laguerre-Radau quadrature formulae prove
to be the most stable in all hydrodynamic regimes.

Numerical validations of the proposed method are carried out for the linear
advection-diffusion and the viscous Burgers’ equation. First, a relatively
large number of basis functions are in the semi-infinite portion of the
domain, in order to evaluate the errors introduced by considering different
bases on either side of the finite/semi-infinite interface. By taking as
reference a standard single-domain DG discretization, spurious reflections are
found to be of negligible entity. A reaction damping term is then introduced
in the semi-infinite layer, in order to simulate the propagation of an
isolated Gaussian initial datum and boundary generated wave train from the
finite subdomain into the semi-infinite subdomain, where they are erased using
a decreasing number of Laguerre basis functions. Few basis functions suffice
to absorb outgoing signals with minimal reflections and low computational
cost, thereby providing an efficient implementation of an absorbing layer.

The outline of the paper is as follows. Section 2 contains the model equation
and outlines the numerical discretization. Stability properties of the
proposed method are analyzed in detail in Section 3, and Section 4 contains
the results of the numerical experiments. The final Section 5 draws
conclusions, discussing possible extensions and future work. The Appendix A
summarizes the analysis of different possible discretizations on the semi-
infinite domain, based on either scaled Laguerre functions or polynomials or
on different choices for the numerical quadrature rules. As already shown in
[7], some of these alternatives are shown to be more problematic in the
advection-dominated case.

##  2 The extended DG discretization approach

We consider as a model problem the one-dimensional nonlinear hyperbolic-
parabolic conservation law for the unknown $c$

$\frac{\partial c}{\partial t}+\frac{\partial f(c)}{\partial
z}=\frac{\partial}{\partial z}\left(\mu(z,t)\frac{\partial c}{\partial
z}\right)+s(c,z,t)$ (1)

for $z\in[0,+\infty)$ and $t\in[0,T].$ We assume that the diffusion
coefficient $\mu$ is a smooth function of its variables and that there are two
positive constants $\mu_{0}$ and $\mu_{1}$ such that

$0<\mu_{0}\leq\mu(z,t)\leq\mu_{1}\qquad\forall z\in[0,+\infty),\quad\forall
t>0.$ (2)

For simplicity, we only consider here Dirichlet boundary conditions

$c(0,t)=g_{0}(t)\ \ \ \ \lim_{z\rightarrow+\infty}c(z,t)=0.$ (3)

We refer to [24] for a comprehensive analysis of other boundary conditions in
the linear advection-diffusion case.

Next, drawing on the approach of [6, 7] we split the domain as
${\mathbb{R}}^{+}=[0,L]\cup[L,+\infty)$ and introduce an extended DG finite
element discretization on ${\mathbb{R}}^{+}$ using a standard polynomial basis
on $[0,L]$ and the scaled Laguerre functions as both basis and test functions
on $[L,+\infty)$. More specifically, on the $[0,L]$ interval a mesh of $N$
non-overlapping elements $K_{m}$ of size $\Delta z_{m}\leq h$ is considered,
such that $[0,L]=\bigcup_{m=1}^{N}K_{m}$. The center of the generic element
$K_{m}$ is denoted by $z_{m}$, while $z_{m\pm 1/2}$ denote its boundary
points. The affine local maps $z=Z_{m}(\xi)=\xi\Delta z_{m}/2+z_{m}$ map the
master element $\hat{K}=[-1,1]$ onto each $K_{m}$. For each non-negative
integer $p$, we then denote by $\mathbb{P}_{p}$ the set of all polynomials of
degree less or equal to $p$ on $\hat{K}$. We also define
$\mathbb{P}_{p}(K_{m})=\left\\{w:w=v\circ{Z}^{-1}_{m},\quad
v\in\mathbb{P}_{p}\right\\}$. For each polynomial degree $p$, the
discontinuous finite element spaces are defined as:

$V^{p}_{h}=\left\\{v\in L^{2}([0,L]):v|_{K_{m}}\in\mathbb{P}_{p}(K_{m})\quad
m=1,\dots,N\right\\}.$ (4)

The bases of $\mathbb{P}_{p}(K_{m})$ are obtained from Legendre polynomials as
follows. First, for $\xi\in\hat{K}$, Legendre polynomials are defined by the
recurrence relation:

$\displaystyle L_{k+1}$ $\displaystyle=\frac{2k+1}{k+1}\xi
L_{k}(\xi)-\frac{k}{k+1}L_{k-1}(\xi),$ $\displaystyle k=1,2,\ldots$ (5)
$\displaystyle L_{0}(\xi)$ $\displaystyle=1,\quad L_{1}(\xi)=\xi.$ (6)

Legendre polynomials form an orthogonal basis for polynomials on $\hat{K}$
since

$\int_{-1}^{1}L_{k}(\xi)L_{l}(\xi)d\xi=\frac{2}{2k+1}\delta_{kl}.$ (7)

For each element $K_{m},\,m=1,\dots,N$ we then denote by
$\phi^{m}_{j}(z),j=0,\dots,p$ the basis and test functions given by

$\phi^{m}_{l}(z)=\sqrt{2l+1}L_{l}\Big{(}2\frac{z-z_{m}}{\Delta z_{m}}\Big{)}.$
(8)

Notice that the normalization is chosen so that

$\int_{z_{m-\frac{1}{2}}}^{z_{m+\frac{1}{2}}}\phi^{m}_{k}(z)\phi^{m}_{l}(z)\,dz=\Delta
z_{m}\delta_{kl}.$ (9)

Therefore, the solution of (1) will be represented on each subinterval $K_{m}$
as

$c(z,t)\approx\sum_{j=0}^{p}c^{(j)}_{m}(t)\phi^{m}_{j}(z),\ \ \ \ z\in K_{m}$
(10)

and standard Gauss-Legendre formulae will be used to discretize the resulting
integrals. For the semi-infinite interval $K_{\infty}=[L,+\infty)$, we
consider the scaled Laguerre functions as modal basis. The possible
alternatives are discussed and analyzed in [7] for the purely hyperbolic case
and in Appendix A of this paper for the hyperbolic-parabolic case. More
specifically, defining scaled Laguerre polynomials on $[0,+\infty)$ by

$\displaystyle(k+1)\mathscr{L}^{\beta}_{k+1}(x)$ $\displaystyle=(2k+1-\beta
x)\mathscr{L}^{\beta}_{k}(x)-k\mathscr{L}^{\beta}_{k-1}(x),$ (11)
$\displaystyle\mathscr{L}^{\beta}_{0}(x)$
$\displaystyle=1,\;\mathscr{L}^{\beta}_{1}(x)=1-\beta x,$ (12)

scaled Laguerre functions are defined by

$\hat{\mathscr{L}}_{k}^{\beta}(x)=e^{-\beta x/2}\mathscr{L}_{k}^{\beta}(x).$
(13)

for the scaling factor $\beta>0$, and are a complete orthogonal system in
$L^{2}(\mathbb{R}^{+})$, such that

$\int_{0}^{+\infty}\hat{\mathscr{L}}_{k}^{\beta}(x)\hat{\mathscr{L}}_{l}^{\beta}(x)dx=\frac{1}{\beta}\delta_{kl}.$
(14)

We then define

$\phi^{\infty}_{j}(z)=\hat{\mathscr{L}}^{\beta}_{j}(z-L),\ \ \ j=0,\dots,q$
(15)

for which the analog of (9) holds

$\int_{L}^{+\infty}\phi^{\infty}_{k}(z)\phi^{\infty}_{l}(z)\,dz=\frac{1}{\beta}\delta_{kl}.$
(16)

and we assume that

$c(z,t)\approx\sum_{j=0}^{q}c^{(j)}_{\infty}(t)\phi^{\infty}_{j}(z)\ \ \ \
z\in K_{\infty}.$ (17)

For the resulting integrals, Gauss-Laguerre-Radau quadrature will be employed,
see [6] for the definition. Approximation (17) amounts to say that the
restriction of the numerical approximation of $c$ to $K_{\infty}$ will be
sought in the linear space $V^{q}_{\infty}$ spanned by the functions defined
in (15).

Therefore, the global finite element space employed in the proposed extended
DG discretization can be identified with $V^{p,q}_{h}=V^{p}_{h}\oplus
V^{q}_{\infty}.$ For $v\in V^{p,q}_{h}$ we can then introduce the jump and
average operators as (see, e.g., [4])

$\displaystyle\llbracket
v(z)\rrbracket=v(z^{-})-v(z^{+}),\qquad\\{v(z)\\}=\frac{1}{2}(v(z^{-})+v(z^{+}))$
(18)

and we remark that for $u,v\in V^{p,q}_{h}$ one has

$\llbracket uv\rrbracket=\left\\{u\right\\}\llbracket
v\rrbracket+\left\\{v\right\\}\llbracket u\rrbracket.$ (19)

The extended DG discretization then involves integration of equation (1)
against a test function $v\in V^{p,q}_{h}$, integrating by parts and imposing
for $m=1,\dots,N$ the appropriate continuity constraints at the interelement
boundaries. Setting

$N_{\infty}=\\{1,2,\dots,N,\infty\\},$ (20)

denoting by $\sigma>0$ the stabilization parameter for the DG approximation of
the parabolic terms, and redefining the jump and average operators at $z=0$ so
as to account for the boundary conditions, we obtain the following weak
extended DG formulation of the problem:

$\displaystyle\sum_{m\in N_{\infty}}\int_{K_{m}}\frac{\partial c}{\partial
t}vdz$ $\displaystyle=-\sum_{m=1}^{N}\hat{f}_{m+1/2}\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket$ $\displaystyle+\sum_{m\in
N_{\infty}}\int_{K_{m}}f(c)v^{\prime}dz-\sum_{m\in
N_{\infty}}\int_{K_{m}}\mu\frac{\partial c}{\partial z}v^{\prime}dz$
$\displaystyle+\sum_{m=0}^{N}\left\\{\mu\left(z_{m+1/2}\right)\frac{\partial
c}{\partial z}\left(z_{m+1/2}\right)\right\\}\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket$
$\displaystyle+\sum_{m=0}^{N}\left\\{\mu\left((z_{m+1/2}\right)v^{\prime}\left(z_{m+1/2}\right)\right\\}\left\llbracket
c\left(z_{m+1/2}\right)\right\rrbracket$
$\displaystyle-\sum_{m=0}^{N}\frac{\sigma}{\Delta z_{m}}\left\llbracket
c\left(z_{m+1/2}\right)\right\rrbracket\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket$ $\displaystyle+\sum_{m\in
N_{\infty}}\int_{K_{m}}svdz$
$\displaystyle+\mu(0,t)v^{\prime}(0)g_{0}(t)+f(g_{0}(t))v(0).$ (21)

At $z=L$, the limit from the left of the approximate solution is computed as
$\sum_{j=0}^{p}c^{(j)}_{N}(t)\phi^{N}_{j}(L)$, while the limit from the right
is computed as $\sum_{j=0}^{q}c^{(j)}_{\infty}(t)$, according to the
approximations (10), (17), respectively. Among the many possible formulations
for the parabolic terms, for definiteness we choose that corresponding to the
Symmetric Interior Penalty Galerkin method (SIPG), see, e.g., [3, 26] and the
review in [20]. Furthermore, the Rusanov numerical flux is employed for the
hyperbolic terms, so that

$\displaystyle\begin{split}\hat{f}_{m+1/2}&=\frac{1}{2}\left[f\left(c_{h,m+1/2}^{+}\right)+f\left(c_{h,m+1/2}^{-}\right)\right]\\\
&-\frac{\Lambda_{m+1/2}}{2}\left(c_{h,m+1/2}^{+}-c_{h,m+1/2}^{-}\right),\end{split}$
(22)

where the time dependency is omitted for simplicity,

$\Lambda_{m+1/2}=\text{max}\left(\left\lvert\frac{df}{dc}\left(c_{h,m+1/2}^{+}\right)\right\rvert,\left\lvert\frac{df}{dc}\left(c_{h,m+1/2}^{-}\right)\right\rvert\right),$
(23)

and
$c_{h,m+1/2}^{+}=c_{h}\left(z_{m+1/2}^{+}\right),\,c_{h,m+1/2}^{-}=c_{h}\left(z_{m+1/2}^{-}\right)$.

We can now define the bilinear form

$a:V^{p,q}_{h}\times V^{p,q}_{h}\times[0,+\infty)\rightarrow\mathbb{R}\rm\ \
as$ (24) $\displaystyle\begin{split}a(w,v,t)=&\sum_{m\in
N_{\infty}}\int_{K_{m}}\mu\frac{\partial w}{\partial z}v^{\prime}dz\\\
&-\sum_{m=0}^{N}\left\\{\mu\left(z_{m+1/2},t\right)\frac{\partial w}{\partial
z}\left(z_{m+1/2}\right)\right\\}\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket\\\
&-\sum_{m=0}^{N}\left\\{\mu\left(z_{m+1/2},t\right)v^{\prime}\left(z_{m+1/2}\right)\right\\}\left\llbracket
w\left(z_{m+1/2}\right)\right\rrbracket\\\
&+\sum_{m=0}^{N}\frac{\sigma}{\Delta z_{m}}\left\llbracket
w\left(z_{m+1/2}\right)\right\rrbracket\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket,\end{split}$ (25)

and the nonlinear function $b:V^{p,q}_{h}\times
V^{p,q}_{h}\rightarrow\mathbb{R}$ as

$\displaystyle\begin{split}b(w,v)&=\sum_{m=1}^{N}f\left(w\left(z_{m+1/2}\right)\right)\left\llbracket
v\left(z_{m+1/2}\right)\right\rrbracket\\\ &-\sum_{m\in
N_{\infty}}\int_{K_{m}}f(w)v^{\prime}dz-\sum_{m\in
N_{\infty}}\int_{K_{m}}s(w)vdz.\end{split}$ (26)

We also introduce $g:V^{p,q}_{h}\times V^{p,q}_{h}\to\mathbb{R}$,
$h:V^{p,q}_{h}\times V^{p,q}_{h}\to\mathbb{R}$ as

$\displaystyle g(w,v,t)$
$\displaystyle=-\mu(L,t)v^{\prime}(L)w(L)+\frac{\sigma}{\Delta z_{N}}v(L)w(L)$
(27) $\displaystyle h(w,v)$ $\displaystyle=-f(w(L))v(L)$ (28)

and the linear operator $L:V^{p,q}_{h}\times[0,+\infty)\rightarrow\mathbb{R}$

$L(v,t)=\mu(0,t)v^{\prime}(0)g_{0}(t)+\frac{\sigma}{\Delta
z_{1}}v(0)g_{0}(t)+f(g_{0}(t))v(0),$ (29)

which is related to the Dirichlet condition at the left endpoint $z=0$. The
extended DG weak formulation can then be written more compactly as follows:

For all $t>0$, find $c_{h}(t)\in V^{p,q}_{h}$ such that, $\forall v\in
V^{p,q}_{h}$,

$\displaystyle\begin{split}\int_{0}^{+\infty}\frac{\partial c_{h}}{\partial
t}vdz=&-a(c_{h},v,t)-b(c_{h},v)\\\
&+L(v,t)+g(c_{h},v,t)+h(c_{h},v).\end{split}$ (30)

Approximating $c_{h}$ using (10) and (17), and taking $v=\phi_{j}^{m}$ and
$v=\phi_{k}^{\infty}$, one obtains a set of equations for the discrete degrees
of freedom $c_{m}^{(j)}\ \ j=0,\dots,p,m=1,\dots,N$ $c_{\infty}^{(k)}\
k=0,\dots,q.$ Collecting these in two time-dependent vectors
$\mathbf{c}_{Dg}\in\mathbb{R}^{N(p+1)}$ and
$\mathbf{c}_{Lg}\in\mathbb{R}^{q+1},$ one obtains the systems

$\displaystyle\begin{split}\frac{d\mathbf{c}_{Dg}}{dt}&=\mathbf{A}_{Dg}\mathbf{c}_{Dg}+\mathbf{A}_{Dg,Lg}\mathbf{c}_{Lg}\\\
&+\mathbf{b}_{Dg}(\mathbf{c}_{Dg})+\mathbf{h}_{Dg}(\mathbf{c}_{Dg},\mathbf{c}_{Lg})+\mathbf{g}_{0}(t)\end{split}$
(31)
$\displaystyle\begin{split}\frac{d\mathbf{c}_{Lg}}{dt}&=\mathbf{A}_{Lg}\mathbf{c}_{Lg}+\mathbf{A}_{Lg,Dg}\mathbf{c}_{Dg}\\\
&+\mathbf{b}_{Lg}(\mathbf{c}_{Lg})+\mathbf{h}_{Lg}(\mathbf{c}_{Dg},\mathbf{c}_{Lg}).\end{split}$
(32)

The time-dependent matrices $\mathbf{A}_{Dg}$ and $\mathbf{A}_{Lg}$ result
from the discretization of the diffusion operator in the interior of the
$[0,L]$ and $[L,+\infty)$ subdomains, respectively. The coupling matrices
$\mathbf{A}_{Dg,Lg}$ and $\mathbf{A}_{Lg,Dg}$ result from the discretization
of the diffusion operator involving discrete degrees of freedom of both
subdomains. The nonlinear functions $\mathbf{b}_{Dg}$ and $\mathbf{b}_{Lg}$
result from the discretization of the hyperbolic part and source terms in the
interior of the $[0,L]$ and $[L,+\infty)$ subdomains, respectively. The term
$\mathbf{g}_{0}(t)$ is associated with boundary conditions at $z=0$, while the
coupling nonlinear functions
$\mathbf{h}_{Dg}(\mathbf{c}_{Dg},\mathbf{c}_{Lg})$ and
$\mathbf{h}_{Lg}(\mathbf{c}_{Dg},\mathbf{c}_{Lg})$ result from the
discretization of the hyperbolic part involving discrete degrees of freedom of
both subdomains.
Next, we define the global unknown vector as

$\mathbf{c}(t)=\left(\mathbf{c}_{Dg}(t),\mathbf{c}_{Lg}(t)\right)^{T}\in\mathbb{R}^{N(p+1)+q+1},$
(33)

and the global vectors

$\displaystyle\mathbf{b}\left(\mathbf{c}(t)\right)=\left(\mathbf{b}_{Dg}\left(\mathbf{c}_{Dg}(t)\right),\mathbf{b}_{Lg}\left(\mathbf{c}_{Lg}(t)\right)\right)^{T}\in\mathbb{R}^{N(p+1)+q+1}$
(34)
$\displaystyle\mathbf{h}(\mathbf{c})=\left(\mathbf{h}_{Dg}\left(\mathbf{c}_{Dg},\mathbf{c}_{Lg}\right),\mathbf{h}_{Lg}\left(\mathbf{c}_{Dg},\mathbf{c}_{Lg}\right)\right)^{T}\in\mathbb{R}^{N(p+1)+q+1}$
(35)
$\displaystyle\mathbf{g}(t)=\left(\mathbf{g}_{Dg}(t),0,\dots,0\right)^{T}\in\mathbb{R}^{N(p+1)+q+1}.$
(36)

Defining the global extended DG matrix

$\mathbf{A}(t)=\begin{pmatrix}\mathbf{A}_{Dg}(t)&\mathbf{A}_{Dg,Lg}(t)\\\ \
\mathbf{A}_{Lg,Dg}(t)&\mathbf{A}_{Lg}(t)\end{pmatrix}\in\mathbb{R}^{(N(p+1)+q+1)\times(N(p+1)+q+1)},$
(37)

the extended DG semi-discrete formulation reads

$\frac{d\mathbf{c}(t)}{dt}=\mathbf{A}(t)\mathbf{c}(t)+\mathbf{b}(\mathbf{c}(t))+\mathbf{h}(\mathbf{c}(t))+\mathbf{g}(t).$
(38)

The matrix $\mathbf{A}$ is the discretization of the diffusion term, the
vector $\mathbf{b}$ is the discretization of the non-linear advective part and
the optional source-reaction term, the vector $\mathbf{h}$ contains the flux
exchange at the interface $z=L$ by means of the application of the Rusanov
flux to the flux function $f$, and the vector $\mathbf{g}$ encodes the
Dirichlet condition at the left endpoint $c(0)=a(t)$. We remark that, because
of the vectors $\mathbf{b}$ and $\mathbf{h}$, problem (38) is non-linear.
However, if the functions $f$ and $s$ in (1) are linear, then
$\mathbf{b}(\mathbf{c}(t))$ and $\mathbf{h}(\mathbf{c}(t))$ can be written as
the product between a matrix and the unknown vector $\mathbf{c}(t)$. In this
case, (38) is a linear system of equations.

The semi-discrete extended DG formulation (38) can then be discretized in time
by any standard method for the numerical solution of ODE systems. In this
paper, we use the Crank-Nicolson method for the linear test problems
considered in Section 4. For the non-linear problems, a second order implicit-
explicit (IMEX) method is used, that is described, e.g., in [10, 14]. The
terms associated with the diffusion process are treated implicitly, as they
can entail otherwise rather restrictive stability constraint on the time step
size, while the terms associated with the hyperbolic conservation law are
treated explicitly.

##  3 Stability analysis

In order to study the numerical stability of the global semi-discrete extended
DG formulation (38), we consider the special case of the linear, constant
coefficient, advection-diffusion equation:

$\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial
z}=\mu\frac{\partial^{2}c}{\partial z^{2}}+s(c,z,t)$ (39)

i.e. equation (1) with $f(c)=uc$, $\mu(z,t)\equiv\mu$,
$\mu,u\in\mathbb{R}^{+}$, and in the homogeneous case the source term $s\equiv
0$.

We study the spectrum of the extended DG matrix $\mathbf{A}$ (37) for the
discretization (38) of (39) as a function of the Péclet number
$Pe=u\mathcal{L}/\mu$, where $\mathcal{L}$ is a reference length scale. More
specifically, we fix $\mathcal{L}=1$, $\beta=1$ for the scaling Laguerre
parameter, $\mu=1,\,u=Pe\mu$. We consider a uniform mesh on the finite
subdomain and empirically study under which conditions on the element size
$\Delta z$ the discretization is stable, in the sense that
$\max{(\Re{[\lambda(\mathbf{A})]})}\leq 0.$ This study is carried out for a
range of values of the maximum DG polynomial degree $p$, number of Laguerre
basis functions $q$, and the values $\sigma=20$ and $\sigma=200$ of the
penalization parameter.

Results of the analysis are reported in Table 1. It can be observed that
larger values of $\sigma$ allow for larger sub-intervals without losing
stability. In addition, critical values of the spacing $\Delta z$ are found to
be independent of the number $q$ of basis functions in the semi-infinite
subdomain at fixed polynomial degree $p$ in the finite subdomain. If the
Péclet number is small enough – approximately less than $100$ – the spectrum
of $\mathbf{A}$ has negative real part regardless of the spacing $\Delta z$,
and a few sub-intervals in the finite subdomain $[0,L]$ are enough to make the
scheme stable. By contrast, large values of $Pe$ require a fine grid in
$[0,L]$ to prevent spurious growth caused by eigenvalues with positive real
part. As a consequence, sufficiently large values of $N$ will be considered in
the following numerical validations.

Table 1: Critical values $\Delta z_{cr}$ of the finite subdomain grid spacing
$\Delta z$ for the stability of $\mathbf{A}$ in the extended DG scheme,
$\mu=1$, $u=Pe\mu$, $\sigma=200$ (left), $\sigma=20$ (right).

$q$ | $p$ | $Pe$
---|---|---
100 | 500 | 1000
180 | 3 | 1/2 | 1/4 | 1/9
2 | 1/2 | 1/8 | 1/20
1 | 1/2 | 1/9 | 1/20
90 | 3 | 1/2 | 1/5 | 1/13
2 | 1/2 | 1/9 | 1/21
1 | 1/2 | 1/10 | 1/21
50 | 3 | 1/2 | 1/6 | 1/10
2 | 1/2 | 1/10 | 1/22
1 | 1/2 | 1/10 | 1/21
20 | 3 | 1/2 | 1/5 | 1/8
2 | 1/2 | 1/10 | 1/22
1 | 1/2 | 1/10 | 1/21

$q$ | $p$ | $Pe$
---|---|---
100 | 500 | 1000
180 | 3 | 1/6 | 1/21 | 1/40
2 | 1/7 | 1/35 | 1/73
1 | 1/6 | 1/30 | 1/61
90 | 3 | 1/5 | 1/20 | 1/34
2 | 1/7 | 1/37 | 1/75
1 | 1/5 | 1/30 | 1/61
50 | 3 | 1/4 | 1/18 | 1/30
2 | 1/7 | 1/37 | 1/75
1 | 1/5 | 1/30 | 1/61
20 | 3 | 1/5 | 1/15 | 1/27
2 | 1/8 | 1/38 | 1/76
1 | 1/7 | 1/30 | 1/61

##  4 Numerical experiments

We present here the results of several numerical tests with the extended DG
approach described in the previous sections.

First, a number of validation tests are carried out, considering both linear
and non-linear model problems. The tests assess the accuracy of the extended
DG scheme by comparing it with a stand-alone, single-domain reference DG
discretization on a wider domain. Errors are computed using a relatively large
number of modes in the semi-infinite subdomain of the extended DG scheme, and
the wider domain for the stand-alone reference run also covers part of the
semi-infinite subdomain.

Next, we add a damping reaction term in the semi-infinite subdomain to
simulate an absorbing layer. We show that the extended DG scheme efficiently
damps signals leaving the finite subdomain with negligible reflections into
the finite region.

The experiments consider the linear advection-diffusion equation with constant
coefficients (39), both in the non-homogeneous case ($s\neq 0$) and the
homogeneous case ($s=0$), and the nonlinear, homogeneous viscous Burgers’
equation with constant viscosity, i.e., equation (1) with $f(c)=c^{2}/2$,
$\mu(z,t)\equiv\mu\in\mathbb{R}^{+}$, and $s=0$.

Errors are computed on the finite region $[0,L]$ using a suitable Gaussian
quadrature rule on the sub-intervals $K_{m}$, whose width is $\Delta z=L/N$
for all $m=1,\dots,M$. In particular, we introduce the discrete norms

$\displaystyle\lVert c_{h}\rVert_{L^{2}}$
$\displaystyle=\sqrt{\sum_{m=1}^{N}\frac{\Delta
z}{2}\sum_{k=1}^{ng}\left[c_{h}\left(\frac{\Delta
z}{2}x_{k}+z_{m}\right)\right]^{2}w_{k}}$ (40) $\displaystyle\lVert
c_{h}\rVert_{L^{1}}$ $\displaystyle=\sum_{m=1}^{N}\frac{\Delta
z}{2}\sum_{k=1}^{ng}\left\lvert c_{h}\left(\frac{\Delta
z}{2}x_{k}+z_{m}\right)\right\rvert w_{k}$ (41) $\displaystyle\lVert
c_{h}\rVert_{L^{\infty}}$
$\displaystyle=\max_{m=1,\dots,N}\max_{k=1,\dots,ng}\left\lvert
c_{h}\left(\frac{\Delta z}{2}x_{k}+z_{m}\right)\right\rvert,$ (42)

where $\\{x_{k}\\}_{k=1}^{ng}$ and $\\{w_{k}\\}_{k=1}^{ng}$ are the Gaussian
nodes and weights on the reference interval $[-1,1]$, with $ng$ the number of
quadrature points. Absolute errors with respect to a reference solution are
defined as

$\mathcal{E}_{L^{r}}=\lVert c_{h}-c_{ref}\rVert_{L^{r}},\quad
r\in{1,2,\infty}$ (43)

where $c_{h}$ and $c_{ref}$ are the numerical and the reference solution,
respectively; the latter may be either the exact solution or a single-domain
DG discretization. In some tests we will be also interested in relative errors
with respect to the reference solution $c_{ref}$ defined as

$\mathcal{E}_{L^{r}}^{rel}=\frac{\lVert c_{h}-c_{ref}\rVert_{L^{r}}}{\lVert
c_{ref}\rVert_{L^{r}}},\quad r\in{1,2,\infty}.$ (44)

###  4.1 Validation of the extended DG scheme coupling strategy

We start by testing the proposed method for the linear advection-diffusion
equation with constant coefficients (39). We first consider the non-
homogeneous case by setting $s\neq 0$, assuming an exact solution
$c(z,t)=ze^{-z}\text{sin}^{2}(z-t)$, and computing the right-hand side
analytically. We run the scheme for a variable number of modes $q$ in the
semi-infinite subdomain, also setting $N=100$, $L=2\,\textrm{m}$, $p=2$,
$\mu=1\,\textrm{m}^{2}/\textrm{s}$, $u=2Pe\mu$, and a final time
$T=10\,\textrm{s}$ with $n=200$ time steps. In order to minimize errors at the
interface, the value of $\beta$ is chosen in such a way that the distance
between the first and the second node in the semi-infinite subdomain matches
the grid spacing $\Delta z$ in the finite subdomain. Relative errors $L^{1}$,
$L^{2}$ and $L^{\infty}$ errors at time $T$ with respect to the exact solution
due to the use of different sets of basis functions in the finite and semi-
infinite subdomains are below $5*10^{-6}$ using at least 20 modes, below
$3*10^{-3}$ using 10 modes, and a few percent using 5 modes, thus displaying
spectral convergence in space (Table 2).

Table 2: Relative $L^{1}$, $L^{2}$ and $L^{\infty}$ errors as functions of the number of modes $q$ of the extended DG scheme with respect to the exact solution, linear non-homogeneous advection-diffusion equation. $\Delta t=0.005\,\textrm{s}$, $\Delta z=0.02\,\textrm{m}$, $C=0.25$. $q$ | $\beta$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---|---
5 | 30 | 5.40E-02 | 5.43E-02 | 7.93E-02
10 | 16 | 2.39E-03 | 2.46E-03 | 3.23E-03
20 | 8 | 3.30E-06 | 4.14E-06 | 2.86E-06
40 | 4 | 3.28E-06 | 4.12E-06 | 2.85E-06
80 | 2 | 3.28E-06 | 4.12E-06 | 2.85E-06

In a second test, we examine the accuracy of the extended DG scheme as a
function of the Péclet number, $Pe$. We choose $N=100,\,\beta=1,\,q=180$ that
match the grid spacing across the interface between the two subdomains and
maximize accuracy in the semi-infinite subdomain. Relative errors with respect
to the exact solution are below $4*10^{-4}$ both in the advection-dominated
and in the diffusion-dominated regimes (Table 3). The accuracy appears to
decrease slightly for large values of $Pe$. We assume this to be related to
the stability issues of the modal Laguerre discretization with Neumann
boundary conditions discussed in the Appendix A.

Table 3: Relative $L^{1}$, $L^{2}$ and $L^{\infty}$ errors as functions of $Pe$ of the extended DG scheme with respect to the exact solution, linear non-homogeneous advection-diffusion equation. $\Delta z=0.01\,\textrm{m}$, $q=180$, $\Delta t=0.05\,\textrm{s}$. $Pe$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---
0.001 | 2.65E-04 | 3.16E-04 | 2.23E-04
10 | 4.69E-05 | 4.63E-05 | 5.15E-05
100 | 4.13E-06 | 4.54E-06 | 4.94E-06
500 | 8.28E-06 | 1.01E-05 | 9.20E-06
1000 | 3.81E-05 | 4.50E-05 | 3.83E-05

Next, we validate the extended DG approach in the homogeneous case. We
consider a Gaussian initial datum:

$c_{0}(z)=\text{exp}\left[-\left(\frac{z-z_{c}}{\sigma_{c}}\right)^{2}\right].$
(45)

The interface is located at $L=10\,\textrm{m}$ and the initial hump is placed
inside the bounded interval $[0,L]$ by choosing $z_{c}=8\,\textrm{m}$. The
velocity is $u=1\,\textrm{m}/\textrm{s}$ and the final time is
$T=4\,\textrm{s}$, so that the peak of the Gaussian crosses the interface, and
the other parameters are $\Delta t=0.02\,\textrm{s}$, $q=40$ modes in the
semi-infinite region, $N=500$ sub-intervals for the DG scheme, so that $\Delta
z=0.02\,\textrm{m}$. As the model evolves, the initial hump expands and its
amplitude decreases because of diffusion (Figure 1).

Figure 1: Numerical solution of the linear homogeneous advection-diffusion
equation. Dashed line: initial datum. Solid line: numerical solution with the
extended DG scheme at $T=4\,\textrm{s}$. $z_{c}=8\,\textrm{m}$,
$\sigma_{c}=1\,\textrm{m}$, $q=40$, $\beta=4$, $\Delta z=0.02\,\textrm{m}$,
$\Delta t=0.02\,\textrm{s}$, $\mu=1\,\textrm{m}^{2}/\textrm{s}$,
$u=1\,\textrm{m}/\textrm{s}$.

Relative errors in the finite subdomain $[0,10\,\textrm{m}]$ are computed for
the extended DG scheme with respect to a single-domain DG solution run on
$[0,50\,\textrm{m}]$. For $q=10$ modes in the semi-infinite subdomain,
relative errors are below a few percent, while for $q=40$ they lower to around
$10^{-10}$ (Table 4).

Table 4: Relative $L^{1}$, $L^{2}$ and $L^{\infty}$ errors in $[0,10\,\textrm{m}]$ of the extended DG discretization with respect to a single-domain DG discretization, linear homogeneous advection-diffusion equation. $z_{c}=8\,\textrm{m}$, $\Delta z=0.02\,\textrm{m}$, $\Delta t=0.02\,\textrm{s}$, $C=1$. $q$ | $\beta$ | $\sigma_{c}$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---|---|---
10 | 16 | 1 | 1.90E-02 | 1.11E-02 | 3.79E-02
2 | 1.97E-02 | 1.13E-02 | 4.10E-02
0.5 | 1.87E-02 | 1.11E-02 | 3.70E-02
40 | 4 | 1 | 4.63E-09 | 2.31E-09 | 5.76E-08
2 | 3.70E-09 | 2.83E-09 | 5.63E-09
0.5 | 2.47E-09 | 1.94E-09 | 2.88E-09

In a final validation test, we consider the case of the viscous Burgers’
equation. The initial datum is the Gaussian profile
$c_{0}(z)=\text{exp}\left(-(z-3)^{2}\right)$, with the interface placed at
$z=3\,\textrm{m}$. As time $t$ evolves, the profile moves rightwards
increasing its steepness – using the viscosity value
$\mu=0.05\,\textrm{m}^{2}/\textrm{s}$ no oscillations are observed until the
final time $T=1\,\textrm{s}$. The extended DG scheme compares well with a
reference solution computed by a single-domain DG discretization on
$[0,10\,\textrm{m}]$ (Figure 2).

Figure 2: Numerical solution at $T=1\,\textrm{s}$ of the viscous Burgers’
equation with Gaussian initial datum
$q_{0}(z)=\text{exp}\left(-(z-3)^{2}\right)$, $N=30$, $q=60$, $\Delta
t=10^{-6}\,\textrm{s}$. Black line: extended DG scheme. Red circles: Single-
domain DG discretization.

SSince the closed form of the solution is not available, we compute the errors
with respect to a stand-alone DG discretization on a larger domain with the
same spacing $\Delta z$. The finite/semi-infinite interface in the extended DG
scheme is placed at $L=3\,\textrm{m}$ and the model is run until
$T=10\,\textrm{s}$, with $\Delta t=10^{-2}\,\textrm{s}$. The cases of $N=15$
and $N=30$ subintervals in the finite subdomain are considered, varying the
number of modes in the semi-infinite subdomain, and, accordingly, the scaling
parameter $\beta$ so that the distance between the first two nodes matches the
grid spacing in $[0,L]$. The stand-alone single-domain DG reference solution
for error computation is obtained on the interval $[0,10\textrm{m}]$. A small
number of Laguerre modes are found to suffice to keep the coupling errors in
the bounded subdomain $[0,L]$ below a few percent (Table 5).

Table 5: Relative $L^{1}$, $L^{2}$ and $L^{\infty}$ errors in $[0,3\,\textrm{m}]$ of the extended DG discretization with respect to a single-domain DG discretization, viscous Burgers’ equation. $T=10\,\textrm{s}$, $\Delta t=10^{-2}\,\textrm{s}$, $N=15$. $N$ | $q$ | $\beta$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---|---|---
15 | 10 | 1.6 | 1.98E-02 | 7.80E-03 | 4.07E-02
20 | 0.85 | 2.51E-02 | 1.04E-02 | 5.03E-02
40 | 0.45 | 2.65E-02 | 1.17E-02 | 5.17E-02
80 | 0.23 | 2.62E-02 | 1.18E-02 | 5.07E-02
30 | 10 | 3.6 | 6.81E-04 | 4.60E-04 | 9.18E-04
30 | 1.2 | 5.62E-04 | 3.82E-04 | 7.59E-04
60 | 0.6 | 6.36E-04 | 4.31E-04 | 7.96E-04
100 | 0.36 | 6.76E-04 | 4.61E-04 | 8.77E-05

###  4.2 Efficiency of the extended DG scheme in absorbing layer tests

The second set of tests assesses the performance of the extended DG scheme in
the absorption of perturbations leaving the finite subdomain when an
artificial damping term $-\gamma c,$ with $\gamma\geq 0$, is added to the
model equations’ right-hand side on the the semi-infinite subdomain
$[L,+\infty)$. As in [6, 7] we choose a sigmoid of the form

$\gamma(z)=\frac{\Delta\gamma}{1+\text{exp}\left(\displaystyle\frac{\alpha
L_{0}-z+L}{\sigma_{D}}\right)},$ (46)

where $\Delta\gamma$ is the sigmoid amplitude, $\alpha\in[0,1]$ the position
of the sigmoid inside the absorbing layer, $L_{0}$ the spatial extension of
the semi-infinite region, i.e. the distance between the first and the last
Gauss-Laguerre-Radau nodes, and $\sigma_{D}$ the sigmoid steepness. As in [6,
7] we set $\alpha=0.3$ $m^{-1}$ and $\sigma_{D}=L_{0}/18$.

In a first experiment, we consider the linear advection-diffusion equation
(39) and analyze the damping of a Gaussian profile defined by (45) and
initially placed inside the finite region $[0,L]$. To this end, we place the
interface at $L=1000\,\textrm{m}$, set the initial data parameters
$z_{c}=750\,\textrm{m}$, $\sigma_{c}=50\,\textrm{m}$, and
$\mu=1\,\textrm{m}^{2}/\textrm{s}$ and $u=1\,\textrm{m}/\textrm{s}$. The crest
moves across the finite region, crosses the interface and is damped in the
semi-infinite region. Spurious reflections into the finite region, measured as
absolute errors of the computed solution in the finite region taking the
absence of perturbation as reference, are below $10^{-3}$ for a range of
values for the semi-infinite Laguerre modes $q$, finite subdomain subintervals
$N$, and number of time steps $n$ – and below $10^{-5}$ for the smallest $q=5$
(Table 6).

Table 6: Absolute $L^{2}$, $L^{1}$ and $L^{\infty}$ residual errors in the finite region for the damping of a Gaussian perturbation with the extended DG scheme, linear advection-diffusion equation. $C=0.33$, $p=2$, $T=500\,\textrm{s}$. $q$ | $N$ | $n$ | $\beta$ | $\mathcal{E}_{2}$ | $\mathcal{E}_{1}$ | $\mathcal{E}_{\infty}$
---|---|---|---|---|---|---
40 | 400 | 600 | 1/28 | 9.51E-05 | 1.34E-04 | 1.03E-04
30 | 1/21 | 5.98E-06 | 1.18E-05 | 6.73E-06
20 | 2/29 | 2.58E-05 | 4.25E-05 | 2.70E-05
10 | 2/15 | 1.85E-06 | 6.84E-06 | 1.28E-06
5 | 1/4 | 1.52E-06 | 6.19E-06 | 8.15E-07
30 | 300 | 450 | 1/28 | 3.55E-04 | 4.99E-04 | 3.53E-04
20 | 1/19 | 2.59E-04 | 3.63E-04 | 2.59E-04
10 | 1/10 | 5.00E-06 | 1.18E-05 | 4.59E-06
5 | 11/60 | 1.67E-06 | 6.57E-06 | 1.03E-06
20 | 250 | 375 | 1/23 | 2.59E-04 | 3.63E-04 | 2.59E-04
10 | 1/12 | 5.00E-06 | 1.18E-05 | 4.59E-06
5 | 1/6 | 1.67E-06 | 6.57E-06 | 1.03E-06
10 | 200 | 300 | 1/15 | 4.56E-05 | 6.83E-05 | 3.84E-05
5 | 1/7 | 1.95E-06 | 7.48E-06 | 1.31E-06

Next, we consider a wave train case, obtained by imposing a Dirichlet boundary
condition $c(0,t)=A\text{sin}(2\pi k/Tt)$ at the left endpoint $z=0$. The
initial condition is $c_{0}=0$. The wave train is generated at $z=0$, crosses
the finite region $[0,L]$ and is damped by the absorbing layer, where we set
$\Delta\gamma=2A$. Numerical parameters are set as $L=500\,\textrm{m}$,
$\mu=1\,\textrm{m}^{2}/\textrm{s}$, $u=1\,\textrm{m}/\textrm{s}$,
$T=5000\,\textrm{s}$, and $n=16000$ time steps. On a range of choices for the
wave number, amplitude, and Laguerre modes $q$, the extended DG scheme absorbs
outgoing perturbations with relative errors computed in $[0,L]$ of less than
$10^{-4}$ for $q=15$, and at most $1.3*10^{-3}$ for $q=5$, with respect to a
reference single-domain DG solution on $[0,2L]$ (Figure 3, Tables 7 and 8).
Results are comparable with those obtained in [7] for the inviscid shallow
water system with a different coupling approach. The efficiency of the tool is
competitive – for $q=5$, less than a hundredth for $N=600$, and less than five
thousandths for $N=1200$, of the computational cost is due to the absorbing
layer.

Figure 3: Damping of a wave train with wavenumber $k=30$ (left), $k=60$ (right), linear advection-diffusion equation. Solid line: extended DG scheme. Red circles: single-domain DG scheme. $A=0.1\,\textrm{m}$, $q=30$, $N=600$, $\beta=0.143$, $T=5000\,\textrm{s}$, $n=16000$. Table 7: Damping of a wave train, linear advection-diffusion equation. Relative $L^{2}$, $L^{1}$ and $L^{\infty}$ errors in $[0,L]$ of the extended DG solution with respect to a single-domain DG solution, $q=15$ modes in the semi-infinite subdomain. $L=500\,\textrm{m}$, $p=1$, $\mu=1\,\textrm{m}^{2}/\textrm{s}$, $u=1\,\textrm{m}/\textrm{s}$, $T=5000\,\textrm{s}$, $n=16000$. $A$ | $k$ | $N$ | $\beta$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---|---|---|---
0.025 | 30 | 600 | 0.286 | 1.67E-06 | 2.25E-05 | 1.30E-07
60 | 1200 | 0.571 | 1.10E-07 | 1.18E-08 | 1.31E-06
0.05 | 30 | 600 | 0.286 | 2.25E-06 | 1.78E-07 | 2.98E-05
60 | 1200 | 0.571 | 2.55E-07 | 2.11E-08 | 3.23E-06
0.1 | 30 | 600 | 0.286 | 2.15E-06 | 1.85E-07 | 2.62E-05
60 | 1200 | 0.571 | 4.74E-07 | 3.92E-08 | 5.99E-06
Table 8: Damping of a wave train, linear advection-diffusion equation. Relative $L^{2}$, $L^{1}$ and $L^{\infty}$ errors in $[0,L]$ of the extended DG solution with respect to a single-domain DG solution, $q=5$ modes in the semi-infinite subdomain. $L=500\,\textrm{m}$, $p=1$, $\mu=1\,\textrm{m}^{2}/\textrm{s}$, $u=1\,\textrm{m}/\textrm{s}$, $T=5000\,\textrm{s}$, $n=16000$. $A$ | $k$ | $N$ | $\beta$ | $\mathcal{E}_{2}^{\textrm{rel}}$ | $\mathcal{E}_{1}^{\textrm{rel}}$ | $\mathcal{E}_{\infty}^{\textrm{rel}}$
---|---|---|---|---|---|---
0.025 | 30 | 600 | 0.74 | 9.23E-05 | 6.74E-06 | 1.30E-03
60 | 1200 | 1.48 | 7.57E-06 | 7.57E-07 | 5.74E-05
0.05 | 30 | 600 | 0.74 | 4.23E-05 | 3.11E-06 | 5.91E-04
60 | 1200 | 1.48 | 1.01E-05 | 1.03E-06 | 7.73E-05
0.01 | 30 | 600 | 0.74 | 3.15E-05 | 2.35E-06 | 4.41E-04
60 | 1200 | 1.48 | 1.36E-05 | 1.44E-06 | 9.93E-05

Finally, we consider the Burgers’ equation. We place the interface at
$L=30\,\textrm{m}$, and center an initial Gaussian profile inside the bounded
region ($z_{c}=25\,\textrm{m}$, $\sigma_{c}=1\,\textrm{m}$). We run the
extended DG scheme until $T=3600\,\textrm{s}$, when most of the initial
perturbation has left the bounded region. Residual errors in the finite region
$[0,L]$ with the extended DG scheme with respect to a single-domain DG
solution on $[0,10L/3]$ are below one percent for as few as $q=5$ modes in the
semi-infinite subdomain (Table 9). Qualitatively equivalent results are
obtained when placing the initial data at the interface,
$z_{c}=30\,\textrm{m}$ (not shown).

Table 9: Damping of a Gaussian profile initially centred at $z_{c}=25\,\textrm{m}$, viscous Burgers’ equation. Relative $L^{2}$, $L^{1}$ and $L^{\infty}$ errors in $[0,L]$ of the extended DG scheme compared with a single-domain DG solution. $T=3600\,\textrm{s}$, $\Delta t=0.1\,\textrm{s}$, $\Delta\gamma=2$, $\sigma_{c}=1\,\textrm{m}$. $q$ | $\beta$ | $\mathcal{E}^{rel}_{2}$ | $\mathcal{E}^{rel}_{1}$ | $\mathcal{E}^{rel}_{\infty}$
---|---|---|---|---
60 | 0.06 | 2.25E-03 | 1.15E-02 | 5.54E-04
40 | 0.09 | 2.24E-03 | 1.14E-02 | 5.50E-04
20 | 0.175 | 2.17E-03 | 1.10E-02 | 5.34E-04
10 | 0.34 | 1.96E-03 | 9.95E-03 | 4.83E-04
5 | 0.68 | 1.54E-03 | 7.76E-03 | 3.84E-04

##  5 Conclusions and perspectives

This paper proposed an extended DG approach for the numerical simulation of
nonlinear hyperbolic-parabolic problems on unbounded domains. Built on earlier
developments of coupled DG-Laguerre discretizations for purely hyperbolic
systems, the scheme models a finite portion of the semi-infinite half-line
using standard Legendre basis functions and the adjacent unbounded portion
using scaled Laguerre basis functions.

Compared to a standard DG discretization, the extended DG scheme only differs
for the presence of two off-diagonal terms in the system matrix, representing
the numerical fluxes at the finite/semi-infinite interface. The new framework
improves on previous endeavours that used bespoke coupling strategies, and
provides a completely seamless coupling approach.

A linear analysis of the extended DG scheme in terms of the Péclet number
identified stability thresholds on the grid spacing in the finite subdomain.
The constraints become more stringent for increasingly advection-dominated
flows and smaller values of the penalization parameter of the DG
discretization in the finite subdomain. The analysis and numerical experiments
used scaled Laguerre basis functions, Gauss-Laguerre-Radau quadrature in the
unbounded subdomain, and Dirichlet boundary conditions. Results covering other
possible options are reported in [24] and are summarized in the Appendix,
corroborating the findings for the standalone Laguerre scheme in a purely
hyperbolic framework [7]. To the best of the authors’ knowledge, a stability
analysis on numerical schemes using different sets of basis functions as the
one presented in this work is not currently available in the literature.

The correctness of the extended DG scheme, particularly regarding the
finite/semi-infinite interface fluxes, was validated in a series of numerical
experiments with the linear homogeneous and non-homogeneous advection-
diffusion equation and the nonlinear viscous Burgers’ equation. By comparison
with a standard single domain implementation, spurious signals due to the
presence of different basis functions are of negligible entity on a range of
spatial resolutions, thereby complementing and strengthening results obtained
with hyperbolic systems in [6, 7].

In tests where the semi-infinite subdomain featured a reactive damping term,
the extended DG scheme displayed compelling performance in efficiently
absorbing outgoing waves in linear and nonlinear models. A very small number
of Laguerre modes, both in absolute terms and as a proportion of the total
computational load, was sufficient to damp single Gaussian signals and wave
trains without spurious phenomena spoiling the simulation in the finite
subdomain. Relative errors with respect to a single-domain DG scheme are small
enough to make the extended DG scheme an interesting technique for the
discretization of fluid dynamics problems on unbounded domains.

The results achieved in this work offer a number of perspectives for future
investigation. First, a similar extended DG approach can be developed coupling
a strong form, nodal DG discretization on the finite domain to the strong
form, nodal approach with scaled Laguerre functions for the semi-infinite
domain, a choice that displayed stability advantages in the large Péclet
number case. The scheme can then be implemented in multiple dimensions, using
tensor product-based discretization approaches on semi-infinite strips or
circular domains, where the problem is discretized using the extended DG
scheme in the vertical or radial direction and a discontinuous Galerkin
approach in the horizontal or azimuthal direction. Such a model may find
applications, for example, in the modelling of the solar corona. The extension
to systems of parabolic equations or to non-linear diffusion may be
considered, such as are found in turbulence modelling. From a more theoretical
perspective, the possibility to prove inf-sup conditions for the extended DG
approach could also be investigated.

## Acknowledgements

This work summarizes and extends results obtained by the first author (F.V.)
in his Master’s Thesis in Mathematical Engineering [24], discussed at
Politecnico di Milano in 2020 and prepared under the supervision of the other
two authors. T.B. and L.B. have been supported by the ESCAPE-2 project of the
Horizon 2020 research and innovation programme (grant agreement No 800897).

## Note

This work has not yet been peer-reviewed and is provided by the contributing
authors as a means to ensure timely dissemination of scholarly and technical
work on a noncommercial basis. Copyright and all rights therein are maintained
by the authors or by other copyright owners. It is understood that all persons
copying this information will adhere to the terms and constraints invoked by
each author’s copyright. This work may not be reposted without explicit
permission of the copyright owner.

##  Appendix A Alternative discretizations on the semi-infinite subdomain

We summarize here the results presented in [24] on the analysis of various
Laguerre-based discretizations of the advection-diffusion equation with
constant coefficients on ${\mathbb{R}}^{+}=[0,+\infty)$. For the purpose of
deriving some discretizations, it can be helpful to reformulate equation (39),
which we report here for convenience,

$\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial
z}=\mu\frac{\partial^{2}c}{\partial z^{2}}$ (47)

as a system of first order equations

$\displaystyle\begin{split}&\frac{\partial c}{\partial t}-\mu\frac{\partial
v}{\partial z}+uv=0\\\ &\frac{\partial c}{\partial z}-v=0.\end{split}$ (48)

We assume that solutions vanish at infinity

$\lim_{z\rightarrow+\infty}c(z,t)=0$ (49)

and that either Dirichlet boundary conditions

$c(0,t)=c_{L}$ (50)

or Neumann boundary conditions

$\frac{\partial c}{\partial z}(0,t)=Dc_{L}$ (51)

are applied at $z=0.$ We require that $\mu>0$ (ellipticity condition) and
$u>0$. In this case, the Dirichlet datum at $z=0$ corresponds to an inflow
boundary condition, which guarantees well-posedness for the hyperbolic part.
We analyze several possible space discretizations, in order to determine which
one shows the best stability properties and can therefore be chosen for the
extended DG scheme in conjunction with the Legendre basis in the finite sub-
domain. As done in [7] for the pure advection problem, we discretize the PDE
system (48) in space, obtaining, after substitution of the discretization of
the second equation in (48) into the first, a system of ordinary differential
equations of the form

$\frac{d\mathbf{c}}{dt}=\mathbf{A}\mathbf{c}+\mathbf{g},$ (52)

where $\mathbf{c}$ is the unknown vector of the expansion of the solution and
$\mathbf{g}$ contains the contribution of boundary conditions at $z=0$, and we
study the eigenvalue structure of the matrix $\mathbf{A}$. The corresponding
discretization scheme is stable if all the eigenvalues have non-positive real
part.
We analyse the following discretizations:

  * •

_Weak form_. We multiply (48) by a test function, integrate by parts and use
either Gauss-Laguerre-Radau (GLR) or Gauss-Laguerre (GL) quadrature rules. Two
different approaches are possible. In a modal approach, entries of the unknown
vector $\mathbf{c}$ are the coefficients of the expansion of the solution in
the orthogonal basis of Laguerre functions or Laguerre polynomials. In a nodal
approach, the basis functions are the Lagrange basis functions associated with
the integration nodes, so that the unknown vector contains the nodal values of
the approximate solution.

Furthermore, the numerical solution can be expanded in a basis of either
scaled Laguerre functions or scaled Laguerre polynomials.

  * •

_Strong form_. In this case we directly discretize the strong formulation (48)
using a collocation approach and GLR quadrature rules. This is the only
practical choice if Dirichlet boundary conditions have to be imposed, because
the GLR nodes include the left endpoint of the semi-infinite subdomain, unlike
the GL nodes.

We now summarize some definitions we need to introduce the different variants
of the matrix $\mathbf{A}$ and vector $\mathbf{g}.$

For discretizations based on Laguerre functions, we define the matrix
$\mathbf{\hat{L}}=\\{\hat{l}_{ij}\\}$ with entries such that

$\hat{l}_{ij}=\begin{cases}1/2&i=j\\\ 1&j<i\\\ 0&j>i.\end{cases}$ (53)

If discretizations based on Laguerre polynomials are considered instead, we
use the matrix $\mathbf{L}=\\{l_{ij}\\}$ defined as

$l_{ij}=\begin{cases}0&i=j\\\ 1&j<i\\\ 0&j>i\end{cases}$ (54)

For nodal discretizations based on the weak form and on scaled Laguerre
functions, we then denote by $z_{j}^{\beta}$ the $j$-th GLR or GL quadrature
node, by $h_{j}^{\beta}(z)$ the associated Lagrangian polynomial, by
$\omega_{i}$ the $i$-th quadrature weight, and by $\hat{d}_{ij}^{\beta}$ the
entries of the GLR or GL differentiation matrix $\mathbf{\hat{D}}_{\beta}$
associated with scaled Laguerre functions, defined as follows:

  * •

GL nodes

$\hat{d}^{\beta}_{ij}=\begin{cases}\dfrac{\hat{\mathscr{L}}_{M}^{\beta}(z_{i}^{\beta})}{(z_{i}^{\beta}-z_{j}^{\beta})\hat{\mathscr{L}}_{M}^{\beta}(z_{j}^{\beta})}&i\neq
j\\\ \\\ -\dfrac{M+2}{2z_{i}^{\beta}}&i=j\\\ \end{cases}$ (55)

  * •

GLR nodes

$\hat{d}^{\beta}_{ij}=\begin{cases}\dfrac{\hat{\mathscr{L}}_{M+1}^{\beta}(z_{i}^{\beta})}{(z_{i}^{\beta}-z_{j}^{\beta})\hat{\mathscr{L}}_{M+1}^{\beta}(z_{j}^{\beta})}&i\neq
j\\\ \\\ 0&i=j\neq 0\\\ \\\ -\beta\dfrac{M+1}{2}&i=j=0.\\\ \end{cases}$ (56)

We also define as $\bm{\hat{\Omega}}_{\beta}$ the diagonal matrix with the
quadrature weights $\hat{\omega}_{i}^{\beta}$ on the diagonal. For a nodal
discretization based on Laguerre polynomials, instead, the differentiation
matrix $\mathbf{D}_{\beta}$ has entries $d^{\beta}_{ij}$ defined as:

  * •

GL nodes

${d}_{ij}^{\beta}=\begin{cases}\dfrac{{\mathscr{L}}_{M}^{\beta}(z_{i}^{\beta})}{(z_{i}^{\beta}-z_{j}^{\beta}){\mathscr{L}}_{M}^{\beta}(z_{j}^{\beta})}&i\neq
j\\\ \\\ \dfrac{\beta z_{i}^{\beta}-M-2}{2z_{i}^{\beta}}&i=j\\\ \end{cases}$
(57)

  * •

GLR nodes

${d}_{ij}^{\beta}=\begin{cases}\dfrac{{\mathscr{L}}_{M+1}^{\beta}(z_{i}^{\beta})}{(z_{i}^{\beta}-z_{j}^{\beta}){\mathscr{L}}_{M+1}^{\beta}(z_{j}^{\beta})}&i\neq
j\\\ \\\ \dfrac{\beta}{2}&i=j\neq 0\\\ \\\ -\beta\dfrac{M}{2}&i=j=0\\\
\end{cases}$ (58)

We also set

$\displaystyle\mathbf{\hat{g}_{1}}=[(\hat{h}_{0}^{\beta})^{\prime\prime}(z_{1}),\dots,(\hat{h}_{0}^{\beta})^{\prime\prime}(z_{M})],\;\mathbf{\hat{g}_{2}}=[(\hat{h}_{0}^{\beta})^{\prime}(z_{1}),\dots,(\hat{h}_{0}^{\beta})^{\prime}(z_{M})],$
$\displaystyle\mathbf{g_{1}}=[(h_{0}^{\beta})^{\prime\prime}(z_{1}),\dots,(h_{0}^{\beta})^{\prime\prime}(z_{M})],\;\mathbf{g_{2}}=[(h_{0}^{\beta})^{\prime}(z_{1}),\dots,(h_{0}^{\beta})^{\prime}(z_{M})],$
$\displaystyle\mathbf{\hat{h}}=[\hat{h}^{\beta}_{0}(0),\dots,\hat{h}^{\beta}_{M}(0)],\;\mathbf{h}=[h^{\beta}_{0}(0),\dots,h^{\beta}_{M}(0)],$
$\displaystyle\mathbf{\widehat{W}}=\bm{\hat{\Omega}}_{\beta}^{-1}\mathbf{\hat{D}}_{\beta}^{T}\bm{\hat{\Omega}}_{\beta},\;\mathbf{W}=\bm{\Omega}_{\beta}^{-1}\mathbf{D}_{\beta}^{T}\bm{\Omega}_{\beta},$
$\displaystyle\mathbf{\hat{r}}=\bm{\hat{\Omega}}_{\beta}^{-1}\mathbf{\hat{h}},\;\mathbf{r}=\bm{\Omega}_{\beta}^{-1}\mathbf{h},$
$\displaystyle\mathbf{e}=[1,\dots,1]^{T}\in\mathbf{R}^{M+1}.$

We also denote by $(\mathbf{\hat{D}}_{\beta})_{M}$ for scaled Laguerre
functions, and by $(\mathbf{D}_{\beta})_{M}$ for scaled Laguerre polynomials,
the matrices obtained from the differentiation matrices
$\mathbf{\hat{D}}_{\beta}$ and $\mathbf{D}_{\beta}$ by removing the first row
and the first column. Finally we denote by $(\mathbf{\hat{D}}_{\beta})_{0}$
for scaled Laguerre functions, and $(\mathbf{D}_{\beta})_{0}$ for scaled
Laguerre polynomials, the matrices obtained from $\mathbf{\hat{D}}_{\beta}$
and $\mathbf{D}_{\beta}$ by replacing the first row with zeros.

The expressions of matrix $\mathbf{A}$ and right-hand side $\mathbf{g}$ for
the derived discretizations are summarized in Table 10 – note the two use of
the matrix $\mathbf{\hat{L}}$ (53) for scaled Laguerre functions and
$\mathbf{L}$ (54) for scaled Laguerre polynomials.

Table 10: Definition of matrix $\mathbf{A}$ and vector $\mathbf{g}$ in the Laguerre discretization of (39) or (48) for several formulations ‘Form’, basis functions ‘BF’ and boundary conditions ‘BC’. ‘Coll’: collocation, ‘Nod’: nodal, ‘Mod’: modal, ‘Dir’: Dirichlet, ‘Neu’: Neumann, ‘LF’: Scaled Laguerre Functions, ‘LP’: Scaled Laguerre Polynomials. See text for symbol definitions. Form | BF | BC | $\mathbf{A}$ | $\mathbf{g}$
---|---|---|---|---
Coll | LF | Dir | $\mu(\mathbf{\hat{D}}_{\beta}^{2})_{M}-u(\mathbf{\hat{D}}_{\beta})_{M}$ | $\mu c_{L}\mathbf{\hat{g}_{1}}-uc_{L}\mathbf{\hat{g}_{2}}$
Coll | LF | Neu | $\mu\mathbf{\hat{D}}_{\beta}(\mathbf{\hat{D}}_{\beta})_{0}-u(\mathbf{\hat{D}}_{\beta})_{0}$ | $\mu Dc_{L}\mathbf{\hat{g}_{2}}-uDc_{L}\mathbf{e_{1}}$
Coll | LP | Dir | $\mu(\mathbf{{D}}_{\beta}^{2})_{M}-u(\mathbf{{D}}_{\beta})_{M}$ | $\mu c_{L}\mathbf{g_{1}}-uc_{L}\mathbf{g_{2}}$
Coll | LP | Neu | $\mu\mathbf{{D}}_{\beta}(\mathbf{{D}}_{\beta})_{0}-u(\mathbf{{D}}_{\beta})_{0}$ | $\mu Dc_{L}\mathbf{g_{2}}-uDc_{L}\mathbf{e_{1}}$
Nod | LF | Dir | $-\mu\mathbf{\hat{D}}_{\beta}\mathbf{\hat{W}}+u\mathbf{\hat{W}}$ | $-\mu c_{L}\mathbf{\hat{D}}_{\beta}\mathbf{\hat{r}}+uc_{L}\mathbf{\hat{r}}$
Nod | LF | Neu | $-\mu\bm{\hat{\Omega}_{\beta}}^{-1}\mathbf{\hat{D}}_{\beta}^{T}\bm{\hat{\Omega}}_{\beta}\mathbf{\hat{D}}_{\beta}-u\mathbf{\hat{D}}_{\beta}$ | $-\mu Dc_{L}\mathbf{\hat{r}}$
Nod | LP | Dir | $-\mu\mathbf{D}_{\beta}\mathbf{W}+\mu\beta\mathbf{D}_{\beta}+u\mathbf{W}-u\beta\mathbf{I}$ | $-\mu c_{L}\mathbf{D}_{\beta}\mathbf{r}+uc_{L}\mathbf{r}$
Nod | LP | Neu | $-\mu\mathbf{W}\mathbf{D}_{\beta}+\mu\beta\mathbf{D}_{\beta}-u\mathbf{D}_{\beta}$ | $-\mu Dc_{L}\mathbf{r}$
Mod | LF | Dir | $-\mu\beta^{2}\mathbf{\hat{L}}^{T}\mathbf{\hat{L}}-u\beta\mathbf{\hat{L}}$ | $\mu\beta^{2}c_{L}\mathbf{\hat{L}}^{T}\mathbf{e}+u\beta c_{L}\mathbf{e}$
Mod | LF | Neu | $-\mu\beta^{2}\mathbf{\hat{L}}\mathbf{\hat{L}}^{T}+u\beta\mathbf{\hat{L}}^{T}$ | $-\mu\beta Dc_{L}\mathbf{e}$
Mod | LP | Dir | $-\mu\beta^{2}\mathbf{L}^{T}(\mathbf{L}+\mathbf{I})-u\beta(\mathbf{L}+\mathbf{I})$ | $u\beta c_{L}\mathbf{e}$
Mod | LP | Neu | $-\mu\beta^{2}(\mathbf{L}+\mathbf{I})\mathbf{L}^{T}+u\beta\mathbf{L}^{T}$ | $-\mu\beta Dc_{L}\mathbf{e}$

As customary for the advection-diffusion problem, the stability property can
be a function of the Péclet number, which is usually defined as
$Pe=u\mathcal{L}/\mu$, where $\mathcal{L}$ is a reference length scale. For
simplicity we choose the length scale $\mathcal{L}=1$, set $\mu=1$ and analyze
the stability of $\mathbf{A}$ for a fixed value of $Pe$; the corresponding
ranges for $\beta$ are shown in Table 11 for both scaled Laguerre functions
and polynomials.

Table 11: Stability of $\mathbf{A}$ as a function of $\beta$: condition under which the largest real part of the eigenvalues is non-positive. $q=50$, $\mu=1$. ‘Neu’: Neumann b.c., ‘Dir’: Dirichlet b.c., ‘LF’: Scaled Laguerre Functions, ‘LP’: Scaled Laguerre Polynomials. |  |  |  LF |  LP
---|---|---|---|---
|  |  | Neu | Dir | Neu | Dir
Strong |  |  | $\forall\beta$ | $\forall\beta$ | $\beta\leq 2.6Pe$ | $\beta\leq 3Pe$
Weak | Nodal | GLR | $\beta\geq 0.58Pe$ | $\forall\beta$ | $0.017Pe\leq\beta\leq 2.83Pe$ | $\beta\leq 3Pe$
GL | $\beta\geq 2Pe$ | $\forall\beta$ | $0.25Pe\leq\beta\leq 2Pe$ | $\beta\leq 8.5Pe$
Modal |  | $\beta\geq 0.58Pe$ | $\forall\beta$ | $0.017Pe\leq\beta\leq 2.83Pe$ | $\beta\leq 3Pe$

It can be observed that only the strong form discretizations based on Laguerre
functions are stable for all boundary conditions and independently of the
value of the Péclet number. Other discretizations based on Laguerre functions
are instead stable under mild conditions on the value of $\beta$ as a function
of the Péclet number. These conditions become problematic only in the very
large Péclet number limit.

In this paper, only the weak form modal discretization based on Laguerre
functions was considered for the extended DG scheme, due to its hierarchical
nature, that allows in principle for an easy (and if necessary, dynamic)
adjustment of the number of basis functions to perform $p-$adaptation. The
strong form nodal discretization based on Laguerre functions seems otherwise
the most robust option and will be further studied as a basis for extended DG
approaches in future work. Discretizations based on Laguerre polynomials are
instead only stable under more restrictive conditions, which also affect the
choice of $\beta$ in the small Péclet number case. These conclusions
complement the results in [7], where the pure advection problem was discussed.
Such an analysis does not seem to have been carried out in the literature, to
the best of the authors’ knowledge.

## References

  * [1] R.A. Akmaev.  Whole atmosphere modeling: Connecting terrestrial and space weather.  Reviews of Geophysics, 49, 2011.
  * [2] D. Appelö and T. Colonius.  A high-order super-grid-scale absorbing layer and its application to linear hyperbolic systems.  Journal of Computational Physics, 228(11):4200–4217, 2009.
  * [3] D.N. Arnold.  An interior penalty finite element method with discontinuous elements.  SIAM Journal of Numerical Analysis, 19:742–760, 1982.
  * [4] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini.  Unified analysis of discontinuous Galerkin methods for elliptic problems.  SIAM Journal of Numerical Analysis, 39:1749–1779, 2002.
  * [5] R.J. Astley.  Infinite elements for wave problems: a review of current formulations and an assessment of accuracy.  International Journal of Numerical Methods in Engineering, 49(7):951–976, 2000.
  * [6] T. Benacchio and L. Bonaventura.  Absorbing boundary conditions: a spectral collocation approach.  International Journal of Numerical Methods in Fluids, 72(9):913–936, 2013.
  * [7] T. Benacchio and L. Bonaventura.  An extension of DG methods for hyperbolic problems to one-dimensional semi-infinite domains.  Applied Mathematics and Computation, 350:266–282, 2019.
  * [8] K. Black.  Spectral elements on infinite domains.  SIAM Journal of Scientific Computing, 19:1667–1681, 1998.
  * [9] L. Bonaventura.  A Semi–implicit Semi–Lagrangian Scheme Using the Height Coordinate for a Nonhydrostatic and Fully Elastic Model of Atmospheric Flows.  Journal of Computational Physics, 158(2):186–213, 2000.
  * [10] L. Bonaventura, E.D. Fernández-Nieto, J. Garres-Díaz, and G. Narbona-Reina.  Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization.  Journal of Computational Physics, 364:209–234, 2017.
  * [11] J.R. Dea.  An experimental adaptation of Higdon-type non-reflecting boundary conditions to linear first-order systems.  Journal of Computational and Applied Mathematics, 235:1354–1366, 2011.
  * [12] B. Engquist and A. Majda.  Absorbing boundary conditions for numerical simulation of waves.  Mathematics of Computation, 31(139):629–651, 1977.
  * [13] K. Gerdes.  A review of infinite element methods for exterior Helmholtz problems.  Journal of Computational Acoustics, 8(1):43–62, 2000.
  * [14] F.X. Giraldo, J.F. Kelly, and E.M. Constantinescu.  Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA).  SIAM Journal on Scientific Computing, 35, 2013.
  * [15] M. Israeli and S.A. Orszag.  Approximation of Radiation Boundary Conditions.  Journal of Computational Physics, 41:115–135, 1981.
  * [16] D.R. Jackson, T.J. Fuller-Rowell, D.J. Griffin, M. J. Griffith, C. W. Kelly, D. R. Marsh, and M. T. Walach.  Future directions for whole atmosphere modeling: Developments in the context of space weather.  Space Weather, 17:1342–1350, 2019.
  * [17] J.B. Klemp and D.R. Durran.  An Upper Boundary Condition Permitting Internal Gravity Wave Radiation in Numerical Mesoscale Models.  Journal of Atmospheric Sciences, 111:430–444, 1983.
  * [18] J.B. Klemp and D.K. Lilly.  Numerical Simulation of Hydrostatic Mountain Waves.  Journal of Atmospheric Sciences, 35:78–107, 1978.
  * [19] P.J. Rasch.  Toward atmospheres without tops: Absorbing upper boundary conditions for numerical models.  Quarterly Journal of the Royal Meteorological Society, 112:1195–1218, 1986.
  * [20] B. Rivière.  Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation.  SIAM, 2008.
  * [21] J. Shen.  Stable and efficient spectral methods in unbounded domains using Laguerre functions.  SIAM Journal on Numerical Analysis, 38:1113–1133, 2001.
  * [22] J. Shen, T. Tang, and L.-L. Wang.  Spectral Methods: Algorithms, Analysis and Applications, volume 40 of Springer Series in Computational Mathematics.  Springer, 2011.
  * [23] J. Shen and L.-L. Wang.  Some recent advances on spectral methods for unbounded domains.  Communications in Computational Physics, 5:195–241, 2009.
  * [24] F. Vismara.  A coupled scheme for the solution of parabolic problems on unbounded domains.  Master’s thesis, Politecnico di Milano, 2020.  Available at: https://www.politesi.polimi.it/handle/10589/166694 (last accessed 7 December 2020).
  * [25] Z.-Q. Wang, B.-Y. Guo, and Y.-N. Wu.  Pseudospectral method using generalized Laguerre functions for singular problems on unbounded domains.  Discrete and Continuous Dynamical Systems Series B, 11(4):1019–1038, 2009.
  * [26] M.F. Wheeler.  An elliptic collocation-finite element method with interior penalties.  SIAM Journal of Numerical Analysis, 15:152–161, 1978.
  * [27] Q. Zhuang, J. Shen, and C. Xu.  A coupled Legendre–Laguerre spectral–element method for the Navier–Stokes equations in unbounded domains.  Journal of Scientific Computing, 42(1):1–22, 2010.


# Size distribution of emitted energies in local load sharing fiber bundles

Subhadeep Roy<EMAIL_ADDRESS>PoreLab, Department of Physics, Norwegian
University of Science and Technology, NO–7491 Trondheim, Norway.  Soumyajyoti
Biswas<EMAIL_ADDRESS>Department of Physics, SRM University - AP,
Andhra Pradesh 522502, India.

###### Abstract

We study the local load sharing fiber bundle model and its energy burst
statistics. While it is known that the avalanche size distribution of the
model is exponential, we numerically show here that the avalanche size ($s$)
and the corresponding energy burst ($E$) in this version of the model have a
non-linear relation ($E\sim s^{\gamma}$). Numerical results indicate that
$\gamma\approx 2.5$ universally for different failure threshold distributions.
With this numerical observation, it is then possible to show that the energy
burst distribution is a power law, with a universal exponent value of
$-(\gamma+1)$.

##  I Introduction

It is well known experimentally that quasistatically stressed disordered
solids produce intermittent response statistics phys_rep , particularly in
terms of acoustic emissions, that show scale-free size distributions. These
intriguing dynamics is seen universally across scales from microscopic
laboratory samples to the geological scale of earthquakes herrmann90 ; bkc96 ;
rmp_2012 ; wiley_book ; lucilla . Empirically, the scale-free size
distribution of breaking progression is known in different communities for
decades. For example, in geoscience, this is known as the Gutenberg-Richter
law, in magnetic domain walls as crackling noise, and so on.

The interests of statistical physicists in this context stems from the
universal nature of the dynamics across length and energy scales. The scale-
free variations of acoustic emissions, waiting time statistics, etc., are
independent of the microscopic details of the underlying systems, which are
very different from each other. Such behavior indicates critical dynamics,
particularly self-organized critical dynamics for the system, where the
universality hypothesis is still applicable, without having to fine-tune a
driving parameter bonamy . Such a phenomenon is therefore open for analysis
with the tools of critical phase transitions, universality and therefore is an
important step towards predictability of imminent failure alava_lasse ;
kun_record ; sci_rep .

As a consequence of the scale-free dynamics and potential applicability of the
universality hypothesis, many generic models were proposed over the years that
reproduce such a scale-free behavior. Such models include the fiber bundle
model, the random fuse model, the Burridge-Knopoff model, and so on wiley_book
; fbm_rmp ; fbm_book . The common underlying feature of these models is that
they are threshold activated, driven, dynamical models. Particularly, for an
external driving parameter crossing a pre-assigned threshold value for a
single unit of these models, that unit is activated and influences the units
in its ‘neighborhood’, which may in-turn get activated and thereby initiating
an ‘avalanche’. As can be guessed, this type of dynamics is often related to
sandpile models of self-organized criticality soc  and indeed such
associations extensively explored in the past soc_fbm .

The two major parameters that influence the nature of the response in such
models are the range of interaction and the strength of the disorder. It was
explored, particularly in the fiber bundle model that for a moderate disorder,
a scale-free avalanche statistics is only recovered for a ‘sufficiently’ long-
range of interaction pre1 ; pre2 ; kun_ptrc . In the random fuse model, where
the interaction range is not parameters to be tuned, it was shown that the
avalanche statistics is not a power-law in the large system size limit sch .
This is in apparent contradiction with the fact that in reality, the
interaction range in disordered elastic samples is not infinite i.e., not a
mean-field-like interaction. However, experiments routinely reveal scale-free
statistics.

One important distinction between the analytical and numerical results of
avalanche dynamics and that of the experiments is that in the former it is the
number of elements failing in an avalanche that is the measurable quantity,
while in the latter it is the energy released in the avalanche. Now, in the
mean-field limit of the fiber bundle model, it is straightforward to show that
the avalanche size and the energy avalanche size are proportional, hence the
two distributions are identical in shape. But this relation is no longer valid
for local load sharing variants. In those cases, therefore, it is crucial to
explore the size distributions of the energy emissions and compare that with
experiments. In this work, we consider the simplest possible variant of the
local load sharing fiber bundle model and analyze the energy avalanche
statistics of that model. We then compare the results with experiments and
also present a plausible argument for its form.

##  II Description of Fiber Bundle Model

After being introduced by Pierce in 1926 Pierce , the fiber bundle model has
been proven to be important yet arguably the simplest model to study failure
processes in disordered solids. A conventional fiber bundle model consists of
a set of linear elastic fibers or Hookean springs, attached between two
parallel plates. The plates are pulled apart by a force $F$, creating a stress
$\sigma=F/L$ on $L$ fibers. Once the stress crosses the breaking threshold of
a particular fiber, chosen from a random distribution, that fiber breaks
irreversibly. The stress of broken fibers is then redistributed either
globally among all surviving fibers (global load sharing or GLS scheme) or
among the surviving nearest neighbors only (local load sharing or LLS scheme).
For the GLS scheme Pierce ; Daniels  no stress concentration occurs anywhere
around the failed fibers as the stress of the failed fibers is shared among
all surviving fibers. On the other hand, in LLS scheme Phoenix ; Smith ;
Newman ; Harlow2 ; Harlow3 ; Smith2 , stress concentration is observed near a
broken patch (set of neighboring broken fibers) and increases with the size of
such patches. After such redistribution, the load per fiber increases
initiating failure of more fibers starting an avalanche. At the end of an
avalanche, either all fibers are broken (suggesting global failure) or the
bundle comes to a stable state with few broken fibers where an increment of
external stress is required to make the model evolve further. The last applied
stress just before global failure is considered to be the nominal stress or
strength $\sigma_{c}$ of the bundle. The fraction of fibers that survive at
$\sigma_{c}$ just before global failure is defined as the critical unbroken
fraction of fibers ($U_{c}$).

Figure 1: Upper Panel: The figure shows the spectrum of avalanche sizes ($s$)
and corresponding energy values ($E$) for GLS FBM with increasing number of
avalanches $k$. We can see that there is a direct correspondence between the
$s$ and $E$ values for a certain $k$. This means, higher $s$ gives higher $E$.
Since in case of GLS FBM, the fibers break in an increasing order of threshold
values, we get, $E(k+1)>E(k)$ if $s(k+1)=s(k)$. Lower Panel: The figure shows
the spectrum of avalanche sizes ($s$) and corresponding energy values ($E$)
for LLS FBM with an increasing number of avalanches $k$. We do not see a
direct correspondence between the $s$ and $E$ values here like GLS FBM. For
example, the red eclipses show the parts where only 1 fiber breaks at each $k$
value but the corresponding $E$ values show many different values without any
particular order as the fibers themselves do not follow any order while
breaking.

##  III Numerical Results

We have studied the fiber bundle model numerically in both mean-field limit
and with local load sharing scheme in one dimension, though the major part of
the paper will deal with the latter only. Numerical simulations are carried
out for system sizes ranging in between $10^{3}$ and $5\times 10^{5}$ and are
averaged over $10^{4}$ configurations. Our motive is to understand the
dynamics of avalanches and corresponding energy bursts emitted during these
avalanches as the model evolves with increasing externally applied stress.
Unless otherwise stated, we will use a uniform distribution ranging from 0 to
1 in order to assign threshold values to individual fibers beyond which it
breaks.

###  III.1 Relation between $s$ and $E$

Figure 1 shows a comparison between different avalanches and energy emitted
during those avalanches for a bundle of size $10^{5}$. The results are
produced for a single configuration. As usual, an avalanche is defined as the
number of fibers broken in-between two consecutive stress increments; $k$ is
the number of such stress increments in this case. While presenting the energy
spectrum and the avalanches we have excluded the final avalanche leading to
global failure.

Figure 2: The figure shows the variation of energy $E$ with avalanche size
(total number of broken fibers) $s$ for both GLS and LLS fiber bundle model.
We observe $E\sim s$, for GLS FBM. On the other hand, for LLS scheme, $E\sim
s^{\gamma}$, with $\gamma\approx 2.5$.

The upper panel of figure 1 shows the results for the GLS fiber bundle model
while on the lower panel, we have shown the results with the local load
sharing (LLS) scheme. Note that the range of $k$ for the LLS model is much
less than the range of $k$ with the GLS scheme. This is understandable since
with the LLS scheme, the model is more unstable due to stress concentration
and a large number of fibers are broken during the final avalanche. The model
evolves with a lesser number of stress increments in this case prior to a
global failure where the average size of the avalanches is smaller compared to
that in the GLS scheme. Now, for an avalanche of size $s$, if $n$ fibers with
threshold values $\tau_{1}$, $\tau_{2}$, $\tau_{3}$, $\cdots$, $\tau_{n}$
break, then the amount of energy emitted during this avalanche will be:

$\displaystyle
E(s)=\displaystyle\frac{1}{2}\displaystyle\sum_{i=1}^{n}\tau_{i}^{2}.$ (1)

This follows from the assumption of linear elastic (stress $\propto$ strain)
behavior of individual fibers up until their individual (brittle) failure
points. With above formalism, for each stress increment $k$, we will obtain an
avalanche $s(k)$ and a corresponding energy burst of magnitude $E(s(k))$.

The energy spectrum follows a particular trend in the case of the GLS scheme.
Since with the GLS scheme the fibers break in the increasing order of their
threshold values, the energy emitted at $k+1$-th load increment will be higher
than the energy emitted at $k$-th increment, even if the avalanche sizes
happen to be the same at $k$ and $k+1$. Due to this, the variations of $s$ and
$E$ with increasing $k$ looks exactly the same, only the values are scaled by
a constant when we transfer from $s$ to $E$. Such correlation between $s$ and
$E$ is not present in the case of the LLS fiber bundle model. In the case of
the LLS scheme, the fibers break due to the interplay between the local stress
profile and the threshold values of the fibers themselves. Due to such
dynamics, the fibers do not break in increasing order of their thresholds.
Then, there might be scenarios where $E(k+1)<E(k)$ when $s(k+1)=s(k)$ or even
$s(k+1)>s(k)$. The red ellipses in the lower panel of figure 1 shows this
absence of correlation between $s(k)$ and $E(s(k))$. For both ellipses, $s=1$
for that period. In spite of that, we see a fluctuation in energy values
without a particular trend. In the following, we will discuss this relation
between $s$ and $E$ in detail.

Figure 3: Distribution of energies for an uniform distribution (0,1) and
system sizes ranging in between $10^{3}$ and $10^{5}$. The results are shown
for GLS FBM. (a) We already know that in the mean-field limit $P(s)\sim
s^{-\beta}$, with $\beta\approx 2.5$. (b) We observe $Q(E)\sim E^{-\alpha}$
where $\alpha\approx 2.5$ as well independent of the system size.

Figure 2 highlights how energy $E$ emitted as a function of average avalanche
size $s$ for a bundle of size $10^{5}$ and configuration $10^{4}$. Results for
both GLS and LLS schemes are shown in the figure. We observe the following
behavior:

$E\sim\left\\{\begin{array}[]{ll}s&\mbox{, \text{for GLS},}\\\
s^{\gamma}&\mbox{, \text{with $\gamma=2.5$ for LLS}.}\\\ \end{array}\right.$
(2)

This behavior can be used to understand the relation between distributions
$P(s)$ of avalanche size $s$ and $Q(E)$ of emitted energies $E$. For this we
simply need to implement a change in variable theorem1  scheme as follows:

$\displaystyle Q(E)\sim
P[s(E)].|s^{\prime}(E)|=P[s(E)].|\displaystyle\frac{ds(E)}{dE}|$ (3)

Change in variable: GLS scheme

In case of GLS scheme, we observe

$\displaystyle E(s)\sim s$ $\displaystyle s(E)\sim E$ (4)

This makes

$\displaystyle s^{\prime}(E)=\displaystyle\frac{ds(E)}{dE}\sim 1$ (5)

We also know that the avalanche size distribution in case of GLS scheme is a
scale free distribution with an exponent 2.5 hh92 .

$\displaystyle P(s)\sim s^{-\beta},\ \ \ \text{with $\beta=2.5$}$ (6)

Then, combining Eq.3, III.1, 5 and 6, we get,

$\displaystyle Q(E)\sim P(E).1\sim E^{-\beta}$ (7)

Figure 4: Distribution of energies for an uniform distribution [0:1] and
system sizes ranging in between $10^{3}$ and $10^{5}$. The results are shown
for LLS FBM. (a) Avalanche size distribution for LLS FBM is an exponential
function: $P(s)\sim e^{-s/s_{0}}$, where $s_{0}$ depends weakly on the system
size. (b) Scale free distribution for energy emitted: $Q(E)\sim E^{-\alpha}$,
with $\alpha\approx 3.5$. Figure 5: (a) The finite size effect of $\alpha(L)$
is shown, which follows a scaling: $[\alpha(\infty)-\alpha(L)]\sim L^{-\eta}$
with $\eta=0.15$, where $\alpha(\infty)$ is the value of the exponent in the
thermodynamic limit. We get $\alpha(\infty)=3.47$. (b) & (c) The least square
fit error and corresponding exponent $\eta$ for the scaling
$[\alpha(\infty)-\alpha(L)]\sim L^{-\eta}$ is given for different values of
$\alpha(\infty)$. We consider the value of $\alpha(\infty)$ and $\eta$ which
produces the minimum error. This same procedure has been followed next while
exploring the same thing for different threshold distributions.

Change in variable: LLS scheme

In case of LLS scheme, we observe

$\displaystyle E(s)\sim s^{\gamma}$ $\displaystyle s(E)\sim E^{-\gamma}$ (8)

This makes

$\displaystyle
s^{\prime}(E)=\displaystyle\frac{ds(E)}{dE}\sim(-\gamma)E^{-(\gamma+1)}$ (9)

where $\gamma=2.5$. We also know that the avalanche size distribution in case
of LLS scheme is an exponential distribution khh97 .

$\displaystyle P(s)\sim e^{-s/s_{0}}$ (10)

Then, combining Eq.3, III.1, 9 and 10, we get,

$\displaystyle Q(E)\sim P(E).\gamma E^{-(\gamma+1)}\sim\gamma
e^{-\displaystyle\frac{E^{-\gamma}}{s_{0}}}E^{-(\gamma+1)}$ (11)

In the limit of high $E$ value, Eq.11 can be simplified as follows

$\displaystyle Q(E)\sim E^{-\alpha}\ \ \ \text{where $\alpha=\gamma+1=3.5$}$
(12)

Above treatment shows that, in case of LLS scheme, in spite of an exponential
distribution for avalanche sizes, the distribution of emitted energy is still
observed to be scale-free.

###  III.2 Distribution of $s$ and $E$: Uniform distribution

Figure 4(a) shows the avalanche size distribution $P(s)$ for a GLS fiber
bundle model with system size ranging from $10^{3}$ to $10^{5}$. This scale-
free decrease of $P(s)$ with $s$ is already known in the literature. We also
observe the same universal exponent 2.5 hh92 . Figure 4(b) shows the
corresponding distribution for the energy emitted. We observe the same scale-
free distribution for the energy as well, with the same exponent 2.5. This
behavior is consistent with Eq.6 and Eq.7 respectively.

Figure 4(a), on the other hand, shows the avalanche size distribution with the
LLS scheme. The distribution is exponential as derived analytically by Kloster
et al.khh97 . The inset of the same results in log scale in order to compare
them with the previous claim by Zhang and Ding zd94 , that $P(s)$ shows a
scale-free behavior with a very high exponent closer to $-4.8$. This claim of
scale-free nature is not substantiated and the exponential form for $P(s)$ is
accepted in the literature.

The distribution of energy in figure 4(b) shows a scale-free distribution, in
spite of the fact that the avalanche size distribution is an exponential
distribution. The exponent of the scale-free distribution is observed to an
increasing function of the size of the bundle

$\displaystyle Q(E)\sim E^{-\alpha(L)}$ (13)

The above behavior is similar to Eq.12, but with a $L$ dependent exponent
instead of a constant value. To compare this $L$ dependent exponent with the
value in Eq.12, we have to study the variation of $\alpha$ in Eq.13 in details
as the size of the bundle is increased. We have discussed this next.

Figure 5 shows the scaling of exponent $\alpha$ in figure 4(b) as the model
approaches the thermodynamic limit. We observe the following scaling,

$\displaystyle\alpha(\infty)-\alpha(L)\sim L^{-\eta}$ (14)

where $\eta=-0.15$ and $\alpha(\infty)$ (= 3.47) has a value close to
$\gamma+1$ (see Eq.12). The fitting and the exponent $\eta$ is calculated from
the minimization of least square fit error. This is shown in figure 5(b) and
(c). We choose a certain value of $\alpha(\infty)$ and fit our numerical
results. This in turn will produce a value of $\eta$ and corresponding least
square fit error. If we repeat this for a number of $\alpha(\infty)$ values,
then we can express the error (see figure 5b) and the exponent $\eta$ (see
figure 5c) as a function of $\alpha(\infty)$. The dotted line in figure 5(a)
corresponds to the value of $\alpha(\infty)$ (= 3.47) and $\eta$ (= 0.15) for
which the least square fit error is minimum.

###  III.3 Universality

So far, we have generated the numerical results where a uniform distribution
from 0 to 1 is used to assign random thresholds to individual fibers. In this
section, we will verify the universality of our results. For this purpose, we
will mainly explore 4 other distributions: (i) linearly increasing from 0 to
1, (ii) linearly decreasing from 0 to 1, (ii) a Weibull distribution with
scale parameter 1 and Weibull modulus 1, and (iii) A power law distribution
from 0 to 1 with exponent 2.0.

Figure 6: (a) Distribution of energies for 4 different threshold
distributions: (a) Linearly increasing [0:1], (b) Linearly decreasing [0:1]
(c) Scale free distribution of exponent 2 between 0 and 1, (d) Weibull
distribution with scale factor 1.0 and Weibull modulus 1.0. The system sizes
varies in between $10^{3}$ and $5\times 10^{5}$. The results are shown for LLS
FBM. We observe a scale free distribution for $E$: $Q(E)\sim E^{-\alpha}$ for
all thresholds. The exponent value $\alpha(L)$ shows finite size effect. The
dotted line shows the slope with highest system size in our simulation. The
exponents are close but less than 3.5. As before, this obeys the following
scaling: $[\alpha(\infty)-\alpha(L)]\sim L^{-\eta}$. The value of $\eta$ for
above mentioned distributions are 0.14, 0.12, 0.13 and 0.11 respectively. We
observe the exponent $\alpha(\infty)$ in the thermodynamic limit to be very
close to 3.5 irrespective of the choice of the threshold distribution. Figure
7: The figure shows the variation of energy $E$ with avalanche size (total
number of broken fibers) $s$ for LLS fiber bundle model. We have repeated the
study for 5 different threshold distributions. We observe for all
distributions, $E\sim s^{\gamma}$, where $\gamma$ has a value 2.5 independent
of the nature of the distribution.

In all these cases, the energy burst size distributions were found to be
scale-free with an exponent value close to $-3.5$ (see figure 6), as is
predicted from Eq. (12). The variation with system size also universal across
these different threshold distributions. These results suggest that the scale-
free nature of the energy burst size distribution in the local load sharing
fiber bundle model is a universal feature.

We have further checked that the relation between an avalanche size and an
energy burst size i..e., $E\sim s^{2.5}$ is valid for all these threshold
distributions, as can be seen from Fig. 7.

##  IV Discussions and conclusions

The local load sharing fiber bundle model is known to be lacking in
reproducing the scale-free avalanche statistics often seen in the experimental
setup of fracturing brittle solids. In all the interpolation schemes between
global (equal) and local load sharing versions of fiber bundles, the avalanche
size distribution $P(s)$ only show a cross-over between the mean-field
($P(s)\sim s^{-\beta}$) and local load sharing ($P(s)\sim e^{-s/s_{0}}$)
limits. The mean-field limit, however, is a rather idealized condition for
modeling real samples.

However, one important distinction between avalanche sizes ($s$) of the fiber
bundle model and what is usually measured in the experiments is that in the
latter case it is the energy burst ($E$) emitted in an avalanche that is
measured. However, that distinction is not at all significant in the mean-
field i.e., the global load-sharing limit of the model, because in that limit
$E\sim s$. However, in the local load sharing version, we numerically find
$E\sim s^{\gamma}$. Given an exponential distribution for the avalanche size
distribution in the local load sharing limit and this numerical observation,
it is possible to show that the size distribution of the energy bursts is
scale-free ($Q(E)\sim E^{-\alpha}$) with $\alpha=\gamma+1$ (see Eq. (12)). We
have then numerically checked that $\gamma\approx 2.5$ for various different
threshold distributions (see Fig. 6) and independently checked that the size
distribution exponent for the energy bursts are close to $-3.5$ (see Figs. 4,
7). Indeed, there are indications in experiments with sandstones that the
avalanche amplitude distribution was exponential while the energy burst
distribution was found to be a power law (see e.g., arma ; expt ). Our results
reproduce the same for the local load sharing fiber bundle model.

In conclusion, the local load sharing fiber bundle model is shown to have a
non-trivial relation between the avalanche size (number of fibers broken) and
the energy burst size (elastic energy released from the broken fibers).
Consequently, the energy burst size distribution is shown to have scale-free
nature, with an exponent value independent of the threshold distributions of
the fibers. Given that experimentally one measures the energy released, these
results indicate that local load sharing fiber bundles can have a significant
role in modeling fracture of brittle solids without having to resort to the
equal load sharing mean-field limit.

The authors thank Prof. Bikas K. Chakrabarti for his valuable comments. This
work was partly supported by the Research Council of Norway through its
Centres of Excellence funding scheme, project number 262644.

## References

  * (1) D. Bonamy, E. Bouchaud, Failure of heterogeneous materials: A dynamic phase transition?, Phys. Rep. 498, 1 (2011).
  * (2) H. Herrmann, S. Roux (eds.), Statistical models for the fracture of disordered media, Elsevier, Amsterdam (1990).
  * (3) B. K. Chakrabarti, L. Benguigui, Statistical physics of fracture and breakdown in disordered systems, Oxford University Press, Oxford (1997).
  * (4) H. Kawamura, T. Hatano, N. Kato, S. Biswas, B. K. Chakrabarti, Statistical physics of fracture, friction, and earthquakes, Rev. Mod. Phys. 84, 839 (2012).
  * (5) S. Biswas, P. Ray, B. K. Chakrabarti, Statistical physics of fracture, breakdown, and earthqauke: Effect of disorder and heterogeneity, Wiley-VCH, Singapore (2015).
  * (6) L. de Arcangelis, C. Godano, J. R. Grasso, E. Lippiello, Statistical physics approach to earthquake occurrence and forecasting, Phys. Rep. 628, 1 (2016).
  * (7) D. Bonamy, S. Santucci, L. Ponson, Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition, Phys. Rev. Lett. 101, 045501 (2008).
  * (8) J. Koivisto, M. Ovaska, A. Miksic, L. Laurson, M. J. Alava, Predicting sample lifetimes in creep fracture of heterogeneous materials, Phys. Rev. E 94, 023002 (2016).
  * (9) V. Kádár, G. Pál, F. Kun, Record statistics of bursts signals the onset of acceleration towards failure, Sci. Rep. 10, 2508 (2020).
  * (10) S. Biswas, D. F. Castellanos, M. Zaiser, Prediction of creep failure time using machine learning, Sci. Rep. 10, 16910 (2020).
  * (11) S. Pradhan, A. Hansen, B. K. Chakrabarti, Failure processes in elastic fiber bundles, Rev. Mod. Phys. 82, 499 (2010).
  * (12) A. Hansen, P. C. Hemmer, S. Pradhan, The fiber bundle model: Modeling failure in materials, Wiley-VCH, Singapore (2015).
  * (13) P. Bak, C. Tang, K. Wiesenfeld, Self-organized criticality: an explanation of $1/f$ noise, Phys. Rev. Lett. 59, 381 (1987).
  * (14) S. Pradhan, B. K. Chakrabarti, Precursors of catastrophe in the Bak-Tang-Wiesenfeld, Manna, and random-fiber-bundle models of failure, Phys. Rev. E 65, 016113 (2001).
  * (15) S. Biswas, S. Roy, P. Ray, Nucleation versus percolation: Scaling criterion for failure in disordered solids, Phys. Rev. E 91, 050105 (2015).
  * (16) S. Roy, S. Biswas, P. Ray, Modes of failure in disordered solids, Phys. Rev. E 96, 063003 (2017).
  * (17) F. Kun, G. Pál, I. Varga, I. G. Main, Effect of disorder on the spatial structure of damage in slowly compressed porous rocks, Phil. Trans. R. Soc. A 377, 20170393 (2019).
  * (18) A. Shekhawat, S. Zapperi, J. P. Sethna, From damage percolation to crack nucleation through finite size criticality, Phys. Rev. Lett. 110, 185505 (2013).
  * (19) F. T. Pierce, Tensile tests for cotton yarns, “The weakest link” theorems on the strength of long and composite specimens, J. Text. Inst. 17, 355 (1926).
  * (20) H. E. Daniels, The statistical theory of strength of bundles of threads, Proc. R. Soc. London, Ser. A 183, 405 (1945).
  * (21) S. L. Phoenix, The asymptotic distribution for the time to failure of a fiber bundle, Adv. Appl. Probab. 11, 153 (1979).
  * (22) R. L. Smith and S. L. Phoenix, Asymptotic distributions for the failure of fibrous materials under series-parallel structure and equal load-sharing, J. Appl. Mech. 48, 75 (1981).
  * (23) W. I. Newman and S. L. Phoenix, Time-dependent fiber bundles with local load sharing, Phys. Rev. E 63, 021507 (2001).
  * (24) D. G. Harlow and S. L. Phoenix, The chain-of-bundles probability model for the strength of fibrous materials I: analysis and conjectures, J. Compos. Mater. 12, 314 (1978).
  * (25) D. G. Harlow and S. L. Phoenix, Probability distributions for the strength of fibrous materials under local load sharing I: two-level failure and edge effects, Adv. Appl. probab. 14, 68 (1982).
  * (26) R. L. Smith, Proc. R. Soc. London, Ser. A 382, 179 (1982).
  * (27) Let’s assume $x$ and $y$ are a continuous functions with p.d.f $f(x)$ and $g(y)$. Also, $y=u(x)$ is a function with inverse, which means it is possible to find $x=v(y)$. In this circumstances, $g(y)$ and $f(x)$ can be expressed as: $g(y)=f[v(y)].|v^{\prime}(y)|$, where $|v^{\prime}(y)|=|dv(y)/dy|$.
  * (28) M. Kloster, A. Hansen, P. C. Hemmer, Burst avalanches in solvable models of fibrous materials, Phys. Rev. E 56, 2615 (1997).
  * (29) P. C. Hammer, A. Hansen, The distribution of simultaneous fiber failures in fiber bundles, ASME J. Appl. Mech. 59, 909 (1992).
  * (30) S. D. Zhang and E. J. Ding, Burst-size distribution in fiber-bundles with local load-sharing, Phys. Lett. A. 193, 425 (1994).
  * (31) S. Pradhan, A. M. Stroisz, E. Fjær, H. K. Lund, E. F. Sønstebø, S. Roy, Fracturing tests on reservoir rocks: Analysis of AE events and radial strain evolution, ARMA 14, 7442 (2014).
  * (32) S. Pradhan, A. M. Stroisz, E. Fjær, J. F. Stenebråten, H. K. Lund, E. V. Sønstebø, Stress-induced fracturing of reservoir rocks: Acoustic monitoring and $\mu$CT image analysis, Rock Mech. Rock Eng. 48, 2529 (2015).