diff --git a/src/shared/SIEM_math_library.F90 b/src/shared/SIEM_math_library.F90 index 07949e470..f1ef6df0a 100644 --- a/src/shared/SIEM_math_library.F90 +++ b/src/shared/SIEM_math_library.F90 @@ -385,6 +385,197 @@ function i_nearless (XVAL) result (I_nl) end function i_nearless + !======================================================= + + subroutine dspherical_unitvect(theta,phi,unitr,unittheta,unitphi) + + ! this function computes unit vectors in spherical coordinates at a point + ! defined by (r,theta,phi) + + implicit none + double precision,intent(in) :: theta,phi + double precision,intent(out) :: unitr(3),unittheta(3),unitphi(3) + double precision :: ctheta,cphi,stheta,sphi + double precision,parameter :: zero = 0.d0 + + ctheta = cos(theta); stheta = sin(theta) + cphi = cos(phi); sphi = sin(phi) + + unitr(1) = stheta*cphi; unitr(2) = stheta*sphi; unitr(3) = ctheta + unittheta(1) = ctheta*cphi; unittheta(2) = ctheta*sphi; unittheta(3) = -stheta + unitphi(1) = -sphi; unitphi(2) = cphi; unitphi(3) = zero + + end subroutine dspherical_unitvect + !======================================================= + + subroutine spherical_unitvect(theta,phi,unitr,unittheta,unitphi) + + ! this function computes unit vectors in spherical coordinates at a point + ! defined by (r,theta,phi) + + implicit none + real(kind=CUSTOM_REAL),intent(in) :: theta,phi + real(kind=CUSTOM_REAL),intent(out) :: unitr(3),unittheta(3),unitphi(3) + real(kind=CUSTOM_REAL) :: ctheta,cphi,stheta,sphi + real(kind=CUSTOM_REAL),parameter :: zero=0.0_CUSTOM_REAL + + ctheta = cos(theta); stheta = sin(theta) + cphi = cos(phi); sphi = sin(phi) + + unitr(1) = stheta*cphi; unitr(2) = stheta*sphi; unitr(3) = ctheta + unittheta(1) = ctheta*cphi; unittheta(2) = ctheta*sphi; unittheta(3) = -stheta + unitphi(1) = -sphi; unitphi(2) = cphi; unitphi(3) = zero + + end subroutine spherical_unitvect + + !======================================================= + + subroutine dlegendreP2_costheta(theta,P2,dP2,d2P2) + + implicit none + double precision,intent(in) :: theta ! radian + double precision,intent(out) :: P2,dP2,d2P2 + double precision :: ctheta + double precision,parameter :: half = 0.5d0,one = 1.d0,two = 2.d0,three = 3.d0 + + ctheta = cos(theta) + + P2 = half*(three*ctheta*ctheta-one) + dP2 = -three*sin(theta)*ctheta + d2P2 = -3*cos(two*theta) + + end subroutine dlegendreP2_costheta + + !======================================================= + + subroutine legendreP2_costheta(theta,P2,dP2,d2P2) + + implicit none + real(kind=CUSTOM_REAL),intent(in) :: theta ! radian + real(kind=CUSTOM_REAL),intent(out) :: P2,dP2,d2P2 + real(kind=CUSTOM_REAL) :: ctheta + real(kind=CUSTOM_REAL),parameter :: half=0.5_CUSTOM_REAL,one=1.0_CUSTOM_REAL,two=2.0_CUSTOM_REAL,three=3.0_CUSTOM_REAL + + ctheta = cos(theta) + + P2 = half*(three*ctheta*ctheta-one) + dP2 = -three*sin(theta)*ctheta + d2P2 = -3*cos(two*theta) + + end subroutine legendreP2_costheta + + !======================================================= + + subroutine compute_g_gradg_elliptical(ndim,r,theta,phi,rho,dotrho,g0,eps,eta,twothirdOmega2,g,gradg) + + implicit none + integer,intent(in) :: ndim + double precision,intent(in) :: r,theta,phi,rho,dotrho,g0,eps,eta,twothirdOmega2 + double precision,intent(out) :: g(ndim) + double precision,optional,intent(out) :: gradg(6) + + double precision,parameter :: two = 2.d0,four = 4.d0,two_third = two/3.d0 + + double precision :: hmat(ndim,ndim) + double precision :: lfac,rinv + double precision :: cotthetaXdP2,doteps,ddoteps,doteta,dotg + double precision :: facP2,P2,dP2,d2P2 + double precision :: unitr(ndim,1),unittheta(ndim,1),unitphi(ndim,1) + double precision :: unitrT(1,ndim),unitthetaT(1,ndim),unitphiT(1,ndim) + + ! compute unit vectors + call dspherical_unitvect(theta,phi,unitr,unittheta,unitphi) + + unitrT(1,:) = unitr(:,1); + unitthetaT(1,:) = unittheta(:,1) + unitphiT(1,:) = unitphi(:,1) + + call dlegendreP2_costheta(theta,P2,dP2,d2P2) + + rinv = 1.d0/r + + dotg = four*rho-two*rinv*g0 !dotg = four_pi_G*rho-two*rinv*g_ss + doteps = eta*eps*rinv !eta*eps/r + ddoteps = 6.d0*rinv*rinv*eps-8.0d0*rho*(doteps+rinv*eps)/g0 !two*four + + !doteta=doteps/eps-0.5_CUSTOM_REAL*r0*rinv*rinv*doteps*doteps+r0*ddoteps/eps + !doteta=doteps/eps-r*doteps*doteps/(eps*eps)+r*ddoteps/eps + doteta = doteps/eps+r*(eps*ddoteps-doteps*doteps)/(eps*eps) + + !cottheta=cos(theta)/sin(theta) + cotthetaXdp2 = -3.d0*cos(theta)*cos(theta) + + ! lfac=g_ss+two_third*(eta*eps*g_ss+four_pi_G*r0*rho*eps-eps*g_ss)*P2-two_third*Omega2*r0 + lfac = g0+two_third*(eta*eps*g0+four*r*rho*eps-eps*g0)*P2-twothirdOmega2*r + + ! compute g + g = -lfac*unitr(:,1)-two_third*eps*g0*dP2*unittheta(:,1) + + if (.not.present(gradg)) return + + ! compute grad g + facP2 = (doteta*eps*g0+eta*doteps*g0+eta*eps*dotg+four*(rho*eps+r*dotrho*eps+r*rho*doteps)-doteps*g0-eps*dotg) + + !hmat=-(dotg+two_third*facP2*P2-two_third*Omega2)*matmul(unitr,unitrT) & + ! -two_third*(doteps*g0+eps+dotg)*dP2*(matmul(unitr,unitthetaT)+matmul(unittheta,unitrT)) & + ! -rinv*(lfac-two_third*Omega2*r+two_third*rinv*eps*g0*d2p2)*matmul(unittheta,unitthetaT) & + ! -rinv*(lfac-two_third*Omega2*r+two_third*rinv*eps*g0*cotthetaXdp2)*matmul(unitphi,unitphiT) + hmat = -(dotg+two_third*facP2*P2-twothirdOmega2)*matmul(unitr,unitrT) & + -two_third*(doteps*g0+eps+dotg)*dP2*(matmul(unitr,unitthetaT)+matmul(unittheta,unitrT)) & + -rinv*(lfac+two_third*rinv*eps*g0*d2p2)*matmul(unittheta,unitthetaT) & + -rinv*(lfac+two_third*rinv*eps*g0*cotthetaXdp2)*matmul(unitphi,unitphiT) + + gradg = (/ hmat(1,1),hmat(2,2),hmat(3,3),hmat(1,2),hmat(1,3),hmat(2,3) /) + + end subroutine compute_g_gradg_elliptical + + !======================================================= + + subroutine compute_g_gradg(ndim,r,theta,phi,rho,g0,g,gradg) + + implicit none + integer,intent(in) :: ndim + double precision,intent(in) :: r,theta,phi,rho,g0 + double precision,intent(out) :: g(ndim) + double precision,optional,intent(out) :: gradg(6) + + double precision,parameter :: two = 2.d0,four = 4.d0 + double precision :: hmat(ndim,ndim) + double precision :: lfac,rinv + double precision :: dotg + double precision :: ctheta,P2 !,dP2,d2P2 + double precision :: unitr(ndim,1),unittheta(ndim,1),unitphi(ndim,1) + double precision :: unitrT(1,ndim),unitthetaT(1,ndim),unitphiT(1,ndim) + + ! compute unit vectors + call dspherical_unitvect(theta,phi,unitr,unittheta,unitphi) + + unitrT(1,:) = unitr(:,1); + unitthetaT(1,:) = unittheta(:,1) + unitphiT(1,:) = unitphi(:,1) + + !call dlegendreP2_costheta(theta,P2,dP2,d2P2) + + ctheta = cos(theta) + P2 = 0.5d0*(3.0d0*ctheta*ctheta-1.0d0) + + rinv = 1.d0/r + dotg = four*rho-two*rinv*g0 !dotg = four_pi_G*rho-two*rinv*g_ss + + lfac = g0 + + ! compute g + g = -lfac*unitr(:,1) + + if (.not.present(gradg)) return + + ! compute grad g + hmat = -dotg*matmul(unitr,unitrT)-rinv*lfac*matmul(unittheta,unitthetaT) & + -rinv*lfac*matmul(unitphi,unitphiT) + + gradg = (/ hmat(1,1),hmat(2,2),hmat(3,3),hmat(1,2),hmat(1,3),hmat(2,3) /) + + end subroutine compute_g_gradg end module siem_math_library diff --git a/src/specfem3D/SIEM_compute_kernels.F90 b/src/specfem3D/SIEM_compute_kernels.F90 index 2de77ad95..b23475edc 100644 --- a/src/specfem3D/SIEM_compute_kernels.F90 +++ b/src/specfem3D/SIEM_compute_kernels.F90 @@ -70,13 +70,14 @@ subroutine calculate_first_gravity_kernel() use specfem_par_full_gravity, only: gdof_cm1, inode_elmt_cm, inode_map_cm, nmir_cm, nnode_cm1, & gknl1, & - gravload1, pgrav1, pgrav_cm1, pgrav_cm, & - is_active_gll,igll_active_on + gravload1, pgrav1, pgrav_cm1, pgrav_cm, neq1, ndscale1, dprecon1, & + is_active_gll,igll_active_on, & + CG_SCALING use siem_math_library_mpi, only: maxvec use siem_poisson, only: compute_grav_kl1_load use siem_solver_petsc, only: petsc_set_vector1, petsc_zero_initialguess1, petsc_solve1 - use siem_solver_mpi, only: interpolate3to5 + use siem_solver_mpi, only: cg_solver3, diagpcg_solver3, interpolate3to5 implicit none @@ -87,13 +88,6 @@ subroutine calculate_first_gravity_kernel() real(kind=CUSTOM_REAL) :: maxload real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: gknl1_cm ! (3,NGLOB_CRUST_MANTLE) - ! safety check - if (POISSON_SOLVER == ISOLVER_BUILTIN) then - !TODO: full gravity builtin solver for kernels - print *,'ERROR: builtin solver not setup for gravity kernels' - call exit_MPI(myrank,'Error builtin solver not setup for gravity kernels') - endif - ! Allocate 1st grav kernel array : allocate(gknl1(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier) if (ier /= 0) stop 'Error allocating gknl1' @@ -117,7 +111,17 @@ subroutine calculate_first_gravity_kernel() maxload = maxvec(abs(gravload1)) if (myrank == 0) print *,' -- Max load: ', maxload - if (POISSON_SOLVER == ISOLVER_PETSC) then + if (POISSON_SOLVER == ISOLVER_BUILTIN) then + ! builtin solver + pgrav1(:) = 0.0_CUSTOM_REAL + if (CG_SCALING) then + gravload1(:) = ndscale1(:) * gravload1(:) + call cg_solver3(myrank,neq1,pgrav1,gravload1) + pgrav1(:) = ndscale1(:) * pgrav1(:) + else + call diagpcg_solver3(myrank,neq1,pgrav1,gravload1,dprecon1) + endif + else ! Petsc solver call petsc_set_vector1(gravload1) @@ -126,14 +130,14 @@ subroutine calculate_first_gravity_kernel() call petsc_zero_initialguess1() pgrav1(:) = 0.0_CUSTOM_REAL - ! Here we use pgrav as the vector the solution is put into, just to + ! Here we use pgrav1 as the vector the solution is put into, just to ! save on allocating another array. Note that we could allocate a - ! separate temporary array but pgrav isnt used after this (apart from + ! separate temporary array but pgrav1 isnt used after this (apart from ! if the forward wavefield is being saved, but it shouldnt be for SIMTYPE 3) call petsc_solve1(pgrav1(1:)) endif - ! Now interpolate this component to the 5GLL setup (Crust mantle only): + ! Now interpolate this component to the 5GLL setup (Crust mantle only) pgrav_cm1(:) = zero pgrav_cm1(:) = pgrav1(gdof_cm1(:)) @@ -207,13 +211,14 @@ subroutine calculate_second_gravity_kernel() use specfem_par_full_gravity, only: gdof_cm1, inode_elmt_cm, inode_map_cm, nmir_cm, nnode_cm1, & gknl2, & - gravload1, pgrav1, pgrav_cm1, pgrav_cm, & - is_active_gll,igll_active_on + gravload1, pgrav1, pgrav_cm1, pgrav_cm, neq1, ndscale1, dprecon1, & + is_active_gll,igll_active_on, & + CG_SCALING use siem_math_library_mpi, only: maxvec use siem_poisson, only: compute_grav_kl2_load use siem_solver_petsc, only: petsc_set_vector1, petsc_zero_initialguess1, petsc_solve1 - use siem_solver_mpi, only: interpolate3to5 + use siem_solver_mpi, only: cg_solver3, diagpcg_solver3, interpolate3to5 implicit none @@ -225,13 +230,6 @@ subroutine calculate_second_gravity_kernel() real(kind=CUSTOM_REAL),dimension(:,:,:), allocatable :: gknl2_cm real(kind=CUSTOM_REAL),dimension(:,:,:,:,:), allocatable :: div_gknl2_cm - ! safety check - if (POISSON_SOLVER == ISOLVER_BUILTIN) then - !TODO: full gravity builtin solver for kernels - print *,'ERROR: builtin solver not setup for gravity kernels' - call exit_MPI(myrank,'Error builtin solver not setup for gravity kernels') - endif - ! Allocate 2nd gravity kernel array allocate(gknl2(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier) if (ier /= 0) stop 'Error allocating gknl2' @@ -257,7 +255,17 @@ subroutine calculate_second_gravity_kernel() maxload = maxvec(abs(gravload1)) if (myrank == 0) print *,' -- Max load: ', maxload - if (POISSON_SOLVER == ISOLVER_PETSC) then + if (POISSON_SOLVER == ISOLVER_BUILTIN) then + ! builtin solver + pgrav1(:) = 0.0_CUSTOM_REAL + if (CG_SCALING) then + gravload1(:) = ndscale1(:) * gravload1(:) + call cg_solver3(myrank,neq1,pgrav1,gravload1) + pgrav1(:) = ndscale1(:) * pgrav1(:) + else + call diagpcg_solver3(myrank,neq1,pgrav1,gravload1,dprecon1) + endif + else ! PETSc solver call petsc_set_vector1(gravload1) @@ -267,10 +275,12 @@ subroutine calculate_second_gravity_kernel() call petsc_solve1(pgrav1(1:)) endif - pgrav_cm(:) = zero + ! Now interpolate this component to the 5GLL setup (Crust mantle only) pgrav_cm1(:) = zero pgrav_cm1(:) = pgrav1(gdof_cm1(:)) + pgrav_cm(:) = zero !initialise before interpolating to be safe + call interpolate3to5(NSPEC_CRUST_MANTLE, NGLOB_CRUST_MANTLE, nnode_cm1, & inode_elmt_cm, nmir_cm, inode_map_cm, is_active_gll, igll_active_on, pgrav_cm1, pgrav_cm) @@ -537,6 +547,7 @@ subroutine SIEM_compute_crust_mantle_kernels() ! Local values for the GLL point sdagtmp(:) = displ_crust_mantle(:,iglob) stmp(:) = b_displ_crust_mantle(:,iglob) + hmatloc(:) = gradg_cm(:, iglob) ! (3) Gravity contraction s \cdot \nabla\nabla\Phi \cdot \s_dagger diff --git a/src/specfem3D/SIEM_index_region.F90 b/src/specfem3D/SIEM_index_region.F90 index 45aa476e6..39b1a72c2 100755 --- a/src/specfem3D/SIEM_index_region.F90 +++ b/src/specfem3D/SIEM_index_region.F90 @@ -570,7 +570,7 @@ subroutine SIEM_get_index_region() print *,'Error: invalid NGLL setting for SIEM indexing' print *,' NGLLX/NGLLY/NGLLZ = ',NGLLX,NGLLY,NGLLZ, ' - all must be equal to 5 for SIEM' print *,' NGLLX_INF/NGLLY_INF/NGLLZ_INF = ',NGLLX_INF,NGLLY_INF,NGLLZ_INF,' - all must be equal to 3 for SIEM' - stop 'SIEM indexing only works for NGLLX==NGLLY==NGLLZ==5 and NGLLX_INF==NGLLY_INF==NGLLZ_INF==3 setting for now!' + stop 'SIEM indexing only works for NGLLX == 5 and NGLLX_INF == 3 setting for now!' endif ! Level-1 solver--------------- diff --git a/src/specfem3D/SIEM_poisson.F90 b/src/specfem3D/SIEM_poisson.F90 index c2ef3f4b1..b83f051aa 100644 --- a/src/specfem3D/SIEM_poisson.F90 +++ b/src/specfem3D/SIEM_poisson.F90 @@ -90,7 +90,8 @@ subroutine poisson_stiffness(iregion,nelmt,nnode,ibool,xstore,ystore,zstore, & !dnx = NGLLX-1; dny = NGLLY-1; dnz = NGLLZ-1 - storekmat(:,:,:) = zero; dprecon(:) = zero + storekmat(:,:,:) = zero + dprecon(:) = zero do i_elmt = 1,nelmt ! suppress fictitious elements in central cube @@ -120,7 +121,11 @@ subroutine poisson_stiffness(iregion,nelmt,nnode,ibool,xstore,ystore,zstore, & deriv = real(matmul(jac,dlagrange_gll(:,i,:)),kind=CUSTOM_REAL) ! use der for gll kmat = kmat + real(matmul(transpose(deriv),deriv)*detjac*gll_weights(i),kind=CUSTOM_REAL) enddo + + ! stiffness matrix storekmat(:,:,i_elmt) = kmat + + ! preconditioner do k = 1,NGLLCUBE dprecon(egdof(k)) = dprecon(egdof(k))+kmat(k,k) enddo @@ -190,7 +195,8 @@ subroutine poisson_stiffnessINF(nelmt,nnode,ibool,xstore,ystore,zstore, & !dnx = NGLLX-1; dny = NGLLY-1; dnz = NGLLZ-2 - storekmat(:,:,:) = zero; dprecon(:) = zero + storekmat(:,:,:) = zero + dprecon(:) = zero do i_elmt = 1,nelmt !ignod=reshape(ibool(1:NGLLX:dnx,1:NGLLY:dny,1:NGLLZ-1:dnz,i_elmt),(/ngnod/)) @@ -225,7 +231,11 @@ subroutine poisson_stiffnessINF(nelmt,nnode,ibool,xstore,ystore,zstore, & deriv = real(matmul(jac,dlagrange_gl(:,i,:)),kind=CUSTOM_REAL) kmat = kmat + real(matmul(transpose(deriv),deriv)*detjac*GLw(i),kind=CUSTOM_REAL) enddo + + ! stiffness matrix storekmat(:,:,i_elmt) = kmat + + ! preconditioner do k = 1,NGLLCUBE dprecon(egdof(k)) = dprecon(egdof(k))+kmat(k,k) enddo @@ -288,12 +298,13 @@ subroutine poisson_stiffness3(iregion,nelmt,nnode,ibool,xstore,ystore,zstore, & !dnx = NGLLX-1; dny = NGLLY-1; dnz = NGLLZ-1 - storekmat(:,:,:) = zero; dprecon(:) = zero + storekmat(:,:,:) = zero + dprecon(:) = zero do i_elmt = 1,nelmt ! suppress fictitious elements in central cube if (iregion == IREGION_INNER_CORE) then - if (idoubling_inner_core(i_elmt) == IFLAG_IN_FICTITIOUS_CUBE)cycle + if (idoubling_inner_core(i_elmt) == IFLAG_IN_FICTITIOUS_CUBE) cycle endif !ignod=reshape(ibool(1:NGLLX:dnx,1:NGLLY:dny,1:NGLLZ:dnz,i_elmt),(/ngnod/)) ! this is wrong!!!! @@ -318,7 +329,11 @@ subroutine poisson_stiffness3(iregion,nelmt,nnode,ibool,xstore,ystore,zstore, & deriv = real(matmul(jac,dlagrange_gll(:,i,:)),kind=CUSTOM_REAL) ! use der for gll kmat = kmat + real(matmul(transpose(deriv),deriv)*detjac*gll_weights(i),kind=CUSTOM_REAL) enddo + + ! stiffness matrix storekmat(:,:,i_elmt) = kmat + + ! preconditioner do k = 1,NGLLCUBE_INF dprecon(egdof(k)) = dprecon(egdof(k))+kmat(k,k) enddo @@ -378,7 +393,8 @@ subroutine poisson_stiffnessINF3(nelmt,nnode,ibool,xstore,ystore,zstore,nnode1, !dnx = NGLLX-1; dny = NGLLY-1; dnz = NGLLZ-2 - storekmat(:,:,:) = zero; dprecon(:) = zero + storekmat(:,:,:) = zero + dprecon(:) = zero do i_elmt = 1,nelmt !ignod=reshape(ibool(1:NGLLX:dnx,1:NGLLY:dny,1:NGLLZ-1:dnz,i_elmt),(/ngnod/)) @@ -413,7 +429,11 @@ subroutine poisson_stiffnessINF3(nelmt,nnode,ibool,xstore,ystore,zstore,nnode1, deriv = real(matmul(jac,dlagrange_gl(:,i,:)),kind=CUSTOM_REAL) kmat = kmat + real(matmul(transpose(deriv),deriv)*detjac*GLw(i),kind=CUSTOM_REAL) enddo + + ! stiffness matrix storekmat(:,:,i_elmt) = kmat + + ! preconditioner do k = 1,NGLLCUBE_INF dprecon(egdof(k)) = dprecon(egdof(k))+kmat(k,k) enddo diff --git a/src/specfem3D/SIEM_prepare_iteration.F90 b/src/specfem3D/SIEM_prepare_iteration.F90 index b91ef3832..10373ded2 100644 --- a/src/specfem3D/SIEM_prepare_iteration.F90 +++ b/src/specfem3D/SIEM_prepare_iteration.F90 @@ -83,7 +83,7 @@ subroutine SIEM_prepare_iteration() ! method in make_gravity.f90 ! compute background gravitational field only once ! WARNING: We have not currently using this for the time marching. - if (NUMBER_OF_THIS_RUN == 1) call compute_background_gravity_SIEM() + if (NUMBER_OF_THIS_RUN == 1) call SIEM_compute_background_gravity() ! initialize for dynamic solver pgrav1(:) = 0.0_CUSTOM_REAL @@ -109,6 +109,68 @@ end subroutine SIEM_prepare_iteration subroutine SIEM_interpolate_gravity() +! interpolates Level-1 forward pgrav1 values onto Level-2 pgrav + + use constants, only: myrank,ADD_TRINF + + use constants_solver, only: NSPEC_CRUST_MANTLE,NSPEC_OUTER_CORE,NSPEC_INNER_CORE, & + NSPEC_TRINFINITE,NSPEC_INFINITE, & + NGLOB_CRUST_MANTLE,NGLOB_OUTER_CORE,NGLOB_INNER_CORE,NGLOB_TRINFINITE,NGLOB_INFINITE + + use specfem_par_full_gravity, only: is_active_gll,igll_active_on, & + inode_elmt_ic,inode_elmt_oc,inode_elmt_cm,inode_elmt_trinf,inode_elmt_inf, & + inode_map_ic,inode_map_oc,inode_map_cm,inode_map_trinf,inode_map_inf, & + nmir_ic,nmir_oc,nmir_cm,nmir_trinf,nmir_inf, & + nnode_ic1,nnode_oc1,nnode_cm1,nnode_trinf1,nnode_inf1, & + gdof_ic1,gdof_oc1,gdof_cm1,gdof_trinf1,gdof_inf1 + + use specfem_par_full_gravity, only: & + pgrav1,pgrav_ic1,pgrav_oc1,pgrav_cm1,pgrav_trinf1,pgrav_inf1, & + pgrav_ic,pgrav_oc,pgrav_cm,pgrav_trinf,pgrav_inf + + use siem_solver_mpi, only: interpolate3to5 + + implicit none + + ! interpolate the gravity values + ! transfer to array for each region (cm, ic, oc) + pgrav_ic1(:) = pgrav1(gdof_ic1(:)) + pgrav_oc1(:) = pgrav1(gdof_oc1(:)) + pgrav_cm1(:) = pgrav1(gdof_cm1(:)) + if (ADD_TRINF) pgrav_trinf1(:) = pgrav1(gdof_trinf1(:)) + pgrav_inf1(:) = pgrav1(gdof_inf1(:)) + + ! interpolate values for each region + call interpolate3to5(NSPEC_INNER_CORE,NGLOB_INNER_CORE,nnode_ic1, & + inode_elmt_ic,nmir_ic,inode_map_ic,is_active_gll,igll_active_on,pgrav_ic1,pgrav_ic) + + call interpolate3to5(NSPEC_OUTER_CORE,NGLOB_OUTER_CORE,nnode_oc1, & + inode_elmt_oc,nmir_oc,inode_map_oc,is_active_gll,igll_active_on,pgrav_oc1,pgrav_oc) + + call interpolate3to5(NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE,nnode_cm1, & + inode_elmt_cm,nmir_cm,inode_map_cm,is_active_gll,igll_active_on,pgrav_cm1,pgrav_cm) + + if (ADD_TRINF) then + call interpolate3to5(NSPEC_TRINFINITE,NGLOB_TRINFINITE,nnode_trinf1, & + inode_elmt_trinf,nmir_trinf,inode_map_trinf,is_active_gll,igll_active_on,pgrav_trinf1,pgrav_trinf) + endif + + call interpolate3to5(NSPEC_INFINITE,NGLOB_INFINITE,nnode_inf1, & + inode_elmt_inf,nmir_inf,inode_map_inf,is_active_gll,igll_active_on,pgrav_inf1,pgrav_inf) + + !debug + if (myrank == 0) print *,' Finished gravity interpolation of loaded (saved) forward wavefield' + + end subroutine SIEM_interpolate_gravity + +! +!------------------------------------------------------------------------------- +! + + subroutine SIEM_interpolate_backward_gravity() + +! interpolates Level-1 backward b_pgrav1 values onto Level-2 b_pgrav + use constants, only: myrank,ADD_TRINF use constants_solver, only: NSPEC_CRUST_MANTLE,NSPEC_OUTER_CORE,NSPEC_INNER_CORE, & @@ -116,13 +178,15 @@ subroutine SIEM_interpolate_gravity() NGLOB_CRUST_MANTLE,NGLOB_OUTER_CORE,NGLOB_INNER_CORE,NGLOB_TRINFINITE,NGLOB_INFINITE use specfem_par_full_gravity, only: is_active_gll,igll_active_on, & - b_pgrav1,b_pgrav_ic1,b_pgrav_oc1,b_pgrav_cm1,b_pgrav_trinf1,b_pgrav_inf1, & - b_pgrav_ic,b_pgrav_oc,b_pgrav_cm,b_pgrav_trinf,b_pgrav_inf, & - gdof_ic1,gdof_oc1,gdof_cm1,gdof_trinf1,gdof_inf1, & inode_elmt_ic,inode_elmt_oc,inode_elmt_cm,inode_elmt_trinf,inode_elmt_inf, & inode_map_ic,inode_map_oc,inode_map_cm,inode_map_trinf,inode_map_inf, & nmir_ic,nmir_oc,nmir_cm,nmir_trinf,nmir_inf, & - nnode_ic1,nnode_oc1,nnode_cm1,nnode_trinf1,nnode_inf1 + nnode_ic1,nnode_oc1,nnode_cm1,nnode_trinf1,nnode_inf1, & + gdof_ic1,gdof_oc1,gdof_cm1,gdof_trinf1,gdof_inf1 + + use specfem_par_full_gravity, only: & + b_pgrav1,b_pgrav_ic1,b_pgrav_oc1,b_pgrav_cm1,b_pgrav_trinf1,b_pgrav_inf1, & + b_pgrav_ic,b_pgrav_oc,b_pgrav_cm,b_pgrav_trinf,b_pgrav_inf use siem_solver_mpi, only: interpolate3to5 @@ -155,9 +219,9 @@ subroutine SIEM_interpolate_gravity() inode_elmt_inf,nmir_inf,inode_map_inf,is_active_gll,igll_active_on,b_pgrav_inf1,b_pgrav_inf) !debug - if (myrank == 0) print *,' Finished gravity interpolation of loaded (saved) wavefield' + if (myrank == 0) print *,' Finished gravity interpolation of loaded (saved) backward wavefield' - end subroutine SIEM_interpolate_gravity + end subroutine SIEM_interpolate_backward_gravity ! !------------------------------------------------------------------------------- @@ -204,11 +268,13 @@ subroutine initialise_inbuilt_solver() if (istat /= 0) stop 'Error allocating ndscale1 array' ndscale1(:) = 1.0_CUSTOM_REAL + ! scaling based on inverted preconditioner values do i = 1,neq1 val = sqrt(abs(dprecon1(i))) if (val /= 0.0_CUSTOM_REAL) ndscale1(i) = 1.0_CUSTOM_REAL / val enddo + ! scales stiffness matrices (storekmat_*) do i_elmt = 1,NSPEC_INNER_CORE inodes1 = inode_elmt_ic1(:,i_elmt) igdof1 = gdof_ic1(inodes1) @@ -267,7 +333,7 @@ subroutine initialise_inbuilt_solver() enddo !debug - if (myrank == 0) print *,'builtin solver: check stiffness',minval(ndscale1),maxval(ndscale1) + if (myrank == 0) print *,'builtin solver: scaling min/max = ',minval(ndscale1),'/',maxval(ndscale1) if (myrank == 0) print * call synchronize_all() @@ -291,7 +357,7 @@ end subroutine initialise_inbuilt_solver !------------------------------------------------------------------------------- ! - subroutine compute_background_gravity_SIEM() + subroutine SIEM_compute_background_gravity() use specfem_par @@ -699,7 +765,7 @@ end subroutine write_ensight_perelementAS call flush_IMAIN() endif - end subroutine compute_background_gravity_SIEM + end subroutine SIEM_compute_background_gravity ! !------------------------------------------------------------------------------- @@ -879,9 +945,11 @@ subroutine SIEM_finalize() endif ! free arrays + ! Level-1 solver if (allocated(storederiv_cm1)) deallocate(storederiv_cm1,storerhojw_cm1) if (allocated(storederiv_ic1)) deallocate(storederiv_ic1,storerhojw_ic1) if (allocated(storederiv_oc1)) deallocate(storederiv_oc1,storerhojw_oc1) + ! Level-2 solver if (allocated(storederiv_cm)) deallocate(storederiv_cm,storerhojw_cm) if (allocated(storederiv_ic)) deallocate(storederiv_ic,storerhojw_ic) if (allocated(storederiv_oc)) deallocate(storederiv_oc,storerhojw_oc) @@ -898,18 +966,21 @@ subroutine SIEM_finalize() if (allocated(ndscale)) deallocate(ndscale) ! gravity perturbation arrays + ! Level-1 solver if (allocated(pgrav1)) deallocate(pgrav1) if (allocated(pgrav_ic1)) deallocate(pgrav_ic1,pgrav_oc1,pgrav_cm1,pgrav_trinf1,pgrav_inf1) - if (allocated(pgrav)) deallocate(pgrav) - if (allocated(pgrav_ic)) deallocate(pgrav_ic,pgrav_oc,pgrav_cm,pgrav_trinf,pgrav_inf) if (allocated(dprecon1)) deallocate(dprecon1,gravload1) if (allocated(b_gravload1)) deallocate(b_gravload1) + ! Level-2 solver + if (allocated(pgrav)) deallocate(pgrav) + if (allocated(pgrav_ic)) deallocate(pgrav_ic,pgrav_oc,pgrav_cm,pgrav_trinf,pgrav_inf) if (allocated(dprecon)) deallocate(dprecon,gravload) - ! sparse petsc solver arrays + ! Level-1 solver if (allocated(l2gdof1)) deallocate(l2gdof1) if (allocated(krow_sparse1)) deallocate(krow_sparse1,kcol_sparse1) if (allocated(kgrow_sparse1)) deallocate(kgrow_sparse1,kgcol_sparse1) + ! Level-2 solver if (allocated(l2gdof)) deallocate(l2gdof) if (allocated(krow_sparse)) deallocate(krow_sparse,kcol_sparse) if (allocated(kgrow_sparse)) deallocate(kgrow_sparse,kgcol_sparse) diff --git a/src/specfem3D/SIEM_prepare_solver.F90 b/src/specfem3D/SIEM_prepare_solver.F90 index e11edb629..f9d5b3fb5 100644 --- a/src/specfem3D/SIEM_prepare_solver.F90 +++ b/src/specfem3D/SIEM_prepare_solver.F90 @@ -100,6 +100,9 @@ subroutine SIEM_prepare_solver() ! get MPI starting time tstart = wtime() + ! indexify regions (determine inode_elmt_*,inode_map_*,nmir_* arrays) + call SIEM_get_index_region() + ! sets the stiffness matrices for Poisson's solver ! calculate dprecon ! allocates gravload, regional pgravs (e.g. pgrav_cm1) @@ -603,26 +606,8 @@ subroutine SIEM_prepare_solver_poisson() endif ! estimated memory size required in MB - ! inode_elmt_cm,.. - sizeval = dble(NGLLCUBE)*dble(NSPEC_CRUST_MANTLE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE)*dble(NSPEC_INNER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE)*dble(NSPEC_OUTER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE)*dble(NSPEC_TRINFINITE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE)*dble(NSPEC_INFINITE)*dble(SIZE_INTEGER) - ! inode_elmt_cm1,.. - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NSPEC_CRUST_MANTLE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NSPEC_INNER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NSPEC_OUTER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NSPEC_TRINFINITE)*dble(SIZE_INTEGER) - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NSPEC_INFINITE)*dble(SIZE_INTEGER) - ! inode_map_cm,.. - sizeval = sizeval + 2.d0*dble(NGLOB_CRUST_MANTLE)*dble(SIZE_INTEGER) - sizeval = sizeval + 2.d0*dble(NGLOB_INNER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + 2.d0*dble(NGLOB_OUTER_CORE)*dble(SIZE_INTEGER) - sizeval = sizeval + 2.d0*dble(NGLOB_TRINFINITE)*dble(SIZE_INTEGER) - sizeval = sizeval + 2.d0*dble(NGLOB_INFINITE)*dble(SIZE_INTEGER) ! storekmat_crust_mantle1,dprecon_crust_mantle1 - sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NGLLCUBE_INF)*dble(NSPEC_CRUST_MANTLE)*dble(CUSTOM_REAL) + sizeval = dble(NGLLCUBE_INF)*dble(NGLLCUBE_INF)*dble(NSPEC_CRUST_MANTLE)*dble(CUSTOM_REAL) sizeval = sizeval + dble(nnode_cm1)*dble(CUSTOM_REAL) ! storekmat_inner_core1,dprecon_inner_core1 sizeval = sizeval + dble(NGLLCUBE_INF)*dble(NGLLCUBE_INF)*dble(NSPEC_INNER_CORE)*dble(CUSTOM_REAL) @@ -672,37 +657,6 @@ subroutine SIEM_prepare_solver_poisson() call flush_IMAIN() endif - ! allocate inode arrays - allocate(inode_elmt_cm(NGLLCUBE,NSPEC_CRUST_MANTLE), & - inode_elmt_ic(NGLLCUBE,NSPEC_INNER_CORE), & - inode_elmt_oc(NGLLCUBE,NSPEC_OUTER_CORE), & - inode_elmt_trinf(NGLLCUBE,NSPEC_TRINFINITE), & - inode_elmt_inf(NGLLCUBE,NSPEC_INFINITE),stat=ier) - if (ier /= 0) stop 'Error allocating inode_elmt_cm,.. arrays' - inode_elmt_cm(:,:) = 0; inode_elmt_ic(:,:) = 0; inode_elmt_oc(:,:) = 0 - inode_elmt_trinf(:,:) = 0; inode_elmt_inf(:,:) = 0 - - allocate(inode_elmt_cm1(NGLLCUBE_INF,NSPEC_CRUST_MANTLE), & - inode_elmt_ic1(NGLLCUBE_INF,NSPEC_INNER_CORE), & - inode_elmt_oc1(NGLLCUBE_INF,NSPEC_OUTER_CORE), & - inode_elmt_trinf1(NGLLCUBE_INF,NSPEC_TRINFINITE), & - inode_elmt_inf1(NGLLCUBE_INF,NSPEC_INFINITE),stat=ier) - if (ier /= 0) stop 'Error allocating inode_elmt_cm1,.. arrays' - inode_elmt_cm1(:,:) = 0; inode_elmt_ic1(:,:) = 0; inode_elmt_oc1(:,:) = 0 - inode_elmt_trinf1(:,:) = 0; inode_elmt_inf1(:,:) = 0 - - allocate(inode_map_ic(2,NGLOB_INNER_CORE), & - inode_map_oc(2,NGLOB_OUTER_CORE), & - inode_map_cm(2,NGLOB_CRUST_MANTLE), & - inode_map_trinf(2,NGLOB_TRINFINITE), & - inode_map_inf(2,NGLOB_INFINITE),stat=ier) - if (ier /= 0) stop 'Error allocating inode_map_ic,.. arrays' - inode_map_ic(:,:) = 0; inode_map_oc(:,:) = 0; inode_map_cm(:,:) = 0 - inode_map_trinf(:,:) = 0; inode_map_inf(:,:) = 0 - - ! indexify regions - call SIEM_get_index_region() - ! Level-1 solver------------------- allocate(storekmat_crust_mantle1(NGLLCUBE_INF,NGLLCUBE_INF,NSPEC_CRUST_MANTLE), & dprecon_crust_mantle1(nnode_cm1)) @@ -743,8 +697,13 @@ subroutine SIEM_prepare_solver_poisson() pgrav_ic1(:) = 0.0_CUSTOM_REAL; pgrav_oc1(:) = 0.0_CUSTOM_REAL; pgrav_cm1(:) = 0.0_CUSTOM_REAL; pgrav_trinf1(:) = 0.0_CUSTOM_REAL; pgrav_inf1(:) = 0.0_CUSTOM_REAL - allocate(dprecon1(0:neq1),gravload1(0:neq1)) - dprecon1(:) = 0.0_CUSTOM_REAL; gravload1(:) = 0.0_CUSTOM_REAL + ! Level-1 solver RHS + allocate(gravload1(0:neq1)) + gravload1(:) = 0.0_CUSTOM_REAL + + ! preconditioner for PCG solver + allocate(dprecon1(0:neq1)) + dprecon1(:) = 0.0_CUSTOM_REAL if (SIMULATION_TYPE == 3) then allocate(b_pgrav_ic(NGLOB_INNER_CORE_ADJOINT), & @@ -798,6 +757,7 @@ subroutine SIEM_prepare_solver_poisson() debug_rho_kl_cm(:,:,:,:,:) = 0.0_CUSTOM_REAL endif + ! setup stiffness matrix (storekmat_*) and preconditioner (dprecon_*) arrays for Level-1 solver ! crust mantle call poisson_stiffness3(IREGION_CRUST_MANTLE,NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE, & ibool_crust_mantle,xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, & @@ -827,7 +787,7 @@ subroutine SIEM_prepare_solver_poisson() call synchronize_all() - ! assemble stiffness matrices + ! assemble preconditioner matrices ! assemble across the MPI processes in a region ! crust_mantle call assemble_MPI_scalar(NPROCTOT_VAL,nnode_cm1,dprecon_crust_mantle1, & @@ -928,6 +888,7 @@ subroutine SIEM_prepare_solver_poisson() call flush_IMAIN() endif + ! Level-2 solver arrays allocate(storekmat_crust_mantle(NGLLCUBE,NGLLCUBE,NSPEC_CRUST_MANTLE), & dprecon_crust_mantle(NGLOB_CRUST_MANTLE)) storekmat_crust_mantle(:,:,:) = 0.0_CUSTOM_REAL; dprecon_crust_mantle(:) = 0.0_CUSTOM_REAL diff --git a/src/specfem3D/SIEM_solve.F90 b/src/specfem3D/SIEM_solve.F90 index 33483479c..f27684ec2 100644 --- a/src/specfem3D/SIEM_solve.F90 +++ b/src/specfem3D/SIEM_solve.F90 @@ -80,18 +80,17 @@ subroutine solve_poisson_equation() use siem_math_library_mpi, only: maxvec - use siem_solver_petsc, only: petsc_set_vector1,petsc_solve1, & - petsc_set_vector,petsc_solve + use siem_solver_petsc, only: petsc_set_vector1, petsc_solve1, & + petsc_set_vector, petsc_solve - use siem_poisson, only: compute_poisson_load,compute_poisson_load3 !,poisson_gravity + use siem_poisson, only: compute_poisson_load, compute_poisson_load3 !,poisson_gravity - use siem_solver_mpi, only: cg_solver,cg_solver3,diagpcg_solver,diagpcg_solver3, interpolate3to5 + use siem_solver_mpi, only: cg_solver, cg_solver3, diagpcg_solver, diagpcg_solver3, interpolate3to5 implicit none ! Local variables real(kind=CUSTOM_REAL) :: maxf,upscale - real(kind=CUSTOM_REAL),parameter :: rone = 1.0_CUSTOM_REAL ! load for Level-1 solver ! Compute the RHS (gravload1) @@ -105,8 +104,8 @@ subroutine solve_poisson_equation() if (POISSON_SOLVER == ISOLVER_BUILTIN) then ! builtin solver maxf = maxvec(abs(gravload1)) - upscale = rone - if (maxf > zero) upscale = rone / maxf + upscale = 1.0_CUSTOM_REAL + if (maxf > 0.0_CUSTOM_REAL) upscale = 1.0_CUSTOM_REAL / maxf if (CG_SCALING) then gravload1(:) = gravload1(:) * ndscale1(:) * upscale @@ -221,6 +220,7 @@ subroutine solve_poisson_equation_backward() gdof_trinf, gdof_trinf1, inode_elmt_trinf, inode_map_trinf, nmir_trinf, & gdof_inf, gdof_inf1, inode_elmt_inf, inode_map_inf, nmir_inf, & nnode_cm1,nnode_ic1,nnode_oc1,nnode_trinf1,nnode_inf1, & + ndscale1,neq1, & is_active_gll,igll_active_on, & CG_SCALING @@ -228,24 +228,35 @@ subroutine solve_poisson_equation_backward() b_pgrav1,b_pgrav_ic1,b_pgrav_oc1,b_pgrav_cm1,b_pgrav_trinf1,b_pgrav_inf1, & b_pgrav,b_pgrav_ic,b_pgrav_oc,b_pgrav_cm,b_pgrav_trinf,b_pgrav_inf - use siem_solver_petsc, only: petsc_set_vector1,petsc_solve1 + use siem_math_library_mpi, only: maxvec + + use siem_solver_petsc, only: petsc_set_vector1, petsc_solve1 use siem_poisson, only: compute_backward_poisson_load3 !,poisson_gravity - use siem_solver_mpi, only: cg_solver,cg_solver3,diagpcg_solver,diagpcg_solver3, interpolate3to5 + use siem_solver_mpi, only: cg_solver, cg_solver3, diagpcg_solver, diagpcg_solver3, interpolate3to5 implicit none - - ! safety check - if (POISSON_SOLVER == ISOLVER_BUILTIN) then - print *,'Error. Adjoint FG only implemented for petsc solver' - stop 'FULL_GRAVITY adjoint FG not fully implemented yet of builtin solver' - endif + ! local parameters + real(kind=CUSTOM_REAL) :: maxf,upscale ! load for Level-1 solver ! Compute the RHS (b_gravload1) call compute_backward_poisson_load3() - if (POISSON_SOLVER == ISOLVER_PETSC) then + if (POISSON_SOLVER == ISOLVER_BUILTIN) then + ! builtin solver + maxf = maxvec(abs(b_gravload1)) + upscale = 1.0_CUSTOM_REAL + if (maxf > 0.0_CUSTOM_REAL) upscale = 1.0_CUSTOM_REAL / maxf + if (CG_SCALING) then + b_gravload1(:) = b_gravload1(:) * ndscale1(:) * upscale + b_pgrav1(1:) = upscale * b_pgrav1(1:) / ndscale1(1:) + endif + call cg_solver3(myrank,neq1,b_pgrav1,b_gravload1) + if (CG_SCALING) then + b_pgrav1(:) = ndscale1(:) * b_pgrav1(:) / upscale + endif + else ! PETSc solver call petsc_set_vector1(b_gravload1) call petsc_solve1(b_pgrav1(1:)) @@ -306,7 +317,7 @@ end subroutine solve_poisson_equation_backward !------------------------------------------------------------------------------- ! - subroutine SIEM_compute_element_add_full_gravity(ispec,nspec,nglob,gravity_rho,deriv_loc,ibool,pgrav,rho_s_H) + subroutine SIEM_solve_element_add_full_gravity(ispec,nspec,nglob,gravity_rho,deriv_loc,ibool,pgrav,rho_s_H) ! routine for crust/mantle and inner core elements to add full gravity contribution to rho_s_H array @@ -397,4 +408,4 @@ subroutine SIEM_compute_element_add_full_gravity(ispec,nspec,nglob,gravity_rho,d enddo enddo - end subroutine SIEM_compute_element_add_full_gravity + end subroutine SIEM_solve_element_add_full_gravity diff --git a/src/specfem3D/SIEM_solver_mpi.F90 b/src/specfem3D/SIEM_solver_mpi.F90 index ad299888a..c3e36d169 100644 --- a/src/specfem3D/SIEM_solver_mpi.F90 +++ b/src/specfem3D/SIEM_solver_mpi.F90 @@ -923,6 +923,7 @@ subroutine scatter_and_assemble3(neq,array,array_g) real(kind=CUSTOM_REAL),intent(in) :: array(0:neq) real(kind=CUSTOM_REAL),intent(out) :: array_g(0:neq) + ! local parameters real(kind=CUSTOM_REAL) :: array_ic(nnode_ic1),array_oc(nnode_oc1),array_cm(nnode_cm1), & array_trinf(nnode_trinf1),array_inf(nnode_inf1) @@ -1000,32 +1001,41 @@ subroutine interpolate3to5(nelmt,nnode,nnode1,inode_elmt,nmir,inode_map,is_activ use siem_gll_library, only: gll_quadrature3inNGLL,zwgljd implicit none - integer,intent(in) :: nelmt,nnode,nnode1,inode_elmt(NGLLCUBE,nelmt),nmir(nnode), & - inode_map(2,nnode),igll_active_on(NGLLCUBE_INF) + integer,intent(in) :: nelmt,nnode,nnode1 + integer,intent(in) :: inode_elmt(NGLLCUBE,nelmt) + integer,intent(in) :: nmir(nnode) + integer,intent(in) :: inode_map(2,nnode) + integer,intent(in) :: igll_active_on(NGLLCUBE_INF) logical,intent(in) :: is_active_gll(NGLLCUBE) + real(kind=CUSTOM_REAL),intent(in) :: x3(nnode1) ! array for 3 GLLX points - real(kind=CUSTOM_REAL),intent(out) :: x5(nnode) ! aray for 5 GLLX points + real(kind=CUSTOM_REAL),intent(out) :: x5(nnode) ! array for 5 GLLX points - double precision :: lagrange_gll3inNGLL(NGLLCUBE,27),xigll(NGLLX),wxgll(NGLLX), & - xigll1(NGLLX_INF),wxgll1(NGLLX_INF) + ! local parameters + double precision :: lagrange_gll3inNGLL(NGLLCUBE,27) + double precision :: xigll(NGLLX),wxgll(NGLLX) + double precision :: xigll1(NGLLX_INF),wxgll1(NGLLX_INF) - integer :: i_node,ielmt,igll,inode1,inodes1(NGLLCUBE_INF) + integer :: i_node,ielmt,igll,inode1 + integer :: inodes1(NGLLCUBE_INF) call zwgljd(xigll1,wxgll1,NGLLX_INF,0.d0,0.d0) call zwgljd(xigll,wxgll,NGLLX,0.d0,0.d0) call gll_quadrature3inNGLL(NDIM,NGLLX,NGLLCUBE,xigll,xigll1,lagrange_gll3inNGLL) - ! inner core - x5 = 0.0_CUSTOM_REAL - do i_node = 1,nnode!NGLOB_INNER_CORE - inode1 = nmir(i_node) + x5(:) = 0.0_CUSTOM_REAL + + do i_node = 1,nnode ielmt = inode_map(1,i_node) - if (ielmt <= 0) then - cycle ! skip fictitious nodes - endif + + ! skip fictitious nodes + if (ielmt <= 0) cycle + igll = inode_map(2,i_node) if (is_active_gll(igll)) then + ! takes value of mirrored node in Level-1 array + inode1 = nmir(i_node) x5(i_node) = x3(inode1) else ! interpolate values diff --git a/src/specfem3D/compute_forces_crust_mantle_Dev.F90 b/src/specfem3D/compute_forces_crust_mantle_Dev.F90 index 3a852b42a..83e04ecba 100644 --- a/src/specfem3D/compute_forces_crust_mantle_Dev.F90 +++ b/src/specfem3D/compute_forces_crust_mantle_Dev.F90 @@ -374,8 +374,8 @@ subroutine compute_forces_crust_mantle_Dev( NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, & if (GRAVITY_VAL) then ! full gravity if (FULL_GRAVITY_VAL .and. .not. DISCARD_GCONTRIB) then - call SIEM_compute_element_add_full_gravity(ispec,NSPEC_CRUST_MANTLE,NGLOB,gravity_rho,deriv(:,:,:,:,ispec),ibool, & - pgrav_crust_mantle,rho_s_H) + call SIEM_solve_element_add_full_gravity(ispec,NSPEC_CRUST_MANTLE,NGLOB,gravity_rho,deriv(:,:,:,:,ispec),ibool, & + pgrav_crust_mantle,rho_s_H) endif #ifdef FORCE_VECTORIZATION diff --git a/src/specfem3D/compute_forces_crust_mantle_noDev.f90 b/src/specfem3D/compute_forces_crust_mantle_noDev.f90 index 08f68092f..60764a8c6 100644 --- a/src/specfem3D/compute_forces_crust_mantle_noDev.f90 +++ b/src/specfem3D/compute_forces_crust_mantle_noDev.f90 @@ -508,8 +508,8 @@ subroutine compute_forces_crust_mantle_noDev(NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, & ! full gravity if (FULL_GRAVITY_VAL .and. .not. DISCARD_GCONTRIB) then - call SIEM_compute_element_add_full_gravity(ispec,NSPEC_CRUST_MANTLE,NGLOB,gravity_rho,deriv_loc,ibool, & - pgrav_crust_mantle,rho_s_H) + call SIEM_solve_element_add_full_gravity(ispec,NSPEC_CRUST_MANTLE,NGLOB,gravity_rho,deriv_loc,ibool, & + pgrav_crust_mantle,rho_s_H) endif do k = 1,NGLLZ diff --git a/src/specfem3D/compute_forces_inner_core_Dev.F90 b/src/specfem3D/compute_forces_inner_core_Dev.F90 index 8387d9f95..926da54a4 100644 --- a/src/specfem3D/compute_forces_inner_core_Dev.F90 +++ b/src/specfem3D/compute_forces_inner_core_Dev.F90 @@ -293,8 +293,8 @@ subroutine compute_forces_inner_core_Dev( NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, & if (GRAVITY_VAL) then ! full gravity if (FULL_GRAVITY_VAL .and. .not. DISCARD_GCONTRIB) then - call SIEM_compute_element_add_full_gravity(ispec,NSPEC_INNER_CORE,NGLOB,gravity_rho,deriv(:,:,:,:,ispec),ibool, & - pgrav_inner_core,rho_s_H) + call SIEM_solve_element_add_full_gravity(ispec,NSPEC_INNER_CORE,NGLOB,gravity_rho,deriv(:,:,:,:,ispec),ibool, & + pgrav_inner_core,rho_s_H) endif #ifdef FORCE_VECTORIZATION diff --git a/src/specfem3D/compute_forces_inner_core_noDev.f90 b/src/specfem3D/compute_forces_inner_core_noDev.f90 index 0ceae322b..6c85cce82 100644 --- a/src/specfem3D/compute_forces_inner_core_noDev.f90 +++ b/src/specfem3D/compute_forces_inner_core_noDev.f90 @@ -434,8 +434,8 @@ subroutine compute_forces_inner_core_noDev( NSPEC_STR_OR_ATT,NGLOB,NSPEC_ATT, & ! full gravity if (FULL_GRAVITY_VAL .and. .not. DISCARD_GCONTRIB) then - call SIEM_compute_element_add_full_gravity(ispec,NSPEC_INNER_CORE,NGLOB,gravity_rho,deriv_loc,ibool, & - pgrav_inner_core,rho_s_H) + call SIEM_solve_element_add_full_gravity(ispec,NSPEC_INNER_CORE,NGLOB,gravity_rho,deriv_loc,ibool, & + pgrav_inner_core,rho_s_H) endif do k = 1,NGLLZ diff --git a/src/specfem3D/finalize_simulation.F90 b/src/specfem3D/finalize_simulation.F90 index d488bdfed..1ee16f937 100644 --- a/src/specfem3D/finalize_simulation.F90 +++ b/src/specfem3D/finalize_simulation.F90 @@ -81,6 +81,11 @@ subroutine finalize_simulation() if (myrank == 0 ) call finish_vtkwindow() endif + ! finalize full gravity + if (FULL_GRAVITY) then + call SIEM_finalize() + endif + ! adios finalizes if (ADIOS_ENABLED) then call finalize_adios() diff --git a/src/specfem3D/iterate_time.F90 b/src/specfem3D/iterate_time.F90 index 0096052ea..48f4d7f99 100644 --- a/src/specfem3D/iterate_time.F90 +++ b/src/specfem3D/iterate_time.F90 @@ -252,13 +252,9 @@ subroutine iterate_time() endif ! full gravity - if (FULL_GRAVITY) then + if (SIMULATION_TYPE == 3 .and. FULL_GRAVITY) then ! calculate the gravity kernels (convolution) using SIEM - if (SIMULATION_TYPE == 3) then - call SIEM_compute_gravity_kernels() - endif - ! finalize - call SIEM_finalize() + call SIEM_compute_gravity_kernels() endif ! close the huge file that contains a dump of all the time steps to disk diff --git a/src/specfem3D/iterate_time_undoatt.F90 b/src/specfem3D/iterate_time_undoatt.F90 index d9e11f418..5c23708ba 100644 --- a/src/specfem3D/iterate_time_undoatt.F90 +++ b/src/specfem3D/iterate_time_undoatt.F90 @@ -628,13 +628,9 @@ subroutine iterate_time_undoatt() endif ! full gravity - if (FULL_GRAVITY) then + if (SIMULATION_TYPE == 3 .and. FULL_GRAVITY) then ! calculate the gravity kernels (convolution) using SIEM - if (SIMULATION_TYPE == 3) then - call SIEM_compute_gravity_kernels() - endif - ! finalize - call SIEM_finalize() + call SIEM_compute_gravity_kernels() endif ! close the huge file that contains a dump of all the time steps to disk diff --git a/src/specfem3D/prepare_gravity.f90 b/src/specfem3D/prepare_gravity.f90 index 8ee5269cf..36d97b936 100644 --- a/src/specfem3D/prepare_gravity.f90 +++ b/src/specfem3D/prepare_gravity.f90 @@ -36,6 +36,8 @@ subroutine prepare_gravity() use specfem_par_innercore use specfem_par_full_gravity + use siem_math_library, only: compute_g_gradg,compute_g_gradg_elliptical + implicit none ! local parameters @@ -57,6 +59,20 @@ subroutine prepare_gravity() integer :: int_radius,idummy,nspl_gravity,iglob,ier integer :: index_fluid + ! full gravity + double precision,dimension(:),allocatable :: g_rad + double precision :: dotrho + double precision :: gvec(NDIM),gradg(6) + ! ellipticity spline + integer :: nspl + double precision,dimension(NR_DENSITY) :: rspl,e_spline,e_spline2,eta_spline,eta_spline2 + ! table + double precision :: eps,eta + double precision,dimension(:),allocatable :: eps_rad,eta_rad,dotrho_rad + ! rotation rate + double precision :: omega,twothirdOmega2 + double precision,parameter :: two_third = 2.d0/3.d0 + ! debugging logical, parameter :: DEBUG = .false. @@ -132,6 +148,25 @@ subroutine prepare_gravity() ! for example: r = RCMB/R_PLANET -> int_radius = RCMB/R_PLANET * NRAD_GRAVITY / range_max enddo + ! for full gravity + if (FULL_GRAVITY_VAL) then + ! gravity kernels need g values and grad(g) + if (SIMULATION_TYPE == 3) then + ! note: this setup might not be needed if one takes gxl,.. and Hxxl,.. directly + ! and store them as g_cm and gradg_cm arrays + allocate(g_rad(NRAD_GRAVITY)) + ! prepare ellipticity splines for epsilon and eta values + if (ELLIPTICITY_VAL) then + allocate(eps_rad(NRAD_GRAVITY),eta_rad(NRAD_GRAVITY),dotrho_rad(NRAD_GRAVITY)) + ! ellipticity + call make_ellipticity(nspl,rspl,e_spline,e_spline2,eta_spline,eta_spline2) + ! rotation rate omega + omega = TWO_PI / (HOURS_PER_DAY * SECONDS_PER_HOUR) ! rotation rate in rad/sec + twothirdOmega2 = two_third * omega*omega / (PI*GRAV*RHOAV) + endif + endif + endif + ! store g, rho and dg/dr=dg using normalized radius in lookup table every 100 m ! get density and velocity from PREM model using dummy doubling flag ! this assumes that the gravity perturbations are small and smooth @@ -169,6 +204,26 @@ subroutine prepare_gravity() minus_deriv_gravity_table(int_radius) = - dg density_table(int_radius) = rho minus_rho_g_over_kappa_fluid(int_radius) = - g / vp**2 ! for fluid: vp**2 = kappa/rho + + ! for full gravity + if (FULL_GRAVITY_VAL) then + ! gravity kernels need g values and grad(g) + if (SIMULATION_TYPE == 3) then + ! note: this setup might not be needed if one takes gxl,.. and Hxxl,.. directly + ! and store them as g_cm and gradg_cm arrays + ! store table data + g_rad(int_radius) = g + if (ELLIPTICITY_VAL) then + ! spline evaluations for epsilon and eta + call spline_evaluation(rspl,e_spline,e_spline2,nspl,radius,eps) + call spline_evaluation(rspl,eta_spline,eta_spline2,nspl,radius,eta) + ! store into table + eps_rad(int_radius) = eps + eta_rad(int_radius) = eta + dotrho_rad(int_radius) = drhodr + endif + endif + endif enddo ! make sure fluid array is only assigned in outer core between 1222 and 3478 km @@ -254,7 +309,30 @@ subroutine prepare_gravity() ! for full gravity contribution in compute forces if (FULL_GRAVITY_VAL) then + ! store factor rho_g_over_kappa gravity_rho_g_over_kappa_outer_core(iglob) = real(fac,kind=CUSTOM_REAL) + + ! gravity kernels need g values and grad(g) + ! outer core kernel not implemented yet... + !if (SIMULATION_TYPE == 3) then + ! ! get g in outer core + ! g = g_rad(int_radius) + ! rho = density_table(int_radius) + ! if (ELLIPTICITY_VAL) then + ! ! w/ ellipticity + ! dotrho = dotrho_rad(int_radius) + ! eps = eps_rad(int_radius) + ! eta = eta_rad(int_radius) + ! ! compute g + ! call compute_g_gradg_elliptical(NDIM,radius,theta,phi,rho,dotrho,g,eps,eta,twothirdOmega2,gvec) + ! else + ! ! no ellipticity + ! ! compute g + ! call compute_g_gradg(NDIM,radius,theta,phi,rho,g,gvec) + ! endif + ! ! in outer core we need only g + ! g_oc(:,iglob) = real(gvec(:),kind=CUSTOM_REAL) + !endif endif enddo @@ -370,6 +448,28 @@ subroutine prepare_gravity() allocate(gravity_rho_crust_mantle(NGLOB_CRUST_MANTLE),stat=ier) if (ier /= 0) stop 'Error allocating gravity rho array for crust/mantle' gravity_rho_crust_mantle(:) = 0._CUSTOM_REAL + + ! gravity kernels need g values and grad(g) + if (SIMULATION_TYPE == 3) then + ! for crust/mantle kernels + allocate(g_cm(NDIM,NGLOB_CRUST_MANTLE),stat=ier) + if (ier /= 0) stop 'Error allocating gravity g_cm array' + g_cm(:,:) = 0._CUSTOM_REAL + allocate(gradg_cm(6,NGLOB_CRUST_MANTLE),stat=ier) + if (ier /= 0) stop 'Error allocating gravity gradg_cm array' + gradg_cm(:,:) = 0._CUSTOM_REAL + ! for outer core kernels - not implemented yet + !allocate(g_oc(NDIM,NGLOB_OUTER_CORE),stat=ier) + !if (ier /= 0) stop 'Error allocating gravity g_oc array' + !g_oc(:,:) = 0._CUSTOM_REAL + ! for inner core kernels - not implemented yet + !allocate(g_ic(NDIM,NGLOB_INNER_CORE),stat=ier) + !if (ier /= 0) stop 'Error allocating gravity g_ic array' + !g_ic(:,:) = 0._CUSTOM_REAL + !allocate(gradg_ic(6,NGLOB_INNER_CORE),stat=ier) + !if (ier /= 0) stop 'Error allocating gravity grad_ic array' + !gradg_ic(:,:) = 0._CUSTOM_REAL + endif endif do iglob = 1,NGLOB_CRUST_MANTLE @@ -444,6 +544,41 @@ subroutine prepare_gravity() ! for full gravity contribution in compute forces if (FULL_GRAVITY_VAL) then gravity_rho_crust_mantle(iglob) = real(rho,kind=CUSTOM_REAL) + + ! gravity kernels need g values and grad(g) + if (SIMULATION_TYPE == 3) then + ! note: this setup might not be needed if one takes gxl,.. and Hxxl,.. directly + ! and store them as g_cm and gradg_cm arrays + ! + ! use pre-computed g vector (for consistency) + !g_cm(:,iglob) = gravity_pre_store_crust_mantle(:,iglob)/rho + !gradg(:) = gravity_H_crust_mantle(:,iglob)/rho + ! + ! or re-compute g & grad(g) + ! + ! get g in crust/mantle + g = g_rad(int_radius) + if (ELLIPTICITY_VAL) then + ! w/ ellipicity + dotrho = dotrho_rad(int_radius) + eps = eps_rad(int_radius) + eta = eta_rad(int_radius) + ! compute g + call compute_g_gradg_elliptical(NDIM,radius,theta,phi,rho,dotrho,g,eps,eta,twothirdOmega2,gvec,gradg) + else + ! no ellipticity + ! compute g + call compute_g_gradg(NDIM,radius,theta,phi,rho,g,gvec,gradg) + endif + ! store g vector & grad(g) + g_cm(:,iglob) = real(gvec(:),kind=CUSTOM_REAL) + gradg_cm(:,iglob) = real(gradg(:),kind=CUSTOM_REAL) + + !debug - compare gvec with gxl,gyl,gzl divided by rho + !print *,'debug: gvec ',gvec(:),'gxl/gyl/gzl',gravity_pre_store_crust_mantle(:,iglob)/rho,'rho',rho + !debug - compare gradg with Hxxl,.. divided by rho + !print *,'debug: gradg ',gradg(:),'Hxxl/..',gravity_H_crust_mantle(:,iglob)/rho,'rho',rho + endif endif enddo else @@ -548,6 +683,29 @@ subroutine prepare_gravity() ! for full gravity contribution in compute forces if (FULL_GRAVITY_VAL) then gravity_rho_inner_core(iglob) = real(rho,kind=CUSTOM_REAL) + + ! gravity kernels need g values and grad(g) + ! inner core kernels not implemented yet... + !if (SIMULATION_TYPE == 3) then + ! ! get g in inner core + ! g = g_rad(int_radius) + ! rho = rho_rad(int_radius) + ! if (ELLIPTICITY_VAL) then + ! ! w/ ellipicity + ! dotrho = dotrho_rad(int_radius) + ! eps = eps_rad(int_radius) + ! eta = eta_rad(int_radius) + ! ! compute g + ! call compute_g_gradg_elliptical(NDIM,radius,theta,phi,rho,dotrho,g,eps,eta,twothirdOmega2,gvec,gradg) + ! else + ! ! no ellipticity + ! ! compute g + ! call compute_g_gradg(NDIM,radius,theta,phi,rho,g,gvec,gradg) + ! endif + ! ! store g vector & grad(g) + ! g_ic(:,iglob) = real(gvec(:),kind=CUSTOM_REAL) + ! gradg_ic(:,iglob) = real(gradg(:),kind=CUSTOM_REAL) + !endif endif enddo else diff --git a/src/specfem3D/read_forward_arrays.F90 b/src/specfem3D/read_forward_arrays.F90 index 857e6903b..763fae92e 100644 --- a/src/specfem3D/read_forward_arrays.F90 +++ b/src/specfem3D/read_forward_arrays.F90 @@ -33,12 +33,15 @@ subroutine read_forward_arrays_startrun() use specfem_par_crustmantle use specfem_par_innercore use specfem_par_outercore + use specfem_par_full_gravity implicit none ! local parameters integer :: ier character(len=MAX_STRING_LEN) outputname + ! full gravity + integer :: neq_read,neq1_read ! checks if anything to do ! undoing attenuation doesn't support the following checkpointing @@ -88,32 +91,49 @@ subroutine read_forward_arrays_startrun() read(IIN) epsilondev_yz_inner_core ! rotation - read(IIN) A_array_rotation - read(IIN) B_array_rotation + if (ROTATION_VAL) then + read(IIN) A_array_rotation + read(IIN) B_array_rotation + endif ! attenuation memory variables - read(IIN) R_xx_crust_mantle - read(IIN) R_yy_crust_mantle - read(IIN) R_xy_crust_mantle - read(IIN) R_xz_crust_mantle - read(IIN) R_yz_crust_mantle - - read(IIN) R_xx_inner_core - read(IIN) R_yy_inner_core - read(IIN) R_xy_inner_core - read(IIN) R_xz_inner_core - read(IIN) R_yz_inner_core + if (ATTENUATION_VAL) then + read(IIN) R_xx_crust_mantle + read(IIN) R_yy_crust_mantle + read(IIN) R_xy_crust_mantle + read(IIN) R_xz_crust_mantle + read(IIN) R_yz_crust_mantle + + read(IIN) R_xx_inner_core + read(IIN) R_yy_inner_core + read(IIN) R_xy_inner_core + read(IIN) R_xz_inner_core + read(IIN) R_yz_inner_core + endif ! full gravity if (FULL_GRAVITY_VAL) then - !TODO: implement full gravity adios read forward - stop 'FULL_GRAVITY read_forward_arrays_startrun() not fully implemented yet' - !read(55) neq1 - !allocate(pgrav1_oldrun(0:neq1)) - !read(55) pgrav1_oldrun + read(IIN) neq_read + read(IIN) neq1_read + ! check if array sizes match + if (neq_read /= neq) then + print *,'Error reading forward array for startrun: rank ',myrank,'has read neq =',neq_read,' - shoud be ',neq + call exit_MPI(myrank,'Invalid forward array neq for startrun') + endif + if (neq1_read /= neq1) then + print *,'Error reading forward array for startrun: rank ',myrank,'has read neq1 =',neq1_read,' - shoud be ',neq1 + call exit_MPI(myrank,'Invalid forward array neq1 for startrun') + endif + read(IIN) pgrav1 endif close(IIN) + + ! full gravity + if (FULL_GRAVITY) then + ! need to interpolate the gravity values + call SIEM_interpolate_gravity() + endif endif endif @@ -131,11 +151,13 @@ subroutine read_forward_arrays() use specfem_par_crustmantle use specfem_par_innercore use specfem_par_outercore + use specfem_par_full_gravity implicit none ! local parameters integer :: ier + integer :: b_neq_read,b_neq1_read character(len=MAX_STRING_LEN) :: outputname ! checks if anything to do @@ -199,13 +221,21 @@ subroutine read_forward_arrays() read(IIN) b_R_yz_inner_core endif + ! full gravity if (FULL_GRAVITY_VAL) then - ! Read in the gravity - !TODO: implement full gravity adios read forward - stop 'FULL_GRAVITY read_forward_arrays() not fully implemented yet' - !read(55) b_neq - !read(55) b_neq1 - !read(55) b_pgrav1 + read(IIN) b_neq_read + read(IIN) b_neq1_read + ! check if array sizes match + if (b_neq_read /= neq) then + print *,'Error reading forward array: rank ',myrank,'has read neq =',b_neq_read,' - shoud be ',neq + call exit_MPI(myrank,'Invalid forward array neq') + endif + if (b_neq1_read /= neq1) then + print *,'Error reading forward array: rank ',myrank,'has read neq1 =',b_neq1_read,' - shoud be ',neq1 + call exit_MPI(myrank,'Invalid forward array neq1') + endif + ! Level-1 solver array + read(IIN) b_pgrav1 endif close(IIN) @@ -214,7 +244,7 @@ subroutine read_forward_arrays() ! full gravity if (FULL_GRAVITY) then ! need to interpolate the gravity values - call SIEM_interpolate_gravity() + call SIEM_interpolate_backward_gravity() endif ! transfers fields onto GPU @@ -271,12 +301,14 @@ subroutine read_forward_arrays_undoatt() use specfem_par_crustmantle use specfem_par_innercore use specfem_par_outercore + use specfem_par_full_gravity implicit none ! local parameters integer :: iteration_on_subset_tmp integer :: ier + integer :: b_neq_read,b_neq1_read character(len=MAX_STRING_LEN) :: outputname ! current subset iteration @@ -327,9 +359,32 @@ subroutine read_forward_arrays_undoatt() read(IIN) b_R_yz_inner_core endif + ! full gravity + if (FULL_GRAVITY_VAL) then + read(IIN) b_neq_read + read(IIN) b_neq1_read + ! check if array sizes match + if (b_neq_read /= neq) then + print *,'Error reading forward array: rank ',myrank,'has read neq =',b_neq_read,' - shoud be ',neq + call exit_MPI(myrank,'Invalid forward array neq') + endif + if (b_neq1_read /= neq1) then + print *,'Error reading forward array: rank ',myrank,'has read neq1 =',b_neq1_read,' - shoud be ',neq1 + call exit_MPI(myrank,'Invalid forward array neq1') + endif + ! Level-1 solver array + read(IIN) b_pgrav1 + endif + close(IIN) endif + ! full gravity + if (FULL_GRAVITY) then + ! need to interpolate the gravity values + call SIEM_interpolate_backward_gravity() + endif + ! transfers fields onto GPU if (GPU_MODE) then ! transfers initialized wavefields to GPU device diff --git a/src/specfem3D/read_forward_arrays_adios.F90 b/src/specfem3D/read_forward_arrays_adios.F90 index c15048c52..4c752d489 100644 --- a/src/specfem3D/read_forward_arrays_adios.F90 +++ b/src/specfem3D/read_forward_arrays_adios.F90 @@ -191,9 +191,9 @@ subroutine read_intermediate_forward_arrays_adios() if (FULL_GRAVITY_VAL) then !TODO: implement full gravity adios read forward stop 'FULL_GRAVITY read_intermediate_forward_arrays_adios() not fully implemented yet' - !read(55) neq1 + !read(IIN) neq1 !allocate(pgrav1_oldrun(0:neq1)) - !read(55) pgrav1_oldrun + !read(IIN) pgrav1_oldrun endif ! closes ADIOS handler to the restart file. @@ -357,9 +357,9 @@ subroutine read_forward_arrays_adios() ! Read in the gravity !TODO: implement full gravity adios read forward stop 'FULL_GRAVITY read_forward_arrays_adios() not fully implemented yet' - !read(55) b_neq - !read(55) b_neq1 - !read(55) b_pgrav1 + !read(IIN) b_neq_read + !read(IIN) b_neq1_read + !read(IIN) b_pgrav1 endif ! closes ADIOS handler to the restart file. @@ -582,6 +582,15 @@ subroutine read_forward_arrays_undoatt_adios(iteration_on_subset_tmp) call delete_adios_selection(sel) enddo + if (FULL_GRAVITY_VAL) then + ! Read in the gravity + !TODO: implement full gravity adios read forward + stop 'FULL_GRAVITY read_forward_arrays_undoatt_adios() not fully implemented yet' + !read(IIN) b_neq_read + !read(IIN) b_neq1_read + !read(IIN) b_pgrav1 + endif + ! Close ADIOS handler to the restart file. ! note: for single file, only close at the very end if (do_close_file) then diff --git a/src/specfem3D/save_forward_arrays.F90 b/src/specfem3D/save_forward_arrays.F90 index 31d939624..2f44a4f52 100644 --- a/src/specfem3D/save_forward_arrays.F90 +++ b/src/specfem3D/save_forward_arrays.F90 @@ -40,6 +40,7 @@ subroutine save_forward_arrays() use specfem_par_crustmantle use specfem_par_innercore use specfem_par_outercore + use specfem_par_full_gravity implicit none @@ -48,7 +49,7 @@ subroutine save_forward_arrays() character(len=MAX_STRING_LEN) outputname ! checks if anything to do - if (UNDO_ATTENUATION ) return + if (UNDO_ATTENUATION) return ! save files to local disk or tape system if restart file if (NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS) then @@ -98,28 +99,32 @@ subroutine save_forward_arrays() write(IOUT) epsilondev_yz_inner_core ! rotation - write(IOUT) A_array_rotation - write(IOUT) B_array_rotation + if (ROTATION_VAL) then + write(IOUT) A_array_rotation + write(IOUT) B_array_rotation + endif ! attenuation memory variables - write(IOUT) R_xx_crust_mantle - write(IOUT) R_yy_crust_mantle - write(IOUT) R_xy_crust_mantle - write(IOUT) R_xz_crust_mantle - write(IOUT) R_yz_crust_mantle - - write(IOUT) R_xx_inner_core - write(IOUT) R_yy_inner_core - write(IOUT) R_xy_inner_core - write(IOUT) R_xz_inner_core - write(IOUT) R_yz_inner_core + if (ATTENUATION_VAL) then + write(IOUT) R_xx_crust_mantle + write(IOUT) R_yy_crust_mantle + write(IOUT) R_xy_crust_mantle + write(IOUT) R_xz_crust_mantle + write(IOUT) R_yz_crust_mantle + + write(IOUT) R_xx_inner_core + write(IOUT) R_yy_inner_core + write(IOUT) R_xy_inner_core + write(IOUT) R_xz_inner_core + write(IOUT) R_yz_inner_core + endif ! full gravity if (FULL_GRAVITY_VAL) then - !TODO: implement full gravity adios save forward - stop 'FULL_GRAVITY save_forward_arrays() not fully implemented yet' - !write(IOUT) neq1 - !write(IOUT pgrav1 + write(IOUT) neq + write(IOUT) neq1 + ! Level-1 solver array + write(IOUT) pgrav1 endif close(IOUT) @@ -188,11 +193,10 @@ subroutine save_forward_arrays() ! full gravity if (FULL_GRAVITY_VAL) then - !TODO: implement full gravity adios save forward - stop 'FULL_GRAVITY save_forward_arrays() not fully implemented yet' - !write(IOUT) neq - !write(IOUT) neq1 - !write(IOUT pgrav1 + write(IOUT) neq + write(IOUT) neq1 + ! Level-1 solver array + write(IOUT) pgrav1 endif close(IOUT) @@ -211,6 +215,7 @@ subroutine save_forward_arrays_undoatt() use specfem_par_crustmantle use specfem_par_innercore use specfem_par_outercore + use specfem_par_full_gravity implicit none @@ -289,11 +294,10 @@ subroutine save_forward_arrays_undoatt() ! full gravity if (FULL_GRAVITY_VAL) then - !TODO: implement full gravity adios save forward - stop 'FULL_GRAVITY save_forward_arrays_undoatt() not fully implemented yet' - !write(IOUT) neq - !write(IOUT) neq1 - !write(IOUT pgrav1 + write(IOUT) neq + write(IOUT) neq1 + ! Level-1 solver array + write(IOUT) pgrav1 endif close(IOUT) diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90 index b03b3679b..7e1f72395 100644 --- a/src/specfem3D/specfem3D_par.F90 +++ b/src/specfem3D/specfem3D_par.F90 @@ -1230,11 +1230,11 @@ module specfem_par_full_gravity ! Level-1 solver arrays real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_pgrav_ic1, b_pgrav_oc1, b_pgrav_cm1, b_pgrav_trinf1, b_pgrav_inf1 real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_pgrav1 - real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_dprecon1, b_gravload1 + real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_gravload1 ! Level-2 solver arrays real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_pgrav_ic, b_pgrav_oc, b_pgrav_cm, b_pgrav_trinf, b_pgrav_inf real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_pgrav - real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_dprecon, b_gravload + real(kind=CUSTOM_REAL),dimension(:),allocatable :: b_gravload ! Level-2 solver ! number of global degrees of freedom