updates SIEM routines (avoids division by zero in CG solver; adds fin…

…alize; updates array allocations; renames parameters, routines and files)

DATA/Par_file

@@ -76,6 +76,9 @@ ATTENUATION                     = .true.
 # full gravity calculation by solving Poisson's equation for gravity potential instead of using a Cowling approximation
 # (must have also GRAVITY flag set to .true. to become active)
 FULL_GRAVITY                    = .false.
+# for full gravity calculation, set to 0 == builtin or 1 == PETSc Poisson solver
+# (the PETSc solver option needs the PETSc library installed; code configuration --with-petsc)
+POISSON_SOLVER                  = 0
 
 # record length in minutes
 RECORD_LENGTH_IN_MINUTES        = 2.5d0

Makefile.in

@@ -551,9 +551,8 @@ SUBDIRS = \
 	tomography/postprocess_sensitivity_kernels \
 	$(EMPTY_MACRO)
 
-ifeq ($(PETSC), yes)
-	SUBDIRS += gindex3D
-endif
+# full gravity tool
+SUBDIRS += gindex3D
 
 #ifeq ($(HAS_GPU),yes)
 #  SUBDIRS := gpu $(SUBDIRS)
@@ -592,11 +591,8 @@ endif
 #	$(EMPTY_MACRO)
 #endif
 
-ifeq ($(PETSC), yes)
-DEFAULT += 	\
-	xgindex3D \
-	$(EMPTY_MACRO)
-endif
+# full gravity tool
+DEFAULT += 	xgindex3D
 
 
 all: default aux movies postprocess tomography
@@ -713,11 +709,9 @@ endif
 	@echo "for unit testing:"
 	@echo "    tests"
 	@echo ""
-ifeq ($(PETSC), yes)
 	@echo "for full-gravity simulations:"
 	@echo "    xgindex3D"
 	@echo ""
-endif
 
 .PHONY: all default clean realclean help tests
 

src/shared/broadcast_computed_parameters.f90

@@ -34,7 +34,7 @@ subroutine broadcast_computed_parameters()
 
   ! local parameters
   ! broadcast parameter arrays
-  integer, parameter :: nparam_i = 49
+  integer, parameter :: nparam_i = 50
   integer, dimension(nparam_i) :: bcast_integer
 
   integer, parameter :: nparam_l = 74
@@ -76,7 +76,8 @@ subroutine broadcast_computed_parameters()
             ATT1,ATT2,ATT3,ATT4,ATT5, &
             GPU_RUNTIME,NUMBER_OF_SIMULTANEOUS_RUNS, &
             MODEL_GLL_TYPE,USER_NSTEP, &
-            NSTEP_STEADY_STATE,NTSTEP_BETWEEN_OUTPUT_SAMPLE /)
+            NSTEP_STEADY_STATE,NTSTEP_BETWEEN_OUTPUT_SAMPLE, &
+            POISSON_SOLVER /)
 
     bcast_logical = (/ &
             TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
@@ -272,6 +273,7 @@ subroutine broadcast_computed_parameters()
     USER_NSTEP = bcast_integer(47)
     NSTEP_STEADY_STATE = bcast_integer(48)
     NTSTEP_BETWEEN_OUTPUT_SAMPLE = bcast_integer(49)
+    POISSON_SOLVER = bcast_integer(50)
 
     ! logicals
     TRANSVERSE_ISOTROPY = bcast_logical(1)

src/shared/read_compute_parameters.f90

@@ -628,6 +628,10 @@ subroutine rcp_SIEM_set_mesh_parameters()
   ! safety check
   if (.not. FULL_GRAVITY) return
 
+  ! checks solver setting
+  if (POISSON_SOLVER /= 0 .and. POISSON_SOLVER /= 1) &
+    stop 'For FULL_GRAVITY calculations, POISSON_SOLVER must be set to either 0 == builtin or 1 == PETSc'
+
   ! start region
   iregion0 = IREGION_CRUST_MANTLE
 

src/shared/read_parameter_file.F90

@@ -102,6 +102,10 @@ subroutine read_parameter_file()
   ! full gravity support
   call read_value_logical(FULL_GRAVITY, 'FULL_GRAVITY', ier)
   if (ier /= 0) stop 'an error occurred while reading the parameter file: FULL_GRAVITY'
+  if (FULL_GRAVITY) then
+    call read_value_integer(POISSON_SOLVER, 'POISSON_SOLVER', ier)
+    if (ier /= 0) stop 'an error occurred while reading the parameter file: POISSON_SOLVER'
+  endif
 
   call read_value_double_precision(RECORD_LENGTH_IN_MINUTES, 'RECORD_LENGTH_IN_MINUTES', ier)
   if (ier /= 0) stop 'an error occurred while reading the parameter file: RECORD_LENGTH_IN_MINUTES'

src/shared/shared_par.f90

@@ -102,6 +102,7 @@ module shared_input_parameters
 
   ! full gravity support
   logical :: FULL_GRAVITY = .false.
+  integer :: POISSON_SOLVER = 0   ! 0 == builtin / 1 == PETSc solver
 
   ! regional mesh cut-off
   logical :: REGIONAL_MESH_CUTOFF

src/specfem3D/SIEM_compute_kernels.F90

@@ -74,7 +74,7 @@ subroutine calculate_first_gravity_kernel()
     is_active_gll,igll_active_on
 
   use siem_math_library_mpi, only: maxvec
-  use siem_poisson, only: compute_grav_kl1_load, compute_poisson_load3
+  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
 
@@ -88,7 +88,7 @@ subroutine calculate_first_gravity_kernel()
   real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: gknl1_cm ! (3,NGLOB_CRUST_MANTLE)
 
   ! safety check
-  if (POISSON_SOLVER == BUILTIN) then
+  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')
@@ -111,18 +111,15 @@ subroutine calculate_first_gravity_kernel()
     if (myrank == 0) print *,'icomponent: ', icomponent
 
     ! Calculate the RHS (gravload1)
-    gravload1(:) = 0.0_CUSTOM_REAL
-
     call compute_grav_kl1_load(icomponent)
 
     !debug
     maxload = maxvec(abs(gravload1))
     if (myrank == 0) print *,'  -- Max load: ', maxload
 
-    if (POISSON_SOLVER == PETSC) then
+    if (POISSON_SOLVER == ISOLVER_PETSC) then
       ! Petsc solver
       call petsc_set_vector1(gravload1)
-      call synchronize_all()
 
       ! Need to zero guess since was previously for poisson eqn
       ! Do we need this - is zero guess set to PETSC_TRUE?
@@ -229,7 +226,7 @@ subroutine calculate_second_gravity_kernel()
   real(kind=CUSTOM_REAL),dimension(:,:,:,:,:), allocatable :: div_gknl2_cm
 
   ! safety check
-  if (POISSON_SOLVER == BUILTIN) then
+  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')
@@ -250,22 +247,19 @@ subroutine calculate_second_gravity_kernel()
   ! Solve the equation 9 times!
   do icomponent = 1,3
     do jcomponent = 1,3
-
-      gravload1(:) = 0.0_CUSTOM_REAL
-
       !debug
       if (myrank == 0) print *,'component ', icomponent, jcomponent
 
+      ! Calculate the RHS (gravload1)
       call compute_grav_kl2_load(icomponent, jcomponent)
 
-      if (POISSON_SOLVER == PETSC) then
+      !debug
+      maxload = maxvec(abs(gravload1))
+      if (myrank == 0) print *,'  -- Max load: ', maxload
+
+      if (POISSON_SOLVER == ISOLVER_PETSC) then
         ! PETSc solver
         call petsc_set_vector1(gravload1)
-        call synchronize_all()
-
-        !debug
-        maxload = maxvec(abs(gravload1))
-        if (myrank == 0) print *,'  -- Max load: ', maxload
 
         call petsc_zero_initialguess1()
 

src/specfem3D/SIEM_index_region.F90

@@ -26,7 +26,7 @@
 !=====================================================================
 
 
-  subroutine get_index_region()
+  subroutine SIEM_get_index_region()
 
   use constants, only: myrank,NGLLX,NGLLY,NGLLZ,NGLLCUBE,IFLAG_IN_FICTITIOUS_CUBE, &
     NGLLX_INF,NGLLY_INF,NGLLZ_INF,NGLLCUBE_INF,ADD_TRINF
@@ -1093,5 +1093,5 @@ subroutine get_index_region()
   deallocate(inode_ic,inode_oc,inode_cm,inode_trinf,inode_inf)
   deallocate(inode_ic1,inode_oc1,inode_cm1,inode_trinf1,inode_inf1)
 
-  end subroutine get_index_region
+  end subroutine SIEM_get_index_region
 

src/specfem3D/SIEM_infinite_element.F90

@@ -64,26 +64,28 @@ subroutine shape_function_infiniteGLHEX8ZW_GLLR(ndim,ngllx,nglly,ngllz,ngll,nip,
   real(kind=kdble),dimension(ngllz) :: gllpz,gllwz,igllpz,igllwz ! GLL points and weights
   real(kind=kdble),dimension(ndim,ngllx) :: lagrange_x,lagrange_dx
   real(kind=kdble),dimension(ndim,2) :: lagrangeINF_x,lagrangeINF_dx
-  integer :: iINF !,ngllxINF(ndim)
+  integer :: iINF
 
+  iINF = 0
   ddir = one
+
   if (iface == 1) then
-    iINF = 2; ddir=-one
+    iINF = 2; ddir = -one
   else if (iface == 2) then
     iINF = 1; ddir = one
   else if (iface == 3) then
     iINF = 2; ddir = one
   else if (iface == 4) then
-    iINF = 1; ddir=-one
+    iINF = 1; ddir = -one
   else if (iface == 5) then
-    iINF = 3; ddir=-one
+    iINF = 3; ddir = -one
   else if (iface == 6) then
     iINF = 3; ddir = one
   endif
 
-  nipx(1)=ngllx
-  nipx(2)=nglly
-  nipx(3)=ngllz
+  nipx(1) = ngllx
+  nipx(2) = nglly
+  nipx(3) = ngllz
 
   ! compute everything in indexed order
   ! get GLL points
@@ -195,19 +197,20 @@ end subroutine shape_function_infiniteGLHEX8ZW_GLLR
 !  real(kind=kdble),dimension(nipx) :: igllpz,igllwz ! GLL points and weights
 !  real(kind=kdble),dimension(ndim,ngllx) :: lagrange_x,lagrange_dx
 !  real(kind=kdble),dimension(ndim,2) :: lagrangeINF_x,lagrangeINF_dx
-!  integer :: iINF !,ngllxINF(ndim)
+!  integer :: iINF
 !
+!  iINF = 0
 !  ddir = one
 !  if (iface == 1) then
-!    iINF = 2; ddir=-one
+!    iINF = 2; ddir = -one
 !  else if (iface == 2) then
 !    iINF = 1; ddir = one
 !  else if (iface == 3) then
 !    iINF = 2; ddir = one
 !  else if (iface == 4) then
-!    iINF = 1; ddir=-one
+!    iINF = 1; ddir = -one
 !  else if (iface == 5) then
-!    iINF = 3; ddir=-one
+!    iINF = 3; ddir = -one
 !  else if (iface == 6) then
 !    iINF = 3; ddir = one
 !  endif
@@ -319,16 +322,19 @@ subroutine get_gnodinfHEX8(ndim,ngllx,nglly,ngllz,nginf,iface,gnodinf)
   ngllxINF(2) = nglly
   ngllxINF(3) = ngllz
 
+  iINF = 0
+  ddir = one
+
   if (iface == 1) then
-    iINF = 2; ddir=-one
+    iINF = 2; ddir = -one
   else if (iface == 2) then
     iINF = 1; ddir = one
   else if (iface == 3) then
     iINF = 2; ddir = one
   else if (iface == 4) then
-    iINF = 1; ddir=-one
+    iINF = 1; ddir = -one
   else if (iface == 5) then
-    iINF = 3; ddir=-one
+    iINF = 3; ddir = -one
   else if (iface == 6) then
     iINF = 3; ddir = one
   endif
@@ -346,7 +352,7 @@ subroutine get_gnodinfHEX8(ndim,ngllx,nglly,ngllz,nginf,iface,gnodinf)
     do j = ngllxINF0(2),ngllxINF(2),inc(2)
       do i = ngllxINF0(1),ngllxINF(1),inc(1)
         inum = inum+1
-        gnodinf(inum) = nglly*ngllx*(k-1)+ngllx*(j-1)+i
+        gnodinf(inum) = nglly*ngllx*(k-1) + ngllx*(j-1)+i
       enddo
     enddo
   enddo