Merge pull request #851 from danielpeter/devel

adds version backward compatibility for Berkeley model; adds speedup modifications to Berkeley model meshing

EXAMPLES/regional_Berkeley/DATA/Par_file

@@ -81,7 +81,7 @@ FULL_GRAVITY                    = .false.
 POISSON_SOLVER                  = 0
 
 # record length in minutes
-RECORD_LENGTH_IN_MINUTES        = 2.5d0
+RECORD_LENGTH_IN_MINUTES        = 5.0d0
 
 #-----------------------------------------------------------
 #
@@ -269,6 +269,17 @@ USE_MONOCHROMATIC_CMT_SOURCE    = .false.
 # print source time function
 PRINT_SOURCE_TIME_FUNCTION      = .true.
 
+## Berkeley source time function
+STF_IS_UCB_HEAVISIDE            = .true.
+# UCB source frequency content (i.e., heaviside function)
+SOURCE_T1                       = 500.d0
+SOURCE_T2                       = 400.d0
+SOURCE_T3                       =  50.d0
+SOURCE_T4                       =  40.d0
+# UCB source time shift
+TAU                             = 500.d0
+
+
 #-----------------------------------------------------------
 #
 # Seismograms
@@ -347,7 +358,7 @@ APPROXIMATE_HESS_KL             = .false.
 # forces transverse isotropy for all mantle elements
 # (default is to use transverse isotropy only between crust and 220)
 # means we allow radial anisotropy throughout the whole crust/mantle region
-USE_FULL_TISO_MANTLE            = .false.
+USE_FULL_TISO_MANTLE            = .true.
 
 # output kernel mask to zero out source region
 # to remove large values near the sources in the sensitivity kernels

setup/constants.h.in

@@ -245,12 +245,6 @@
 !!
 !!-----------------------------------------------------------
 
-! source-time function is a UCB-style filtered heaviside
-  logical, parameter :: STF_IS_UCB_HEAVISIDE = .false.
-
-! Save mirror files
-  logical, parameter :: SAVE_MIRRORS = .false.
-
 ! Path to A3d model
 ! Make sure this line ends on "/"
   character (len=*), parameter :: A3d_folder = "DATA/SEMUCB_A3d/"
@@ -266,13 +260,19 @@
 ! pathname of the topography file (un-smoothed)
   character (len=*), parameter :: PATHNAME_TOPO_FILE_BERKELEY = trim(A3d_folder)//"ETOPO5_1x1_filtre.dat"
 
-!! uncomment this to use Berkeley topography as default - and comment out corresponding parameters above !!
+!! uncomment this only to use Berkeley topography as default for all other Earth models as well
+!! note: Berkeley topography setting is used by default for SEMUCB model (see in get_model_parameters.F90 file).
+!!       uncommenting the lines below here (- and commenting out corresponding parameters in the topo/bathy section above)
+!!       will force the Berkeley topo to be used as the default topography for all Earth models.
 !! Topography defaults to Berkeley
 !  integer, parameter :: EARTH_NX_BATHY = NX_BATHY_BERKELEY
 !  integer, parameter :: EARTH_NY_BATHY = NY_BATHY_BERKELEY
 !  double precision, parameter :: EARTH_RESOLUTION_TOPO_FILE = RESOLUTION_TOPO_FILE_BERKELEY
 !  character (len=*), parameter :: EARTH_PATHNAME_TOPO_FILE = PATHNAME_TOPO_FILE_BERKELEY
 
+! Save mirror files - not implemented yet
+!  logical, parameter :: SAVE_MIRRORS = .false.
+
 
 !!-----------------------------------------------------------
 !!

src/gpu/compute_stacey_acoustic_gpu.c

@@ -76,7 +76,7 @@ void FC_FUNC_ (compute_stacey_acoustic_gpu,
     size_t local_work_size[2];
     cl_uint idx = 0;
 
-    //daniel debug
+    // debug
     //clCheck (clFinish (mocl.command_queue));
     //printf ("rank %d: stacey a %i, %i save %i num blocks x/y= %i %i nglob %i nspec2D %i nspec %i\n",
     //          mp->myrank,interface_type,num_abs_boundary_faces,mp->save_forward,num_blocks_x,num_blocks_y,
@@ -104,7 +104,7 @@ void FC_FUNC_ (compute_stacey_acoustic_gpu,
     clCheck (clEnqueueNDRangeKernel (mocl.command_queue, mocl.kernels.compute_stacey_acoustic_kernel, 2, NULL,
                                      global_work_size, local_work_size, 0, NULL, NULL));
 
-    //daniel debug
+    // debug
     //clCheck (clFinish (mocl.command_queue));
     //printf ("rank %d: stacey b %i, %i \n", mp->myrank,interface_type,num_abs_boundary_faces);
     //fflush (stdout);

src/meshfem3D/add_topography.f90

@@ -31,6 +31,10 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo)
   use shared_parameters, only: REGIONAL_MESH_CUTOFF,REGIONAL_MESH_CUTOFF_DEPTH,USE_LOCAL_MESH,ELLIPTICITY
   use meshfem_par, only: R220,NX_BATHY,NY_BATHY,R_PLANET
 
+  ! for old version Berkeley compatibility
+  use constants, only: USE_OLD_VERSION_FORMAT,ICRUST_BERKELEY,THREE_D_MODEL_BERKELEY
+  use meshfem_models_par, only: THREE_D_MODEL,REFERENCE_CRUSTAL_MODEL
+
   implicit none
 
   double precision,dimension(NGNOD), intent(inout) :: xelm,yelm,zelm
@@ -45,6 +49,13 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo)
 
   integer :: ia
 
+  ! for compatibility
+  double precision :: theta,phi
+  double precision :: vpvc,vphc,vsvc,vshc,etac,rhoc
+  double precision :: moho
+  double precision :: rmoho
+  logical :: found_crust,elem_in_crust,moho_only
+
   ! we loop on all the points of the element
   do ia = 1,NGNOD
 
@@ -76,6 +87,29 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo)
     endif
     gamma = (r - rbottom) / (R_UNIT_SPHERE - rbottom)
 
+    ! old version compatility
+    if (USE_OLD_VERSION_FORMAT) then
+      ! Berkeley model
+      if (THREE_D_MODEL == THREE_D_MODEL_BERKELEY .and. &
+          REFERENCE_CRUSTAL_MODEL == ICRUST_BERKELEY) then
+        ! convert lat/lon to theta/phi
+        call latlon_2_thetaphi_dble(lat,lon,theta,phi)
+
+        ! gets smoothed moho depth
+        elem_in_crust = .true.
+        moho_only = .true.
+        call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust,elem_in_crust,moho_only)
+
+        rmoho = R_UNIT_SPHERE - moho
+
+        ! if point above moho then move points, otherwise skip
+        if (r <= rmoho) cycle
+
+        ! adjust gamma stretching to moho boundary distance
+        gamma = (r - rmoho) / (R_UNIT_SPHERE - rmoho)
+      endif
+    endif
+
     ! add elevation to all the points of that element
     ! also make sure gamma makes sense
     if (gamma < -0.02 .or. gamma > 1.02) then
@@ -107,7 +141,6 @@ subroutine add_topography_gll(xstore,ystore,zstore,ispec,nspec,ibathy_topo)
   use shared_parameters, only: REGIONAL_MESH_CUTOFF,REGIONAL_MESH_CUTOFF_DEPTH,USE_LOCAL_MESH,ELLIPTICITY
   use meshfem_par, only: R220,NX_BATHY,NY_BATHY
 
-
   implicit none
 
   ! input parameters

src/meshfem3D/compute_element_properties.f90

@@ -188,7 +188,7 @@ subroutine compute_element_properties(ispec,iregion_code,idoubling,ipass, &
   endif ! IREGION_CRUST_MANTLE
 
   ! sets element tiso flag
-  call compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,ispec,nspec,idoubling)
+  call compute_element_tiso_flag(elem_is_tiso,elem_in_crust,elem_in_mantle,iregion_code,ispec,nspec,idoubling)
 
   ! stores as element flags
   ispec_is_tiso(ispec) = elem_is_tiso
@@ -473,14 +473,14 @@ end subroutine compute_element_GLL_locations
 !
 
 
-  subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,ispec,nspec,idoubling)
+  subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_crust,elem_in_mantle,iregion_code,ispec,nspec,idoubling)
 
 ! sets transverse isotropic flag for elements in crust/mantle
 
   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,REFERENCE_MODEL_SEMUCB, &
-    THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE, &
+    THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE,THREE_D_MODEL_BERKELEY, &
     USE_OLD_VERSION_FORMAT
 
   use meshfem_models_par, only: &
@@ -489,7 +489,7 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
   implicit none
 
   logical,intent(out) :: elem_is_tiso
-  logical,intent(in) :: elem_in_mantle
+  logical,intent(in) :: elem_in_crust,elem_in_mantle
   integer,intent(in) :: iregion_code,ispec
   integer,intent(in) :: nspec
   integer,dimension(nspec),intent(in) :: idoubling
@@ -501,12 +501,11 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
   elem_is_tiso = .false.
 
   ! checks if anything to do
-  if (.not. TRANSVERSE_ISOTROPY) return
   if (iregion_code /= IREGION_CRUST_MANTLE) return
+  if (.not. TRANSVERSE_ISOTROPY) return
 
   ! user output
   if (is_first_call) then
-    is_first_call = .false.
     if (myrank == 0) then
       ! only output once
       write(IMAIN,*) '  setting tiso flags in mantle model'
@@ -516,6 +515,8 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
         write(IMAIN,*) '    using formatting from old versions (7.0 to 8.0)'
       call flush_IMAIN()
     endif
+    ! turns off flag to show output only once
+    is_first_call = .false.
   endif
 
   ! transverse isotropic models
@@ -526,6 +527,10 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
     ! note: this will increase the computation time by ~ 45 %
     if (USE_OLD_VERSION_FORMAT) then
       if (elem_in_mantle) elem_is_tiso = .true.
+      ! adds special case for Berkeley SEMUCB model
+      if (THREE_D_MODEL == THREE_D_MODEL_BERKELEY) then
+        if (elem_in_crust) elem_is_tiso = .true.
+      endif
     else
       if (idoubling(ispec) == IFLAG_MANTLE_NORMAL &
           .or. idoubling(ispec) == IFLAG_670_220 &
@@ -557,8 +562,7 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
       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
+          .or. (idoubling(ispec) == IFLAG_CRUST .and. elem_in_mantle)) then
         elem_is_tiso = .true.
       endif
     else
@@ -572,14 +576,30 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
     endif
 
   case (REFERENCE_MODEL_SEMUCB)
-    ! SEMUCB - allows tiso for full mantle & crust
-    ! (same as USE_FULL_TISO_MANTLE option)
-    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.
+    ! SEMUCB
+    if (USE_OLD_VERSION_FORMAT) then
+      ! default from version 6.0
+      if (idoubling(ispec) == IFLAG_220_80 &
+          .or. idoubling(ispec) == IFLAG_80_MOHO) then
+        ! models use only transverse isotropy between moho and 220 km depth
+        elem_is_tiso = .true.
+      endif
+      ! or additional setting for SEMUCB was to use USE_FULL_TISO_MANTLE
+      !
+      ! note: this sets tiso for elements fully in mantle or fully in crust
+      !       however, this will still have tiso == .false. for elements with the moho inside going through the element,
+      !       having corner nodes both in mantle and crust (above and below actual moho)
+      !if (elem_in_crust) elem_is_tiso = .true.
+      !if (elem_in_mantle) elem_is_tiso = .true.
+    else
+      ! allows tiso for full mantle & crust
+      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.
+      endif
     endif
 
   case default
@@ -639,6 +659,11 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
     ! note: THREE_D_MODEL_SGLOBE_ISO
     !       sgloberani_iso model based on PREM, it will have tiso already set from crust down to 220
 
+  case (THREE_D_MODEL_BERKELEY)
+    ! SEMUCB
+    ! double-check to use tiso for elements fully in crust
+    if (elem_in_crust) elem_is_tiso = .true.
+
   case default
     ! nothing special to add
     continue

src/meshfem3D/create_mass_matrices.f90

@@ -543,19 +543,24 @@ subroutine create_mass_matrices_ocean_load(ibool,xstore,ystore,zstore,NSPEC2D_TO
   ! checks if anything to do
   if (.not. OCEANS) return
 
-  ! user output
-  if (myrank == 0) then
-    write(IMAIN,*) '     updates mass matrix with ocean load'
-    call flush_IMAIN()
-  endif
-
   ! initializes
   do_ocean_load = .false.
 
   ! note: new version:
   !       for 3D Earth with topography, compute local height of oceans
   if (TOPOGRAPHY) do_ocean_load = .true.
 
+  ! user output
+  if (myrank == 0) then
+    write(IMAIN,*) '     updates mass matrix with ocean load'
+    if (do_ocean_load) then
+      write(IMAIN,*) '       using minimum ocean thickness: ',MINIMUM_THICKNESS_3D_OCEANS,'(m)'
+    else
+      write(IMAIN,*) '       using 1D ocean PREM thickness: ',THICKNESS_OCEANS_PREM * EARTH_R,'(m)'
+    endif
+    call flush_IMAIN()
+  endif
+
   ! create ocean load mass matrix for degrees of freedom at ocean bottom
   rmass_ocean_load(:) = 0._CUSTOM_REAL
 

src/meshfem3D/meshfem3D_models.F90

@@ -924,7 +924,8 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
         case (THREE_D_MODEL_BERKELEY)
           ! 3D Berkeley Model SEMUCB
           ! note: passes r/theta/phi (geocentric coordinates)
-          call model_berkeley_shsv(r_used,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL)
+          call model_berkeley_shsv(r_used,theta,phi,vpv,vph,vsv,vsh,rho, &
+                                   dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL)
 
           ! updates velocities from reference model
           if (TRANSVERSE_ISOTROPY) then
@@ -1625,9 +1626,9 @@ subroutine meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc, &
       ! Berkeley crustal model
       ! note: passes r/theta/phi (geocentric coordinates)
       !       Berkeley crustal model is referencing geocentric positions in a spherical Earth frame
-      call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust)
+      call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust,elem_in_crust,moho_only)
       ! old version - isotropic crustal velocities
-      !call model_berkeley_crust(r,theta,phi,vpc,vsc,rhoc,moho,found_crust)
+      !call model_berkeley_crust(r,theta,phi,vpc,vsc,rhoc,moho,found_crust,elem_in_crust,moho_only)
       !vpvc = vpc
       !vphc = vpc
       !vsvc = vsc