2D seismic wave propagation using the distributional finite-difference method: further developments and potential for global seismology

We present a time-domain distributional finite-difference scheme based on the Lebedev staggered grid for the numerical simulation of wave propagation in acoustic and elastic media. The central aspect of the proposed method is the representation of the stresses and displacements with different sets of B-splines functions organized according to the stag-2


INTRODUCTION
Over time, seismology has become increasingly dependent on high-performance computing resources to compute the seismic wavefield, and efficient algorithms for solving the wave equation have been developed.Different methods have been adopted within the earthquake seismology community and in exploration geophysics.
Many earthquake seismologists has adopted the Spectral Element Method (SEM), beginning with pioneering modeling of the crust (Tape et al. 2009) and uppermost mantle (Fichtner 2010) at local and regional scales, since then extended to many regions of the world (e.g.Magnoni et al. (2022)).Notably, the latest global tomographic models that unveil the structure of our planet's deep interior typically required millions of CPU hours for their construction (French & Romanowicz 2014;Bozdag et al. 2016;Lei et al. 2020).Indeed, until the advent of the SEM, seismic waveform modeling of the Earth's mantle at the global scale relied on the computation of synthetics using normal mode perturbation theory (Woodhouse & Dziewonski 1984;Li & Romanowicz 1996;Mégnin & Romanowicz 2000), the validity of which required strong assumptions on the strength and smoothness of the 3D heterogeneity.
The introduction of the SEM to global seismology (Komatitsch & Vilotte 1998;Komatitsch & Tromp 1999) opened up a new era, providing a comparatively efficient method for the accurate numerical computation of the entire seismic wavefield, in the presence of topography and potentially arbitrary lateral variations of 3D structure.
Given the computational challenges associated with full 3D simulations, where the computational Distributional finite-difference modeling 3 cost increases with the fourth power of frequency, several approaches have been introduced to approximate the wavefield under specific conditions.For example, when assuming an axisymmetric or partly axisymmetric model is a good enough approximation to compute the wavefield generated by localized structures.This led to the development of C-SEM, a "coupled mode-SEM" code (Capdeville et al. 2003a) in which SEM computations in the mantle were coupled through a DtN operator to fast 1D normal mode computations in the core, which was used in the development of the first SEM-based global radially anisotropic shear velocity models of the upper mantle (Lekić & Romanowicz 2011;French et al. 2013) and of the whole mantle (French & Romanowicz 2014).Capdeville et al. (2003b) extended the C-SEM method to an annulus of SEM coupled above and below with 1D modes, enabling the modeling of ultra-low velocity zones at the base of the mantle in 3D, using diffracted S waves at periods down to 8 s (Cottaar & Romanowicz 2012) or even shorter (To et al. 2005).Hypothesizing that it is more important to accurately represent the heterogeneity in the vertical plane containing the source and the receiver, Nissen-Meyer et al. (2014) proposed an efficient global 2.5D solver, AxiSEM, in which cylindrical symmetry is assumed around each source.The requirement of cylindrical symmetry in AxiSEM was later released in AxiSEM3D (Leng et al. 2019), which allows for a smooth variation of the structure in the direction perpendicular to the great circle plane.AxiSEM3D has already received much attention in the global seismological community, in particular for the modeling of complex structures in D" (Li et al. 2022).Meanwhile, SPECFEM3D GLOBE (Tromp et al. 2008), which makes no approximations on the structure, is open-source and has been significantly optimized over time, has now become the standard for full seismic wavefield computations at regional and global scales and has been used in particular for the development of global radially anisotropic shear velocity models GLAD-M15 (Bozdag et al. 2016) and GLAD-M25 (Lei et al. 2020).Other groups have developed their own SEM codes, notably in the framework of the CSEM (Collaborative Seismic Earth Model, e.g.Fichtner et al. 2018) or in the framework of regional and crustal seismic imaging (Trinh et al. 2019;Lu et al. 2018).These codes are however not open source.Although very higher-order spectral element methods, such as the 12th-order method, are more computationally efficient than the classical 4th-order SEM (Lyu et al. 2020), the overall computational efficiency of SEM decreases significantly when its order exceeds 20 due to the extremely small available time step, which limits the usefulness of super-high-order (≥ 20th-order) SEM for global-scale simulations.
Still, there is a need for more efficient methods to reach the higher frequencies required for the study of complex structures in the deep Earth, as the computational burden of the SEM is prohibitive for global waveform modeling at periods shorter than 10 s.In this study, we investigate an alternative to SEM for modeling seismic wave propagation at the global scale called the distributional finitedifference method (DFDM) introduced in Masson (2022).The proposed algorithm shall allow for do-main decomposition, medium heterogeneity, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grids, free surfaces, and fluid-solid interfaces.
SEM is the favored approach to model wave propagation at the planet's scale; on the other hand, the finite-difference method (FDM) is predominant in exploration geophysics and earthquake ground motion simulations in local surface sedimentary structures (e.g.Virieux et al. 2011;Moczo et al. 2014Moczo et al. , 2021)).This apparent discrepancy may be due to communities of practice, software availability, and, most likely, the nature of the seismic data and the model geometry.Active seismic imaging relies primarily on body wave analysis with a focus on reflections, while surface waves represent an essential part of the data processed in seismology and global tomography.Thus, SEM, especially suited for surface wave modeling, is preferred in seismology.FDM needs to be adapted to model surface waves accurately.Otherwise, it is remarkably simple and efficient: an attractive feature in exploration geophysics.Fluid-solid interfaces also need to be modeled accurately and influence the method to be chosen (De Basabe & Sen 2015).There are, of course, numerous exceptions to the general statements above.FDM has been successfully used to model wave propagation at the global and regional scale (Igel & Weber 1995;Igel et al. 1995;Chaljub & Tarantola 1997;Jahnke et al. 2005;Kawai et al. 2006;Li et al. 2014).Summation-by-part finite-difference schemes, named SBP-SAT approaches, which naturally account for the free surface, have reached recently an improved accuracy (Sjögreen & Petersson 2012;Petersson & Sjögreen 2015).Such an advanced SBP-SAT software is openly available to the seismology community (SW4, Petersson & Sjogreen 2017).SEM has also been employed for exploration and geotechnical applications (Zheng et al. 2014).The literature in FDM is abundant, and we refer the reader to review papers (Virieux et al. 2011;Moczo et al. 2014Moczo et al. , 2021)), and the references therein for a comprehensive review.
We focus on the staggered-grid finite-difference approach, which is relevant in this study.We highlight the pioneering work of Yee (1966) in electromagnetism, Madariaga (1976) in seismic source modeling, and Virieux (1986) in seismic wave propagation.In staggered-grid finite-difference approaches, the field variables are discretized on different grids, leading to centered finite-difference operators with improved accuracy for spatial derivatives of fields.Further, the staggered-grid approach effectively annihilates the spurious (non-physical) parasitic modes that are problematic in FDM (see, e.g., New et al. 1998;Dovgilovich & Sofronov 2015).Thus, this staggered-grid FDM is praised for its simplicity and computational efficiency and is widely used for the 3D modeling of seismic wave propagation (Graves 1993;Olsen 1995;Pitarka et al. 1998;Cotton et al. 1998;Matsushima 1998;Wald & Graves 1998;Kristek et al. 1999).The former Virieux algorithm (Virieux 1986) permits limited anisotropy (i.e.orthotropic) and can accommodate orthogonal curvilinear coordinates.Several staggered-grid FDM schemes have been introduced to account for general anisotropy and non-Distributional finite-difference modeling 5 orthogonal curvilinear grids, notably, partly-staggered grid scheme (e.g., Saenger & Bohlen 2004) and Lebedev scheme (e.g., Lisitsa & Vishnevskiy 2010;de la Puente et al. 2014).Even though they do not strictly fall in the staggered-grid FDM category, one should also mention the collocated approaches relying on the MacCormack scheme where left and right FD operators are employed to prevent the parasitic modes from appearing (Zhang & Chen 2006;Zhu et al. 2009;Lan & Zhang 2011;Sun et al. 2016).Bernth & Chapman (2011) emphasizes that the Lebedev scheme is more efficient and accurate than the partly staggered grid scheme.The stability and grid dispersion of the staggered-grid FDM is analyzed carefully in Moczo et al. (2000).Kristek et al. (2002) and Bohlen & Saenger (2006) addressed the difficulty of accurately modeling the free surface using staggered-grid FDM.Nonetheless, several approaches have been introduced to address this issue (Lombard et al. 2008;Zhang et al. 2012a,b;de la Puente et al. 2014;Gao et al. 2015;Sun et al. 2016;Shragge & Konuk 2020;Zang et al. 2021).An aspect of special importance in the finite-difference modeling in both seismic exploration and ground motion simulations is the implementation of material interfaces.It is not trivial in the case of a general shape and position in the uniform grid.The problem has been recently addressed by Mittet (2017), Jiang & Zhang (2021), Koene et al. (2022) and Moczo et al. (2022).In this work, we build up the Lebedev approach to construct our DFDM algorithm.
Besides SEM, numerous methods that decompose the computational domain in multiple elements exist.In particular, the classical finite-element method (FEM), (Lysmer & Drake 1972;Smith 1975;Marfurt 1984).Boundary-element methods (BEM) (Papageorgiou & Pei 1998;Ba & Gao 2017;Bouchon & Coutant 1994), finite-volume methods (FVM) (Dormy & Tarantola 1995;Dumbser et al. 2007) and Discontinuous Galerkin methods (DGM) (Käser et al. 2008;De Basabe et al. 2008;de la Puente et al. 2008;Étienne et al. 2010;Warburton 2013;Bonnasse-Gahot et al. 2018).Mixed FEM approaches following partially the staggered-grid strategy have been developed (Nédélec 1980;Bécache et al. 2002).Over the last decade, Isogeometric Analysis (IGA) (Hughes et al. 2005;Auricchio et al. 2010) has garnered significant attention within the field of structural engineering.This approach discretize Partial Differential Equations (PDEs) using the same type of functions employed in Computer Aided Design (CAD): Non-Uniform Rational B-Splines (NURBS).These NURBS not only provide a precise representation of the geometric aspects of the computational domain but also serve as the basis for defining the solution space for PDEs.Element-based approaches allow to mesh complex geological features and are well suited to model free surfaces accurately and discontinuities thanks to their weak formalism.However, they lead to several numerical challenges.They are more difficult to implement than FDM, often more expensive in computational time, and their accuracy depends on the quality of the meshing (Virieux et al. 2011).One should also mention pseudo-spectral methods (PSM) that are highly accurate but less flexible for representing media with complex geometries (Fornberg 1987;Carcione 1996).Masson (2022) introduced a distributional finite-difference method (DFDM) that combines some essential features of FDM, FEM, SEM, IGA, and DGM.DFDM decomposes the space into multiple elements with arbitrary sizes.Thus, a volumetric mesh may be employed to describe the model better whenever needed.Within the elements, DFDM has an algorithmic structure similar to FDM and relies on compactly supported operators, which ensures efficiency.Our algorithm relies on B-spline bases similarly to IGA with two important differences.First, DFDM does not require designing optimal quadrature rules (e.g.Zou et al. 2022) for achieving efficiency.The mass and stiffness matrices are assembled independently prior to the computations, while the element shapes and heterogeneity are accounted for approximately using weights associated with virtual grid points.This approach is wellsuited for geophysical applications where the exact structure of the medium in which waves propagate is rarely known precisely.Second, in DFDM, field variables are represented using different basis functions (similarly to staggered grid FDM).This has the dual benefit of reducing the number of degrees of freedom in the problem and eliminating parasitic modes that are notoriously problematic in wave propagation modeling (Virieux 1986;Lisitsa & Vishnevskiy 2010).Finally, in DFDM, the numerical wavefield is discontinuous between neighboring elements, both the displacement and the stress can be defined explicitly at the boundaries, similarly to DGM.
The algorithm proposed by Masson (2022) relies on a MacCormack type of strategy to remove the parasitic modes; all the displacements are represented in a first basis (with the basis function skewed to the left), and all the stresses are represented in a second basis (with basis functions skewed to the right).This algorithm also shares the structure of the partly-staggered-grid schemes with all the displacement components located at one grid position and all the stress components located at another grid position (e.g., Saenger & Bohlen 2004).Masson & Virieux (2023) proposed a DFDM scheme analogous to the classic algorithm from Virieux (1986) where field variables are represented in different bases organized according to the staggered grid: specific differential operators switching consistently from one basis to the other one when estimating a field and its spatial derivative are designed as we shall also see in this work (this operation is not performed in mixed FEM approaches).They showed that the DFDM algorithm is more accurate and faster than its FDM counterpart.Interestingly, the staggered approach permits the computation of the partial derivatives by applying 1D operators in a single spatial direction.
In contrast, the method of Masson (2022) requires using the 1D operators in the two spatial directions.
The switch between bases is still performed in both strategies.The two latter approaches also differ in the time scheme.Masson (2022) employs a centered second-order time integration scheme with a displacement formulation and Masson & Virieux (2023) employ a first-order leap-frog strategy.Last, the two studies account for material heterogeneity in different manners.Masson (2022), includes material heterogeneity using a mass lumping technique and Masson & Virieux (2023) discretizes the material properties on a virtual grid associated with the basis functions.For clarity and to outline the similarities with FDM, Masson & Virieux (2023) only considered a single isotropic element.Nonetheless, their algorithm applies to limited anisotropy (i.e., orthotropic material) and adapts to orthogonal curvilinear grid (e.g., Xie et al. 2002).
In this study, we expand the Masson & Virieux (2023) algorithm to the Lebedev grid, which consists of two Virieux grids.It is needed in seismology to model wave propagation in arbitrary anisotropic media and to work with non-orthogonal curvilinear coordinates.Further, it allows the meshing of complex geological features and accurately model interfaces, particularly free surfaces and fluid-solid interfaces with topography.A key aspect of this paper is the evaluation of the proposed algorithm for modeling seismic wave propagation through the Earth at the global scale.This manuscript is organized as follows.The first section starts with the second-order elastodynamic system of equations.We apply the non-orthogonal curvilinear coordinate transformation.Then, we discretize the field variables using the basis functions organized according to the Lebedev grid.
Thereafter, we briefly recall how a dual basis can be employed to approximate functions and their derivatives.We proceed with the derivation of the one-dimensional DFDM operators employed for approximating the partial derivatives.The 1D operators are then extended to 2D using a tensor product.The treatment of boundary conditions and interfaces is discussed together with the introduction of the 2D operators.Finally, we detail our 2D time-domain algorithm that models wave propagation.
In the second section, we present numerical examples demonstrating that the proposed numerical algorithm can accurately and efficiently model wave propagation through the Earth at the global scale.
We compare the results obtained with DFDM to those obtained using SEM.DFDM is significantly more accurate than SEM when using a similar setup.We finish with concluding remarks and future perspectives.

Governing equations
We consider wave propagation in a two-dimensional linear anisotropic elastic medium.The proposed algorithm relies on the second-order elastodynamic system of equations. where σ ij are the component of the stress tensor, ρ is the density, C ij can be inferred from elastic parameters of the linear Hooke law, ∂ t is the partial derivative with respect to time, ∂ i is the partial derivative with respect to spatial direction i.

Curvilinear transformation
To allow for curvilinear meshes that follow the topography of the interfaces, we apply a smooth, locally invertible coordinate transformation where x and y are the cartesian coordinates in the physical space and x ′ and, y ′ are the coordinates in the computational space.We illustrate such a transformation in Fig. 1.Using the transformation in Equation 2, we can write Equation 1 in the conservative form (Thompson et al. 1982;Appelö & Petersson 2009;Dovgilovich & Sofronov 2015) Distributional finite-difference modeling 9 ) ) ) Computational space

Physical space
Figure 1.Illustration of the geometrical mapping between the physical space coordinates (x, y) and the com- where and |J| is the determinant of the Jacobian matrix J.We need to transform the linear system above into a discrete linear system to obtain a numerical solution.

Discrete representation of the field variables
We discretize the field variables in Equation 3 following the approach introduced in Masson ( 2022) and Masson & Virieux (2023).We represent the wavefield using orthogonal basis functions that are staggered according to the so-called Lebedev grid for the following reasons: • Orthonormal bases maintain the self-adjoint antisymmetric nature of the wave equation after discretization (Masson 2022).This is essential to achieve numerical reciprocity and ensure the stability of the numerical scheme (subject to a CFL condition).
• Staggered basis functions ensure that the spurious(non-physical) parasitic modes associated with the centered finite-difference operators are effectively removed (see, e.g., New et al. 1998).
In two dimensions, the Lebedev grid consists of two Here, the upper indices indicate the basis in which the field variable are represented in the two spatial directions.For example, u 12 is represented using the basis B 1 in direction x and the basis B 2 in direction y.Each set of variables is organized following the staggered grid (Virieux 1986).For the first Virieux grid, we have and for the second Virieux grid we have where the Bkl ij (k = 1, 2) are the two-dimensional orthonormal B-spline basis functions.They are constructed from the one-dimensional basis functions Bk x ij (k = 1, 2) and Bl y ij (l = 1, 2) using a tensorial product.The orthonormal bases Bkl ij (k = 1, 2) may be constructed from arbitrary basis functions B kl ij (k = 1, 2) that we will refer to as the canonical bases.In this study, the canonical bases are taken as the two-dimensional B-spline bases detailed in Appendix A and pictured in Figure 2.

Approximation using dual bases
Before constructing the DFD operators, we show that the orthonormal bases employed for representing the field variables can be used to approximate functions and their derivatives.The orthonormal bases Bk (k=1,2) defined in Appendix C are self-dual.Thus, they provide us with a direct way to find the approximation f k (θ) of an arbitrary function f (θ), where θ is a spatial coordinate (see, e.g., Woźny 2013).We have where the approximation f k (θ) is represented in the basis Bk and the expansion coefficients f k i are obtained directly using the integrals In Equation 6, the function f k (θ) is an optimal approximation of f (θ) in the least square sense; it Notably, the approximation error vanishes and we obtain an exact result (i.e., ||f − f k || = 0) when the function f (θ)=f k (θ) can be represented in the basis Bk .We can approximate the derivative f ′ (θ) of the function f (θ) with a similar approach.We Distributional finite-difference modeling 13 have where f ′ k i are the expansion coefficients representing the approximation f ′ k (θ) of f ′ (θ) in the basis Bk .Let us notice that this basis for the derivative might not be the one used for approximating the function itself.By substituting f (θ) with f ′ (θ) in Equation 7and using integration by parts, we obtain the expansion coefficients where ( B) ′ are the derivatives of B. Equations 7 and 9 are employed to approximate the derivatives of the field variables and thus construct the DFD operators.Equations 7 and 9 are central in DFDM; they allow one to compute the value/derivative of a function in a basis different from the one in which they are represented.It is similar to what is done in staggered-grid FDM, where we take wavefield values on one grid and compute their derivative on the corresponding shifted grid.This is the key point regarding the staggered-grid strategy we promote.

Distributional finite-difference operators in 1D
We now construct the one-dimensional DFD operators that are the basic bricks of our algorithm.Following (Masson & Virieux 2023), they should act on the field variables expressed in the first basis and return their derivatives' representations in the second basis.We assimilate the function f (θ) introduced in the previous section to a discrete field variable f l (θ) represented in the orthogonal basis Bl .We assume that the discrete representation f l of f l (θ) is known.We can express f l (θ) and its where f l j are the expansion coefficients stored in the vector f l , Bl j are the basis functions, and ( Bl j ) ′ are their derivatives.Keeping in mind that the upper indices k and l indicate the basis in which the variables are represented; our objective is to obtain the approximations expressed in Bk of the function and its derivative in Equation 10.By substituting Bk i (θ) and f l (θ) in Equation ( 7) with their expressions in Equations (C.2) and (10a) we obtain the linear system where The mass matrix M kl and its Cholesky factorization L n are defined in Appendix B and Appendix C.
The operator T kl can be thought of as a basis transform; it acts on the vector f l representing the function f l (θ) in the basis Bl and returns the vector in the basis Bk .
By substituting Bk i (θ) and f ′ l (θ) in eq. ( 9) with their expressions in eqs.(C.2) and ( 10b) we obtain the linear system associated with the derivative where and The stiffness matrix K lk and the vectors p k and q k are defined in Appendix B. When the boundary values f l (θ − ) and f l (θ + ) are equal to the first and last elements of the vector f l , the linear system in Equation 13becomes and by definition, the first-order distributional finite-difference operator is The operator D kl acts on the vector f l representing the function f l (θ) in the basis Bl .It returns the Notice that the factorization employed here is different (but equivalent) to that employed in Masson (2022) where the DFD operators have an adjoint form and partly encompass the boundary values.
In this paper, the self-adjoint form of the global or assembled linear system is recovered by using an appropriate averaging of the boundary values between neighboring elements, as detailed in the next section.Notice that all the matrices involved in calculating the DFD operators in Equations 12 and 14 Distributional finite-difference modeling 15 are band-diagonal matrices because the B-spline canonical basis functions have compact support.It ensures the efficiency of the proposed algorithm.

Distributional finite-difference operators in 2D
We construct the 2D-DFD operators needed to compute the partial derivative of a function using a tensor product, i.e., by applying the 1D-DFD operators in the spatial dimensions x ′ and y ′ as detailed in Masson (2022).
Let f kl (x ′ , y ′ ) be one of the field variables in Equations 4 and 5.It is represented in the basis Bkl through its expansion coefficients F kl mn stored in the matrix F kl .Throughout this paper, we assume that the indices (i,j) of the arrays correspond to (row, column) and are associated with directions (x,y).
We have where the 2D basis functions Bkl mn (x ′ , y ′ ) = Bk xm (x ′ ) Bl yn (y ′ ) are constructed from the 1D basis functions Bk xm (x ′ ) and Bl yn (y ′ ).According to the staggered grid, one should represent the partial derivative ∂ ∂x f kl (x ′ , y ′ ) in the basis Bil (with l ∈ [1, 2] and i = 3 − l) and the partial derivative ∂ ∂y f kl (x ′ , y ′ ) in the basis Bkj (with k ∈ [1, 2] and j = 3 − k).We have where the expansion coefficients D il Xmn and D kj Ymn stored in the matrices D il X and D kj Y represent the approximations ∂ ∂x f kl (x ′ , y ′ ) and ∂ ∂y f kl (x ′ , y ′ ) in the bases Bil and Bkj .Notice that the function f kl (x ′ , y ′ ) and its partial derivatives are represented in the same 1D basis perpendicular to the partial derivative's direction.Following the approach in Masson (2022), the expansion coefficients representing the partial derivatives can be computed using where In the domain Ω, the internal boundary values in Equation 21shall be set to impose the boundary conditions or computed from the internal expansion coefficients Here, the brackets [] Ω indicate that the expansion coefficients are taken inside the current domain Ω.
To model a free surface, for example, one can take boundary values of the traction equal to zero when evaluating the partial derivatives of the stresses and use the internal values when evaluating the partial derivatives of the displacements.Similarly, for a domain with fixed boundaries, one can take boundary values of displacement equal to zero when evaluating the partial derivatives of displacements and use the internal values in Equation 22 when evaluating the partial derivatives of the stresses.
We average the boundary values between neighboring domains to connect multiple domains or elements.This continuity condition is needed to ensure numerical reciprocity and the stability of the numerical scheme (Masson 2022).Notice however, that the continuity of the wavefield is not satisfied numerically and the values to the left and right of an interface may be different.The partial derivatives in Equation 21 are evaluated in the domain Ω using of the current domain Ω.This operation can be thought of as an approximate basis transform.We have where the quantities within the square brackets [] Ω and [] Ω n are those associated with the current Ω and the neighboring Ω n domains.The matrices L n θ are defined in Equation C.1.The matrix T nn θ in Equation 24a is similar to the zeroth-order DFD operator in Equation 12 but acts on the boundary values taken from the neighboring domains [] Ω n .Accordingly, the matrix M nn θ in Equation 24b is similar to the mass matrix in Equation B.1 but the second basis functions B n θj (θ) are taken from the neighboring domain Ω n .When the wavefield is represented in the same basis along the interface between two elements (i.e., for conformal interfaces), the matrices T nn θ in Equation 23 may be omitted because T nn θ = I.To obtain a stable algorithm that satisfies numerical reciprocity, the boundary values must be uniquely defined between neighboring domains (Masson 2022).We must apply the same weighting scheme in neighboring elements.For example, one enforces the following relations between the five domains pictured in Figure 1: Here, [α a ] Ω denotes the weight associated with the face a of the domain Ω a .The continuity conditions in Equation 25 ensure that the global differential operators (that approximate the partial derivatives in all domains taken together) will be adjoint after applying the boundary conditions (Masson 2022).
Thus far, we derived the DFD operators approximating the partial derivatives in Equation 21 independently of the wave equation.We can use these operators as such to solve the elastodynamic system in Equations 3.

2D algorithm
Our algorithm is obtained by substituting the partial derivatives in Equations 3 with the DFD operators in Equation 21and by substituting the physical quantities by their discrete counterparts.In this work, we adopt the virtual point approach introduced in Masson & Virieux (2023).The material properties and quantities associated with the curvilinear mapping are simply evaluated at the virtual grid points, as detailed in Appendix D. We approximate the second-time derivative using a centered second-order finite-difference scheme.Thus, Equations 3 are solved using a recursive procedure.For clarity, we give the update expressions for the variables associated with the first Virieux grid in Equations 4.
The update expressions for the second Virieux grid in Equation 5 are similar and are obtained by performing the index substitutions (1 → 2; 2 → 1) in Equations (29-31) below.At each step, one first evaluates the stresses in each domain at time t using: The stresses are then transformed according to the curvilinear mapping using: In Equations 29-31, most of the computational burden lies in the evaluation of the partial derivatives D kl X () and D kl Y () using Equation 21.We give detailed expressions for the partial derivatives needed to evaluate Equations 29 and 31 in Appendix E.
In order to obtain a stable scheme, the coefficients α in Equation 23 must be complementary when evaluating the partial derivatives D kl X () and D kl Y () in Equations 29 and 31.If the coefficients 31.This rather simple strategy simply enforces numerical reciprocity and ensures that the discrete wave equation remains self-adjoint (Masson 2022).
Notice that more elaborated yet related approaches have been introduced in DGM (e.g., Grote et al. 2007).

VERIFICATION
To evaluate the accuracy of the proposed algorithm, we model 2D seismic wave propagation in the isotropic 1D radial ak135 Earth model (Kennett et al. 1995) using the DFDM algorithm.While our model is isotropic, the numerical tests demonstrate our algorithm's capability to effectively handle anisotropic media.This is because the coefficient representing the curvilinear mapping and the anisotropic moduli are the same, leading to comparable apparent numerical anisotropy (i.e., when assuming undeformed elements).We do not consider a cylindrical symmetry as in Nissen-Meyer et al.
(2014); nonetheless, such a symmetry can be envisioned and may be relevant for future applications.
We compare the results to those obtained using the SEM method.We present the setups for our DFDM and the SEM benchmark simulations and then discuss our results.

Example Setup
We use a simplified isotropic version of the ak135 Earth model with no attenuation and a single-layer crust.Inside the crustal layer, we use a linear velocity profile.It interpolates the values of the elastic properties of the ak135 model from the surface to the base of the lower crust.
To perform the DFDM simulation, we decompose the model into multiple elements, following a standard cube sphere approach (Komatitsch & Tromp 1999;Nissen-Meyer et al. 2008), as shown in Figure3b.There are 8 layers of elements in the numerical mesh.They match the principal interfaces of the ak135 model, that is, the MOHO at a depth of 35 Km, the transition zone between the upper and the lower mantle at a depth of 660 Km, the core-mantle boundary (CMB) at a depth of 2891.5 Km and the inner core boundary (ICB) at a depth of 5153.5 Km.To achieve a relatively homogeneous number of grid points per minimum wavelength and improve efficiency, we use two sub-layers in the 1.000 0.946 0.982 0.995 Table 1.Summary of the parameters employed in the different simulations, the measured errors, and CPU times.
The SEM-Ref simulation is the reference simulation with the highest accuracy (i.e., at least 11 grid points per wavelength).The SEM-1 simulation has a standard setup (with at least 5.5 grid points per wavelength).
The SEM-2 simulation uses a finer mesh than SEM-1 and is more accurate (with at least 6.9 grid points per wavelength).The DFDM simulation uses a coarse mesh (with at least 3 grid points per wavelength); it is more accurate and faster than SEM-1 and SEM-2.cases, the mesh matches the principal discontinuities, i.e., the MOHO, the 670 discontinuity, the CMB, and the ICB.In the SEM simulations, the mesh has two doubling layers, and the number of elements is adjusted as a function of the desired accuracy.The first doubling layer is just under the 670 Km discontinuity, and the second one is close to the middle of the outer core, between the CMB and the ICB.For the DFDM simulation, the number of elements is independent of the desired accuracy.We adjust the accuracy by varying the number of grid points within the elements.Because DFDM allows for non-conformal interfaces, the number of grid points is adjusted independently inside individual elements.The numbers of elements employed for the SEM-Ref, SEM-1, SEM-2 and DFDM are indicated in Table 1.
lower mantle, the outer core, and the inner core.In the horizontal direction, the mesh consists of 12 slices (i.e. 3 slices connected to each face of the central cube).Thus, the numerical mesh contains 93 elements in total.We set the number of grid points inside the elements so that there are at least 3 grid points per wavelength (using the grid spacing associated with the virtual grid points in Appendix D).
We shall see that this gives an accuracy comparable to a sufficiently accurate SEM simulation.The enhanced accuracy of DFDM is attributed to different factors.First, as opposed to SEM, no approximation is used in the mass matrix, which gives the differential operators a spectral nature.Second, the numerical solution is highly continuous, thanks to the B-spline bases.Because fewer elements are needed, this reduces the number of interfaces at which the numerical wavefield is discontinuous.Last, mixed conditions are employed to share the boundary values between neighboring elements.Because DFDM allows for non-conformal interfaces, the number of grid points in different directions does not need to be the same in neighboring elements.Finally, we assign the model parameters by evaluating the ak135 model at the virtual grid points in Equation D.1.
For the simulations, we consider a seismic event located at the north pole, at an intermediate depth Distributional finite-difference modeling 23 of 50 km, to mitigate the amplitudes of the surface and body waves.We represent the seismic source using the moment tensor point source in Equation D.3 with . We take the source time function as the Ricker wavelet with peak frequency f 0 = 2.5 × 10 −2 Hz and an estimated maximum frequency f max=6.25×10−2 Hz.
We recorded the wavefield at 20 receivers at the surface of the Earth.The receivers represented by triangles in Figure 4 are equally spaced and buried at a depth of 10 m.Their epicentral distances vary from 9 to 180 degrees.

Comparison with SEM
To assess the accuracy of the DFDM algorithm, we compare the DFDM solution to SEM solutions.
The SEM method is well established and known to produce sufficiently accurate results for the setup considered.It is very close to the analytical when a sufficient number of grid points per wavelength is employed.We ran three SEM simulations, we used the SPECFEM2D current source (Komatitsch et al. 2012) published under the GPL3 license.The first SEM simulation we will refer to as SEM-Ref is our reference simulation.It is highly accurate and relies on a fine mesh with at least 11 grid points per wavelength.The second SEM simulation, which we will refer to as SEM-1, employs a relatively standard setup with at least 5.5 grid points per wavelength.The third simulation, which we will refer to as SEM-2, is more accurate than SEM-1 and employs a finer mesh with at least 6.9 grid points per wavelength.
Similarly to the DFDM simulation, the SEM simulations mesh pictured in Figure 3a matches the principal interfaces of the ak135 model.The three SEM meshes employed for the simulations have two doubling layers to obtain a relatively homogeneous number of grid points per minimum wavelength.
The doubling layer is located just under the 670 Km discontinuity.The second one is close to the middle of the outer core, between the CMB and the ICB.We give the total number of elements for SEM-Ref, SEM-1, and SEM-2 in Table 1.They are adjusted to achieve the desired accuracy with a minimum number of grid points per wavelength.Notably, the SEM-1 mesh has a single layer of elements inside the crust, while SEM-Ref and SEM-2 use two layers of elements.Thus, the crust is modeled much more accurately in the SEM-Ref and SEM-2 simulations.In the SEM simulations, we assign the model parameters by evaluating the ak135 model at the grid points.We use the same time step in all simulations.It is small enough to minimize the numerical dispersion associated with the time integration scheme.It is close to the stable time step of reference simulation SEM-Ref, which has the most restrictive stability condition.The stable time step in DFDM depends on the basis functions employed and may be adjusted by varying the spacing between the basis functions toward the sides of the elements.With the setup we employ here, the stability condition is slightly more restrictive for DFDM than for SEM.However, because DFDM requires less grid points per wavelength to achieve a certain accuracy, the stable time step is larger in the DFDM simulations.

Results and discussion
We first compare the results obtained using DFDM to those obtained for the reference solution SEM-Ref, which is assumed to be nearly perfectly accurate.In Figure 4 we present snapshots showing the wavefield obtained for DFDM and SEM-Ref.Even though DFDM uses as few as 3 points per wavelength, we observe no visible difference between the two wavefields, which highlights the high accuracy of DFDM when using coarse meshes (Masson 2022;Masson & Virieux 2023).Notice the P to S and S to P wave conversions at the solid-fluid interfaces are accurately modeled by DFDM, which uses the same scheme in the entire domain, while SEM uses a specific solid-fluid coupling approach.This is due to the specific choice of the boundary coefficients α in Equation 23, as detailed in Appendix E. Generally, when dealing with interconnected purely solid or fluid elements, the choice of alpha values is relatively arbitrary.However, in the case of solid-fluid coupling interfaces, the value of alpha is fixed.On the one hand, we use the internal value of the displacement in the solid elements (i.e.α = 1 ).On the other hand, we take the internal value of the pressure in the fluid elements (i.e.α = 0).This specific choice enforces the continuity of the pressure and gives the most accurate results.
In Figures 5a and 5b we superimpose the seismograms simulated by DFDM and SEM-Ref at the 20 receivers.There is no noticeable difference between the seismograms obtained using DFDM and SEM-Ref, both for the body waves and for the strong dispersive surface waves.This is valid for all epicentral distances, which is remarkable and confirms the high accuracy of DFDM.
To quantify the error, we measured the following quantities for our DFDM simulation: where u and u ref are the vectors containing the displacement components at all times at all receivers for the considered and the reference solutions.
• The standard deviation of the error normalized by that of the reference solution where σ denotes the standard deviation.
where Cov denotes the covariance.
We show values of the errors in Table 1.These errors may appear relatively large, but the measured Pearson correlation coefficient ρ = 0.995 indicates an excellent correlation between the DFDM and the SEM-Ref simulations (ρ = 1 corresponds to a perfect correlation).In the DFDM simulation, the maximum difference and its standard deviation are 18% and 10%.
To investigate the relative accuracy of DFDM with respect to SEM when a standard setup is used, we compare the differences between the seismograms obtained using DFDM, SEM-1, SEM-2, and SEM-Ref in Figures 5c and 5c.We observe that SEM-1 has the most prominent error, while DFDM has the smallest error.Notice that the error is principally associated with the surface waves that are notoriously difficult to model accurately.The surface waves are primarily sensitive to the velocity in the crust, and we know that DFDM is more accurate than SEM when using a similar number of points per wavelength (Masson 2022).Thus, we anticipate that DFDM is more accurate than SEM-1 because DFDM and SEM-1 employ one layer of elements in the crust with 5 grid points in the vertical direction.What is more remarkable is the superior accuracy of DFDM with respect to SEM-2, because SEM-2 uses two layers of elements inside the crust (which corresponds to 9 grid points in the vertical direction).In contrast, DFDM uses only one layer of elements (with 5 grid points in the vertical direction).Further, DFDM uses only 3.5 points per wavelength in the crust in the horizontal direction, while SEM-2 uses almost 7.These observations show that DFDM is remarkably accurate and adapted to model surface waves.Such waves are notoriously challenging in global and regional seismological modeling.
To better quantify the numerical error, we also measured the error, its standard deviation, and the Pearson correlation coefficient for SEM-1 and SEM-2.We show values of the errors in Table 1.In SEM-1, we use the recommended setup for SEM, i.e., at least 5 points per wavelength (Komatitsch & Vilotte 1998)), and we expect no significant numerical dispersion.We find a maximum error of 33%, and its standard deviation is 19%.Thus the typical error is twice as large in SEM-1 as in DFDM.
The measured Pearson correlation coefficient ρ = 0.946 also indicates a degraded correlation with the reference solution.It is remarkable given that the grid size (i.e., the number of grid points) in SEM-1 is about 3.5 times larger than in DFDM.For SEM-2, where there are at least 6.9 points per wavelength, we find a maximum error of 19%, a standard deviation is 12%, and the Pearson correlation coefficient is ρ = 0.982.These values show that the error in SEM-2 is slightly larger than in DFDM.

Distributional finite-difference modeling 29
We consider, however, that these two simulations have comparable accuracy.Looking at the grid dimensions in DFDM and SEM-2 suggests that SEM requires roughly 5 times more grid points to achieve an accuracy similar to DFDM in 2D.Because DFDM requires more computational effort than SEM to update the wavefield, comparing the CPU time measured for the different simulations is interesting.The values given in Table 1 are only indicative because the specfem2D code used for the SEM simulation is a well-established, and optimized software while our DFDM code is still at an early development stage.Nonetheless, we observe that the CPU time for DFDM is slightly shorter than for SEM-1 and significantly smaller than for SEM-2.Thus DFDM appears substantially faster than SEM to achieve comparable accuracy.

DISCUSSION AND OUTLOOK
We showed that DFDM is significantly more accurate than SEM in numerical simulations of the global wave propagation.Through measurements of the CPU time, we observed that the gain in accuracy achieved with DFDM also translates into better computational efficiency.It is especially striking that we can model the surface waves with fewer points per wavelength than in SEM.
For further applications and to better estimate the net gain in computational efficiency with respect to SEM, our algorithm needs to be expanded from two to three dimensions.This is the subject of a forthcoming paper (Lyu et al., 2023a, in prep.).Lateral variations in the velocity structure and topography will also need to be accounted for in real applications.In this study, we limited our numerical tests to an axisymmetric model without topography for an easier comparison with SEM.This is because material heterogeneity is accounted for using slightly different approximations in SEM and DFDM, and obtaining a trustworthy reference solution in heterogeneous media is delicate.Nonetheless, Masson & Virieux (2022) demonstrated that DFDM can accurately account for heterogeneity within the elements, including sharp interfaces.Furthermore, Masson (2022) showed that DFDM is in agreement with SEM in the case of strong topographic variations.Accounting for the 3D structure of the Earth's crust is especially challenging in global seismology.Our numerical examples demonstrate that DFDM allows us to employ non-conformal meshes to efficiently handle the strong contrast between the crust and the upper mantle.This is more delicate in SEM, where doubling layers or similar approaches are needed to refine the numerical grid with respect to depth or as a function of velocity fluctuations.Furthermore, DFDM also allows us to use small elements as SEM, and similar strategies may be adopted.Thus, DFDM brings additional flexibility to accurately account for structures with complicated geometries.In practice, a 3D crust may be implemented by varying the upper element thickness to accommodate Moho topography; additional elements may be added to account for sharp lateral variations.
The data underlying this article will be shared on reasonable request to the corresponding author (yder masson@berkeley.edu).The SPECFEM software can be obtained at the following address https: //geodynamics.org/resources/specfem2d and is described in Komatitsch & Vilotte (1998).we obtain the basis functions of the B-spline basis β p with polynomial order p recursively using the Cox-de Boor formula In our examples, we employ pairs B 1 /B 2 of B-spline bases where In this case, the basis functions are shifted similarly to the staggered finite-difference grid, as illustrated in Figure A1a.
The two-dimensional canonical B-spline bases are constructed from the one-dimensional bases ) and B l y (y ′ ) (l = 1, 2) defined in direction x ′ and y ′ using a tensor product.We have Such two-dimensional canonical basis functions are illustrated in Figure 2a.

APPENDIX B: DEFINITIONS: MASS, STIFFNESS AND BOUNDARY MATRICES
We define the matrices and vectors needed to construct the orthonormal basis functions and the DFD operators.We build those matrices and vectors from an arbitrary pair of canonical bases B 1 /B 2 .In this paper, we employ staggered B-spline bases, as detailed in A and illustrated in Figure A1.

Distributional finite-difference modeling 39
The mass matrix is: Here, B k θi (θ) and B l θj (θ) are the 1D canonical basis functions presented in Appendix A. The subscript θ def = x ′ or y ′ denotes the direction of the computational domain x ′ or y ′ along which M θ is computed.
and [y ′ − , y ′ + ].The superscripts kl indicate the bases from which M kl θ is computed ( B 1 or B 2 in Appendix A ).When the basis functions have compact support, the mass matrix is band diagonal, ensuring the computations' efficiency.Further, if the two bases are identical (i.e., k = l), then the mass matrix M kk θ is symmetric.It is employed to construct the orthonormal bases introduced in the next section.
The stiffness matrix is The boundary matrix Q θ in eq.(B.2) is: Notice that Equation B.3 is specific to the B-spline bases considered in this paper because all the basis functions are equal to zero at the extremities of the interval considered, except the first one and the last one.More general expressions are available in (Masson 2022).

APPENDIX C: ORTHONORMAL BASES: ORTHONORMAL B-SPLINE BASES
To construct the orthonormal bases Bk θ employed to represent the field variables, we start by computing the mass matrices M kk in eq.(B.1) associated with the canonical bases B k (for k = 1, 2).Then, we compute the Cholesky factorizations L k θ of the mass matrices.we have and we define the basis functions Bk x i (x ′ ) of the orthonormal bases Bk x (x ′ ) as the linear combinations x , v where Bk x (x ′ ) (k = 1, 2) and Bl y (y ′ ) (l = 1, 2) are the one-dimensional bases defined in directions x ′ and y ′ .

APPENDIX D: DISCRETIZATION OF THE MATERIAL PARAMETERS USING VIRTUAL GRIDS
To discretize the material parameters, we adopt the approach introduced in Masson & Virieux (2023).
We define virtual grids where each basis function is associated with a unique grid point.The grid coordinates are computed from the 1D basis in the equation C.2; we have The coordinates (x k i , ŷl j ) define the virtual grids G kl associated with the basis Bkl in each domain.They can be employed to assign the values to the material parameter arrays.In the proposed algorithm, we use the following: where ab = xx, xy, yx, yy (D.2b)

MOMENT TENSOR POINT SOURCE REPRESENTATION
We represent the seismic source using the moment tensor point source where s(t) is the source time function, M is the moment tensor, δ is the delta function and x s is the source position vector.Practically, we implement the seismic source by adding the components of the moment tensor to the stresses in Equation 29.We have where we take the delta function as  following: 1) the weights are complementary to those employed for the displacement (this is needed to obtain a stable scheme); 2) If the axes of the computational domain aren't oriented similarly in the two domains, the components of the transformed stress need to be rotated.
The partial derivatives of the stress components needed to evaluate the right-hand-side of Equation 31 are: where the 2D arrays S ij θθ contain the expansion coefficients representing the stresses inside the current domain Ω and the 1D arrays s i θθ () contain the expansion coefficients representing the stresses at the boundaries of the current domain Ω.We compute the boundary values in Equation E.4 using: Finally, when the parametric coordinates don't have the same orientation in the two domains, one needs to rotate the boundary values of the transformed stresses.For example, using where J is the Jacobian matrix that aligns the parametric coordinates in the two domains.
Following Masson (2022), we set the coefficients α θ ′ ± in Equations E.2 and E.5 depending on the nature of the interfaces.We have: • α θ ′ ± = 1 at all interfaces surrounding the computational domain to model a free surface.• α θ ′ ± = 1/2 at all interfaces between domains that are both solid/elastic or fluid/acoustic.• α θ ′ ± = 1 inside the solid/elastic domain and α θ ′ ± = 0 inside the fluid/acoustic domain, at all solidfluid interfaces.

Figure 2 .
Figure 2. (a) Representation of the canonical B-spline bases employed to construct the orthonormal bases that represent the field variables in the proposed DFD algorithm.The stresses σ 11 xx , σ 11 yy and σ 11 xy are represented using the basis pictured in gray.The stresses σ 22 xx , σ 22 yy and σ 22 xy are represented using the basis pictured in red.The displacements u 12 x and u 12 y are represented using the basis pictured in blue.The displacements u 21 x and u 21 y are represented using the basis pictured in green.(b) Representation of the analogous finite-difference staggered grid.The complete grid to the left is called the Lebedev grid.It decomposes into the two staggered grids to the right (i.e., as in Virieux (1986)).
kl contains the internal expansion coefficients representing the function f kl (x ′ , y ′ ) inside the domain Ω.The indices kl indicates the basis B kl in which f kl (x ′ , y ′ ) is represented inside Ω. f j (x ′ ± ) contains the expansion coefficients representing the boundary value f kl (x ′ ± , y ′ ) along the faces separating the domains Ω x± and Ω in Figure 1.The index i indicates the basis B j y in which f kl (x ′ ± , y ′ ) is represented inside Ω. f i (y ′ ± ) contains the expansion coefficients representing the boundary value f kl (x ′ , y ′ ± ) along the faces separating the domains Ω y± and Ω in Figure 1.The index j indicates the basis B i x in which f kl (x ′ , y ′ ± ) is represented inside Ω. D ij θ , b i θ are the 1D-DFD operators introduced in Section 2.5.
) Here, the brackets' subscripts indicate the domain from which the boundary values are taken.[] Ω corresponds to the current domain, and [] Ω x ′ − , [] Ω x ′ + , [] Ω y ′ − and [] Ω y ′ + correspond to the neighboring domains to the left, right, bottom and top, respectively, as pictured in Figure 1.The weighted average coefficients α x ′ − , α x ′ + , α y ′ − and α y ′ + are associated with the left, right, bottom and top faces of the domain Ω, respectively.The matrices T nn θ in Equation 23 are associated with the faces.They transform the boundary values taken from the neighboring domains Ω n and return their representation in the basis

Figure 3 .
Figure 3. Numerical mesh employed for the SEM simulations (left) and the DFDM simulation (right).In both

Figure 5 Figure 6 .
Figure 5. a) Seismograms (vertical component of the displacement vector) simulated at the 20 receivers represented by triangles in Figure 4.The solid black and the dashed red lines show the seismogram modeled in the reference SEM-Ref and the DFDM simulations.b) Magnified view of the seismogram recorded at an epicentral distance of 90 degrees inside the black box in a).c) The errors of the seismograms obtained by the SEM-1, SEM-2 and DFDM simulations with respect to seismograms obtained by the reference SEM-Ref simulation.

Figure A1 .
Figure A1.Illustration of the basis functions B 1 and B 2 and, their orthogonal counterparts B1 and B2 for p = 5.
basis functions B θ .
denotes the elements of the inverse of the matrix L k θ .The basis function in eq.(C.2) are orthonormal in the sense that Mkk θ def = ( Mkk θ ) ij = ∫ θ+ θ− Bk θi (θ) Bk θj (θ)dθ = δ ij def = I, (C.3) where δ ij is the Kronecker Delta and I is the identity matrix.This result is easily verified by inserting Equation C.2 in Equation C.3 and using Equation B.1.Notice that we chose the Cholesky factorization in Equation C.1 for its simplicity and efficiency.However, any factorization of the form M k θ = A k θ ⋅ A k θ T may be employed to construct the proposed algorithm.Equation C.3 shows that the orthonormal bases Bk θ (θ) are their own dual bases, by definition.Finally, we construct the two-dimensional orthonormal B-spline bases employed to represent the discrete field variables (v 21 kl s ) ij = ∬ ϕ kl ij (x)δ(x − x s ) dx dy = ϕ kl ij (x s ).(D.7)In the examples, the source time function employed is the Ricker wavelet in Equation 33.

Figure
Figure A2.(a) Workflow for computing and exchanging the boundary values between two neighboring domains Ω 1 and Ω 2 .Here, we consider the case where the right face of domain Ω 1 coincides with the left face of domain Ω 2 .In both domains, we first evaluate the displacement U ij x at the boundary using Equation E.3.We then send the boundary values [u j x (x ′ )] Ω 1 to the neighboring domain.During the exchange, the transformation matrix T ij y is applied to obtain the boundary values in the basis of the receiving domain.Finally, the internal and external boundary values are weighted (according to the weights on the arrows) and summed to obtain the averaged boundary value in Equation E.2.(b) is similar to (a) but for the stresses.The steps are the same except the 5d) where the internal boundary values [s i θθ ()] Ω are taken from the current domain Ω and the external Distributional finite-difference modeling 45 boundary values [s i θθ ()] Ω θ ′ ± are taken from the neighboring domains.The matrices T ij θ perform basis transformations to represent the external boundary values in the same basis as the displacement inside the current domain Ω.The weight α θ ′± controls whether the boundary values s i θθ () are taken inside (α θ ′ ± = 0), outside (α θ ′ ± = 1)or in between (0 < α θ ′ ± < 1) the current and the neighboring domains.In Equation E.5, the boundary values are computed in all domains Ω using: xx, yy, xy; j = 1, 2; k = 3 − j) xx, yy, xy; j = 1, 2; k = 3 − j) xx, yy, xy; i = 1, 2; l = 3 − i) (E.6c)[s i θθ (y ′ + )] Ω = [d i T x ⋅ S il θθ ]Ω (θ = xx, yy, xy; i = 1, 2; l = 3 − i) (E.6d) senting the components of the moment tensor point source M ij in the basis Bkl , Notice that the displacements at the two previous time step needs to be stored in memory while the stresses are only used as a intermediate step to update the displacements.The arrays above are as follow:Σ kl ij are 2D arrays with dimensions (N k x × N l y )containing the values or expansion coefficients representing the stresses σ ij (t) in the basis Bkl ,U kl i are 2D arrays with dimensions (N k x × N l y )containing the values or expansion coefficients representing the displacements u i (t) in the basis Bkl , ′ (f kl (x ′ , y ′ )) in the basis Bkl , ′ (f kl (x ′ , y ′ )) in the basis Bkl , M kl ij (t) are 2D arrays with dimensions (N k x × N l y ) containing the values or expansion coefficients reprej ) evaluated at the virtual grid points (x k i , ŷl j ) associated with the basis Bkl , ○ denotes the element-wise or the Hadamard product, ⊘ denotes the element-wise or the Hadamard division.
j [s kl xy ] i,j [s kl xy ] i,j [s kl yy ] i,j xx ] i,j [s kl xy ] i,j [s kl xy ] i,j [s kl yy ] i,j