Skip to content

adds compatibility flag (for older version comparisons); updates adjacency array setup #848

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Oct 7, 2024
Merged

adds compatibility flag (for older version comparisons); updates adjacency array setup #848

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fa91265
adds flag USE_OLD_VERSION_FORMAT to constants.h.in for compatibility …
danielpeter Oct 7, 2024
6e76675
updates usage of flag USE_1D_REFERENCE from constants.h to affect all…
danielpeter Oct 7, 2024
5a25fa1
updates adjacency array setup for point location searches (using node…
danielpeter Oct 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 12 additions & 1 deletion setup/constants.h.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -564,6 +564,17 @@
! be careful when changing this flag when computing classical kernel
integer, parameter :: NTSTEP_BETWEEN_COMPUTE_KERNELS = 1

!!-----------------------------------------------------------
!!
!! new version compatibility
!!
!!-----------------------------------------------------------
! for comparison against older versions
! in version 7.0: tiso was constrained to below moho, new version allows also in crust
! in version 8.0: conversions between geodetic & geocentric latitudes were always applied,
! and additional differences in 1-chunk centering & source positioning
logical, parameter :: USE_OLD_VERSION_FORMAT = .false.


!!-----------------------------------------------------------
!!
Expand Down Expand Up @@ -1024,7 +1035,7 @@
! parameters for GLL model (used for iterative model inversions)
character(len=*), parameter :: PATHNAME_GLL_modeldir = 'DATA/GLL/'

! to create a reference model based on 1D_REF but with 3D crust and 410/660 topography
! to create a reference model based on the 1D reference model but with 3D crust and 410/660 topography
logical,parameter :: USE_1D_REFERENCE = .false.

!!-- uncomment for using PREM as reference model (used in CEM inversion)
Expand Down
74 changes: 52 additions & 22 deletions src/meshfem3D/compute_element_properties.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -480,7 +480,8 @@
use constants, only: IMAIN,myrank, &
IFLAG_CRUST,IFLAG_220_80,IFLAG_80_MOHO,IFLAG_670_220,IFLAG_MANTLE_NORMAL,IREGION_CRUST_MANTLE, &
REFERENCE_MODEL_1DREF,REFERENCE_MODEL_1DREF, &
THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE
THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE, &
USE_OLD_VERSION_FORMAT

use meshfem_models_par, only: &
TRANSVERSE_ISOTROPY,USE_FULL_TISO_MANTLE,REFERENCE_1D_MODEL,THREE_D_MODEL
Expand Down Expand Up @@ -511,6 +512,8 @@
write(IMAIN,*) ' setting tiso flags in mantle model'
if (USE_FULL_TISO_MANTLE) &
write(IMAIN,*) ' using fully transverse isotopic mantle'
if (USE_OLD_VERSION_FORMAT) &
write(IMAIN,*) ' using formatting from old versions (7.0 to 8.0)'
call flush_IMAIN()
endif
endif
Expand All @@ -521,12 +524,16 @@
if (USE_FULL_TISO_MANTLE) then
! all elements below the actual moho will be used for transverse isotropy
! note: this will increase the computation time by ~ 45 %
if (idoubling(ispec) == IFLAG_MANTLE_NORMAL &
.or. idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
if (elem_in_mantle) elem_is_tiso = .true.
else
if (idoubling(ispec) == IFLAG_MANTLE_NORMAL &
.or. idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.

Check warning on line 535 in src/meshfem3D/compute_element_properties.f90

View check run for this annotation

Codecov / codecov/patch

src/meshfem3D/compute_element_properties.f90#L534-L535 

Added lines #L534 - L535 were not covered by tests
endif
endif

! all done
Expand All @@ -545,26 +552,49 @@
! THREE_D_MODEL_S29EA
! THREE_D_MODEL_GLL
! which show significant transverse isotropy also below 220km depth
! assigns TI to elements in crust and mantle down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
! assigns TI to elements in mantle elements just below actual moho down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. (idoubling(ispec) == IFLAG_CRUST &
.and. elem_in_mantle) ) then
elem_is_tiso = .true.
endif
else
! assigns TI to elements in crust and mantle down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
endif
endif

case default
! default reference models
! for example, PREM assigns transverse isotropy between Moho and 220km
! assigns TI to elements in crust and mantle elements (down to 220),
! to allow for tiso in crust and below actual moho (especially for oceanic crusts);
! the crustal models will decide if model parameters are tiso or iso
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
! assigns TI to elements in mantle elements just below actual moho
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. (idoubling(ispec) == IFLAG_CRUST &
.and. elem_in_mantle )) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
endif
else
! assigns TI to elements in crust and mantle elements (down to 220),
! to allow for tiso in crust and below actual moho (especially for oceanic crusts);
! the crustal models will decide if model parameters are tiso or iso
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
endif
endif
end select

Expand Down
94 changes: 66 additions & 28 deletions src/meshfem3D/meshfem3D_models.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -703,7 +703,10 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
!
! note: in general, models here make use of perturbation values with respect to their
! corresponding 1-D reference models
if (MODEL_3D_MANTLE_PERTUBATIONS .and. r_prem > RCMB/R_PLANET .and. .not. suppress_mantle_extension) then
if (MODEL_3D_MANTLE_PERTUBATIONS &
.and. r_prem > RCMB/R_PLANET &
.and. .not. suppress_mantle_extension &
.and. .not. USE_1D_REFERENCE) then

! extend 3-D mantle model above the Moho to the surface before adding the crust
if (r_prem > RCMB/R_PLANET .and. r_prem < RMOHO/R_PLANET) then
Expand Down Expand Up @@ -749,6 +752,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case(THREE_D_MODEL_MANTLE_SH)
! full_sh model
! gets lat/lon coordinates from geocentric theta/phi
lat = (PI/2.0d0-theta)*180.0d0/PI
lon = phi*180.0d0/PI
if (lon > 180.0d0) lon = lon - 360.0d0
Expand Down Expand Up @@ -794,20 +798,12 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
case (THREE_D_MODEL_S362ANI,THREE_D_MODEL_S362WMANI, &
THREE_D_MODEL_S362ANI_PREM,THREE_D_MODEL_S29EA)
! 3D Harvard models s362ani, s362wmani, s362ani_prem and s2.9ea
! gets colat/lon coordinates from geocentric theta/phi
xcolat = sngl(theta*180.0d0/PI)
xlon = sngl(phi*180.0d0/PI)
xrad = sngl(r_used*R_PLANET_KM)
call model_s362ani_subshsv(xcolat,xlon,xrad,xdvsh,xdvsv,xdvph,xdvpv)

! to use speed values from the 1D reference model but with 3D mesh variations
if (USE_1D_REFERENCE) then
! sets all 3D variations in the mantle to zero
xdvpv = 0.d0
xdvph = 0.d0
xdvsv = 0.d0
xdvsh = 0.d0
endif

if (TRANSVERSE_ISOTROPY) then
! tiso perturbation
vpv = vpv*(1.0d0+dble(xdvpv))
Expand Down Expand Up @@ -852,7 +848,6 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case (THREE_D_MODEL_SGLOBE,THREE_D_MODEL_SGLOBE_ISO)
! 3D SGLOBE-rani model (Chang)

! normally mantle perturbations are taken from 24.4km (R_MOHO) up.
! we need to add the if statement for sgloberani_iso or sgloberani_aniso to take from 50km up:
if (r_prem > RCMB/R_PLANET .and. r_prem < 6321000.d0/R_PLANET) then
Expand All @@ -861,7 +856,6 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
! this will then "extend the mantle up to the surface" from 50km depth
r_used = 6321000.d0/R_PLANET
endif

call mantle_sglobe(r_used,theta,phi,vsv,vsh,dvsv,dvsh,dvp,drho)

! updates velocities from reference model
Expand Down Expand Up @@ -902,6 +896,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case (THREE_D_MODEL_SPIRAL)
! SPiRaL v1.4 anisotropic model
! gets lat/lon coordinates from geocentric theta/phi
lat = (PI/2.0d0-theta)*180.0d0/PI
lon = phi*180.0d0/PI
! input: lat in range [-90,90]; lon in range [-180,180]
Expand Down Expand Up @@ -938,7 +933,8 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
! (based on density variations) on top of reference 3D model
if (HETEROGEN_3D_MANTLE &
.and. .not. suppress_mantle_extension &
.and. THREE_D_MODEL /= THREE_D_MODEL_HETEROGEN_PREM) then
.and. THREE_D_MODEL /= THREE_D_MODEL_HETEROGEN_PREM &
.and. .not. USE_1D_REFERENCE) then
! gets spherical coordinates of actual point location
r_used = r
! adds hetergeneous perturbations (isotropic)
Expand Down Expand Up @@ -1199,7 +1195,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &

! returns velocities and density for points in 3D crustal region

use shared_parameters, only: ELLIPTICITY
use meshfem_models_par

implicit none
Expand All @@ -1218,7 +1213,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &
double precision,intent(inout) :: moho,sediment

! local parameters
double precision :: lat,lon
double precision :: vpvc,vphc,vsvc,vshc,etac
double precision :: vpc,vsc,rhoc !vpc_eu
double precision :: c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
Expand All @@ -1231,13 +1225,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &
! for point radius smaller than deepest possible crust radius (~80 km depth)
if (r < R_DEEPEST_CRUST) return

! converts geocentric colatitude/lon (theta/phi) to geographic latitude/lon (lat/lon)
! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
call thetaphi_2_geographic_latlon_dble(theta,phi,lat,lon,ELLIPTICITY)

! double-check lon in range [-180,180]
if (lon > 180.0d0 ) lon = lon - 360.0d0

!---
!
! ADD YOUR MODEL HERE
Expand All @@ -1260,14 +1247,14 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &

if (.not. is_inside_region) then
! uses default crust outside of model region
call meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
call meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)
endif

case default
! default crust
call meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
call meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)

Expand Down Expand Up @@ -1358,20 +1345,21 @@ end subroutine meshfem3D_models_get3Dcrust_val
!


subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
subroutine meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc, &
moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)

! returns velocity/density for default crust

use constants, only: USE_OLD_VERSION_FORMAT
use meshfem_models_par

implicit none

! lat/lon - in degrees (range lat/lon = [-90,90] / [-180,180]
! radius r - normalized by globe radius [0,1.x]
double precision,intent(in) :: lat,lon,r
! theta/phi - geocentric coordinates
double precision,intent(in) :: r,theta,phi

double precision,intent(out) :: vpvc,vphc,vsvc,vshc,etac,rhoc
double precision,intent(out) :: moho,sediment
Expand All @@ -1380,6 +1368,7 @@ subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
double precision,intent(out) :: c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c
! local parameters
double precision :: lat,lon
! for isotropic crust
double precision :: vpc,vsc
double precision :: vpc_area,vsc_area,rhoc_area,moho_area,sediment_area
Expand Down Expand Up @@ -1420,7 +1409,56 @@ subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
! ADD YOUR MODEL HERE
!
!---
! lat/lon range: [-90,90] / [-180,180]

! note: Here we assume that the crustal models are defined with respect to geographic lat/lon positions.
! This is different to mantle models, where we assume that those are referenced to a spherical Earth.
!
! Geocentric and geographic colatitude (theta) are only identical for a spherical Earth.
! Depending on the Par_file 'ELLIPTICITY' flag however, we will stretch out the spherical mesh after assigining
! the model velocities.
!
! Given we want to assign the crustal velocities at point (x/y/z geocentric) which then corresponds to the
! coordinates (lat/lon geographic) in the final, stretched out mesh, we need to calculate the (lat/lon) position
! from (theta/phi) with the ellipticity correction.
!
! again, in other words (see comment in moho_stretching.f90):
! "
! the mesh here by default starts as a spherical mesh. the position (x/y/z) is thus a geocentric position
! and we want to assign the velocity & depth of the resulting geographic (lat/lon) position after stretching.
!
! we would have two options to achieve this:
! (a) we convert first all the crustal model (lat/lon) positions to the geocentric (lat'/lon') ones
! and then search for the corresponding geocentric (lat'/lon') position here.
! (b) we calculate the geographic (lat/lon) for this geocentric point location (x/y/z) and then
! search for the geographic (lat/lon) position in the original crustal model positions.
! the implementation here follows option (b) to leave the crustal model in its original form.
!
! thus, having the geocentric (x/y/z) position and its geocentric (lat'/lon') coordinate, we account here for the
! ellipticity stretching afterwards to find the proper resulting geographic (lat/lon) coordinate
! for which the crustal model is defined.
! in case the mesh stays spherical (no ELLIPTICITY), the geocentric (lat'/lon') coordinates are identical
! to geographic (lat/lon) coordinates and no correction is applied.
!
! this assumes that the crustal models are given for geographic (lat/lon) positions and
! not corrected geocentric (lat'/lon') positions.
! "
!
! TODO: we might want to re-visit this assumption.

if (USE_OLD_VERSION_FORMAT) then
! lat/lon from geocentric position w/out geodetic correction
lat = (PI_OVER_TWO - theta) * RADIANS_TO_DEGREES
lon = phi * RADIANS_TO_DEGREES
else
! converts geocentric colatitude/lon (theta/phi) to geographic latitude/lon (lat/lon)
! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
call thetaphi_2_geographic_latlon_dble(theta,phi,lat,lon,ELLIPTICITY)
endif

! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
! double-check lon in range [-180,180]
if (lon > 180.d0) lon = lon - 360.0d0
if (lon < -180.d0) lon = lon + 360.0d0

select case (REFERENCE_CRUSTAL_MODEL)

Expand Down
Loading