Fix issues in vector3d side exchange and 3d mpi gpu comms

Author: fluidnumerics-joe

CMakeLists.txt

@@ -68,7 +68,8 @@ ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90)
 # Fortran compiler flags
 if( "${CMAKE_Fortran_COMPILER_ID}" STREQUAL "GNU" )
     set( CMAKE_Fortran_FLAGS "${CMAKE_FORTRAN_FLAGS} -cpp -ffree-line-length-512" )
-    set( CMAKE_Fortran_FLAGS_DEBUG "-g -O0 -C -Wall -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow" )
+    set( CMAKE_Fortran_FLAGS_DEBUG "-g -O0 -C -Wall -fbounds-check -fbacktrace" )
+    #set( CMAKE_Fortran_FLAGS_DEBUG "-g -O0 -C -Wall -fbounds-check -fbacktrace -ffpe-trap=invalid,zero,overflow" )
     set( CMAKE_Fortran_FLAGS_COVERAGE "${CMAKE_Fortran_FLAGS_DEBUG} --coverage")
     set( CMAKE_Fortran_FLAGS_PROFILE "-pg -O3")
     set( CMAKE_Fortran_FLAGS_RELEASE "-O3" )

src/SELF_DomainDecomposition_t.f90

@@ -66,7 +66,7 @@ subroutine Init_DomainDecomposition_t(this,enableMPI)
 #undef __FUNC__
 #define __FUNC__ "Init_DomainDecomposition_t"
     implicit none
-    class(DomainDecomposition_t),intent(out) :: this
+    class(DomainDecomposition_t),intent(inout) :: this
     logical,intent(in) :: enableMPI
     ! Local
     integer       :: ierror
@@ -81,9 +81,9 @@ subroutine Init_DomainDecomposition_t(this,enableMPI)
     if(enableMPI) then
       this%mpiComm = MPI_COMM_WORLD
       print*,__FILE__," : Initializing MPI"
-      call MPI_INIT(ierror)
-      call MPI_COMM_RANK(this%mpiComm,this%rankId,ierror)
-      call MPI_COMM_SIZE(this%mpiComm,this%nRanks,ierror)
+      call mpi_init(ierror)
+      call mpi_comm_rank(this%mpiComm,this%rankId,ierror)
+      call mpi_comm_size(this%mpiComm,this%nRanks,ierror)
       print*,__FILE__," : Rank ",this%rankId+1,"/",this%nRanks," checking in."
     else
       print*,__FILE__," : MPI not initialized. No domain decomposition used."

src/SELF_MappedScalar_3D_t.f90

@@ -334,6 +334,7 @@ subroutine SideExchange_MappedScalar3D_t(this,mesh)
         if(r2 == rankId) then
 
           e2 = e2Global-offset
+          
           if(flip == 0) then
 
             do j = 1,this%interp%N+1

src/SELF_MappedVector_3D_t.f90

@@ -316,7 +316,7 @@ subroutine SideExchange_MappedVector3D_t(this,mesh)
     ! Local
     integer :: e1,e2,s1,s2,e2Global
     integer :: flip
-    integer :: i1,i2,j1,j2,ivar
+    integer :: i,i2,j,j2,ivar
     integer :: r2
     integer :: rankId,offset
     integer :: idir
@@ -333,7 +333,6 @@ subroutine SideExchange_MappedVector3D_t(this,mesh)
     do concurrent(s1=1:6,e1=1:mesh%nElem,ivar=1:this%nvar,idir=1:3)
 
       e2Global = mesh%sideInfo(3,s1,e1)
-      e2 = e2Global-offset
       s2 = mesh%sideInfo(4,s1,e1)/10
       flip = mesh%sideInfo(4,s1,e1)-s2*10
 
@@ -343,58 +342,86 @@ subroutine SideExchange_MappedVector3D_t(this,mesh)
 
         if(r2 == rankId) then
 
+          e2 = e2Global-offset
+
           if(flip == 0) then
 
-            do j1 = 1,this%interp%N+1
-              do i1 = 1,this%interp%N+1
-                this%extBoundary(i1,j1,s1,e1,ivar,idir) = &
-                  this%boundary(i1,j1,s2,e2,ivar,idir)
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                this%extBoundary(i,j,s1,e1,ivar,idir) = &
+                  this%boundary(i,j,s2,e2,ivar,idir)
               enddo
             enddo
 
           else if(flip == 1) then
 
-            do j1 = 1,this%interp%N+1
-              do i1 = 1,this%interp%N+1
-
-                i2 = j1
-                j2 = this%interp%N+2-i1
-                this%extBoundary(i1,j1,s1,e1,ivar,idir) = &
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = this%interp%N+2-i
+                j2 = j
+                this%extBoundary(i,j,s1,e1,ivar,idir) = &
                   this%boundary(i2,j2,s2,e2,ivar,idir)
-
               enddo
             enddo
 
           else if(flip == 2) then
 
-            do j1 = 1,this%interp%N+1
-              do i1 = 1,this%interp%N+1
-                i2 = this%interp%N+2-i1
-                j2 = this%interp%N+2-j1
-                this%extBoundary(i1,j1,s1,e1,ivar,idir) = &
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = this%interp%N+2-i
+                j2 = this%interp%N+2-j
+                this%extBoundary(i,j,s1,e1,ivar,idir) = &
                   this%boundary(i2,j2,s2,e2,ivar,idir)
               enddo
             enddo
 
           else if(flip == 3) then
 
-            do j1 = 1,this%interp%N+1
-              do i1 = 1,this%interp%N+1
-                i2 = this%interp%N+2-j1
-                j2 = i1
-                this%extBoundary(i1,j1,s1,e1,ivar,idir) = &
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = i
+                j2 = this%interp%N+2-j
+                this%extBoundary(i,j,s1,e1,ivar,idir) = &
                   this%boundary(i2,j2,s2,e2,ivar,idir)
               enddo
             enddo
 
           else if(flip == 4) then
 
-            do j1 = 1,this%interp%N+1
-              do i1 = 1,this%interp%N+1
-                i2 = j1
-                j2 = i1
-                this%extBoundary(i1,j1,s1,e1,ivar,idir) = &
-                  this%boundary(i2,j2,s2,e2,ivar,idir)
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                this%extBoundary(i,j,s1,e1,ivar,idir) = &
+                  this%boundary(j,i,s2,e2,ivar,idir)
+              enddo
+            enddo
+
+          else if(flip == 5) then
+
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = this%interp%N+2-j
+                j2 = i
+                this%extBoundary(i,j,s1,e1,ivar,idir) = this%boundary(i2,j2,s2,e2,ivar,idir)
+              enddo
+            enddo
+
+          else if(flip == 6) then
+
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = this%interp%N+2-j
+                j2 = this%interp%N+2-i
+                this%extBoundary(i,j,s1,e1,ivar,idir) = this%boundary(i2,j2,s2,e2,ivar,idir)
+              enddo
+            enddo
+
+          else if(flip == 7) then
+
+            do j = 1,this%interp%N+1
+              do i = 1,this%interp%N+1
+                i2 = j
+                j2 = this%interp%N+2-i
+                this%extBoundary(i,j,s1,e1,ivar,idir) = this%boundary(i2,j2,s2,e2,ivar,idir)
               enddo
             enddo
 

src/gpu/SELF_DomainDecomposition.f90

@@ -37,64 +37,14 @@ module SELF_DomainDecomposition
 
   contains
 
-    procedure :: Init => Init_DomainDecomposition
     procedure :: Free => Free_DomainDecomposition
 
     procedure :: SetElemToRank => SetElemToRank_DomainDecomposition
 
   endtype DomainDecomposition
 
-  interface
-    function check_gpu_aware_support() bind(c,name="check_gpu_aware_support")
-      use iso_c_binding
-      integer(c_int) :: check_gpu_aware_support
-    endfunction check_gpu_aware_support
-  endinterface
 contains
 
-  subroutine Init_DomainDecomposition(this,enableMPI)
-    implicit none
-    class(DomainDecomposition),intent(out) :: this
-    logical,intent(in) :: enableMPI
-    ! Local
-    integer       :: ierror
-    integer(c_int) :: gpuaware
-
-    this%mpiComm = 0
-    this%mpiPrec = prec
-    this%rankId = 0
-    this%nRanks = 1
-    this%nElem = 0
-    this%mpiEnabled = enableMPI
-
-    if(enableMPI) then
-      this%mpiComm = MPI_COMM_WORLD
-      print*,__FILE__," : Initializing MPI"
-      call MPI_INIT(ierror)
-
-      if(check_gpu_aware_support() == 0) then
-        print*,__FILE__" : Error! GPU Aware support is not detected. Stopping."
-        call MPI_FINALIZE(ierror)
-        stop 1
-      endif
-
-      call MPI_COMM_RANK(this%mpiComm,this%rankId,ierror)
-      call MPI_COMM_SIZE(this%mpiComm,this%nRanks,ierror)
-      print*,__FILE__," : Rank ",this%rankId+1,"/",this%nRanks," checking in."
-    else
-      print*,__FILE__," : MPI not initialized. No domain decomposition used."
-    endif
-
-    if(prec == real32) then
-      this%mpiPrec = MPI_FLOAT
-    else
-      this%mpiPrec = MPI_DOUBLE
-    endif
-
-    allocate(this%offsetElem(1:this%nRanks+1))
-
-  endsubroutine Init_DomainDecomposition
-
   subroutine Free_DomainDecomposition(this)
     implicit none
     class(DomainDecomposition),intent(inout) :: this