Pressure Seismogram #138
Pressure Seismogram #138
Conversation
|
Can one of the admins verify this patch? |
2 similar comments
|
Can one of the admins verify this patch? |
|
Can one of the admins verify this patch? |
src/writer/seismogram.cpp
Outdated
| case specfem::enums::seismogram::type::pressure: | ||
| filename = { this->output_folder + "/" + network_name + station_name + | ||
| "BX" + ".semp" }; // TODO get the right name | ||
| break; |
There was a problem hiding this comment.
Relevant to this checklist.
There was a problem hiding this comment.
The naming scheme is based off SPECFEM2D seismogram naming scheme. https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/write_seismograms.F90#L302
| //bulk modulus | ||
| const type_real kappa = [&]() -> type_real { | ||
| if constexpr (MediumTag == specfem::element::medium_tag::acoustic){ | ||
| return point_properties.kappa; | ||
| }else if constexpr (MediumTag == specfem::element::medium_tag::elastic){ | ||
| //we may also need an isotropic if statement? | ||
| return point_properties.lambdaplus2mu; | ||
| }else{ | ||
| static_assert(false, "include/domain/impl/receivers/kernel.tpp: unknown medium to retrieve kappa!"); | ||
| return 0; | ||
| } | ||
| }(); |
There was a problem hiding this comment.
These may or may not be the right properties for the bulk modulus. Please verify.
There was a problem hiding this comment.
|
I think your formulation is different then what has been previously used in SPECFEM2D. For elastic domains - https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/compute_pressure.f90#L398 For acoustic domains - https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/compute_pressure.f90#L569 I think we should use the same formulations. |
| @@ -84,6 +84,10 @@ void specfem::domain::impl::kernels::receiver_kernel< | |||
| NGLL, DimensionType, MediumTag, specfem::kokkos::DevScratchSpace, | |||
| Kokkos::MemoryTraits<Kokkos::Unmanaged>, true, true, true, false, using_simd>; | |||
|
|
|||
| using Aux2ComponentFieldType = specfem::element::field< | |||
There was a problem hiding this comment.
I think what you might want to use is using Aux2ComponentFieldType = typename ElementFieldType::ViewType. This should directly give you access to the underlying Kokkos::View used to store the data
There was a problem hiding this comment.
Can we remove Aux2ComponentFieldType if its not getting used?
| #pragma unroll | ||
| #endif | ||
| for (int l = 0; l < specfem::dimension::dimension<DimensionType>::dim; l++) { | ||
| aux_field.displacement(iz,ix,l) = sv_receiver_field(l); |
There was a problem hiding this comment.
With direct access to Kokkos::View as stated in the comment above. This could be aux_field(iz, ix, l)
| @@ -186,10 +194,110 @@ void specfem::domain::impl::kernels::receiver_kernel< | |||
| point_properties, | |||
There was a problem hiding this comment.
I think it would be ideal if we can pass seismogram type to this function. I realize for pressure we need the seismo type within get_field. With your current formulation this seems not possible, since we need the divergence. Can you look at the SPECFEM2D formulation https://github.com/SPECFEM/specfem2d/blob/master/src/specfem2D/compute_pressure.f90 .
With that formulation we should be able to compute pressure using only gradients
src/writer/seismogram.cpp
Outdated
| case specfem::enums::seismogram::type::pressure: | ||
| filename = { this->output_folder + "/" + network_name + station_name + | ||
| "BX" + ".semp" }; // TODO get the right name | ||
| break; |
There was a problem hiding this comment.
The naming scheme is based off SPECFEM2D seismogram naming scheme. https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/write_seismograms.F90#L302
|
I've added seismogram-typing to As for unit testing, I was thinking about just having a set of (SPECFEM++ parameter-file, SPECFEM2D seismogram-output-files) pairs to compare against. The current seismogram tests don't have a ground truth for pressure. |
Yes I think thats the best way of doing it right now. Where we compare seismograms from the 2 versions of the software. We already do that within displacement tests https://github.com/PrincetonUniversity/SPECFEMPP/blob/78365e29dca1ad784f336bdefdbd5b370585fb69/tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp#L283C1-L328C1. In that case I compare displacement seismograms, but you could easily do the same by adding pressure seismograms in |
| } | ||
|
|
||
| // -kappa * (dsxdx + dszdz) | ||
| receiver_field(0) = -properties.lambdaplus2mu * ( |
There was a problem hiding this comment.
Seems like your computation for pressure field is different from the original formulation. In the original formulation pressure = trace of stress tensor (https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/compute_pressure.f90#L398)
|
|
||
| // -kappa * (dsxdx + dszdz) | ||
| receiver_field(0) = -properties.lambdaplus2mu * ( | ||
| (dsx_dxi * partial_derivatives.xix + |
There was a problem hiding this comment.
Also seems like I havent given access to kappa within the point properties type. lambdaplus2mu != kappa (
).| // tells us that fhe first index is the derivative index | ||
| #define du(i,j) (i==0? (j==0? dsx_dx:dsz_dx):(j==0? dsx_dz:dsz_dz)) | ||
|
|
||
| // WET code... should this be refactored with |
There was a problem hiding this comment.
I agree. At later point, we should just call the function from impl/elements. I think I need to do some refactoring before that. For now this implementation looks good.
|
|
||
| //https://specfem2d-kokkos.readthedocs.io/en/devel/api/datatypes/field_derivatives/point.html#_CPPv4N7specfem5point17field_derivatives2duE | ||
| // tells us that fhe first index is the derivative index | ||
| #define du(i,j) (i==0? (j==0? dsx_dx:dsz_dx):(j==0? dsx_dz:dsz_dz)) |
There was a problem hiding this comment.
This LGTM, we should eventually define this as field_Derivatives. I think some refactor is needed before then. For now this is fine.
| #undef du | ||
| // https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/compute_pressure.f90#L398 | ||
| // pressure = - trace(sigma) / 3 | ||
| receiver_field(0) = - (sigma_xx + sigma_zz)/3.0; |
There was a problem hiding this comment.
You missed sigma_yy which in 2D is lambdaplus2mu * (du(0,0) + du(1,1))
There was a problem hiding this comment.
Thank you for catching that! I was planning on looking for a lost factor of 0.75 after the meeting. This might be it.
int-ptr-ptr
force-pushed
the
pressure_seismo
branch
from
23e66bd
to
15b6437
Compare
|
The It seems like serial test 3 has the wrong displacement, but serial tests 1 and 2 are good for all 4 seismogram types. It does not seem like the provenance parfile generates the same simulation. Serial tests 5-7 have acceleration relative errors ~ 5e-2, and there are some stations that give all zeros for pressure Fortran-side, so rel-errors are failing out as |
It might not be the best to test all seismograms in all simulations then. Especially if the seismogram values are very small. I think that might be the issue with acceleration. I am not sure what issue you are seeing with serial test 3, unless the reference Fortran seismograms have changed. The implementation in this PR hasnt modified how displacement seismograms should be computed. If you selectively wanted to add pressure seimogram tests, I would recommend adding it to serial test 1 and 2 since it'll test both elastic and acoustic implementation individually. |
|
Also I am used git commit ID |
|
To add/remove seismgrams to the tests, it suffices to modify the parameter Is the specfem parameter documentation the only part of the docs that has seismograms? I added a table here. Is that sufficient, or is there somewhere else that needs changing? |
Description
This PR attempts to add a pressure seismogram option.
Issue Number
#137
Checklist
Developers' List
Before certain checks in the maintainers' list, the PR should be feature complete (ex.
CMakeListsdoes not need an update at first, but a new feature in the PR creates a new.tppfile). These should be verified maintainer-side:f"{network}{station}BX.semp"as a placeholder. (link to comment)Maintainers' List
Please make sure to check developer documentation on specfem docs.