From 68a8f50de342bcde8e5bb8be3fd47dac2a45f05a Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Fri, 29 Nov 2024 15:26:15 +0100 Subject: [PATCH] moves mass matrix allocations to create_mass_matrices() routine; updates usage of nglob_xy (for absorbing boundary cases) --- src/meshfem3D/create_mass_matrices.f90 | 78 +++++++++--- src/meshfem3D/create_regions_mesh.F90 | 117 ++++-------------- src/meshfem3D/save_arrays_solver.f90 | 8 +- src/shared/memory_eval.f90 | 8 +- src/shared/save_header_file.F90 | 8 +- ...te_forces_viscoelastic_calling_routine.F90 | 6 +- src/specfem3D/prepare_timerun.F90 | 1 + 7 files changed, 100 insertions(+), 126 deletions(-) diff --git a/src/meshfem3D/create_mass_matrices.f90 b/src/meshfem3D/create_mass_matrices.f90 index 4604f7b99..d3675c849 100644 --- 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 index ac94b8bff..d537ca921 100644 --- 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/save_arrays_solver.f90 b/src/meshfem3D/save_arrays_solver.f90 index a311db5f5..0379b62e8 100644 --- 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) diff --git a/src/shared/memory_eval.f90 b/src/shared/memory_eval.f90 index 4854585e7..c4001629c 100644 --- a/src/shared/memory_eval.f90 +++ b/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) diff --git a/src/shared/save_header_file.F90 b/src/shared/save_header_file.F90 index b7c341d0e..8768156b0 100644 --- a/src/shared/save_header_file.F90 +++ b/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,*) diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index 85a13f822..4632d2294 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/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 diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90 index 78adbdc34..47a736d7a 100644 --- a/src/specfem3D/prepare_timerun.F90 +++ b/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