Stable discontinuous staggered grid in the 4 th -order

1. TXY= TYX TYZ= TZY TZX= TXZ TXX , TYY , TZZ U V W European Geosciences Union, General Assembly 2010 Vienna, Austria, 02 – 07 May 2010 Stable discontinuous staggered grid in the 4th-order finite-difference modeling of seismic ground motionKristek1,2, J., Moczo1,2, P., Galis1,2, M. 1 - Comenius University Bratislava, Slovakia 2 - Slovak Academy of Sciences, Bratislava, Slovakia www.nuquake.eu reasonable to use a discontinuous spatial grid with a finer part, with a grid spacing h,covering the upper part of the model, and a coarser part, with a grid spacing H>h, covering the lower part of the model. The total number of grid points in the discontinuous spatial grid can be significantly smaller than that in a uniform grid. This simple idea led modelers to implement discontinuous grids in their numerical modeling. A finer spatial grid near the free surface can be also useful if the free-surface topography is to be included. A number of algorithms to include discontinuous spatial grid have beendeveloped. They mainly differ in the allowed grid ratio H/h and the way they interpolate values at the missing grid positions in the coarser grid. In general, one could guess, the larger is the grid ratioH/h, the larger is the possibility of inaccuracy and, mainly, instability with an increasing number of time steps due to a numerical noise that is generated at the contact of the finer and coarser grids. At the same time, it is obvious that the noise generation and possible instability cannot simply correlate with the grid ratio – they may strongly depend on the material and wavefield. The possibility of instability is only rarely explicitly addressed. The very obvious problem is that of the missing grid points. This is the reason for interpolation. The other, and apparently not so obvious, problem is how to update particle velocity and stress at those grid points of the coarser grid which coincide with the grid points of the finer grid. Abstract We present an algorithm of a stable discontinuous staggered-grid. The grid ratio H/h can be an arbitrary odd number. The key feature of the algorithm is the application of the Lanczos downsampling filter. We demonstrate the stability and accuracy of the algorithm for large number of time steps. Kristek, Moczo, Galis, 2010. Stable discontinuous staggered grid in the 4th-order finite-difference modeling of seismic motion, Submitted to GJI Introduction If the minimum wave speed in an upper part of a computational model is smaller than that in a lower part of the model it may be The velocity-stress discontinuous staggered grid The spatial grid cell in the 3Dvelocity-stress staggered grid Interior 4th-order velocity-stressstaggered-grid schemeon the finer grid, grid spacing h Interior 4th-order velocity-stressstaggered-grid schemeon coarser grid, grid spacing H Overlapping Zone Interior 2nd-order velocity-stressstaggered-grid schemeon the finer grid, grid spacing h Values obtained usingbilinear interpolationin the horizontal plane for the finer grid:interior 4th-order scheme for the coarser grid:values obtained using the Lanczos downsampling filter Boundary of the Coarser Gridin the Overlapping Zone Top: A vertical grid plane of the spatial velocity-stress discontinuous staggered grid. Only the grid plane with positions of the normal stress-tensor components in both the finer and coarser grids is shown. The vertical grid plane with the xy-stress-tensor components in both grids is analogical. There are also vertical grid planes with the finer grid only. Figure shows the simplest possible configuration of the discontinuous grid, that is, the case with H/h = 3. The grid ratio, however, can be an arbitrary odd number:H/h = 3, 5, 7, … Boundary of the Finer Gridin the Overlapping Zone The theoretically minimum wavelength that can propagate in the finer grid is2h. In the same medium, the theoretically minimum wavelength that can propagate in the coarser grid is 2H. Consequently, the wavelengths larger than 2h but smaller than 2H cannot propagate from the finer grid into the coarser grid. Therefore, the finer-grid field values at the red-circumscribed green positions have to be downsampled before they are used to update the field values at the internal grid positions of the coarser grid. The overlapping zone is placed in the medium with the larger S-wave speed. If the finer grid sufficiently oversamplesthe chosen minimum wavelengthin the faster medium, the 2nd-order scheme may be applied. The application of the 2nd-order scheme near the finer-grid boundary significantly reduces the spatial extent of interpolation.

2. Source parameters: Single vertical force acting at free surface Time function – Gabor signal where fp = 1 Hz,  = 1.5,  = 0, tS = 1.4 s 0.0001 0 -0.0001 30 s (15000 Time Levels) 40 s (20000 Time Levels) 50 s (25000 Time Levels) without Lanczos filtration with Lanczos filtration 77 s 35 000th time level 0.2 88 s 40 000th time level 0 99 s 45 000th time level -0.2 European Geosciences Union, General Assembly 2010 Vienna, Austria, 02 – 07 May 2010 Note on the non-reflecting boundariesin the overlapping zone Lanczos downsampling filter (e.g., Duchon 1979, Turkowski and Gabriel 1990) The Lanczos filter can be expressed as Interpolation Some of grid positions required by the bilinear interpolation are not available near a grid boundary. In such cases it is necessary to apply an extrapolation. where A is a scaling factor determined by condition Lanczos downsampling filter Application of the Lanczos filter at a grid position requires n neighboring finer-grid positionsfrom each side. Near the grid boundarythese positions are not available. At such points we do not apply the Lanczos downsampling filter. In other words, the direct finer-grid values are used for updating values at the coarser-grid points. We use PML. Therefore a potential effect (if any) is probably sufficiently eliminated. and Then the filtered value of a field variable Φ at the grid position ( i, j ) of the finer grid is obtained as The filtration is applied in a horizontal grid plane.This is reasonable because the filtered valuesonly enter the FD approximationsto spatial derivatives in the vertical direction. Numerical examples Time historiesof the vertical component of the particle velocity at the free surface Total energy in the grid (PML is not included) Example 1 – homogeneous halfspace Medium parameters : P-wave speed (at 1 Hz) = 5000 m/s S-wave speed (at 1 Hz) = 2600 m/s density = 2600 kg/m3 P-wave quality factor = 400 S-wave quality factor = 200 Snapshots of the vertical component at the free surface Discontinuous grid parameters: n = H/h = 3 finer grid h = 30 m dimensions: 337 x 481 x 19 coarser grid H = 90 m dimensions: 113 x 161 x 107.5 Example 2 – Grenoble valley, France Thestructurally complex model of the deep sedimentary valley was a subject of the international numerical exercise. The model and results are described in detail in the article by Chaljub et al. (2010). We used discontinuous grid with n = H/h = 5. NS component of the particle velocity at the free surface NS component of the particle velocity at the free surface (receiver position indicated by red triangle in figure right) Conclusions We developed an algorithm of the discontinuous gridfor the 3D 4th-order VS staggered-grid finite-difference modeling The key feature of the algorithmis the application of the Lanczos downsampling filter The ratio between the grid spacings of the coarser and finer gridscan bean arbitrary odd number The algorithm allowslarge numbers of time levels without instability The algorithm is directly applicable alsoto the displacement-stress staggered-grid finite-difference scheme The concept of the Lanczos downsampling filter is general and robust – its positive effect should not be dependenton a particular algorithm of the discontinuous grid References Chaljub, Moczo, Tsuno, Bard, Kristek, Käser, Stupazzini, Kristekova, 2010, BSSA, in press. Douchon, 1979, .J. Applied Meteorology 18 Turkowski, Gabriel, 1990, In Graphics Gems I, 147–165, Academic Press Acknowledgements This work was supported in part by the Slovak Research and Development Agency under the contract No. APVV-0435-07 (project OPTIMODE). We also gratefully acknowledge the funding by the European Union through the Initial Training Network QUEST (grant agreement 238007), a Marie Curie Action within the "People" Programme.