Efficient hyperbolic-parabolic models on multi-dimensional unbounded domains using an extended DG approach

We introduce an extended discontinuous Galerkin discretization of hyperbolic-parabolic problems on multidimensional semi-infinite domains. Building on previous work on the one-dimensional case, we split the strip-shaped computational domain into a bounded region, discretized by means of discontinuous finite elements using Legendre basis functions, and an unbounded subdomain, where scaled Laguerre functions are used as a basis. Numerical fluxes at the interface allow for a seamless coupling of the two regions. The resulting coupling strategy is shown to produce accurate numerical solutions in tests on both linear and non-linear scalar and vectorial model problems. In addition, an efficient absorbing layer can be simulated in the semi-infinite part of the domain in order to damp outgoing signals with negligible spurious reflections at the interface. By tuning the scaling parameter of the Laguerre basis functions, the extended DG scheme simulates transient dynamics over large spatial scales with a substantial reduction in computational cost at a given accuracy level compared to standard single-domain discontinuous finite element techniques.


INTRODUCTION
Accurate computational representation of physical phenomena on unbounded domains remains a challenge for numerical methods for transient differential problems, with fields of interest ranging from cell growth modelling 1 to upper atmosphere dynamics 2 and space weather 3,4 .In the quest to make the inherently infinite problem tractable, numerical treatments based on absorbing (also called sponge) layers are usually preferred to analytical approaches relying on non-reflecting boundary conditions 5,6,7,8 (though see also 9 for a recent study challenging that convention).The former class of methods typically relies on dividing the computational domain into a bounded region of interest to the physical phenomena relevant to the problem at hand (e.g., the troposphere and lower stratosphere in currently operational atmospheric models 10,11 ) and an unbounded buffer region for damping the outgoing signals.As the interface between the two subdomains is essentially artificial, the numerical treatment of the interface and the buffer region need to cause reflections in the bounded region that are as low-amplitude as possible.An optimal setup with absorbing boundary conditions maximizes the extent of the region of interest while keeping enough computational points in the buffer zone so as to minimize reflections due to the artificial boundary and the buffer zone itself.To that end, fine-tuning of the sponge layer parameters is usually required, making it challenging to obtain efficient numerical tools agnostic to the properties of the outgoing waves.We refer to comprehensive reviews, e.g., 12,13,14,15 and references therein for a complete overview of the open boundary conditions issue.
In order to discretize differential problems on unbounded domains, spectral collocation methods 16 based on scaled basis functions 17,18,19,20,21,22 were used in unbounded regions coupled with finite volume 23 and discontinuous Galerkin (DG) finite element 24 methods on finite regions for one-dimensional linear and nonlinear hyperbolic problems.The coupled setup enabled accurate simulation of wavelike phenomena and economical damping of outgoing signals with minimal reflections in the region of interest using a small number of scaled Laguerre basis functions and tuning their scaling parameter in the semi-infinite part.The methodology was then broadened to include the diffusive case 25 within a seamless, extended DG (XDG) modal approach using Legendre basis functions in the bounded region, scaled spectral Laguerre basis functions in the unbounded region, and flux-based coupling of the two basis function sets at the finite/semi-infinite interface.A formal analysis confirmed the stability and the potential of the XDG framework in providing reliable and efficient numerical solutions in very large regions.
The present study extends the previous one-dimensional proof-of-concept models 23,24,25 in two unexplored directions: • First, the XDG framework is implemented in two spatial dimensions and used to simulate linear and nonlinear problems on semi-infinite strip-shaped domains.The seamlessness of the approach is evident in the expressions for many of the integrals appearing in the discrete formulation in the finite and semi-infinite regions.Those expressions are formally identical up to replacing the finite-domain grid size with the inverse of the Laguerre scaling parameter; • Second, time-to-solution values with the two-dimensional XDG model are compared one-to-one with a single-domain DG framework that uses Legendre basis functions and an inhomogeneous grid with nodes coinciding with those of the extended scheme in the unbounded part.For several choices of model parameters, the new model provides numerical solutions with several times lower wallclock times at a given accuracy level, with the efficiency advantage growing with the number of Laguerre modes and thus larger simulated portions of the semi-infinite region.
We remark that alternative methods exist for high-order simulations on unbounded domains, but available methodologies are limited to the linear case without like-for-like efficiency comparisons 26 , concern elliptic problems 27 , or lack formal stability analyses 28 .Recent studies 29,20 introduced adaptive methods based on time-evolving distributions of Laguerre and Hermite collocation points for simulations on unbounded domains and applied them within spectrally adaptive neural network models 30 .The rest of the paper is structured as follows.Section 2 describes the problem definition and geometric setup, while section3 provide details on the 2D XDG model and its discrete and algebraic formulation.Section 4 reports the results obtained in (a) convergence tests with a standalone DG-Laguerre 2D formulation; (b) validation tests for the coupling strategy of the extended scheme for the linear advection-diffusion, nonlinear Burgers', and nonlinear shallow water equations; (c) performance tests in the semi-infinite regions; and (d) absorbing layer tests with an added reaction term modelling the damping of single-wave and wave train-like signals.Section 5 concludes the paper.

PROBLEM DEFINITION AND GEOMETRY DISCRETIZATION
In a two-dimensional semi-infinite strip Ω = [0,   ]×ℝ + ≡ [0, +∞), we consider the following conservation law with damping: where Ω  denotes the part of the boundary where a Dirichlet condition is imposed, while a Neumann condition is assigned to Ω  .For the sake of simplicity, in the rest of this section we will focus on the transport-diffusion case, i.e.: (, , ) = −∇ +  (2)   where  is a 2 × 2 symmetric and positive definite matrix uniformly bounded above and below, that is, there exist two positive constants  0 and  1 such that and we will assume that  is constant and diagonal: If   ,   > 0, then  is symmetric and positive definite, and the boundedness condition is automatically satisfied with  0 = min(  ,   ) and  1 = max(  ,   ).The function  denotes a generic transport term which may also depend on the spatial variables.Thus  = (, , ) = [ 1 (, , ),  2 (, , )] ⊤ but we will omit the dependence on  and  in the following expressions.We allow  to be non-linear in  and discontinuous in  and .
where  is a sufficiently regular trial/test space.We split the domain as and we discretize the problem using an XDG scheme 25 .
In two dimensions, XDG is a coupled scheme consisting of: • a standard DG scheme in both directions in the rectangle • a "DG in "-"Laguerre in " scheme with scaled Laguerre basis functions in the semi-infinite region More specifically, we introduce in Ω − the computational grid Then, we consider the discontinuous finite element space: and choose as basis of ℙ   ,  (   ,  ) the normalized Legendre basis where   is the -th Legendre polynomial.Hence, the solution of (1) will be represented in each interval    ,  as Next, the semi-infinite region is discretized using the computational grid +∞).Scaled Laguerre functions are chosen as basis functions in the -direction 23,24,25,19 .Defining we can represent the solution as Note that the same notation  with different superindices is used to denote both types of basis functions of the XDG approach in order to remark the unitary character of the discretization using Legendre-based and Laguerre-based elements.
We then sum over all intervals    ,  to obtain where Γ ℎ is the set of edges, being Γ  the set of internal edges, and Γ  and Γ  the edges on Ω where a Dirichlet or a Neumann boundary condition is imposed, respectively.We also allow periodic boundary conditions to be assigned on a subset of the boundary -in this case we will simply treat the corresponding edges as internal edges.The set Γ  is further subdivided as to distinguish edges internal to Ω − , edges internal to Ω − and edges on the interface between the two regions (Figure 1).We analyze the second and fourth integral in (13), i.e., the integrals over the edges.Starting from the diffusive term, we split the sum as We impose the Neumann boundary condition on Γ  , so that the second sum on the right hand side becomes Regarding the first sum on the right hand side, we first notice that where the jump and average operators across an edge  are defined as where the elements separated by  are denoted as Ω  and Ω  and  is oriented from Ω  to Ω  .If  ∈ Γ  is a boundary edge, we extend the definitions of jump and average as If  is regular, we have that [[∇ ⋅ ]] = 0 and [[]] = 0 on every  ∈ Γ  ∪ Γ  .Using this property, the Dirichlet boundary condition on Γ  on equation (18), and the definition of jumps and averages on a boundary face, we obtain for any choice of the penalization parameters  and .|| denotes the length of the interface .
We can now insert (21) into (13) to obtain the following expanded weak formulation: for all test functions  ∈  , find  ∈  such that

ALGEBRAIC FORMULATION
We can now move from the continuous to the discrete setting.We first introduce the spaces are not uniquely defined at the internal edges.We then replace it by q( ℎ , ) ℎ|  , where q is the numerical flux (see Appendix A for details).
The Galerkin formulation associated with (1)- (22) In order to obtain a discrete formulation, the representations ( 9) or (11) for  ℎ are replaced in the above expression, choosing The resulting integrals are then simplified using the orthogonality of the Legendre and Laguerre basis functions, see Appendix A for details on the computation of the integrals.
In order to cast the spatially discrete problem in algebraic form, we first define the unknown vector as  ∈ ℝ   , with ) .We also define   =     (  +1)(  +1) and   =   (  +1)( +1) as the number of degrees of freedom in the bounded and unbounded region, respectively.We order the unknown degrees of freedom as follows.We first consider the degrees of freedom in Ω − and define the local vector ,  ,  (0,1) which contains the unknowns corresponding to    ,  .The global vector for the DG-DG discretization in the finite region is then Concerning the DG-Laguerre discretization in the semi-infinite region, we order the unknowns as ]  ∈ ℝ   . ( The global unknown vector for the XDG formulation will then be  = [ − ,  − ]  ∈ ℝ   .The algebraic formulation consists of the following components: • The equation for the time derivative of the two unknown vectors is (see also equation (A2)): • The terms referring to the degrees of freedom in the bounded region result in the algebraic representation  −  − , where  − ∈ ℝ   ×  represents the DG-DG discretization in Ω − (see expressions (A12),(A13),(A17),(A18),(A19) and the first terms of equation (A16)).
• The coupling terms are expressed algebraically as  ,12  − and  ,21  − (see the second term and the third term of equation ( A16)).The matrices  𝐶,12 ∈ ℝ   ×  and  ,21 ∈ ℝ   ×  (30)   are the coupling matrices, which are non-zero only in the last   (  + 1)(  + 1) columns and the last   (  + 1)(  + 1) rows, respectively, corresponding to the degrees of freedom in the   DG-DG elements adjacent to the interface.
• The transport term can be assembled in the vector () ∈ ℝ   (see expressions (A23) and (A24) and Appendix C).
The algebraic formulation of the problem is then where In the linear case, the transport term reduces to  , where the matrix  has the same block structure as .
The spatially discretized problem can be solved in time by means of a suitable time integration scheme.In the linear case, equation (31) reads where  =  +  + .The -method can then be employed: given  0 and  ∈ [0, 1], the solution at time step  is given by In the non-linear case we rewrite (31) as then separate the linear and non-linear terms by defining   (, ) = − −1 ( + ) and   (, ) = − −1 (() − ).Problem (35) is discretized in time using a 3-stage IMEX-ARK method 31 .
The algebraic formulation readily extends to the vectorial case.The theoretical details are omitted here for conciseness, section 4.2.3 contains numerical results with the shallow water equations.

NUMERICAL EXPERIMENTS
In this section, several numerical tests of the XDG scheme are carried out.First, we perform a convergence analysis of the standalone unbounded two-dimensional section of the XDG scheme, where the  direction is discretized with Legendre basis functions and the  direction is discretized with scaled Laguerre basis functions.Therefore, we set   = 0 and compute the convergence rates of both the spatial discretization (-convergence and -convergence) and the temporal discretization.Next, the coupling strategy developed in Section 3 between the finite Ω − and the semi-infinite Ω − portions of the 2D strip in the XDG model is validated both on the linear 2D advection-diffusion equation and on non-linear problems -the 2D Burgers' equation and the 2D shallow water equations.As in previous work 23,24,25 , the quality of the coupling is evaluated by computing the errors of the XDG approach with respect to a reference DG discretization on a larger domain in benchmarks with Gaussian signals.Next, we compare the performance of both the 2D XDG scheme and a DG discretization that uses a non-uniform grid in the semi-infinite direction.For both schemes, we compute wallclock times and errors with respect to an exact solution in tests with the advection-diffusion equation and Gaussian initial data using a large number of Laguerre modes and small values of the scaling parameter so as to simulate a larger portion of the unbounded region.Finally, an absorbing layer is implemented in the semi-infinite region by means of a sigmoidal reaction term, and tests are run with wave trains to check whether the obtained model is able to efficiently damp outgoing signals without reflections in the finite part of the two-dimensional strip.
Errors are computed using suitable Gaussian quadrature rules.Specifically, we define the discrete norms where =1 are Gaussian nodes and weights in the reference interval [−1, 1].Then, we define the absolute and relative errors between the numerical solution  ℎ and a reference solution   as In all the following tests we employ a uniform grid in time with   time steps, so that Δ =  ∕  .We compute absolute errors at the final time and we estimate the rates of convergence as follows.Assume that   (ℎ) ≈ ℎ   , where  is a constant and so that The Courant number values in the  and  direction are computed as   =   Δ  ∕Δ and   =   Δ  ∕Δ; in the semiinfinite region, Δ is taken as the distance between the first two Laguerre nodes and   = 1.In all experiments on the linear advection-diffusion equation we employ the Crank-Nicolson method in time, that is, the -method with  = 1∕2, while the Burgers' equation and the non-linear shallow water equations are solved in time by means of a 3-stage IMEX-ARK method 31 .
Concerning the DG penalization method, we focus on NIPG ( = 1), which is stable for any choice of  ≥ 0; we therefore set  = 0 in the following.For the sake of conciseness of the main text, detailed results are placed in the Appendix, see the tables in Appendix D. The model was implemented in MATLAB ® Version 9.9.0.1524771 (R2020b Update 2) on a workstation with processor CPU Intel ® Core TM i7-9750H CPU @ 2.60GHz with 6 cores and 32GB RAM.Timings reported in the performance tests in section 4.3 are therefore to be taken as preliminary.An implementation with a lower-level programming language will allow more accurate performance measurements via batch submissions to compute-only nodes on HPC systems -this is left for future work.

Standalone DG-Laguerre discretization -Convergence tests
The first set of numerical experiments of the two-dimensional model focuses on a single-domain discretization of the semiinfinite strip, with a Laguerre discretization in the  direction and a DG scheme in the  direction (  = 0).A convergence analysis is carried out on the linear undamped non-homogeneous 2D advection-diffusion equation, i.e., equation (1) with We impose the exact solution: with  0 =   ∕2 and  0 =   ∕10, and compute the right-hand side accordingly.We impose homogeneous Dirichlet boundary conditions on the lower boundary  = 0 and periodic conditions on the right and left boundaries  = 0 and  =   .Absolute errors are computed at the final time  with respect to the exact solution.
First, we evaluate spatial convergence in the  direction.To this end, a large number of intervals are used for the DG discretization in the  direction,   = 1000 with   = 3.In order to highlight the convergence in , we also set a small time step, Δ = 5 × 10 −4 , and we run the simulation for   = 100 time steps.We choose  ∈ {5, 10, 20, 30, 35}; the corresponding values of the Courant number in the -direction are   ∈ {8.10 × 10 −3 , 1.40 × 10 −2 , 2.86 × 10 −2 , 4.22 × 10 −2 , 4.90 × 10 −2 }, while   = 1.5.The small values of the Courant number in the  direction are due to the time step Δ being much smaller than the grid spacing Δ.This is not dictated by stability reasons, as the implicit time integration scheme is stable for values of   = (1) (this is confirmed by numerical tests not reported here).Our intent is rather to highlight the convergence in  by making the discretization error in time negligible.Experimental convergence is exponential in the number of Laguerre modes  (straight lines in the semi-logarithmic plot, Figure 2 left, see also table D1 in Appendix D), in line with expectations with a spectral discretization.In Figure 2 center, the coefficient of the exponential (solid black line) was chosen as 0.4 for benchmarking purposes.Behaviour for larger  values (not shown) confirms the super-polynomial convergence.
Next, we study the convergence of the DG discretization in the  direction.The NIPG scheme is stable for any choice of  ≥ 0, but it is suboptimal when the polynomial degree is even, in agreement with the theory 32 .Therefore, we get quadratic convergence for   = 2 as well as   = 1, but convergence of order 4 is recovered when   = 3 (Figure 2 center for   = 3, see also Table D2 in appendix D for the results for all   = 1, 2, 3).The time step and the number of time steps are the same as the previous test, (Δ = 5 × 10 −4 ,   = 100).

Validation of the coupling strategy in the 2D XDG scheme
In the second set of tests, we validate the accuracy of the coupling between the finite part and the semi-infinite part of the 2D strip, comparing the results obtained with the 2D XDG scheme and those obtained with a standalone DG discretization on a larger domain 23,24,25 , still considering the undamped case  = 0. We remark that the current implementation of the XDG model enables the use of polynomial degree up to 4 for the finite part, as shown in the convergence tests and later in the efficiency tests of section 4.3.However, since a high polynomial degree is not needed to validate the coupling strategy, we restrict ourselves to   =   = 1 in the coupling validation tests.

Linear advection-diffusion equation
We start with the advection-diffusion equation with Gaussian initial data placed inside the finite part.More specifically, we set   =   = 10 ,   =   = 1,   = 50,   = 500,  = 4 ,   = 200, corresponding to Δ = 0.2 , Δ = 0.02 , Δ = 0.02 .The physical parameters are   =   = 0.1  2 ∕,   = 0.5 ∕,   = 1 ∕, so that the Courant number in the -direction is   = 1.A homogeneous Dirichlet condition is imposed on the lower boundary, with periodic conditions on the vertical boundaries.The initial data is the Gaussian profile with  = 1 ,  0 = 5  and  0 = 8 .As the simulation progresses, the initial profile is transported and diffused and crosses the finite/semi-infinite interface.We test the coupling approach for different choices of   =   and .For a fixed value of , the scaling parameter  should be chosen carefully in order to represent the semi-infinite part of the domain taking care of two aspects.On the one hand, the spacing between the first two Laguerre nodes should be comparable to the size of the last DG interval for the sake of information exchange at the interface.On the other hand, no perturbation leaving the finite domain should reach the last Laguerre node, so  should be tuned in such a way that the semi-infinite part of the domain is large enough for the problem at hand.The solution computed by the 2D XDG model (Figure 3) with the combinations (, ) ∈ {(10, 4), (40, 6)} is compared to a reference 2D single-domain DG discretization on [0,   ] × [0, 2  ] with the same spacing in the -direction as the finite region of the extended scheme, that is,  ′  = 1000.Relative errors of the XDG scheme with respect to the reference are at most around a few percent (Table D4 in Appendix D), in line with previous published work 25 .It is also worth noting that, with the current set of parameters, the number of total entries of the matrix of the XDG scheme is almost 4 times smaller than the corresponding matrix for the standalone DG scheme.For  = 10 and  = 40 Laguerre basis functions, the matrix of the XDG scheme has about half as many non-zero entries as the matrix for the standalone DG scheme.
Relative and absolute errors in the finite domain due to the coupling are small and limited to a narrow region near the interface (Figure 4 right, see also Table D4 in Appendix D).

Non-linear shallow water equations
The two-dimensional non-linear shallow water equations are 33 where ℎ = ℎ(, , ) is the water depth,  = (, , ) and  = (, , ) the velocities in the  and  direction, respectively, and  is the acceleration of gravity.Introducing the vectors

ℎ(𝑥, 𝑧, 𝑡) ℎ(𝑥, 𝑧, 𝑡)𝑢(𝑥, 𝑧, 𝑡) ℎ(𝑥, 𝑧, 𝑡)𝑣(𝑥, 𝑧, 𝑡)
and we obtain the system We solve the problem in the semi-infinite strip Ω = [0, 10 ] × [0, +∞).To that end, we set   = 50 and   = 1; we place the interface at   = 8  with   = 40 and   = 1.We then test the accuracy of the XDG scheme for different choices of  and : the value of  is selected so that the extension of the semi-infinite part is approximately the same in all tests.The final time is  = 0.75  with   = 750 time steps.We write depth and velocities as ℎ(, , ) =  + h(, , ), (, , ) =  + ũ(, , ) and (, , ) =  + ṽ(, , ), where h, ũ and ṽ are perturbations of the reference configuration, taken to be  = 10 ,  = 1 ∕,  = 2 ∕.The initial condition is homogeneous for both ũ and ṽ.We consider for h(, , 0) a Gaussian profile as in (42) (see Figure 3 left) with  = 1 ,  0 =   ∕2 = 5 ,   =   = 1 ; in the outgoing case, the initial profile is placed inside the finite region, with  0 = 5 , while in the incoming case it is centered beyond the interface, at  0 = 10 .The solution of the XDG scheme is compared to a reference single-domain DG solution on  ′  = 3  = 24  with the same grid spacing,  ′  = 3  = 120 23 .The boundary conditions are periodic in the  direction and an outflow condition is imposed on the lower boundary  = 0.With these setup choices, in both the incoming case and the outgoing case, the initial data is horizontally propagated across the periodic boundaries.
The XDG scheme simulates wave motion at the correct speed both in the outgoing case and in the incoming case (Figure 5 for  = 50).Relative errors compared to the reference single-domain DG solution are mostly far below one percent, and at most a few percent (Tables 1 and 2), thus further validating the coupling strategy in the case of hyperbolic systems.

Extended DG vs. single domain DG performance comparison in large semi-infinite region
Next, we run a performance comparison between the XDG scheme and a single-domain DG discretization in wave simulation in the semi-infinite region.In these efficiency comparison tests, for simplicity we consider the undamped linear advection-diffusion equation with constant coefficients,  = 0 and Gaussian initial profile (42).The exact solution in ℝ 2 at a generic time  is   ∞ () 10 1 5.24e-05 3.90e-04 1.01e-03 1.33e-02 2.36e-04 1.73e-03 30 3 7.02e-07 8.62e-06 4.58e-05 5.92e-04 3.82e-06 3.93e-05 50 5 2.08e-07 1.44e-06 1.11e-05 1.04e-04 9.63e-07 7.22e-06 Table 1 Shallow water equations with Gaussian initial data, outgoing case.  = 10 ,   = 8 ,   = 50,   = 40,   =   = 1,  = 0.75 ,   = 750,  = 10 ,  = 1 ∕,  = 2 ∕.Relative errors in the finite region with respect to a single-domain DG discretization on  ′  = 3  with  ′  = 3  .See also Figure 5, top panel.In order to choose the value of  for the simulations, we carry out an error analysis in the XDG scheme as  varies, for a fixed value of .Relative errors in the semi-infinite domain are computed at the final time with respect to the exact solution (Figure 6 left).We also report in Figure 6 (right) the extension of the semi-infinite domain and the distance between the first two Laguerre nodes for each choice of .The analysis suggests that  should be neither too large nor too small.In the former case, the spatial extension of the semi-infinite region is not large enough to accurately represent the solution beyond the interface, while in the latter case the first two nodes are too far apart, causing spurious signals due to the coupling.If  is in those critical regions, the errors exhibit a noticeable increase.In addition, by fine-tuning the choice of  the accuracy of the numerical solution can be improved by over one order of magnitude, with no effect on the total computational cost (Figure 6 left).On the basis of this analysis, for this set of tests we select the value  = 5, which allows the representation of a large portion of the semi-infinite domain while keeping the errors as small as possible.

𝑞(𝑥, 𝑧, 𝑡)
We remark that techniques exist 20 that further enhance the accuracy of Laguerre spectral discretizations by adjusting the extension of the semi-infinite domain through an adaptive choice of the  and  parameters.While such a feature would be of relatively little interest in the context of damped runs with a small number of Laguerre basis functions such as those considered in the previous sections, adaptive algorithms 20 could be of great use where higher accuracy throughout the simulation might be required.The exploration of such techniques in the framework of the XDG model is left for future work.With the choice  = 5, errors obtained with the 2D XDG scheme are low even using a polynomial degree   = 1 in the finite region (black lines in Figure 7).By contrast, the 2D single-domain DG discretization on a non-uniform grid in  achieves errors lower than one percent only using   = 3 (red and blue lines in Figure 7), at more than 4 times the computational cost of the XDG model.  = 1 Finally, we choose a smaller value of  to simulate large spatial scales in the  direction; more precisely, we set  = 0.05 and we move the interface at   = 100 .In this way, the last Laguerre node (i.e., the extension of the unbounded part of the XDG scheme) is placed at   = 633.7 ,   = 2120.9 and   = 3652.4 for  = 30, 50, 70, respectively.With these parameter choices, the numerical solution obtained with the XDG scheme is at least as accurate as, and much less computationally costly than, the single-domain DG solution, with a significant efficiency gain (Figure 8 and Table 3).For  = 30 points in the semiinfinite region, the single-domain DG scheme runs four times slower than the XDG scheme at given accuracy level.For  = 50, the XDG scheme is five times more efficient.  = 1   3 Comparison between the best relative errors obtained with the XDG scheme (XDG) and the single-domain DG scheme (DG) in the case  = 0.05 and corresponding computational times  in seconds.Polynomial degree   = 1 is used for the XDG and   = 3 for the single-domain DG discretization.The efficiency gain  is defined as the ratio between the computational time of the single-domain DG run and that of the XDG run.

Efficiency of the XDG scheme in absorbing layer tests
We now implement an absorbing layer in the semi-infinite region by setting  in (1) as the sigmoid defined by Here Δ is the sigmoid amplitude,  ∈ [0, 1] the position of the sigmoid inside the absorbing layer,  0 the spatial extension of the semi-infinite region and   the sigmoid steepness.We consider the 2D nonlinear shallow water equations ( 43) and impose a Dirichlet boundary condition sinusoidal in time to simulate a train of waves originating at the lower boundary   = 0 23,24 : The   , the value of the scaling parameter  is chosen so that the distance between the first two Laguerre modes is approximately constant as  varies.First, to assess the efficiency of the absorbing layer, we compare the numerical solution of the 2D XDG scheme with a reference single-domain DG solution on [0,   ] × [0, 5  ] using a uniform grid with the same grid spacing Δ as in the finite portion of the XDG scheme, and we compute errors at the final time.The absorbing layer, active in the semi-infinite part of the XDG setup, is able to damp outgoing perturbations (central panels of Figures 9 and 10).Spurious reflections in the finite region due to the use of two different sets of basis functions in the 2D XDG scheme result in relative errors in the finite region of a few percent at most (Table 4 and right panels of Figures 9 and 10).Moreover, by suitably tuning the scaling parameter , errors can be kept under control even with a small number of modes in the unbounded region.Table 4 Non-linear shallow water equations with an absorbing layer in the semi-infinite region, wave train test, XDG vs. uniform-grid DG.Relative errors of the XDG model at the final time with respect to a single-domain DG discretization in [0, 12 ] × [0, 10 ] with  ′  = 5  .(The grid spacing Δ in the single-domain DG discretization is the same as that used in the finite region of the XDG scheme.)Second, to gauge the damping performance of the XDG scheme compared to a standard approach, we analyze the results of the single-domain DG discretization on a uniform grid with respect to a single-domain DG discretization on a non-uniform grid in the -direction, where the endpoints of the intervals coincide with the Laguerre nodes employed in the semi-infinite region for the XDG discretization.The relative errors in the finite region [0, 12 ] × [0, 2 ] (Table 5) are (10 −2 ) or less, and comparable with those in Table 4, showing that no accuracy penalty is incurred by using the XDG scheme as an absorbing layer.
Although the final time for the tests in this section was kept short because of limited computational resources, additional numerical tests not shown here confirmed that the errors actually level off at the values shown in Tables 4 and 5 when running until later times.In addition, the values did not vary significantly when refining the spatial or temporal spacing, thereby confirming that the errors are chiefly due to the presence of the absorbing layer and the choice of related parameters.While the two-dimensional results confirm the viability of the XDG scheme 23,24,25 , we expect the advantages of the XDG scheme to be more evident in problems with more complex dynamics.An in-depth exploration of the properties of the XDG scheme in damping internal waves arising as solutions of fluid flow equations in stratified media and non-flat bottom boundaries is left for future work.

CONCLUSION
We proposed an extended discontinuous Galerkin numerical scheme for the discretization of hyperbolic and parabolic problems on two-dimensional semi-infinite domains.The method deals with the unbounded computational domain by splitting it into a finite region, where Legendre polynomials are used as basis and test functions, and an unbounded region, discretized using scaled Laguerre functions.The two subdomains are seamlessly coupled by means of numerical fluxes at the interface.For several integral terms within the discrete formulation, the expressions for the integrals was found to be formally identical except for replacing the finite subdomain grid spacing with the inverse Laguerre scaling parameter.Formal accuracy of the algorithm was retrieved in convergence experiments, which confirmed the expected theoretical behaviour.No accuracy degradation issues were encountered linked to the tensor-product structure of the multi-dimensional discretization, obviating the need for approaches such as the hyperbolic cross approximation developed in 34 .
Coupling validation tests were carried out on the 2D linear advection-diffusion equation, the 2D Burgers' equation, and the 2D non-linear shallow water equations using Gaussian initial data.Comparing the XDG solution with a standard single-domain DG discretization, we showed that spurious oscillations due to the coupling are negligible.Numerical experiments showed that the interface between the two regions is transparent to both outgoing and incoming signals, with minimal spurious reflections into the bounded domain.
Numerical evidence shows that the choice of the scaling parameter  affects the accuracy of the solution in the semi-infinite part of the domain.This is an important feature of the method because the optimal selection of  would allow one to achieve smaller errors with the same number of spectral modes .While this work only reports the best results obtained by tuning , a rigorous criterion for its selection remains an important open problem and shall be the subject of future work.
In addition, by suitable tuning of the scaling parameter of Laguerre basis functions, the XDG scheme proved to be a highly efficient tool to accurately represent dynamics over arbitrarily large spatial scales.Achieving comparable accuracy in the semiinfinite region with a single-domain DG discretization required higher polynomial degrees and more than four times as much computational time in performance comparison tests.
By means of a sigmoidal reaction term, an absorbing layer was then implemented in the semi-infinite region.Exponential convergence provided by the Laguerre basis functions enabled the scheme to efficiently damp outgoing wave trains with a small number of modes.Reflections into the finite domain measured as errors with respect to a single-domain DG discretization on a uniform grid were of comparable amplitude to those obtained by using a standard DG approach on a non-uniform grid with nodes coinciding with the Laguerre nodes in the semi-infinite region.
Numerical evidence suggests that the proposed XDG approach can be a competitive tool to discretize transient dynamics efficiently over arbitrarily large computational domains.Inspired by these results, several research directions can be explored in future work.First, a natural extension would be the use of polar coordinates to simulate semi-infinite circular domains.In that setup, more complex nonlinear magnetohydrodynamics models could be considered in order to gauge the potential of the XDG scheme to model phenomena spanning very large length scales such as coronal mass ejections.Second, with the addition of bottom topography, the scheme could be used to discretize compressible Euler equations first in vertical slice two-dimensional domains and then in full-fledged three-dimensional implementations.Established idealized orography benchmarks 35,10 could test the method's effectiveness as an efficient damping layer.If the techniques proposed in our work achieve the same accuracy at a lower computational cost than currently employed approaches, related savings can be reinvested in the region of interest to achieve e.g. higher accuracy or longer lead times in forecasts.
Finally, increasing resolution in the unbounded region would enable comparison with extensions of weather forecast models to the middle and upper atmosphere and, in perspective, use of the method in space weather models 36,4,37 .Future applications of the method could also benefit from a recently proposed adaptive capability 29,20 whereby the scaling parameter can vary within simulations, thus enabling the scheme to simulate multiple phenomena developing over very different spatial and temporal scales.
the first integral can easily be written as The fifth integral is ( The sum over  runs from 0 to   if   = 1, … ,   , from 0 to  if   = ∞.Next, we consider the edge integrals, i.e. the second, third, fourth and sixth integral in (25); to this end we introduce the quantities for   = 1, … ,   , which represent the value of basis functions and their derivatives at the left and right endpoint of each subintervals.Notice that, on a uniform grid, these quantities do not depend on   and   .We start from the second, third and fourth term in (25).For the sake of exposition, in the following we specialize the formulation by considering horizontal edges aligned with the  direction and vertical edges aligned with the  direction.Let  ∈ Γ and the coupling between the two discretizations is described by the matrices   (see Appendix B).
To conclude the analysis of the second, third and fourth terms of (25) we consider Dirichlet edges  ∈ Γ  , where we need to apply the definition of jumps and averages on boundary edges.If  belongs to the left boundary  = 0, then  = [−1, 0] and we find for the second integral and for the third, respectively.

B LOCAL AND COUPLING MATRICES
We report here the expressions for the local matrices in the XDG-Laguerre formulations.For vertical edges  ∈ Γ where we used the properties  ∞  (  ) = 1 and ( ∞  ) ′ (  ) = −(1∕2 + ).

C UNKNOWNS AND TRIAL FUNCTIONS IN FLUX INTEGRAL
In the sixth integral of (25), i.e. the flux integral, expressions for the discrete unknowns and test functions read as follows.

Figure 5
Figure 5 Non-linear two-dimensional shallow water equations, Gaussian initial data.Contour plot (contour interval 0.05) of the numerical solution of the XDG scheme at the final time  = 0.75  in the outgoing case (above) and incoming case (below) for Gaussian initial data.Depth ℎ(, ,  ) (left), horizontal velocity (, ,  ) (center) and vertical velocity (, ,  ) (right). = 50.The last Laguerre mode is placed at 44.52 , only the region [0, 10 ]×[0, 20 ] is shown here.The black cross represents the center of the initial profile, while the red dashed line denotes the finite/semi-infinite interface.
Left: relative errors in the semi-infinite region as a function of  for  = 50.Right: position of the last Laguerre node (  ) and distance between the first two Laguerre nodes () for the different choices of .

Table 5
Non-linear shallow water equations with an absorbing layer in the semi-infinite region, wave train test, non uniform-grid vs. uniform-grid DG.Relative errors of the single-domain DG model on a non-uniform grid at the final time with respect to a single-domain DG discretization in [0, 12 ] × [0, 10 ] with  ′  = 5  .
, under the assumption that  is diagonal, ∇ = [  ∕,   ∕]).In Ω − (  = 1, … ,   ) this gives ,∞ replacing    ,  ,    ,∞ replacing    ,  , and ∫ Ω(ℎ ) ⋅ ∇ ℎ   = − ∫ ,   and a similar expression for   = ∞ with modifications as above.Analogous computations on the right boundary  =   , where  = [1, 0], lead to Finally, we consider the right hand side of (25), which we report here for convenience: computed using suitable quadrature rules.The second and third integrals take different forms depending on whether the edge  belongs to the bottom, right or left boundary.If   = 1, … ,   we have, respectively, If instead   = ∞, we only consider the right and the left boundary, for which we obtain for   = 1, … ,   ; the expression for   = ∞ is identical as long as one replaces Δ with 1∕.