Merge pull request xcompact3d#26 from semi-h/feature

Add divergence and gradient

Author: semi-h

src/backend.f90

@@ -30,7 +30,11 @@ module m_base_backend
       procedure(tds_solve), deferred :: tds_solve
       procedure(transposer), deferred :: trans_x2y
       procedure(transposer), deferred :: trans_x2z
+      procedure(trans_d2d), deferred :: trans_y2z
+      procedure(trans_d2d), deferred :: trans_z2y
+      procedure(trans_d2d), deferred :: trans_y2x
       procedure(sum9into3), deferred :: sum_yzintox
+      procedure(vecadd), deferred :: vecadd
       procedure(get_fields), deferred :: get_fields
       procedure(set_fields), deferred :: set_fields
       procedure(alloc_tdsops), deferred :: alloc_tdsops
@@ -90,6 +94,20 @@ subroutine transposer(self, u_, v_, w_, u, v, w)
          class(field_t), intent(inout) :: u_, v_, w_
          class(field_t), intent(in) :: u, v, w
       end subroutine transposer
+
+      subroutine trans_d2d(self, u_, u)
+         !! transposer subroutines are straightforward, they rearrange
+         !! data into our specialist data structure so that regardless
+         !! of the direction tridiagonal systems are solved efficiently
+         !! and fast.
+         import :: base_backend_t
+         import :: field_t
+         implicit none
+
+         class(base_backend_t) :: self
+         class(field_t), intent(inout) :: u_
+         class(field_t), intent(in) :: u
+      end subroutine trans_d2d
    end interface
 
    abstract interface
@@ -106,6 +124,22 @@ subroutine sum9into3(self, du, dv, dw, du_y, dv_y, dw_y, du_z, dv_z, dw_z)
       end subroutine sum9into3
    end interface
 
+   abstract interface
+      subroutine vecadd(self, a, x, b, y)
+         !! adds two vectors together: y = a*x + b*y
+         import :: base_backend_t
+         import :: dp
+         import :: field_t
+         implicit none
+
+         class(base_backend_t) :: self
+         real(dp), intent(in) :: a
+         class(field_t), intent(in) :: x
+         real(dp), intent(in) :: b
+         class(field_t), intent(inout) :: y
+      end subroutine vecadd
+   end interface
+
    abstract interface
       subroutine get_fields(self, u_out, v_out, w_out, u, v, w)
          !! copy the specialist data structure from device or host back
@@ -136,7 +170,8 @@ end subroutine set_fields
    end interface
 
    abstract interface
-      subroutine alloc_tdsops(self, tdsops, n, dx, operation, scheme)
+      subroutine alloc_tdsops(self, tdsops, n, dx, operation, scheme, n_halo, &
+                              from_to, bc_start, bc_end, sym, c_nu, nu0_nu)
          import :: base_backend_t
          import :: dp
          import :: tdsops_t
@@ -147,6 +182,10 @@ subroutine alloc_tdsops(self, tdsops, n, dx, operation, scheme)
          integer, intent(in) :: n
          real(dp), intent(in) :: dx
          character(*), intent(in) :: operation, scheme
+         integer, optional, intent(in) :: n_halo
+         character(*), optional, intent(in) :: from_to, bc_start, bc_end
+         logical, optional, intent(in) :: sym
+         real(dp), optional, intent(in) :: c_nu, nu0_nu
       end subroutine alloc_tdsops
    end interface
 

src/cuda/backend.f90

@@ -34,7 +34,11 @@ module m_cuda_backend
       procedure :: tds_solve => tds_solve_cuda
       procedure :: trans_x2y => trans_x2y_cuda
       procedure :: trans_x2z => trans_x2z_cuda
+      procedure :: trans_y2z => trans_y2z_cuda
+      procedure :: trans_z2y => trans_z2y_cuda
+      procedure :: trans_y2x => trans_y2x_cuda
       procedure :: sum_yzintox => sum_yzintox_cuda
+      procedure :: vecadd => vecadd_cuda
       procedure :: set_fields => set_fields_cuda
       procedure :: get_fields => get_fields_cuda
       procedure :: transeq_cuda_dist
@@ -101,20 +105,28 @@ function init(globs, allocator) result(backend)
 
    end function init
 
-   subroutine alloc_cuda_tdsops(self, tdsops, n, dx, operation, scheme)
+   subroutine alloc_cuda_tdsops( &
+      self, tdsops, n, dx, operation, scheme, &
+      n_halo, from_to, bc_start, bc_end, sym, c_nu, nu0_nu &
+      )
       implicit none
 
       class(cuda_backend_t) :: self
       class(tdsops_t), allocatable, intent(inout) :: tdsops
       integer, intent(in) :: n
       real(dp), intent(in) :: dx
       character(*), intent(in) :: operation, scheme
+      integer, optional, intent(in) :: n_halo
+      character(*), optional, intent(in) :: from_to, bc_start, bc_end
+      logical, optional, intent(in) :: sym
+      real(dp), optional, intent(in) :: c_nu, nu0_nu
 
       allocate(cuda_tdsops_t :: tdsops)
 
       select type (tdsops)
       type is (cuda_tdsops_t)
-         tdsops = cuda_tdsops_t(n, dx, operation, scheme)
+         tdsops = cuda_tdsops_t(n, dx, operation, scheme, n_halo, from_to, &
+                                bc_start, bc_end, sym, c_nu, nu0_nu)
       end select
 
    end subroutine alloc_cuda_tdsops
@@ -421,6 +433,33 @@ subroutine trans_x2z_cuda(self, u_z, v_z, w_z, u, v, w)
 
    end subroutine trans_x2z_cuda
 
+   subroutine trans_y2z_cuda(self, u_z, u_y)
+      implicit none
+
+      class(cuda_backend_t) :: self
+      class(field_t), intent(inout) :: u_z
+      class(field_t), intent(in) :: u_y
+
+   end subroutine trans_y2z_cuda
+
+   subroutine trans_z2y_cuda(self, u_y, u_z)
+      implicit none
+
+      class(cuda_backend_t) :: self
+      class(field_t), intent(inout) :: u_y
+      class(field_t), intent(in) :: u_z
+
+   end subroutine trans_z2y_cuda
+
+   subroutine trans_y2x_cuda(self, u_x, u_y)
+      implicit none
+
+      class(cuda_backend_t) :: self
+      class(field_t), intent(inout) :: u_x
+      class(field_t), intent(in) :: u_y
+
+   end subroutine trans_y2x_cuda
+
    subroutine sum_yzintox_cuda(self, du, dv, dw, &
                                du_y, dv_y, dw_y, du_z, dv_z, dw_z)
       implicit none
@@ -431,6 +470,17 @@ subroutine sum_yzintox_cuda(self, du, dv, dw, &
 
    end subroutine sum_yzintox_cuda
 
+   subroutine vecadd_cuda(self, a, x, b, y)
+      implicit none
+
+      class(cuda_backend_t) :: self
+      real(dp), intent(in) :: a
+      class(field_t), intent(in) :: x
+      real(dp), intent(in) :: b
+      class(field_t), intent(inout) :: y
+
+   end subroutine vecadd_cuda
+
    subroutine copy_into_buffers(u_send_s_dev, u_send_e_dev, u_dev, n)
       implicit none
 

src/omp/backend.f90

@@ -26,7 +26,11 @@ module m_omp_backend
       procedure :: tds_solve => tds_solve_omp
       procedure :: trans_x2y => trans_x2y_omp
       procedure :: trans_x2z => trans_x2z_omp
+      procedure :: trans_y2z => trans_y2z_omp
+      procedure :: trans_z2y => trans_z2y_omp
+      procedure :: trans_y2x => trans_y2x_omp
       procedure :: sum_yzintox => sum_yzintox_omp
+      procedure :: vecadd => vecadd_omp
       procedure :: set_fields => set_fields_omp
       procedure :: get_fields => get_fields_omp
    end type omp_backend_t
@@ -83,20 +87,28 @@ function init(globs, allocator) result(backend)
 
    end function init
 
-   subroutine alloc_omp_tdsops(self, tdsops, n, dx, operation, scheme)
+   subroutine alloc_omp_tdsops( &
+      self, tdsops, n, dx, operation, scheme, &
+      n_halo, from_to, bc_start, bc_end, sym, c_nu, nu0_nu &
+      )
       implicit none
 
       class(omp_backend_t) :: self
       class(tdsops_t), allocatable, intent(inout) :: tdsops
       integer, intent(in) :: n
       real(dp), intent(in) :: dx
       character(*), intent(in) :: operation, scheme
+      integer, optional, intent(in) :: n_halo
+      character(*), optional, intent(in) :: from_to, bc_start, bc_end
+      logical, optional, intent(in) :: sym
+      real(dp), optional, intent(in) :: c_nu, nu0_nu
 
       allocate(tdsops_t :: tdsops)
 
       select type (tdsops)
       type is (tdsops_t)
-         tdsops = tdsops_t(n, dx, operation, scheme)
+         tdsops = tdsops_t(n, dx, operation, scheme, n_halo, from_to, &
+                           bc_start, bc_end, sym, c_nu, nu0_nu)
       end select
 
    end subroutine alloc_omp_tdsops
@@ -170,6 +182,33 @@ subroutine trans_x2z_omp(self, u_, v_, w_, u, v, w)
 
    end subroutine trans_x2z_omp
 
+   subroutine trans_y2z_omp(self, u_, u)
+      implicit none
+
+      class(omp_backend_t) :: self
+      class(field_t), intent(inout) :: u_
+      class(field_t), intent(in) :: u
+
+   end subroutine trans_y2z_omp
+
+   subroutine trans_z2y_omp(self, u_, u)
+      implicit none
+
+      class(omp_backend_t) :: self
+      class(field_t), intent(inout) :: u_
+      class(field_t), intent(in) :: u
+
+   end subroutine trans_z2y_omp
+
+   subroutine trans_y2x_omp(self, u_, u)
+      implicit none
+
+      class(omp_backend_t) :: self
+      class(field_t), intent(inout) :: u_
+      class(field_t), intent(in) :: u
+
+   end subroutine trans_y2x_omp
+
    subroutine sum_yzintox_omp(self, du, dv, dw, &
                                du_y, dv_y, dw_y, du_z, dv_z, dw_z)
       implicit none
@@ -180,6 +219,17 @@ subroutine sum_yzintox_omp(self, du, dv, dw, &
 
    end subroutine sum_yzintox_omp
 
+   subroutine vecadd_omp(self, a, x, b, y)
+      implicit none
+
+      class(omp_backend_t) :: self
+      real(dp), intent(in) :: a
+      class(field_t), intent(in) :: x
+      real(dp), intent(in) :: b
+      class(field_t), intent(inout) :: y
+
+   end subroutine vecadd_omp
+
    subroutine copy_into_buffers(u_send_s, u_send_e, u, n)
       implicit none