Updates to geoclaw bouss calling of PETSc to support shared memory

src/2d/bouss/amr2.f90

@@ -127,18 +127,18 @@ program amr2
 
     ! Local variables
     integer :: i, iaux, mw, level
-    integer :: ierr
+    integer :: ierr ! error flag for PETSc calls
     integer :: ndim, nvar, naux, mcapa1, mindim, dimensional_split
     integer :: nstart, nsteps, nv1, nx, ny, lentotsave, num_gauge_SAVE
     integer :: omp_get_max_threads, maxthreads
     real(kind=8) :: time, ratmet, cut, dtinit, dt_max
-    logical :: vtime, rest, output_t0    
+    logical :: vtime, rest, output_t0
 
     ! Timing variables
     integer(kind=8) :: clock_start, clock_finish, clock_rate, ttotal, count_max
     real(kind=8) :: ttotalcpu
     integer(kind=8) :: tick_clock_finish
-    integer, parameter :: timing_unit = 48
+    integer :: timing_unit = 0
     character(len=512) :: timing_line, timing_substr
     character(len=*), parameter :: timing_base_name = "timing."
     character(len=*), parameter :: timing_header_format =                      &
@@ -160,10 +160,12 @@ program amr2
     character(len=*), parameter :: parmfile = 'fort.parameters'
 
 #ifdef HAVE_PETSC
-!   needs to be very first line in main
-      call PetscInitialize(ierr)
-      CHKERRA(ierr);
+    !   needs to be very first line in main
+    PetscCallA(PetscOptionsSetValue(PETSC_NULL_OPTIONS,'-mpi_linear_solver_server','',ierr))
+    PetscCallA(PetscOptionsSetValue(PETSC_NULL_OPTIONS,'-mpi_linear_solver_server_view','',ierr))
+    PetscCallA(PetscInitialize(ierr))
 #endif
+    timing_unit = 48
 
     ! Open parameter and debug files
     open(dbugunit,file=dbugfile,status='unknown',form='formatted')
@@ -599,7 +601,7 @@ program amr2
         time = t0
         nstart = 0
     endif
-
+    call PetscViewerASCIIStdoutSetFileUnit(outunit,ierr)
 
     write(parmunit,*) ' '
     write(parmunit,*) '--------------------------------------------'
@@ -860,7 +862,6 @@ program amr2
     write(*,format_string)
     write(*,*)
     write(timing_unit,*)
-    close(timing_unit)
 
     if (isolver.eq. 2) then
 !     if (countIterRef == 0) then
@@ -928,12 +929,16 @@ program amr2
 
     write(outunit,"(//,' ------  end of AMRCLAW integration --------  ')")
 
-    ! Close output and debug files.
-    close(outunit)
+    ! Close debug file.
     close(dbugunit)
+    close(outunit)
 
 #ifdef HAVE_PETSC
-      call PetscFinalize(ierr)
+    call PetscViewerASCIIStdoutSetFileUnit(timing_unit,ierr)
+    call PetscFinalize(ierr)
 #endif
 
-end program amr2
+    if (timing_unit .gt. 0) then
+       close(timing_unit)
+    endif
+    end program amr2

src/2d/bouss/bouss_module.f90

@@ -49,8 +49,8 @@ module bouss_module
        integer :: numColsTot
 
        ! for CRS format
-       integer, allocatable, dimension(:) :: rowPtr, cols
-       real(kind=8), allocatable, dimension(:) :: vals 
+       integer, pointer :: rowPtr(:), cols(:)
+       real(kind=8), pointer :: vals(:)
 
        ! intfCountc are the # equations (cells)  added for you as a coarse grid
        ! these are the equations under a fine grid on the first interior cells
@@ -166,7 +166,7 @@ subroutine set_bouss(rest,time,naux)
         crs = .true.
     endif
 
-    ! crs = .false.  ! uncomment this line to force CRS with SGN for testing
+    !crs = .false.  ! uncomment this line to force CRS with SGN for testing
 
     !------------------------------------------
     if (rest) then

src/2d/bouss/buildSparseMatrixMScoo.f90

@@ -11,7 +11,7 @@ subroutine buildSparseMatrixMScoo(rhs,nvar,naux,levelBouss,numBoussCells)
     integer, intent(in) :: nvar, naux, levelBouss, numBoussCells
 
     !! pass in numBoussCells for this level so can dimension this array
-    real(kind=8) :: rhs(0:2*numBoussCells) 
+    real(kind=8) :: rhs(2*numBoussCells) 
 
     type(matrix_patchIndex), pointer :: mi
     type(matrix_levInfo),  pointer :: minfo

src/2d/bouss/buildSparseMatrixSGNcoo.f90

@@ -16,7 +16,7 @@ subroutine buildSparseMatrixSGNcoo(q,qold,aux,soln,rhs,                    &
     integer, intent(in) :: nvar, naux, levelBouss, numBoussCells
 
     !! pass in numBoussCells for this level so can dimension these arrays
-    real(kind=8) :: soln(0:2*numBoussCells), rhs(0:2*numBoussCells) 
+    real(kind=8) :: soln(2*numBoussCells), rhs(2*numBoussCells) 
 
     type(matrix_patchIndex), pointer :: mi
     type(matrix_levInfo),  pointer :: minfo

src/2d/bouss/buildSparseMatrixSGNcrs.f90

@@ -17,7 +17,7 @@ subroutine buildSparseMatrixSGNcrs(q,qold,aux,soln,rhs,rowPtr,cols,vals,   &
     integer, intent(inout) :: cols(0:24*numBoussCells)
 
     !! pass in numBoussCells for this level so can dimension these arrays
-    real(kind=8) :: soln(0:2*numBoussCells), rhs(0:2*numBoussCells) 
+    real(kind=8) :: soln(0:2*numBoussCells-1), rhs(0:2*numBoussCells-1) 
 
     type(matrix_patchIndex), pointer :: mi
     type(matrix_levInfo),  pointer :: minfo

src/2d/bouss/implicit_update_bouss_2Calls.f90

@@ -12,15 +12,20 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
     use amr_module
     use topo_module, only: topo_finalized
     use bouss_module
-
+#include <petsc/finclude/petscvec.h>
+    use petscvec
+
     implicit none
 
     integer, intent(in) :: nvar, naux, levelBouss, numBoussCells
     logical, intent(in) :: doUpdate
 
     !! pass in numBoussCells for this level so can dimension this array
-    real(kind=8) :: soln(0:2*numBoussCells)
-    real(kind=8) :: rhs(0:2*numBoussCells)
+    Vec :: v_soln,v_rhs
+    integer :: ierr
+
+    real(kind=8),pointer :: soln(:)
+    real(kind=8),pointer :: rhs(:)
 
     type(matrix_patchIndex), pointer :: mi
     type(matrix_levInfo),  pointer :: minfo
@@ -38,6 +43,11 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
 #ifdef WHERE_AM_I
     write(*,*) "starting implicit_update for level ",levelBouss
 #endif
+    PetscCallA(VecCreateSeq(PETSC_COMM_SELF,2*numBoussCells,v_rhs,ierr))
+    PetscCallA(VecCreateSeq(PETSC_COMM_SELF,2*numBoussCells,v_soln,ierr))
+    PetscCallA(VecGetArrayF90(v_rhs,rhs,ierr))
+    PetscCallA(VecGetArrayF90(v_soln,soln,ierr))
+
     dt = possk(levelBouss)
     nD = numBoussCells ! shorthand for size of matrix blocks D1 to D4
 
@@ -102,7 +112,8 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
       else ! CRS format
          minfo%vals = 0
          minfo%cols = -1  ! to recognize if not overwritten with real col indices
-        call prepBuildSparseMatrixSGNcrs(soln,rhs,nvar,naux,levelBouss,numBoussCells)
+         call prepBuildSparseMatrixSGNcrs(soln,rhs,nvar,naux,levelBouss,numBoussCells)
+         CHKMEMQ
       endif
     endif
 
@@ -141,7 +152,12 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
        if (minfo%numColsTot .gt. 0 .or. minfo%matrix_nelt .gt. 0) then
           call system_clock(clock_startLinSolve,clock_rate)
           call cpu_time(cpu_startLinSolve)
-          call petsc_driver(soln,rhs,levelBouss,numBoussCells,time,topo_finalized)
+          PetscCallA(VecRestoreArrayF90(v_rhs,rhs,ierr))
+          PetscCallA(VecRestoreArrayF90(v_soln,soln,ierr))
+          CHKMEMQ;
+          call petsc_driver(v_soln,v_rhs,levelBouss,numBoussCells,time,topo_finalized)
+          PetscCallA(VecGetArrayF90(v_rhs,rhs,ierr))
+          PetscCallA(VecGetArrayF90(v_soln,soln,ierr))
           call system_clock(clock_finishLinSolve,clock_rate)
           call cpu_time(cpu_finishLinSolve)
           !write(89,*)" level ",levelBouss,"  rhs   soln"
@@ -201,6 +217,7 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
         !endif
     end do        
 !$OMP END PARALLEL DO
+    CHKMEMQ
 
     !! next loop to copy from either other grids at same level and to
     !! update ghost cell momenta with newly copied psi
@@ -238,6 +255,10 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time)
     write(*,*) "ending   implicit_update for level ",levelBouss
 #endif
 
+    PetscCallA(VecRestoreArrayF90(v_rhs,rhs,ierr))
+    PetscCallA(VecRestoreArrayF90(v_soln,soln,ierr))
+    PetscCallA(VecDestroy(v_rhs,ierr))
+    PetscCallA(VecDestroy(v_soln,ierr))
 
 contains
 
@@ -255,8 +276,6 @@ integer pure function iaddaux(iaux,i,j)
         iaddaux = locaux + iaux-1 + naux*((j-1+nghost)*mitot+i+nghost-1)
     end function iaddaux
 
-
-
 end subroutine implicit_update
 
 subroutine solnUpdate(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
@@ -279,7 +298,7 @@ subroutine solnUpdate(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
 
     integer, intent(in) :: nvar,nghost,mitot,mjtot,nD,mptr,nb,levelBouss
     logical, intent(in) :: doUpdate
-    real(kind=8), intent(in) :: soln(0:2*nD), dt, rhs(0:2*nD)
+    real(kind=8), intent(in) :: soln(2*nD), dt, rhs(2*nD)
 
     real(kind=8), intent(inout) :: valnew(nvar,mitot,mjtot)
     real(kind=8), intent(inout) :: valOther(nvar,mitot,mjtot)
@@ -415,7 +434,7 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
 
     integer, intent(in) :: nvar,naux,nghost,mitot,mjtot,nD,mptr,nb,levelBouss
     logical, intent(in) :: doUpdate
-    real(kind=8), intent(in) :: soln(0:2*nD), dt, rhs(0:2*nD)
+    real(kind=8), intent(in) :: soln(2*nD), dt, rhs(2*nD)
 
     real(kind=8), intent(inout) :: valnew(nvar,mitot,mjtot)
     real(kind=8), intent(inout) :: valOther(nvar,mitot,mjtot)
@@ -505,9 +524,9 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
                          (grav/alpha * etay(i,j) - soln(k_ij+nD))
                  else  ! CRS format,  different mapping
                   valnew(2,i,j) = valnew(2,i,j) + dt * valnew(1,i,j)*     &
-                         (grav/alpha * etax(i,j) - soln(2*(k_ij-1)))
+                         (grav/alpha * etax(i,j) - soln(2*(k_ij-1)+1))
                   valnew(3,i,j) = valnew(3,i,j) + dt * valnew(1,i,j)*     &
-                         (grav/alpha * etay(i,j) - soln(2*k_ij-1))
+                         (grav/alpha * etay(i,j) - soln(2*k_ij-1+1))
                  endif
 
                 ! save update in components 4,5 for Dirichlet BC next iteration
@@ -521,8 +540,8 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
                 else
                     !valOther(4,i,j) = soln(2*k_ij-1)
                     !valOther(5,i,j) = soln(2*k_ij)
-                    valOther(4,i,j) = soln(2*(k_ij-1))
-                    valOther(5,i,j) = soln(2*k_ij-1)
+                    valOther(4,i,j) = soln(2*(k_ij-1)+1)
+                    valOther(5,i,j) = soln(2*k_ij-1+1)
                  endif
             else  ! not a bouss cell. set to zero
                valOther(4,i,j) = 0.d0
@@ -610,8 +629,8 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, &
                 valnew(4,i,j) = soln(k_ij)
                 valnew(5,i,j) = soln(k_ij+nD)
               else
-                valnew(4,i,j) = soln(2*(k_ij-1))
-                valnew(5,i,j) = soln(2*k_ij-1)
+                valnew(4,i,j) = soln(2*(k_ij-1)+1)
+                valnew(5,i,j) = soln(2*k_ij-1+1)
               endif
             else
               valnew(4,i,j) = 0.d0

src/2d/bouss/petsc_driver.f90

@@ -1,23 +1,18 @@
-subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalized)
+subroutine petsc_driver(solution,rhs,levelBouss,numBoussCells,time,topo_finalized)
 
     ! use petsc to solve sparse linear system
     ! called this way so can allocate storage of correct size
 
     use bouss_module
 #ifdef HAVE_PETSC
 #include <petsc/finclude/petscksp.h>
-!!#include "petscmat.h"
     use petscksp
     implicit none
 
     integer, intent(in) :: levelBouss, numBoussCells
-    real(kind=8), intent(inout) :: rhs_geo(0:2*numBoussCells)
     real(kind=8), intent(in) :: time
     logical, intent(in) :: topo_finalized
 
-    !! pass in numBoussCells for this level so can dimension this array
-    real(kind=8), intent(out) :: soln(0:2*numBoussCells)
-
     type(matrix_levInfo),  pointer :: minfo
 
     integer :: numCores, omp_get_max_threads,i,nelt
@@ -87,6 +82,7 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize
             end do
           else 
             ! using CRS sparse matrix format
+            CHKMEMQ;
             call MatCreateSeqAijWithArrays(PETSC_COMM_SELF,2*numBoussCells,  & 
                        2*numBoussCells,minfo%rowPtr,minfo%cols,  &
                        minfo%vals,Jr(levelBouss),ierr)
@@ -136,17 +132,15 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize
       ! if you want to use previous soln as initial guess,
       ! set soln to old soln here
 
-      if (.not. crs) then
-        ! switch to 0-based indexing for COO triplet format
-        ! (already 0-based for CRS)
-        do i = 1, 2*numBousscells
-           rhs_geo(i-1) = rhs_geo(i) 
-        end do
-      endif 
+!      if (.not. crs) then
+!        ! switch to 0-based indexing for COO triplet format
+!        ! (already 0-based for CRS)
+!        do i = 1, 2*numBousscells
+!           rhs_geo(i-1) = rhs_geo(i) 
+!        end do
+!      endif 
 
-      call VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,rhs_geo,rhs,ierr)
-      CHKERRA(ierr)
-      call VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,soln,solution,ierr)
+
       CHKERRA(ierr)
 
       ! save matrices for debugging. comment out to turn off
@@ -167,23 +161,17 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize
       !endif
 
       ! check if need to adjust soln back to 1 indexing
-      if (.not. crs) then ! bump both  back up
-        do i = 1, 2*numBoussCells
-          soln(2*numBoussCells-i+1) = soln(2*numBoussCells-i)
-          rhs_geo(2*numBoussCells-i+1) = rhs_geo(2*numBoussCells-i) 
-        end do
-      endif
+!      if (.not. crs) then ! bump both  back up
+!        do i = 1, 2*numBoussCells
+!          soln(2*numBoussCells-i+1) = soln(2*numBoussCells-i)
+!          rhs_geo(2*numBoussCells-i+1) = rhs_geo(2*numBoussCells-i) 
+!        end do
+!      endif
 
       !! if itnum > itmax then
       !! call KSPGetResidualNorm(ksp(levelBouss),resmax,ierr)
       CHKERRA(ierr)
       ! petsc already puts solution into my array soln so just return
-
-       ! create and destroy each step since fast
-       call VecDestroy(rhs,ierr)
-       CHKERRA(ierr)
-       call VecDestroy(solution,ierr)
-       CHKERRA(ierr)
 #endif
 
 end subroutine petsc_driver

src/2d/bouss/prepBuildSparseMatrixSGNcoo.f90

@@ -11,7 +11,7 @@ subroutine prepBuildSparseMatrixSGNcoo(soln,rhs,nvar,naux,levelBouss,numBoussCel
     integer, intent(in) :: nvar, naux, levelBouss, numBoussCells
 
     !! pass in numBoussCells for this level so can dimension this array
-    real(kind=8) :: soln(0:2*numBoussCells), rhs(0:2*numBoussCells) 
+    real(kind=8) :: soln(2*numBoussCells), rhs(2*numBoussCells) 
     real(kind=8) fms
 
     type(matrix_patchIndex), pointer :: mi

src/2d/bouss/prepBuildSparseMatrixSGNcrs.f90

@@ -11,7 +11,7 @@ subroutine prepBuildSparseMatrixSGNcrs(soln,rhs,nvar,naux,levelBouss,numBoussCel
     integer, intent(in) :: nvar, naux, levelBouss, numBoussCells
 
     !! pass in numBoussCells for this level so can dimension this array
-    real(kind=8) :: soln(0:2*numBoussCells), rhs(0:2*numBoussCells) 
+    real(kind=8) :: soln(0:2*numBoussCells-1), rhs(0:2*numBoussCells-1) 
 
     type(matrix_patchIndex), pointer :: mi
     type(matrix_levInfo),  pointer :: minfo