updates SIEM gravity setup; adds read/save forward arrays; adds backw…

…ard field solve for builtin solver

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
 

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

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---------------