Merge pull request #857 from danielpeter/devel

updates rmass allocations; updates checks and Berkeley model testing
SPECFEM · Dec 2, 2024 · a6ca022 · a6ca022
2 parents 19dfaf2 + 4208388
commit a6ca022
Show file tree

Hide file tree

Showing 12 changed files with 152 additions and 173 deletions.
diff --git a/README.md b/README.md
@@ -62,6 +62,11 @@ involved in the community and keep them in the specfem3d github wiki:
 [specfem3d wiki](https://github.com/SPECFEM/specfem3d/wiki)
 
 
+## Our contributors :sparkles:
+
+[![SPECFEM3D_GLOBE contributors](https://contrib.rocks/image?repo=SPECFEM/specfem3d_globe)](https://github.com/SPECFEM/specfem3d_globe/graphs/contributors)
+
+
 ## Computational Infrastructure for Geodynamics (CIG)
 
 SPECFEM3D_GLOBE is part of the software that is hosted by the Computational Infrastructure for Geodynamics (CIG). It is available on the CIG website [here (SPECFEM3D_GLOBE)](https://geodynamics.org/resources/specfem3d).

diff --git a/src/meshfem3D/create_mass_matrices.f90 b/src/meshfem3D/create_mass_matrices.f90
@@ -56,8 +56,8 @@ subroutine create_mass_matrices(idoubling,ibool, &
   use regions_mesh_par2, only: &
     xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore, &
     gammaxstore,gammaystore,gammazstore,rhostore,kappavstore, &
-    rmassx,rmassy,rmassz,b_rmassx,b_rmassy, &
-    nglob_xy
+    rmassx,rmassy,rmassz,b_rmassx,b_rmassy,rmass_ocean_load, &
+    nglob_xy,nglob_oceans
 
   implicit none
 
@@ -75,10 +75,61 @@ subroutine create_mass_matrices(idoubling,ibool, &
   ! local parameters
   double precision :: weight
   real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+  integer :: ispec,i,j,k,iglob,ier
 
-  integer :: ispec,i,j,k,iglob
+  ! initializes
+  nglob_xy = 0      ! for rmassx/rmassy
+  nglob_oceans = 0  ! for ocean load mass matrix
+
+  ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+  ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+  if (NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
+    ! crust/mantle region uses different rmassx/rmassy/rmassz for absorbing boundary
+    if (iregion_code == IREGION_CRUST_MANTLE) nglob_xy = nglob
+  endif
+
+  if (ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION) then
+    ! rotation in solid domains uses different rmassx/rmassy/rmassz
+    if (iregion_code == IREGION_CRUST_MANTLE .or. &
+        iregion_code == IREGION_INNER_CORE) nglob_xy = nglob
+  endif
+
+  ! ocean load mass matrix
+  if (OCEANS) then
+    ! ocean load only applies in crust/mantle region
+    if (iregion_code == IREGION_CRUST_MANTLE) nglob_oceans = nglob
+  endif
+
+  ! allocates mass matrices in this slice (will be fully assembled in the solver)
+  allocate(rmassz(nglob),stat=ier)
+  if (ier /= 0) stop 'Error in allocate 22'
+  rmassz(:) = 0.0_CUSTOM_REAL
+
+  allocate(rmassx(nglob_xy), &
+           rmassy(nglob_xy), &
+           stat=ier)
+  if (ier /= 0) stop 'Error in allocate 21'
+  rmassx(:) = 0.0_CUSTOM_REAL
+  rmassy(:) = 0.0_CUSTOM_REAL
+
+  allocate(b_rmassx(nglob_xy), &
+           b_rmassy(nglob_xy),stat=ier)
+  if (ier /= 0) stop 'Error in allocate b_21'
+  b_rmassx(:) = 0.0_CUSTOM_REAL
+  b_rmassy(:) = 0.0_CUSTOM_REAL
+
+  ! allocates ocean load mass matrix as well if oceans
+  allocate(rmass_ocean_load(nglob_oceans),stat=ier)
+  if (ier /= 0) stop 'Error in allocate 22'
+  rmass_ocean_load(:) = 0.0_CUSTOM_REAL
 
-  ! initializes matrices
+  ! checks if anything to do
+  if (iregion_code == IREGION_TRINFINITE .or. iregion_code == IREGION_INFINITE) return
+
+  ! check if anything to do (must have region domain points)
+  if (nglob == 0) return
+
+  ! setup matrices
   !
   ! in the case of Stacey boundary conditions, add C*delta/2 contribution to the mass matrix
   ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
@@ -89,18 +140,7 @@ subroutine create_mass_matrices(idoubling,ibool, &
   !
   ! Now also handle EXACT_MASS_MATRIX_FOR_ROTATION, which requires similar corrections
 
-  rmassx(:) = 0._CUSTOM_REAL
-  rmassy(:) = 0._CUSTOM_REAL
-  rmassz(:) = 0._CUSTOM_REAL
-
-  b_rmassx(:) = 0._CUSTOM_REAL
-  b_rmassy(:) = 0._CUSTOM_REAL
-
-  ! checks if anything to do
-  if (iregion_code == IREGION_TRINFINITE .or. iregion_code == IREGION_INFINITE) return
-
   ! first create the main standard mass matrix with no corrections
-
 ! openmp mesher
 !$OMP PARALLEL DEFAULT(SHARED) &
 !$OMP PRIVATE(ispec,i,j,k,iglob,weight, &
@@ -427,6 +467,14 @@ subroutine create_mass_matrices_Stacey(ibool,iregion_code)
   ! adds contributions to mass matrix to stabilize Stacey conditions
   select case (iregion_code)
   case (IREGION_CRUST_MANTLE)
+
+    ! compatibility with old versions (6.x, 7.x)
+    if (USE_OLD_VERSION_FORMAT) then
+      ! note: this will override the rotation contributions added in case before this routine call
+      rmassx(:) = rmassz(:)
+      rmassy(:) = rmassz(:)
+    endif
+
     do iface = 1,num_abs_boundary_faces
 
       ispec = abs_boundary_ispec(iface)

diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
@@ -55,22 +55,19 @@ subroutine create_regions_mesh(npointot, &
     HDF5_ENABLED
 
   use meshfem_par, only: &
-    myrank,nspec,nglob,iregion_code, &
+    myrank,nspec,iregion_code, &
     ibool,idoubling,xstore,ystore,zstore, &
     xstore_glob,ystore_glob,zstore_glob
 
   use meshfem_par, only: &
     NCHUNKS,SAVE_MESH_FILES,ABSORBING_CONDITIONS,LOCAL_PATH, &
     ADIOS_FOR_ARRAYS_SOLVER,ADIOS_FOR_SOLVER_MESHFILES, &
-    ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,GRAVITY_INTEGRALS, &
+    GRAVITY_INTEGRALS, &
     FULL_GRAVITY, &
     NGLOB1D_RADIAL_CORNER, &
     NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
     volume_total,Earth_mass_total,Earth_center_of_mass_x_total,Earth_center_of_mass_y_total,Earth_center_of_mass_z_total
 
-  use meshfem_models_par, only: &
-    OCEANS
-
 #ifdef USE_CEM
   use meshfem_models_par, only: CEM_REQUEST
 #endif
@@ -96,8 +93,8 @@ subroutine create_regions_mesh(npointot, &
     mu0store,Gc_prime_store,Gs_prime_store, &
     rho_vp,rho_vs, &
     Qmu_store,tau_e_store, &
-    nglob_xy,rmassx,rmassy,rmassz,b_rmassx,b_rmassy, &
-    nglob_oceans,rmass_ocean_load, &
+    rmassx,rmassy,rmassz,b_rmassx,b_rmassy, &
+    rmass_ocean_load, &
     iMPIcut_xi,iMPIcut_eta, &
     ispec_is_tiso
 
@@ -129,9 +126,6 @@ subroutine create_regions_mesh(npointot, &
   ! now perform two passes in this part to be able to save memory
   integer,intent(in) :: ipass
 
-  ! local parameters
-  integer :: ier
-
   ! user output
   if (myrank == 0) then
     write(IMAIN,*)
@@ -374,71 +368,6 @@ subroutine create_regions_mesh(npointot, &
       write(IMAIN,*) '  ...creating mass matrix'
       call flush_IMAIN()
     endif
-
-    ! allocates mass matrices in this slice (will be fully assembled in the solver)
-    !
-    ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
-    ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
-    ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
-    !
-    ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
-    ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
-    if (NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
-      select case (iregion_code)
-      case (IREGION_CRUST_MANTLE)
-        nglob_xy = nglob
-      case (IREGION_INNER_CORE, IREGION_OUTER_CORE)
-        nglob_xy = 1
-      case (IREGION_TRINFINITE, IREGION_INFINITE)
-        nglob_xy = 1
-      case default
-        call exit_mpi(myrank,'Invalid region code for nglob_xy')
-      end select
-    else
-       nglob_xy = 1
-    endif
-
-    if (ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION) then
-      select case (iregion_code)
-      case (IREGION_CRUST_MANTLE,IREGION_INNER_CORE)
-         nglob_xy = nglob
-      case (IREGION_OUTER_CORE)
-         nglob_xy = 1
-      case (IREGION_TRINFINITE, IREGION_INFINITE)
-        nglob_xy = 1
-      case default
-        call exit_mpi(myrank,'Invalid region code for nglob_xy with EXACT_MASS_MATRIX_FOR_ROTATION')
-      end select
-    endif
-
-    allocate(rmassx(nglob_xy), &
-             rmassy(nglob_xy), &
-             stat=ier)
-    if (ier /= 0) stop 'Error in allocate 21'
-    rmassx(:) = 0.0_CUSTOM_REAL
-    rmassy(:) = 0.0_CUSTOM_REAL
-
-    allocate(b_rmassx(nglob_xy), &
-             b_rmassy(nglob_xy),stat=ier)
-    if (ier /= 0) stop 'Error in allocate b_21'
-    b_rmassx(:) = 0.0_CUSTOM_REAL
-    b_rmassy(:) = 0.0_CUSTOM_REAL
-
-    allocate(rmassz(nglob),stat=ier)
-    if (ier /= 0) stop 'Error in allocate 22'
-    rmassz(:) = 0.0_CUSTOM_REAL
-
-    ! allocates ocean load mass matrix as well if oceans
-    if (OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) then
-      nglob_oceans = nglob
-    else
-      ! allocate dummy array if no oceans
-      nglob_oceans = 1
-    endif
-    allocate(rmass_ocean_load(nglob_oceans),stat=ier)
-    if (ier /= 0) stop 'Error in allocate 22'
-    rmass_ocean_load(:) = 0.0_CUSTOM_REAL
-
     ! creating mass matrices in this slice (will be fully assembled in the solver)
     ! note: for Stacey boundaries, needs indexing nimin,.. filled in the first pass
     call create_mass_matrices(idoubling,ibool, &
@@ -522,25 +451,6 @@ subroutine create_regions_mesh(npointot, &
 
     endif ! .not. GRAVITY_INTEGRALS
 
-    ! synchronizes processes before clean up memory
-    call synchronize_all()
-
-    ! frees memory
-    deallocate(rmassx,rmassy,rmassz)
-    deallocate(b_rmassx,b_rmassy)
-    deallocate(rmass_ocean_load)
-    ! frees allocated mesh memory
-    deallocate(xstore_glob,ystore_glob,zstore_glob)
-    ! frees MPI arrays memory
-    call crm_free_MPI_arrays()
-    ! free Stacey helper arrays (not needed anymore)
-    if (allocated(nimin)) deallocate(nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta)
-    if (allocated(abs_boundary_ispec)) then
-      deallocate(abs_boundary_ispec,abs_boundary_npoin,abs_boundary_ijk)
-      deallocate(abs_boundary_jacobian2Dw)
-      deallocate(abs_boundary_normal)
-    endif
-
     ! compute volume, bottom and top area of that part of the slice, and then the total
     call compute_volumes_and_areas(NCHUNKS,iregion_code,nspec,wxgll,wygll,wzgll, &
                                    xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
@@ -569,6 +479,25 @@ subroutine create_regions_mesh(npointot, &
 
     endif
 
+    ! synchronizes processes before clean up memory
+    call synchronize_all()
+
+    ! frees memory
+    deallocate(rmassx,rmassy,rmassz)
+    deallocate(b_rmassx,b_rmassy)
+    deallocate(rmass_ocean_load)
+    ! frees allocated mesh memory
+    deallocate(xstore_glob,ystore_glob,zstore_glob)
+    ! frees MPI arrays memory
+    call crm_free_MPI_arrays()
+    ! free Stacey helper arrays (not needed anymore)
+    if (allocated(nimin)) deallocate(nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta)
+    if (allocated(abs_boundary_ispec)) then
+      deallocate(abs_boundary_ispec,abs_boundary_npoin,abs_boundary_ijk)
+      deallocate(abs_boundary_jacobian2Dw)
+      deallocate(abs_boundary_normal)
+    endif
+
   case default
     stop 'there cannot be more than two passes in mesh creation'
 

diff --git a/src/meshfem3D/model_crust_berkeley.f90 b/src/meshfem3D/model_crust_berkeley.f90
@@ -148,7 +148,7 @@ subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust,elem_in_c
 
   depth = (1.d0 - x) * EARTH_R_KM
 
-  call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth)
+  call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth,moho_only)
 
   ! using crustal values
   if (USE_OLD_VERSION_FORMAT) then
@@ -226,7 +226,7 @@ subroutine model_berkeley_crust_aniso(x,theta,phi,vpv,vph,vsv,vsh,eta_aniso,rho,
 
   depth = (1.d0 - x) * EARTH_R_KM
 
-  call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth)
+  call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth,moho_only)
 
   ! using crustal values
   if (USE_OLD_VERSION_FORMAT) then
@@ -266,14 +266,15 @@ end subroutine model_berkeley_crust_aniso
 !--------------------------------------------------------------------------------------------------
 !
 
-  subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth)
+  subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth,moho_only)
 
   use model_crust_berkeley_par
 
   implicit none
 
   double precision,intent(in) :: theta,phi,z
   double precision,intent(out) :: rho,vp,vsv,vsh,moho_depth
+  logical,intent(in) :: moho_only
 
   ! local parameters
   ! 4-th order GLL positions
@@ -290,13 +291,22 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth)
   double precision,external :: moho_filtre
   double precision,external :: lagrange
 
+  ! initialize
+  rho = 0.d0
+  vp = 0.d0
+  vsv = 0.d0
+  vsh = 0.d0
+
   ! get moho depth
   moho_depth = moho1D_depth - moho_filtre(theta,phi)
 
   !debug
   !print *,"debug: [get_crust_val_csem] Moho depth:",moho_depth, &
   !        "moho1D_depth:",moho1D_depth,"moho_filtre:",moho_filtre(theta,phi)
 
+  ! check if anything further to do or return only moho
+  if (moho_only) return
+
   !
   ! horizontal interpolation for all registered depths
   !

diff --git a/src/meshfem3D/save_arrays_solver.f90 b/src/meshfem3D/save_arrays_solver.f90
@@ -773,13 +773,13 @@ subroutine save_arrays_mesh_parameters()
   !
   ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
   ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be fictitious / unused
+  NGLOB_XY_CM = 0
+  NGLOB_XY_IC = 0
 
   if (NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
-     NGLOB_XY_CM = NGLOB_REGIONS(IREGION_CRUST_MANTLE)
-  else
-     NGLOB_XY_CM = 0
+    NGLOB_XY_CM = NGLOB_REGIONS(IREGION_CRUST_MANTLE)
   endif
-  NGLOB_XY_IC = 0
+
   if (ROTATION .and. EXACT_MASS_MATRIX_FOR_ROTATION) then
     NGLOB_XY_CM = NGLOB_REGIONS(IREGION_CRUST_MANTLE)
     NGLOB_XY_IC = NGLOB_REGIONS(IREGION_INNER_CORE)