Support PETSc >= 3.22

examples/structured/test.py

@@ -1,10 +1,14 @@
 # this is a simple script to test the IDWarp module
+import os
 from mpi4py import MPI
 from idwarp import USMesh
 
+# Get the path of this file
+baseDir = os.path.dirname(os.path.abspath(__file__))
+
 options = {
     "fileType": "CGNS",
-    "gridFile": "../../input_files/o_mesh.cgns",
+    "gridFile": os.path.join(baseDir,"../../input_files/o_mesh.cgns"),
     "aExp": 3.0,
     "bExp": 5.0,
     "LdefFact": 100.0,

src/IO/readStructuredCGNS.F90

@@ -38,7 +38,7 @@ subroutine readStructuredCGNS(cg)
         ! Always the first base
         base = 1_intType
 
-        call cg_nzones_f(cg, base, nZones, ierr);
+        call cg_nzones_f(cg, base, nZones, ierr)
         if (ierr .eq. CG_ERROR) call cg_error_exit_f
         allocate (zones(nZones))
 

src/IO/writeCGNS.F90

@@ -4,6 +4,7 @@ subroutine writeCGNS(cgns_file)
     use communication
     use CGNSGrid
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Input Arguments
@@ -30,7 +31,7 @@ subroutine writeCGNS(cgns_file)
     ! Only do writing on root proc:
     rootProc: if (myid == 0) then
 
-        call VecGetArrayF90(XvLocal, xx, ierr)
+        call VecGetArrayCompat(XvLocal, xx, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
         ! Open and get the number of zones:
@@ -79,7 +80,7 @@ subroutine writeCGNS(cgns_file)
         call cg_close_f(cg, ierr)
         if (ierr .eq. CG_ERROR) call cg_error_exit_f
 
-        call VecRestoreArrayF90(XvLocal, xx, ierr)
+        call VecRestoreArrayCompat(XvLocal, xx, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
     end if rootProc

src/IO/writePlot3d.F90

@@ -4,6 +4,7 @@ subroutine writePlot3d(plot3d_file)
     use communication
     use CGNSGrid
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Input Arguments
@@ -29,7 +30,7 @@ subroutine writePlot3d(plot3d_file)
     ! Only do writing on root proc:
     rootProc: if (myid == 0) then
 
-        call VecGetArrayF90(XvLocal, xx, ierr)
+        call VecGetArrayCompat(XvLocal, xx, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
         ! Open the new file
@@ -79,7 +80,7 @@ subroutine writePlot3d(plot3d_file)
             deallocate (coorX, coorY, coorZ)
         end do
 
-        call VecRestoreArrayF90(XvLocal, xx, ierr)
+        call VecRestoreArrayCompat(XvLocal, xx, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
         ! Close the output file

src/modules/Makefile

@@ -20,6 +20,7 @@ FF90_OBJECTS_3 = gridInput.o\
 		 communication.o\
 
 FF90_OBJECTS_4 = gridData.o\
+		 petscCompat.o\
 
 FF90_OBJECTS_5 = kd_tree.o\
 		 kd_tree2.o\
@@ -30,6 +31,7 @@ FILES_TO_COMPLEXIFY =	precision.F90\
 			gridInput.f90\
 			communication.f90\
 			gridData.F90\
+			petscCompat.F90\
 			kd_tree.F90\
 			cgnsGrid.F90\
 			kd_tree2.f90\

src/modules/gridData.F90

@@ -54,7 +54,7 @@ module gridData
 
     real(kind=realType), dimension(:), allocatable :: d2wall
     real(kind=realType), dimension(:), allocatable :: denominator, denominator0
-    real(kind=realType), dimension(:, :), allocatable :: numerator
+    real(kind=realType), dimension(:, :), allocatable, target :: numerator
 
     ! Symmetry Information
     integer(kind=intType) :: nLoop

src/modules/kd_tree.F90

@@ -1552,6 +1552,7 @@ subroutine initNodalProperties(tp)
         use gridInput
         use gridData
         use communication
+        use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
         implicit none
 
         Type(tree_master_record), Pointer :: tp
@@ -1565,14 +1566,14 @@ subroutine initNodalProperties(tp)
         call EChk(ierr, __FILE__, __LINE__)
 
         ! Extract a pointer from Xs to use in the main routine
-        call VecGetArrayF90(XsLocal, XsPtr, ierr)
+        call VecGetArrayCompat(XsLocal, XsPtr, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
         call computeNodalProperties(tp, .True.)
         call determineCorners(tp)
 
         ! Restore all the arrays
-        call VecRestoreArrayF90(XsLocal, XsPtr, ierr)
+        call VecRestoreArrayCompat(XsLocal, XsPtr, ierr)
         call EChk(ierr, __FILE__, __LINE__)
 
     end subroutine initNodalProperties

src/modules/petscCompat.F90

@@ -0,0 +1,131 @@
+module petscCompat
+! Include PETSc version macros
+#include <petsc/finclude/petsc.h>
+    ! Module for PETSc compatibility functions
+    ! This module provides compatibility for various subroutines in PETSc
+    ! This is because of the changes in the API between different versions of PETSc
+    ! This module is used to ensure that the code works with different versions of PETSc
+
+    use petscvec
+    use constants
+    use communication
+    implicit none
+    save
+
+    ! Internal static array needed for PETSc <3.23 VecGetOwnershipRanges
+    integer(kind=intType), allocatable, target :: array(:)
+    logical :: arrayAllocated = .false.
+
+contains
+
+    subroutine VecGetArrayCompat(v, array, ierr)
+        ! PETSc compatibility routine to get array from vector
+        implicit none
+        Vec :: v
+        real(kind=realType), pointer :: array(:)
+        integer(kind=intType) :: ierr
+
+#if PETSC_VERSION_LT(3,20,0)
+        call VecGetArrayF90(v, array, ierr)
+#else
+        call VecGetArray(v, array, ierr)
+#endif
+
+    end subroutine VecGetArrayCompat
+
+    subroutine VecRestoreArrayCompat(v, array, ierr)
+        ! PETSc compatibility routine to restore array
+        implicit none
+        Vec :: v
+        real(kind=realType), pointer :: array(:)
+        integer(kind=intType) :: ierr
+
+#if PETSC_VERSION_LT(3,20,0)
+        call VecRestoreArrayF90(v, array, ierr)
+#else
+        call VecRestoreArray(v, array, ierr)
+#endif
+
+    end subroutine VecRestoreArrayCompat
+
+    subroutine VecCreateMPIWithArrayCompat(comm, bs, n, nGlobal, array, v, ierr)
+        ! PETSc compatibility routine to create a MPI vector with an array.
+        ! In some cases array is not supplied so it supports optional array argument.
+        implicit none
+        integer(kind=intType), intent(in) :: comm, bs
+        integer(kind=intType), intent(in) :: n, nGlobal
+        real(kind=realType), dimension(:), intent(in), optional :: array
+        Vec :: v
+        integer(kind=intType) :: ierr
+
+#if PETSC_VERSION_GE(3,22,0)
+        ! PETSc >= 3.22, need to use PETSC_NULL_SCALAR_ARRAY
+        if (present(array)) then
+            call VecCreateMPIWithArray(comm, bs, n, nGlobal, array, v, ierr)
+        else
+            call VecCreateMPIWithArray(comm, bs, n, nGlobal, PETSC_NULL_SCALAR_ARRAY, v, ierr)
+        end if
+#else
+        ! Older PETSc, use PETSC_NULL_SCALAR for older versions
+        if (present(array)) then
+            call VecCreateMPIWithArray(comm, bs, n, nGlobal, array, v, ierr)
+        else
+            call VecCreateMPIWithArray(comm, bs, n, nGlobal, PETSC_NULL_SCALAR, v, ierr)
+        end if
+#endif
+
+    end subroutine VecCreateMPIWithArrayCompat
+
+    subroutine VecGetOwnershipRangesCompat(v, ptr, ierr)
+        ! PETSc compatibility routine to get ownership ranges
+        Vec :: v
+        integer(kind=intType), pointer :: ptr(:)
+        integer(kind=intType), intent(out) :: ierr
+
+#if PETSC_VERSION_GE(3,23,0)
+        ! PETSc >= 3.23, directly get pointer
+        call VecGetOwnershipRanges(v, ptr, ierr)
+#else
+        ! Older PETSc, allocate array that gets passed to PETSc and set pointer
+        integer(kind=intType), allocatable, target, save :: array(:)
+
+        ! Check if we already have an allocated array
+        if (allocated(array)) then
+            deallocate (array)
+        end if
+        ! Use zero-based indexing to match to rank
+        allocate (array(0:nProc))
+        array = zero
+        arrayAllocated = .true.
+
+        ! Get the ownership ranges into the static array and set the pointer
+        call VecGetOwnershipRanges(v, array, ierr)
+        ptr => array
+#endif
+
+    end subroutine VecGetOwnershipRangesCompat
+
+    subroutine VecRestoreOwnershipRangesCompat(v, ptr, ierr)
+        ! PETSc compatibility routine to restore ownership ranges pointer
+        Vec :: v
+        integer(kind=intType), pointer :: ptr(:)
+        integer(kind=intType), intent(out) :: ierr
+
+#if PETSC_VERSION_GE(3,23,0)
+        ! PETSc >= 3.23, just restore ranges using the pointer
+        call VecRestoreOwnershipRanges(v, ptr, ierr)
+#else
+        ! Older PETSc, deallocate the internal static array and nullify pointer
+        if (arrayAllocated) then
+            if (allocated(array)) then
+                deallocate (array)
+            end if
+            nullify (ptr)
+            arrayAllocated = .false.
+        end if
+        ierr = 0
+#endif
+
+    end subroutine VecRestoreOwnershipRangesCompat
+
+end module petscCompat

src/warp/getCommonVolumeCoordinates.F90

@@ -3,6 +3,7 @@
 subroutine getCommonVolumeCoordinates(gridNodes, nDOF)
 
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Subroutine variables
@@ -20,20 +21,21 @@ subroutine getCommonVolumeCoordinates(gridNodes, nDOF)
                        INSERT_VALUES, SCATTER_REVERSE, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
-    call VecGetArrayF90(commonGridVec, xx, ierr)
+    call VecGetArrayCompat(commonGridVec, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Perform the actual copy
     gridNodes = xx
 
-    call VecRestoreArrayF90(Xv, xx, ierr)
+    call VecRestoreArrayCompat(Xv, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
 end subroutine getCommonVolumeCoordinates
 
 subroutine setCommonVolumeCoordinates(gridNodes, nDOF)
 
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Subroutine variables
@@ -42,13 +44,13 @@ subroutine setCommonVolumeCoordinates(gridNodes, nDOF)
     real(kind=realType), pointer, dimension(:) :: xx
     integer(kind=intType) :: ierr
 
-    call VecGetArrayF90(commonGridVec, xx, ierr)
+    call VecGetArrayCompat(commonGridVec, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Perform actual copy
     xx = gridNodes
 
-    call VecRestoreArrayF90(commonGridVec, xx, ierr)
+    call VecRestoreArrayCompat(commonGridVec, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Now do the vecScatter: warp_to_common in REVERSE

src/warp/getVolumeCoordinates.F90

@@ -2,6 +2,7 @@
 subroutine getVolumeCoordinates(gridNodes, nDOF)
 
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Subroutine variables
@@ -11,20 +12,21 @@ subroutine getVolumeCoordinates(gridNodes, nDOF)
     integer(kind=intType) :: ierr
 
     ! No error checking here!
-    call VecGetArrayF90(Xv, xx, ierr)
+    call VecGetArrayCompat(Xv, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Perform copy
     gridNodes = xx
 
-    call VecRestoreArrayF90(Xv, xx, ierr)
+    call VecRestoreArrayCompat(Xv, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
 end subroutine getVolumeCoordinates
 
 subroutine setVolumeCoordinates(gridNodes, nDOF)
 
     use gridData
+    use petscCompat, only: VecGetArrayCompat, VecRestoreArrayCompat
     implicit none
 
     ! Subroutine variables
@@ -34,13 +36,13 @@ subroutine setVolumeCoordinates(gridNodes, nDOF)
     integer(kind=intType) :: ierr
 
     ! No error checking here!
-    call VecGetArrayF90(Xv, xx, ierr)
+    call VecGetArrayCompat(Xv, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Perform copy
     xx = gridNodes
 
-    call VecRestoreArrayF90(Xv, xx, ierr)
+    call VecRestoreArrayCompat(Xv, xx, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
 end subroutine setVolumeCoordinates