Closed jwbishop closed 1 year ago
EtienneBachmann
commented
1 year ago Hi, if you think the bug may be related to inaccuracies in the receiver localization, you can try to turn the flag "DO_BRUTE_FORCE_POINT_SEARCH" to .true. in the setup/constants.h, recompile and rerun. A few years ago, a kd-tree algorithm was introduced to speed up that point localization task, and became the default option. However, I noticed this was producing sometimes inaccurate results.
danielpeter
commented
1 year ago hi Jordan,
there are two main changes with respect to the older version you compare against:
pressure is now also output for elastic elements. the older version had no pressure output for elastic domains, thus pressure seismograms for receivers in the elastic domain would be just zeroed out. the new version now computes pressure with the following recipe, as commented in the source code file compute_interpolated_dva.f90 :
! for an elastic element:
!
! isostatic stress or pressure: p = - 1/3 trace(sigma)
! (corresponds to hydrostatic pressure in fluids)
!
! from L. S. Bennethum, Compressibility Moduli for Porous Materials Incorporating Volume Fraction,
! J. Engrg. Mech., vol. 132(11), p. 1205-1214 (2006), below equation (5):
! for a 3D isotropic solid, pressure is defined in terms of the trace of the stress tensor as
! p = -1/3 (t11 + t22 + t33) where t is the Cauchy stress tensor.
!
! to compute pressure in 3D in an elastic solid, one uses: pressure = - trace(sigma) / 3
!
! sigma_ij = lambda delta_ij trace(epsilon) + 2 mu epsilon_ij
! = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_ij
!
! sigma_xx = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_xx
! sigma_yy = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_yy
! sigma_zz = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_zz
!
! pressure = - trace(sigma) / 3 = - (lambda + 2/3 mu) trace(epsilon) = - kappa * trace(epsilon)
!
! this routines limits the pressure computations to: non-anisotropic, non-attenuation casethis explains the difference in peak absolute pressure for receivers in the elastic halfspace. note, this is not a bug but a new feature in the new version :)
the station locations in STATIONS_FILTERED: the formatting of the file output values has changed (due to the handling of Cray outputting to ascii in a weird way). I would assume, these differences might just be rounding errors introduced when going from e22.10 (old version) to an f24.12 (new version, based on your X-range 652'000 to 662'000). so your difference of +/- 0.03 is probably within rounding errors for the large X-range, whereas the Y- and Z-range are smaller and less prone to numerical rounding and this updated output format.
so from my point of view, you're all good. the simulation is doing the right thing, only the output formats changed a bit.
rmodrak
commented
1 year ago Hi Jordan, Daniel, Etienne,
Very helpful, thanks for this issue and for the useful explanations. Below is an attempt at context, knowing that Jordan has a background primarily in acoustics, but has been working recently with seismologists.
Like other linearized wave solvers, SPECFEM acoustic and elastic simulations ignore the initial stress state. For acoustic media, there is little confusion, because it is common in that context to speak about pressure as the deviation from some background state. For elastic media, however, I think it more typical to speak about pressure as the absolute "confining pressure", without reference to some initial state, which differs I think from how SPECFEM works and from what Daniel describes above.
Please correct me if I'm wrong in any details-- regardless I think the bottom line is that careful explanations of terminology could be warranted if any of us were to write something like this up.
Best regards, Ryan
On Fri, Nov 11, 2022 at 6:21 AM daniel peter @.***> wrote:
hi Jordan,
there are two main changes with respect to the older version you compare against:
pressure is now also output for elastic elements. the older version had no pressure output for elastic domains, thus pressure seismograms for receivers in the elastic domain would be just zeroed out. the new version now computes pressure with the following recipe, as commented in the source code file compute_interpolated_dva.f90 :
! for an elastic element: ! ! isostatic stress or pressure: p = - 1/3 trace(sigma) ! (corresponds to hydrostatic pressure in fluids) ! ! from L. S. Bennethum, Compressibility Moduli for Porous Materials Incorporating Volume Fraction, ! J. Engrg. Mech., vol. 132(11), p. 1205-1214 (2006), below equation (5): ! for a 3D isotropic solid, pressure is defined in terms of the trace of the stress tensor as ! p = -1/3 (t11 + t22 + t33) where t is the Cauchy stress tensor. ! ! to compute pressure in 3D in an elastic solid, one uses: pressure = - trace(sigma) / 3 ! ! sigma_ij = lambda delta_ij trace(epsilon) + 2 mu epsilon_ij ! = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_ij ! ! sigma_xx = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_xx ! sigma_yy = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_yy ! sigma_zz = lambda (epsilon_xx + epsilon_yy + epsilon_zz) + 2 mu epsilon_zz ! ! pressure = - trace(sigma) / 3 = - (lambda + 2/3 mu) trace(epsilon) = - kappa * trace(epsilon) ! ! this routines limits the pressure computations to: non-anisotropic, non-attenuation case
this explains the difference in peak absolute pressure for receivers in the elastic halfspace. note, this is not a bug but a new feature in the new version :)
- the station locations in STATION_FILTERED: the formatting of the file output values has changed (due to the handling of Cray outputting to ascii in a weird way). I would assume, these differences might just be rounding errors introduced when going from e22.10 (old version) to an f24.12 (new version, based on your X-range 652'000 to 662'000). so your difference of +/- 0.03 is probably within rounding errors for the large X-range, whereas the Y- and Z-range are smaller and less prone to numerical rounding and this updated output format.
so from my point of view, you're all good. the simulation is doing the right thing, only the output formats changed a bit.
— Reply to this email directly, view it on GitHub https://github.com/SPECFEM/specfem3d/issues/1557#issuecomment-1311689396, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABCGSSVEFYGYJUPIZMPIWODWHZB5JANCNFSM6AAAAAAR5APGYI . You are receiving this because you are subscribed to this thread.Message ID: @.***>
danielpeter
commented
1 year ago hi Ryan,
that's an excellent remark. the pressure output in SPECFEM is only the pressure induced by the (elastic) deformation due to the wave displacement. the confining pressure in elastic domains, or initial hydrostatic pressure in acoustic domains, is ignored here.
thus, in-situ stress comparisons would need different handling. somewhat related, we would also need to change wave propagation under induced stresses, as mentioned in Jeroen & Jeannot's paper. it would be great to include this in future versions... :)
jwbishop
commented
1 year ago Hi All,
Thanks for the many helpful comments. I have since confirmed that we do get a waveform match between the version of SPECFEM3D mentioned above and the current devel branch.
There is still a small offset for many of the receiver locations relative to the previous locations when using the brute force search instead of the k-d tree, but this difference appears neglible for all but a few receivers, which may now be on the interface.
In relation to conversation #1559, we propose adding two small examples to reproduce the acoustic to seismic coupling examples (Figures 4 and 8) from our (recent paper)[https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021JB023142]. This example would complement the case from waterlayered_halfspace by using fluid values for air and also implement a non-planar interface between fluid and elastic media.
Hi All,
@ammcpherson, @carltape, and I are trying to reproduce results from the published paper by Bishop et al. 2022, which uses SPECFEM3D Cartesian to model infrasound waves coupling into the Earth over topography.
We think we may have found a bug in how SPECFEM3D now handles coupling between acoustic and elastic media. This issue may be connected to #1522.
The original paper was written using modeling from the following SPECFEM3D version: commit 8926d3d3354b9fde5fc3b356189b8ddc4c3a30e5 tree 75a3e6c8c57314cb4cb8ec87eb94807cc1f7a7ec parent 19ebf5511b8d2fb6aed269bc2aba203e7346d850 author Daniel Peter daniel.peter22@gmail.com 1568991976 +0300 committer Daniel Peter daniel.peter22@gmail.com 1568991976 +0300
Compiled on an Intel Fortran compiler: 2018.5.274-GCC-5.4.0-2.26
For comparison, we're using the version from Git commit 456b86423d89d8ddd695ef70b590d66c6683c68a From Date: Sat Oct 22 20:12:37 2022 +0200.
For a halfspace case (such as Figure 4 from the paper), we have receivers along a transect. The peak absolute pressure (max(abs(pressure_signal))) for receivers at different burial heights (Z) is shown below for both versions:
In both figures, the CMTSOLUTION location is noted with the red, dotted line and the black, dotted line marks the interface level between the acoustic medium (top) and the elastic medium (below). The pressure should be zero below the dotted line for all receivers. We see this on the left side figure, but the new version from October 2022 shows nonzero pressure values in the elastic medium.
This comparison uses Case 4 in the manuscript, so there is an enormous impedance contrast between acoustic and elastic media.
We see a similar issue when the acoustic-elastic interface has terrain (such as Figure 9 of the above paper; Case 7).
An example pressure waveform in the elastic domain is shown below. Our (external) source time function is an error function, so the resulting pressure waves (in acoustic media) have the form of a Gaussian derivative.
It's not clear to us if this is connected or not, but the location of the receivers in STATIONS_FILTERED also changed between the releases. The original STATIONS are the same, but the new release has an offset "X" location compared the original filtered stations. The Y and Z coordinates are the same. Notably, the difference between the receiver locations is maximal at the center of the simulation domain. For our simulations, the domain decomposition is such that there are 9 slices along xi and 8 slices along eta.
Thank you.