moves mass matrix allocations to create_mass_matrices() routine; updates usage of nglob_xy (for absorbing boundary cases)

Author: danielpeter

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)

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'
 

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)

src/shared/memory_eval.f90

@@ -555,12 +555,12 @@ subroutine memory_eval(NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
   !
   ! 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)

src/shared/save_header_file.F90

@@ -745,18 +745,18 @@ subroutine save_header_file(NSPEC_REGIONS,NGLOB_REGIONS,NPROC,NPROCTOT, &
   !
   ! 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)
   endif
+
   write(IOUT,*) 'integer, parameter :: NGLOB_XY_CM = ',NGLOB_XY_CM
   write(IOUT,*) 'integer, parameter :: NGLOB_XY_IC = ',NGLOB_XY_IC
   write(IOUT,*)

src/specfem3D/compute_forces_viscoelastic_calling_routine.F90

@@ -850,9 +850,7 @@ subroutine compute_forces_inner_core( NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, &
 ! (left in this file to let compiler decide about inlining)
 
   use constants_solver, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,USE_DEVILLE_PRODUCTS_VAL, &
-    ATT5_VAL,N_SLS,NSPEC_INNER_CORE_STRAIN_ONLY,NCHUNKS_VAL
-
-  use shared_parameters, only: ABSORBING_CONDITIONS
+    ATT5_VAL,N_SLS,NSPEC_INNER_CORE_STRAIN_ONLY
 
   ! note: passes sum_terms array as subroutine argument which will help for performance (better than use-statement)
   use specfem_par_innercore, only: sum_terms_inner_core,factor_common_inner_core,NSPEC_INNER_CORE
@@ -891,8 +889,6 @@ subroutine compute_forces_inner_core( NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, &
   real(kind=CUSTOM_REAL), dimension(NGLOB),intent(in) :: pgrav_inner_core
 
   ! checks if anything to
-  ! no need for inner core, absorbing boundaries placed at outer core bottom surface
-  if (NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) return
   ! for regional mesh cut-offs, there are no inner core elements
   if (NSPEC_INNER_CORE == 0) return
 

src/specfem3D/prepare_timerun.F90

@@ -336,6 +336,7 @@ subroutine prepare_timerun_mass_matrices()
      endif
   endif
   rmassz_inner_core = 1._CUSTOM_REAL / rmassz_inner_core
+
   ! outer core
   rmass_outer_core = 1._CUSTOM_REAL / rmass_outer_core