updates storage and usage of rhostore

Author: danielpeter

src/generate_databases/save_arrays_solver.F90

@@ -143,8 +143,11 @@ subroutine save_arrays_solver_mesh()
   write(IOUT) gammazstore
   write(IOUT) jacobianstore
 
+  ! these arrays are needed for acoustic and elastic simulations (and ..)
+  ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic, ..)
   write(IOUT) kappastore
   write(IOUT) mustore
+  write(IOUT) rhostore
 
   write(IOUT) ispec_is_acoustic
   write(IOUT) ispec_is_elastic
@@ -159,10 +162,6 @@ subroutine save_arrays_solver_mesh()
     write(IOUT) rmass_acoustic
   endif
 
-  ! this array is needed for acoustic simulations but also for elastic simulations with CPML,
-  ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-  write(IOUT) rhostore
-
   ! elastic
   if (ELASTIC_SIMULATION) then
     write(IOUT) rmass

src/generate_databases/save_arrays_solver_hdf5.F90

@@ -394,8 +394,11 @@ subroutine save_arrays_solver_mesh_hdf5()
 
     dset_name = "kappastore" ! 4 r (/0,0,0,offset_nspec/)
     call h5_create_dataset_gen(dset_name,(/NGLLX,NGLLY,NGLLZ,sum(offset_nspec(:))/), 4, CUSTOM_REAL)
-    dset_name = "mustore" ! 4 r (/0,0,0,offset_nspec/)
+    dset_name = "mustore"    ! 4 r (/0,0,0,offset_nspec/)
     call h5_create_dataset_gen(dset_name,(/NGLLX,NGLLY,NGLLZ,sum(offset_nspec(:))/), 4, CUSTOM_REAL)
+    dset_name = "rhostore"   ! 4 r (/0,0,0,offset_nspec/)
+    call h5_create_dataset_gen(dset_name,(/NGLLX,NGLLY,NGLLZ,sum(offset_nspec(:))/), 4, CUSTOM_REAL)
+
     dset_name = "ispec_is_acoustic" ! 1 l (/0,0,0,offset_nspec/)
     call h5_create_dataset_gen(dset_name,(/sum(offset_nspec(:))/), 1, 0)
     dset_name = "ispec_is_elastic" ! 1 l (/0,0,0,offset_nspec/)
@@ -409,11 +412,6 @@ subroutine save_arrays_solver_mesh_hdf5()
       call h5_create_dataset_gen(dset_name,(/sum(offset_nglob(:))/), 1, CUSTOM_REAL)
     endif
 
-    ! this array is needed for acoustic simulations but also for elastic simulations with CPML,
-    ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-    dset_name = "rhostore" ! 4 r (/0,0,0,offset_nspec/)
-    call h5_create_dataset_gen(dset_name,(/NGLLX,NGLLY,NGLLZ,sum(offset_nspec(:))/), 4, CUSTOM_REAL)
-
     ! elastic
     if (ELASTIC_SIMULATION) then
       dset_name = "rmass" ! 1 r (/offset_nglob/)
@@ -940,10 +938,14 @@ subroutine save_arrays_solver_mesh_hdf5()
   call h5_write_dataset_collect_hyperslab(dset_name,gammazstore,(/0,0,0,sum(offset_nspec_irregular(0:myrank-1))/),H5_COL)
   dset_name = "jacobianstore" ! 4 r (/0,0,0,offset_nspec_irregular = offset_nspec/)
   call h5_write_dataset_collect_hyperslab(dset_name,jacobianstore,(/0,0,0,sum(offset_nspec_irregular(0:myrank-1))/),H5_COL)
+
   dset_name = "kappastore" ! 4 r (/0,0,0,offset_nspec/)
   call h5_write_dataset_collect_hyperslab(dset_name,kappastore,(/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
   dset_name = "mustore" ! 4 r (/0,0,0,offset_nspec/)
   call h5_write_dataset_collect_hyperslab(dset_name,mustore,(/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
+  dset_name = "rhostore" ! 4 r (/0,0,0,offset_nspec/)
+  call h5_write_dataset_collect_hyperslab(dset_name, rhostore, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
+
   dset_name = "ispec_is_acoustic" ! 1 l (/0,0,0,offset_nspec/)
   call h5_write_dataset_collect_hyperslab(dset_name, ispec_is_acoustic, (/sum(offset_nspec(0:myrank-1))/),H5_COL)
   dset_name = "ispec_is_elastic" ! 1 l (/0,0,0,offset_nspec/)
@@ -957,11 +959,6 @@ subroutine save_arrays_solver_mesh_hdf5()
     call h5_write_dataset_collect_hyperslab(dset_name, rmass_acoustic,(/sum(offset_nglob(0:myrank-1))/),H5_COL)
   endif
 
-  ! this array is needed for acoustic simulations but also for elastic simulations with CPML,
-  ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-  dset_name = "rhostore" ! 4 r (/0,0,0,offset_nspec/)
-  call h5_write_dataset_collect_hyperslab(dset_name, rhostore, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
-
   ! elastic
   if (ELASTIC_SIMULATION) then
     dset_name = "rmass" ! 1 r (/offset_nglob/)

src/generate_databases/setup_color_perm.f90

@@ -631,12 +631,12 @@ subroutine setup_permutation(myrank,nspec,nglob,ibool,ANISOTROPY,perm,SAVE_MESH_
 
   call permute_elements_real(kappastore,temp_array_real,perm,nspec)
   call permute_elements_real(mustore,temp_array_real,perm,nspec)
+  call permute_elements_real(rhostore,temp_array_real,perm,nspec)
 
   if (ATTENUATION) then
      call permute_elements_real(qmu_attenuation_store,temp_array_real,perm,nspec)
      call permute_elements_real(qkappa_attenuation_store,temp_array_real,perm,nspec)
   endif
-  call permute_elements_real(rhostore,temp_array_real,perm,nspec)
 
   ! elastic arrays
   if (ELASTIC_SIMULATION) then

src/shared/define_mass_matrices.f90

@@ -418,9 +418,9 @@ subroutine define_mass_matrices_Stacey_elastic(nglob,nspec,DT,ibool,rho_vp,rho_v
         ny = abs_boundary_normal(2,igll,iface)
         nz = abs_boundary_normal(3,igll,iface)
 
-        vn = deltatover2*(nx+ny+nz)
+        vn = deltatover2 * (nx+ny+nz)
 
-        ! C*deltat/2 contributions
+        ! C * DT/2 contributions
         tx = rho_vp(i,j,k,ispec)*vn*nx + rho_vs(i,j,k,ispec)*(deltatover2-vn*nx)
         ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(deltatover2-vn*ny)
         tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(deltatover2-vn*nz)
@@ -506,7 +506,7 @@ subroutine define_mass_matrices_Stacey_acoustic(nglob,nspec,DT,ibool,rho_vp, &
         k = abs_boundary_ijk(3,igll,iface)
 
         ! C * DT/2 contribution
-        sn = deltatover2/rho_vp(i,j,k,ispec)
+        sn = deltatover2 / rho_vp(i,j,k,ispec)
 
         ! gets associated, weighted jacobian
         jacobianw = abs_boundary_jacobian2Dw(igll,iface)

src/specfem3D/compute_boundary_kernel.f90

@@ -58,8 +58,7 @@ subroutine compute_boundary_kernel()
 
         ! computes contribution from top element
         call compute_boundary_kernel_elem(kernel_moho_top, &
-                                          mustore(i,j,k,ispec_top), &
-                                          kappastore(i,j,k,ispec_top),rho_vs(i,j,k,ispec_top), &
+                                          mustore(i,j,k,ispec_top),kappastore(i,j,k,ispec_top),rhostore(i,j,k,ispec_top), &
                                           accel(:,iglob_top),b_displ(:,iglob_top), &
                                           dsdx_top(:,:,i,j,k,ispec2D),b_dsdx_top(:,:,i,j,k,ispec2D), &
                                           normal_moho_top(:,igll,ispec2D) )
@@ -76,12 +75,11 @@ subroutine compute_boundary_kernel()
           ! iglob_top == iglob_bot!
 
           ! computes contribution from bottom element
-          call compute_boundary_kernel_elem( kernel_moho_bot, &
-                      mustore(i,j,k,ispec_bot), &
-                      kappastore(i,j,k,ispec_bot),rho_vs(i,j,k,ispec_bot), &
-                      accel(:,iglob_bot),b_displ(:,iglob_bot), &
-                      dsdx_bot(:,:,i,j,k,ispec2D),b_dsdx_bot(:,:,i,j,k,ispec2D), &
-                      normal_moho_bot(:,jgll,ispec2D) )
+          call compute_boundary_kernel_elem(kernel_moho_bot, &
+                                            mustore(i,j,k,ispec_bot),kappastore(i,j,k,ispec_bot),rhostore(i,j,k,ispec_bot), &
+                                            accel(:,iglob_bot),b_displ(:,iglob_bot), &
+                                            dsdx_bot(:,:,i,j,k,ispec2D),b_dsdx_bot(:,:,i,j,k,ispec2D), &
+                                            normal_moho_bot(:,jgll,ispec2D) )
 
           ! note: kernel point position: indices given by ijk_moho_top(:,igll,ispec2D)
           moho_kl(igll,ispec2D) = moho_kl(igll,ispec2D) &
@@ -114,12 +112,11 @@ subroutine compute_boundary_kernel()
           iglob_top = ibool(i,j,k,ispec_top)
 
           ! computes contribution from top element
-          call compute_boundary_kernel_elem( kernel_moho_top, &
-                    mustore(i,j,k,ispec_top), &
-                    kappastore(i,j,k,ispec_top),rho_vs(i,j,k,ispec_top), &
-                    accel(:,iglob_top),b_displ(:,iglob_top), &
-                    dsdx_top(:,:,i,j,k,ispec2D),b_dsdx_top(:,:,i,j,k,ispec2D), &
-                    normal_moho_top(:,igll,ispec2D) )
+          call compute_boundary_kernel_elem(kernel_moho_top, &
+                                            mustore(i,j,k,ispec_top),kappastore(i,j,k,ispec_top),rhostore(i,j,k,ispec_top), &
+                                            accel(:,iglob_top),b_displ(:,iglob_top), &
+                                            dsdx_top(:,:,i,j,k,ispec2D),b_dsdx_top(:,:,i,j,k,ispec2D), &
+                                            normal_moho_top(:,igll,ispec2D) )
 
           ! note: kernel point position igll: indices given by ijk_moho_top(:,igll,ispec2D)
           moho_kl(igll,ispec2D) = moho_kl(igll,ispec2D) + kernel_moho_top * deltat
@@ -131,12 +128,11 @@ subroutine compute_boundary_kernel()
           iglob_bot = ibool(i,j,k,ispec_bot)
 
           ! computes contribution from bottom element
-          call compute_boundary_kernel_elem( kernel_moho_bot, &
-                    mustore(i,j,k,ispec_bot), &
-                    kappastore(i,j,k,ispec_bot),rho_vs(i,j,k,ispec_bot), &
-                    accel(:,iglob_bot),b_displ(:,iglob_bot), &
-                    dsdx_bot(:,:,i,j,k,ispec2D),b_dsdx_bot(:,:,i,j,k,ispec2D), &
-                    normal_moho_bot(:,igll,ispec2D) )
+          call compute_boundary_kernel_elem(kernel_moho_bot, &
+                                            mustore(i,j,k,ispec_bot),kappastore(i,j,k,ispec_bot),rhostore(i,j,k,ispec_bot), &
+                                            accel(:,iglob_bot),b_displ(:,iglob_bot), &
+                                            dsdx_bot(:,:,i,j,k,ispec2D),b_dsdx_bot(:,:,i,j,k,ispec2D), &
+                                            normal_moho_bot(:,igll,ispec2D) )
 
           ! note: kernel point position igll: indices given by ijk_moho_bot(:,igll,ispec2D)
           moho_kl(igll,ispec2D) = moho_kl(igll,ispec2D) - kernel_moho_bot * deltat
@@ -153,24 +149,27 @@ end subroutine compute_boundary_kernel
 !-------------------------------------------------------------------------------------------------
 !
 
-subroutine compute_boundary_kernel_elem(kernel, mul, kappal, rho_vsl, &
+subroutine compute_boundary_kernel_elem(kernel, mul, kappal, rhol, &
                                         accel, b_displ, ds, b_ds, norm)
 
 ! compute the boundary kernel contribution from one side of the boundary
 ! see e.g.: Tromp et al. (2005), eq. (25), or Liu & Tromp (2008), eq. (65)
 
-  use constants
+  use constants, only: CUSTOM_REAL,NDIM,ONE
 
   implicit none
 
-  real(kind=CUSTOM_REAL)  kernel, mul, kappal, rho_vsl
-  real(kind=CUSTOM_REAL) :: accel(NDIM), b_displ(NDIM), ds(NDIM,NDIM), b_ds(NDIM,NDIM), norm(NDIM)
+  real(kind=CUSTOM_REAL), intent(inout) ::  kernel
+  real(kind=CUSTOM_REAL), intent(in) ::  mul, kappal, rhol
+  real(kind=CUSTOM_REAL), intent(in) :: accel(NDIM), b_displ(NDIM)
+  real(kind=CUSTOM_REAL), intent(in) :: ds(NDIM,NDIM), b_ds(NDIM,NDIM), norm(NDIM)
 
+  ! local parameters
   real(kind=CUSTOM_REAL) :: eps3, eps(NDIM,NDIM), epsdev(NDIM,NDIM), normal(NDIM,1)
   real(kind=CUSTOM_REAL) :: b_eps3, b_eps(NDIM,NDIM), b_epsdev(NDIM,NDIM)
-  real(kind=CUSTOM_REAL) :: temp1(NDIM,NDIM), rhol, kl(1,1), one_matrix(1,1)
-
+  real(kind=CUSTOM_REAL) :: temp1(NDIM,NDIM), kl(1,1), one_matrix(1,1)
 
+  ! initializes
   normal(:,1) = norm
   one_matrix(1,1) = ONE
 
@@ -194,7 +193,6 @@ subroutine compute_boundary_kernel_elem(kernel, mul, kappal, rho_vsl, &
   epsdev(2,2) = eps(2,2) - eps3 / 3
   epsdev(3,3) = eps(3,3) - eps3 / 3
 
-
   ! backward/reconstructed-forward strain (epsilon) trace
   b_eps3 = b_ds(1,1) + b_ds(2,2) + b_ds(3,3)
 
@@ -218,9 +216,6 @@ subroutine compute_boundary_kernel_elem(kernel, mul, kappal, rho_vsl, &
   ! matrix multiplication
   temp1 = matmul(epsdev,b_epsdev)
 
-  ! density value
-  rhol = rho_vsl ** 2 / mul
-
   ! isotropic kernel value
   ! see e.g.: Tromp et al. (2005), eq. (25), or Liu & Tromp 2008, eq. (65)
   kl = ( rhol * dot_product(accel(:), b_displ(:)) + kappal * eps3 * b_eps3 &

src/specfem3D/finalize_simulation.f90

@@ -142,8 +142,7 @@ subroutine finalize_simulation_cleanup()
   use fault_solver_kinematic, only: SIMULATION_TYPE_KIN,BC_KINFLT_free
 
   !! solving wavefield discontinuity problem with non-split-node scheme
-  use wavefield_discontinuity_solver, only: &
-               finalize_wavefield_discontinuity
+  use wavefield_discontinuity_solver, only: finalize_wavefield_discontinuity
 
   implicit none
 
@@ -233,6 +232,56 @@ subroutine finalize_simulation_cleanup()
   deallocate(kappastore,mustore,rhostore)
   deallocate(ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic)
 
+  if (ELASTIC_SIMULATION) then
+    if (allocated(rho_vp)) deallocate(rho_vp)
+    if (allocated(rho_vs)) deallocate(rho_vs)
+  endif
+
+  ! kernels
+  if (SIMULATION_TYPE == 3) then
+    ! acoustic
+    if (allocated(rho_ac_kl)) deallocate(rho_ac_kl)
+    if (allocated(kappa_ac_kl)) deallocate(kappa_ac_kl)
+    if (allocated(rhop_ac_kl)) deallocate(rhop_ac_kl)
+    if (allocated(alpha_ac_kl)) deallocate(alpha_ac_kl)
+    ! elastic
+    if (allocated(rho_kl)) deallocate(rho_kl)
+    if (allocated(mu_kl)) deallocate(mu_kl)
+    if (allocated(kappa_kl)) deallocate(kappa_kl)
+    if (allocated(cijkl_kl)) deallocate(cijkl_kl)
+    ! poroelastic
+    if (allocated(rhot_kl)) deallocate(rhot_kl)
+    if (allocated(rhof_kl)) deallocate(rhof_kl)
+    if (allocated(sm_kl)) deallocate(sm_kl)
+    if (allocated(eta_kl)) deallocate(eta_kl)
+    if (allocated(mufr_kl)) deallocate(mufr_kl)
+    if (allocated(B_kl)) deallocate(B_kl)
+    if (allocated(C_kl)) deallocate(C_kl)
+    if (allocated(M_kl)) deallocate(M_kl)
+    if (allocated(rhob_kl)) deallocate(rhob_kl)
+    if (allocated(rhofb_kl)) deallocate(rhofb_kl)
+    if (allocated(Bb_kl)) deallocate(Bb_kl)
+    if (allocated(Cb_kl)) deallocate(Cb_kl)
+    if (allocated(Mb_kl)) deallocate(Mb_kl)
+    if (allocated(mufrb_kl)) deallocate(mufrb_kl)
+    if (allocated(phi_kl)) deallocate(phi_kl)
+    if (allocated(rhobb_kl)) deallocate(rhobb_kl)
+    if (allocated(rhofbb_kl)) deallocate(rhofbb_kl)
+    if (allocated(phib_kl)) deallocate(phib_kl)
+    if (allocated(ratio_kl)) deallocate(ratio_kl)
+    if (allocated(cs_kl)) deallocate(cs_kl)
+    if (allocated(cpI_kl)) deallocate(cpI_kl)
+    if (allocated(cpII_kl)) deallocate(cpII_kl)
+    ! hessians
+    if (allocated(hess_ac_kl)) deallocate(hess_ac_kl)
+    if (allocated(hess_rho_ac_kl)) deallocate(hess_rho_ac_kl)
+    if (allocated(hess_kappa_ac_kl)) deallocate(hess_kappa_ac_kl)
+    if (allocated(hess_kl)) deallocate(hess_kl)
+    if (allocated(hess_rho_kl)) deallocate(hess_rho_kl)
+    if (allocated(hess_kappa_kl)) deallocate(hess_kappa_kl)
+    if (allocated(hess_mu_kl)) deallocate(hess_mu_kl)
+  endif
+
   ! faults
   if (SIMULATION_TYPE_DYN) call BC_DYNFLT_free()
   if (SIMULATION_TYPE_KIN) call BC_KINFLT_free()

src/specfem3D/gravity_perturbation.f90

@@ -61,14 +61,13 @@ subroutine gravity_init()
 
   use constants, only: IMAIN,IIN_G,myrank,NGLLX,NGLLY,NGLLZ
 
-  use specfem_par, only: NGLOB_AB, NSTEP, NSPEC_AB, mustore, &
+  use specfem_par, only: NGLOB_AB, NSTEP, NSPEC_AB, &
+       rhostore, &
        xstore, ystore, zstore, &
        xigll, yigll, zigll, &
        wxgll, wygll, wzgll, &
        NGNOD, ibool
 
-  use specfem_par_elastic, only: rho_vs
-
   implicit none
 
   ! local parameters
@@ -174,7 +173,7 @@ subroutine gravity_init()
 
   do ispec = 1,NSPEC_AB
 
-    rho_elem = rho_vs(:,:,:,ispec)*rho_vs(:,:,:,ispec)/mustore(:,:,:,ispec)
+    rho_elem = rhostore(:,:,:,ispec)
 
     do ia = 1,NGNOD
       iglob = ibool(iax(ia),iay(ia),iaz(ia),ispec)

src/specfem3D/read_mesh_databases.F90

@@ -182,6 +182,7 @@ subroutine read_mesh_databases()
 
     read(IIN) kappastore
     read(IIN) mustore
+    read(IIN) rhostore
 
     read(IIN) ispec_is_acoustic
     read(IIN) ispec_is_elastic
@@ -220,6 +221,7 @@ subroutine read_mesh_databases()
 
   call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore,kind=4))
   call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore,kind=4))
+  call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4))
 
   call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic,kind=4))
   call bcast_all_l_for_database(ispec_is_elastic(1), size(ispec_is_elastic,kind=4))
@@ -289,11 +291,6 @@ subroutine read_mesh_databases()
     rmassz_acoustic(:) = 0.0_CUSTOM_REAL
   endif
 
-! rho array is needed for acoustic simulations but also for elastic simulations with CPML,
-! read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-  if (I_should_read_the_database) read(IIN) rhostore
-  call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4))
-
   ! elastic simulation
   if (ELASTIC_SIMULATION) then
     ! checks

src/specfem3D/read_mesh_databases_hdf5.F90

@@ -265,6 +265,8 @@ subroutine read_mesh_databases_hdf5()
     call h5_read_dataset_collect_hyperslab(dsetname,kappastore,(/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
     dsetname = "mustore"
     call h5_read_dataset_collect_hyperslab(dsetname,mustore,(/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
+    dsetname = "rhostore"
+    call h5_read_dataset_collect_hyperslab(dsetname, rhostore, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
   endif
 
 ! note: the size(..) function returns either integer(kind=4) or integer(kind=8)
@@ -295,6 +297,7 @@ subroutine read_mesh_databases_hdf5()
 
   call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore,kind=4))
   call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore,kind=4))
+  call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4))
 
   if (size(ispec_is_acoustic) > 0) &
     call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic,kind=4))
@@ -369,15 +372,6 @@ subroutine read_mesh_databases_hdf5()
     rmassz_acoustic(:) = 0.0_CUSTOM_REAL
   endif
 
-! rho array is needed for acoustic simulations but also for elastic simulations with CPML,
-! read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-  if (I_should_read_the_database) then
-    dsetname = "rhostore"
-    !call h5_read_dataset_p(dsetname, rhostore)
-    call h5_read_dataset_collect_hyperslab(dsetname, rhostore, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL)
-  endif
-  call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4))
-
   ! elastic simulation
   if (ELASTIC_SIMULATION) then
     ! checks