From 1fac7990966a79a56890a9f30fe497f7514ef54e Mon Sep 17 00:00:00 2001 From: Semih Akkurt Date: Fri, 20 Dec 2024 15:13:39 +0000 Subject: [PATCH] Add infrastructure for 010 BCs. --- src/cuda/poisson_fft.f90 | 56 ++++++++++++++++++++++++++++++++++++++-- src/omp/poisson_fft.f90 | 27 +++++++++++++++++++ src/poisson_fft.f90 | 55 +++++++++++++++++++++++++++++---------- src/solver.f90 | 2 +- 4 files changed, 124 insertions(+), 16 deletions(-) diff --git a/src/cuda/poisson_fft.f90 b/src/cuda/poisson_fft.f90 index 6a23b6e7..f958eab6 100644 --- a/src/cuda/poisson_fft.f90 +++ b/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<<>>( & !& + 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 diff --git a/src/omp/poisson_fft.f90 b/src/omp/poisson_fft.f90 index 4b7fefb2..5e5a54ec 100644 --- a/src/omp/poisson_fft.f90 +++ b/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 diff --git a/src/poisson_fft.f90 b/src/poisson_fft.f90 index f9af64dc..1423c373 100644 --- a/src/poisson_fft.f90 +++ b/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 !! diff --git a/src/solver.f90 b/src/solver.f90 index eb3855fa..edb8b300 100644 --- a/src/solver.f90 +++ b/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)