Add infrastructure for 010 BCs.

src/cuda/poisson_fft.f90

@@ -13,7 +13,7 @@ module m_cuda_poisson_fft
   use m_tdsops, only: dirps_t
 
   use m_cuda_allocator, only: cuda_field_t
-  use m_cuda_spectral, only: process_spectral_000
+  use m_cuda_spectral, only: process_spectral_000, process_spectral_010
 
   implicit none
 
@@ -35,7 +35,9 @@ module m_cuda_poisson_fft
     procedure :: fft_forward => fft_forward_cuda
     procedure :: fft_backward => fft_backward_cuda
     procedure :: fft_postprocess_000 => fft_postprocess_000_cuda
-
+    procedure :: fft_postprocess_010 => fft_postprocess_010_cuda
+    procedure :: enforce_periodicity_y => enforce_periodicity_y_cuda
+    procedure :: undo_periodicity_y => undo_periodicity_y_cuda
   end type cuda_poisson_fft_t
 
   interface cuda_poisson_fft_t
@@ -213,4 +215,54 @@ subroutine fft_postprocess_000_cuda(self)
 
   end subroutine fft_postprocess_000_cuda
 
+  subroutine fft_postprocess_010_cuda(self)
+    implicit none
+
+    class(cuda_poisson_fft_t) :: self
+
+    type(cudaXtDesc), pointer :: descriptor
+
+    complex(dp), device, dimension(:, :, :), pointer :: c_dev
+    type(dim3) :: blocks, threads
+    integer :: tsize
+
+    ! obtain a pointer to descriptor so that we can carry out postprocessing
+    call c_f_pointer(self%xtdesc%descriptor, descriptor)
+    call c_f_pointer(descriptor%data(1), c_dev, &
+                     [self%nx_spec, self%ny_spec, self%nz_spec])
+
+    ! tsize is different than SZ, because here we work on a 3D Cartesian
+    ! data structure, and free to specify any suitable thread/block size.
+    tsize = 16
+    blocks = dim3((self%nx_spec - 1)/tsize + 1, self%nz_spec, 1)
+    threads = dim3(tsize, 1, 1)
+
+    ! Postprocess div_u in spectral space
+    call process_spectral_010<<<blocks, threads>>>( & !&
+      c_dev, self%waves_dev, self%nx_spec, self%ny_spec, self%y_sp_st, &
+      self%nx_glob, self%ny_glob, self%nz_glob, &
+      self%ax_dev, self%bx_dev, self%ay_dev, self%by_dev, &
+      self%az_dev, self%bz_dev &
+      )
+
+  end subroutine fft_postprocess_010_cuda
+
+  subroutine enforce_periodicity_y_cuda(self, f_out, f_in)
+    implicit none
+
+    class(cuda_poisson_fft_t) :: self
+    class(field_t), intent(inout) :: f_out
+    class(field_t), intent(in) :: f_in
+
+  end subroutine enforce_periodicity_y_cuda
+
+  subroutine undo_periodicity_y_cuda(self, f_out, f_in)
+    implicit none
+
+    class(cuda_poisson_fft_t) :: self
+    class(field_t), intent(inout) :: f_out
+    class(field_t), intent(in) :: f_in
+
+  end subroutine undo_periodicity_y_cuda
+
 end module m_cuda_poisson_fft

src/omp/poisson_fft.f90

@@ -15,6 +15,9 @@ module m_omp_poisson_fft
     procedure :: fft_forward => fft_forward_omp
     procedure :: fft_backward => fft_backward_omp
     procedure :: fft_postprocess_000 => fft_postprocess_000_omp
+    procedure :: fft_postprocess_010 => fft_postprocess_010_omp
+    procedure :: enforce_periodicity_y => enforce_periodicity_y_omp
+    procedure :: undo_periodicity_y => undo_periodicity_y_omp
   end type omp_poisson_fft_t
 
   interface omp_poisson_fft_t
@@ -57,4 +60,28 @@ subroutine fft_postprocess_000_omp(self)
     class(omp_poisson_fft_t) :: self
   end subroutine fft_postprocess_000_omp
 
+  subroutine fft_postprocess_010_omp(self)
+    implicit none
+
+    class(omp_poisson_fft_t) :: self
+  end subroutine fft_postprocess_010_omp
+
+  subroutine enforce_periodicity_y_omp(self, f_out, f_in)
+    implicit none
+
+    class(omp_poisson_fft_t) :: self
+    class(field_t), intent(inout) :: f_out
+    class(field_t), intent(in) :: f_in
+
+  end subroutine enforce_periodicity_y_omp
+
+  subroutine undo_periodicity_y_omp(self, f_out, f_in)
+    implicit none
+
+    class(omp_poisson_fft_t) :: self
+    class(field_t), intent(inout) :: f_out
+    class(field_t), intent(in) :: f_in
+
+  end subroutine undo_periodicity_y_omp
+
 end module m_omp_poisson_fft

src/poisson_fft.f90

@@ -30,6 +30,9 @@ module m_poisson_fft
     procedure(fft_forward), deferred :: fft_forward
     procedure(fft_backward), deferred :: fft_backward
     procedure(fft_postprocess), deferred :: fft_postprocess_000
+    procedure(fft_postprocess), deferred :: fft_postprocess_010
+    procedure(field_process), deferred :: enforce_periodicity_y
+    procedure(field_process), deferred :: undo_periodicity_y
     procedure :: solve_poisson
     procedure :: base_init
     procedure :: waves_set
@@ -63,13 +66,22 @@ end subroutine fft_postprocess
   end interface
 
   abstract interface
-    subroutine poisson_xxx(self, f)
+    subroutine poisson_xxx(self, f, temp)
       import :: poisson_fft_t
       import :: field_t
 
       class(poisson_fft_t) :: self
-      class(field_t), intent(inout) :: f
+      class(field_t), intent(inout) :: f, temp
     end subroutine poisson_xxx
+
+    subroutine field_process(self, f_out, f_in)
+      import :: poisson_fft_t
+      import :: field_t
+
+      class(poisson_fft_t) :: self
+      class(field_t), intent(inout) :: f_out
+      class(field_t), intent(in) :: f_in
+    end subroutine field_process
   end interface
 
 contains
@@ -114,36 +126,53 @@ subroutine base_init(self, mesh, xdirps, ydirps, zdirps)
     ! waves_set requires some of the preprocessed tdsops variables.
     call self%waves_set(mesh%geo, xdirps, ydirps, zdirps)
 
-    self%poisson => poisson_000
+    ! use correct procedure based on BCs
+    if (self%periodic_x .and. self%periodic_y .and. self%periodic_z) then
+      self%poisson => poisson_000
+    else if (self%periodic_x .and. (.not. self%periodic_y) &
+             .and. (self%periodic_z)) then
+      self%poisson => poisson_010
+    end if
   end subroutine base_init
 
-  subroutine solve_poisson(self, f)
+  subroutine solve_poisson(self, f, temp)
     implicit none
 
     class(poisson_fft_t) :: self
-    class(field_t), intent(inout) :: f
+    class(field_t), intent(inout) :: f, temp
 
-    call self%poisson(f)
+    call self%poisson(f, temp)
 
   end subroutine solve_poisson
 
-  subroutine poisson_000(self, f)
+  subroutine poisson_000(self, f, temp)
     implicit none
 
     class(poisson_fft_t) :: self
-    class(field_t), intent(inout) :: f
+    class(field_t), intent(inout) :: f, temp
 
-    ! call forward FFT
     call self%fft_forward(f)
-
-    ! postprocess
     call self%fft_postprocess_000
-
-    ! call backward FFT
     call self%fft_backward(f)
 
   end subroutine poisson_000
 
+  subroutine poisson_010(self, f, temp)
+    implicit none
+
+    class(poisson_fft_t) :: self
+    class(field_t), intent(inout) :: f, temp
+
+    call self%enforce_periodicity_y(temp, f)
+
+    call self%fft_forward(temp)
+    call self%fft_postprocess_010
+    call self%fft_backward(temp)
+
+    call self%undo_periodicity_y(f, temp)
+
+  end subroutine poisson_010
+
   subroutine waves_set(self, geo, xdirps, ydirps, zdirps)
     !! Spectral equivalence constants
     !!

src/solver.f90

@@ -352,7 +352,7 @@ subroutine poisson_fft(self, pressure, div_u)
     call self%backend%reorder(p_temp, div_u, RDR_Z2C)
 
     ! solve poisson equation with FFT based approach
-    call self%backend%poisson_fft%solve_poisson(p_temp)
+    call self%backend%poisson_fft%solve_poisson(p_temp, pressure)
 
     ! reorder back to our specialist data structure from 3D Cartesian
     call self%backend%reorder(pressure, p_temp, RDR_C2Z)