Pressure solver which uses a tridiagonal algorithm operating on the pressure terms in Fourier space. It uses the pencil FFT module for 3D FFTs in pencil decomposition. These use FFTW for the actual FFT kernel. More...

Functions/Subroutines
type(component_descriptor_type) function, public	fftsolver_get_descriptor ()
	Descriptor of this component for registration. More...

subroutine	initialisation_callback (current_state)
	This initialisation callback sets up the pencil fft module, allocates data for the fourier space pressure term (might be different size than p) and populates cos x and y. More...

subroutine	timestep_callback (current_state)
	Timestep call back, which will transform to Fourier space, do a tridiagonal solve and then back into time space. More...

subroutine	finalisation_callback (current_state)
	Called at MONC finalisation, will call to the pencil fft module to clean itself up and free the pressure term. More...

subroutine	complete_psrce_calculation (current_state, y_halo_size, x_halo_size)
	Completes the psrce calculation by waiting on all outstanding psrce communications to complete and then combine the received values with the P field for U and V. More...

subroutine	tridiagonal_solver (current_state, pressure_term, fourier_space_sizes)
	The tridiagonal solver which runs in Fourier space on the pressure terms. Note that because we are going from n reals to n/2+1 complex numbers (and into their reals for this solver) then there might be more data in the pressure term than in p. Therefore use the space sizes which are provided as an argument to this procedure rather than the current state local grid dimensions. More...

subroutine	copy_p_to_halo_buffer (current_state, neighbour_description, dim, source_index, pid_location, current_page, source_data)
	Copies the p field data to halo buffers for a specific process in a dimension and halo cell. More...

subroutine	copy_halo_buffer_to_p (current_state, neighbour_description, dim, target_index, neighbour_location, current_page, source_data)

subroutine	perform_local_data_copy_for_p (current_state, halo_depth, involve_corners, source_data)
	Does local data copying for P variable halo swap. More...

Variables
real(kind=default_precision), dimension(:), allocatable	cos_x

real(kind=default_precision), dimension(:), allocatable	cos_y

real(kind=default_precision)	pi

real(kind=default_precision), dimension(:,:,:), allocatable	pt

integer, dimension(3)	fourier_space_sizes

type(halo_communication_type), save	halo_swap_state

Detailed Description

Pressure solver which uses a tridiagonal algorithm operating on the pressure terms in Fourier space. It uses the pencil FFT module for 3D FFTs in pencil decomposition. These use FFTW for the actual FFT kernel.

Function/Subroutine Documentation

◆ complete_psrce_calculation()

subroutine fftsolver_mod::complete_psrce_calculation	(	type(model_state_type), intent(inout), target	current_state,
		integer, intent(in)	y_halo_size,
		integer, intent(in)	x_halo_size
	)

private

Completes the psrce calculation by waiting on all outstanding psrce communications to complete and then combine the received values with the P field for U and V.

Parameters

current_state	The current model state
y_halo_size	The halo size in the Y dimension
x_halo_size	The halo size in the X dimension

Definition at line 145 of file fftsolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     integer, intent(in) :: y_halo_size, x_halo_size
 
     integer :: ierr, combined_handles(2), i, j, k
 
     combined_handles(1)=current_state%psrce_x_hs_recv_request
     combined_handles(2)=current_state%psrce_y_hs_recv_request
     call mpi_waitall(2, combined_handles, mpi_statuses_ignore, ierr)
 
     do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
       do k=2,current_state%local_grid%size(z_index)
 #ifdef U_ACTIVE
         current_state%p%data(k,j,x_halo_size+1)=current_state%p%data(k,j,x_halo_size+1)-&
                current_state%psrce_recv_buffer_x(k-1,j-x_halo_size)
 #endif
 #ifdef V_ACTIVE
         if (j .gt. y_halo_size+1) current_state%p%data(k, j, x_halo_size+1)=current_state%p%data(k, j, x_halo_size+1)-&
              current_state%global_grid%configuration%horizontal%cy * current_state%sv%data(k, j-1, x_halo_size+1)
 #endif
       end do
     end do
       
 #ifdef V_ACTIVE
     do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
       do k=2,current_state%local_grid%size(z_index)
         current_state%p%data(k,y_halo_size+1,i)=current_state%p%data(k,y_halo_size+1,i)-&
              current_state%psrce_recv_buffer_y(k-1,i-y_halo_size)
       end do
     end do
 #endif
 
     combined_handles(1)=current_state%psrce_x_hs_send_request
     combined_handles(2)=current_state%psrce_y_hs_send_request
     call mpi_waitall(2, combined_handles, mpi_statuses_ignore, ierr)

Here is the caller graph for this function:

◆ copy_halo_buffer_to_p()

subroutine fftsolver_mod::copy_halo_buffer_to_p	(	type(model_state_type), intent(inout)	current_state,
		type(neighbour_description_type), intent(inout)	neighbour_description,
		integer, intent(in)	dim,
		integer, intent(in)	target_index,
		integer, intent(in)	neighbour_location,
		integer, dimension(:), intent(inout)	current_page,
		type(field_data_wrapper_type), dimension(:), intent(in), optional	source_data
	)

private

Definition at line 269 of file fftsolver.F90.

     type(model_state_type), intent(inout) :: current_state
     integer, intent(in) :: dim, target_index, neighbour_location
     integer, intent(inout) :: current_page(:)
     type(neighbour_description_type), intent(inout) :: neighbour_description
     type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
     call copy_buffer_to_field(current_state%local_grid, neighbour_description%recv_halo_buffer, current_state%p%data, &
          dim, target_index, current_page(neighbour_location))
 
     current_page(neighbour_location)=current_page(neighbour_location)+1

Here is the caller graph for this function:

◆ copy_p_to_halo_buffer()

subroutine fftsolver_mod::copy_p_to_halo_buffer	(	type(model_state_type), intent(inout)	current_state,
		type(neighbour_description_type), intent(inout)	neighbour_description,
		integer, intent(in)	dim,
		integer, intent(in)	source_index,
		integer, intent(in)	pid_location,
		integer, dimension(:), intent(inout)	current_page,
		type(field_data_wrapper_type), dimension(:), intent(in), optional	source_data
	)

private

Copies the p field data to halo buffers for a specific process in a dimension and halo cell.

Parameters

current_state	The current model state
neighbour_descriptions	Description of the neighbour halo swapping status
dim	Dimension to copy from
source_index	The source index of the dimension we are reading from in the prognostic field
pid_location	Location of the neighbouring process in the local stored data structures
current_page	The current (next) buffer page to copy into
source_data	Optional source data which is read from

Definition at line 250 of file fftsolver.F90.

     type(model_state_type), intent(inout) :: current_state
     integer, intent(in) :: dim, pid_location, source_index
     integer, intent(inout) :: current_page(:)
     type(neighbour_description_type), intent(inout) :: neighbour_description
     type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
     call copy_field_to_buffer(current_state%local_grid, neighbour_description%send_halo_buffer, current_state%p%data, &
          dim, source_index, current_page(pid_location))
 
     current_page(pid_location)=current_page(pid_location)+1

Here is the caller graph for this function:

◆ fftsolver_get_descriptor()

type(component_descriptor_type) function, public fftsolver_mod::fftsolver_get_descriptor ( )

Descriptor of this component for registration.

Returns: The fft solver component descriptor

Definition at line 33 of file fftsolver.F90.

     fftsolver_get_descriptor%name="fftsolver"
     fftsolver_get_descriptor%version=0.1
     fftsolver_get_descriptor%initialisation=>initialisation_callback
     fftsolver_get_descriptor%timestep=>timestep_callback
     fftsolver_get_descriptor%finalisation=>finalisation_callback

Here is the call graph for this function:

Here is the caller graph for this function:

◆ finalisation_callback()

subroutine fftsolver_mod::finalisation_callback ( type(model_state_type), intent(inout), target current_state )

private

Called at MONC finalisation, will call to the pencil fft module to clean itself up and free the pressure term.

Parameters

current_state The current model state

Definition at line 131 of file fftsolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
 
     call finalise_pencil_fft(current_state%parallel%monc_communicator)
     deallocate(pt)
     if (allocated(cos_x)) deallocate(cos_x)
     if (allocated(cos_y)) deallocate(cos_y)

Here is the caller graph for this function:

◆ initialisation_callback()

subroutine fftsolver_mod::initialisation_callback ( type(model_state_type), intent(inout), target current_state )

private

This initialisation callback sets up the pencil fft module, allocates data for the fourier space pressure term (might be different size than p) and populates cos x and y.

Definition at line 43 of file fftsolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
 
     integer :: i, my_y_start, my_x_start
 
     if (.not. is_component_enabled(current_state%options_database, "diverr")) then
       call log_master_log(log_error, "The FFT solver component requires the diverr component to be enabled")
     end if 
 
     pi=4.0_default_precision*atan(1.0_default_precision)
 
     fourier_space_sizes=initialise_pencil_fft(current_state, my_y_start, my_x_start)
 
     call init_halo_communication(current_state, get_single_field_per_halo_cell, halo_swap_state, 1, .false.)
 
     allocate(pt(fourier_space_sizes(z_index), fourier_space_sizes(y_index), fourier_space_sizes(x_index)))
 
 #ifdef U_ACTIVE
     allocate(cos_x(fourier_space_sizes(x_index)))
     do i=1, fourier_space_sizes(x_index)
       cos_x(i)=(2.0_default_precision/(current_state%global_grid%configuration%horizontal%dx*&
            current_state%global_grid%configuration%horizontal%dx))&
            *(cos(2.0_default_precision*pi*real(((i-1)+(my_x_start-1))/2, kind=default_precision)/&
            real(current_state%global_grid%size(X_INDEX), kind=default_precision))-1.0_default_precision)
     end do    
 #endif
 
 #ifdef V_ACTIVE
     allocate(cos_y(fourier_space_sizes(y_index)))
     do i=1, fourier_space_sizes(y_index)
       cos_y(i)=(2.0_default_precision/(current_state%global_grid%configuration%horizontal%dy*&
            current_state%global_grid%configuration%horizontal%dy))&
            *(cos(2.0_default_precision*pi*real(((i-1)+(my_y_start-1))/2, kind=default_precision)/&
            real(current_state%global_grid%size(Y_INDEX), kind=default_precision))-1.0_default_precision)
     end do    
 #endif
     current_state%psrce_x_hs_send_request=mpi_request_null
     current_state%psrce_y_hs_send_request=mpi_request_null
     current_state%psrce_x_hs_recv_request=mpi_request_null
     current_state%psrce_y_hs_recv_request=mpi_request_null 

Here is the caller graph for this function:

◆ perform_local_data_copy_for_p()

subroutine fftsolver_mod::perform_local_data_copy_for_p	(	type(model_state_type), intent(inout)	current_state,
		integer, intent(in)	halo_depth,
		logical, intent(in)	involve_corners,
		type(field_data_wrapper_type), dimension(:), intent(in), optional	source_data
	)

private

Does local data copying for P variable halo swap.

Parameters

current_state	The current model state_mod
source_data	Optional source data which is written into

Definition at line 285 of file fftsolver.F90.

     type(model_state_type), intent(inout) :: current_state
     integer, intent(in) :: halo_depth
     logical, intent(in) :: involve_corners
     type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
     call perform_local_data_copy_for_field(current_state%p%data, current_state%local_grid, &
          current_state%parallel%my_rank, halo_depth, involve_corners)

Here is the caller graph for this function:

◆ timestep_callback()

subroutine fftsolver_mod::timestep_callback ( type(model_state_type), intent(inout), target current_state )

private

Timestep call back, which will transform to Fourier space, do a tridiagonal solve and then back into time space.

Parameters

current_state The current model state

Definition at line 87 of file fftsolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
 
     integer :: start_loc(3), end_loc(3), i
 
     do i=1,3
       start_loc(i)=current_state%local_grid%local_domain_start_index(i)
       end_loc(i)=current_state%local_grid%local_domain_end_index(i)
     end do
 
     call complete_psrce_calculation(current_state, current_state%local_grid%halo_size(y_index), &
          current_state%local_grid%halo_size(x_index))
 
 #ifdef FFT_TEST_MODE
     current_state%p%data=real(current_state%parallel%my_rank, kind=8) + 1.d0
 #endif
 
     call perform_forward_3dfft(current_state, current_state%p%data(start_loc(z_index):end_loc(z_index), &
          start_loc(y_index):end_loc(y_index), start_loc(x_index):end_loc(x_index)), pt)
 
 #ifndef FFT_TEST_MODE
     call tridiagonal_solver(current_state, pt, fourier_space_sizes)
 #endif
     
     ! Here it is complex space, distributed as needed. The other option is real space but with +1 for the last process in Y
     ! as require that complex number - not sure the best approach. We can extract the values here from complex anyway so
     ! worth trying that first. Operating on the size of pt rather than grid sizes. Downside of real space xform is then we have +1
     ! in the forwards transformation    
     call perform_backwards_3dfft(current_state, pt, current_state%p%data(start_loc(z_index):end_loc(z_index), &
          start_loc(y_index):end_loc(y_index), start_loc(x_index):end_loc(x_index)))   
 
 #ifdef FFT_TEST_MODE
     if (current_state%parallel%my_rank == 1) then
       write(*,*) current_state%parallel%my_rank, current_state%p%data(:,3,3)
     end if
 #endif
 
     call blocking_halo_swap(current_state, halo_swap_state, copy_p_to_halo_buffer, &
          perform_local_data_copy_for_p, copy_halo_buffer_to_p)
     

Here is the call graph for this function:

Here is the caller graph for this function:

◆ tridiagonal_solver()

subroutine fftsolver_mod::tridiagonal_solver	(	type(model_state_type), intent(inout), target	current_state,
		real(kind=default_precision), dimension(:,:,:), intent(inout)	pressure_term,
		integer, dimension(3), intent(in)	fourier_space_sizes
	)

private

The tridiagonal solver which runs in Fourier space on the pressure terms. Note that because we are going from n reals to n/2+1 complex numbers (and into their reals for this solver) then there might be more data in the pressure term than in p. Therefore use the space sizes which are provided as an argument to this procedure rather than the current state local grid dimensions.

Parameters

current_state	The current model state
pressure_term	The pressure term in Fourier space
fourier_space_sizes	The size of the pressure term in each dimension

Definition at line 189 of file fftsolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     real(kind=DEFAULT_PRECISION), dimension(:,:,:), intent(inout) :: pressure_term
     integer, intent(in) :: fourier_space_sizes(3)
 
     integer :: i, j, k, j_start
     real(kind=DEFAULT_PRECISION) :: b(fourier_space_sizes(z_index)), b1(fourier_space_sizes(z_index)),   &
         s(fourier_space_sizes(z_index)), s1(fourier_space_sizes(z_index)), cij
     
     do i=1,fourier_space_sizes(x_index)
       j_start=merge(3, 1, current_state%local_grid%start(y_index) ==1 .and. current_state%local_grid%start(x_index) == 1 &
            .and. i .lt. 3)
       do j=j_start,fourier_space_sizes(y_index)
 
         cij=cos_x(i)+cos_y(j)        
         b(2)=cij-current_state%global_grid%configuration%vertical%cza(2)
         b1(2)=1.0_default_precision/b(2)
         s1(2)=pressure_term(2,j,i)
 
         do k=3,fourier_space_sizes(z_index)
            b(k)=cij+current_state%global_grid%configuration%vertical%czg(k)
            b1(k)=1.0_default_precision/(b(k)-current_state%global_grid%configuration%vertical%czh(k)*b1(k-1))
            s1(k)=pressure_term(k,j,i)-current_state%global_grid%configuration%vertical%czb(k)*s1(k-1)*b1(k-1)
         end do
 
         pressure_term(fourier_space_sizes(z_index),j,i)=s1(fourier_space_sizes(z_index))* b1(fourier_space_sizes(z_index))
 
         do k=fourier_space_sizes(z_index)-1,2,-1
           pressure_term(k,j,i)=(s1(k)-current_state%global_grid%configuration%vertical%cza(k)*&
                pressure_term(k+1,j,i))*b1(k)
         end do                
       end do
     end do    
 
     ! Handle the zero wavenumber, which will only be on processes where the X and Y  dimension starts at 1
     if (current_state%local_grid%start(y_index) == 1 .and. current_state%local_grid%start(x_index) == 1) then
       do i=1,2
         do j=1,2
           s(2)=pressure_term(2,j,i)
           pressure_term(1,j,i)=0.0_default_precision
           pressure_term(2,j,i)=0.0_default_precision
           do k=3, fourier_space_sizes(z_index)
             s(k)=pressure_term(k,j,i)
             pressure_term(k,j,i)=(s(k-1)-current_state%global_grid%configuration%vertical%czg(k-1)*pressure_term(k-1,j,i)-&
                  current_state%global_grid%configuration%vertical%czb(k-1)*pressure_term(k-2,j,i))/&
                  current_state%global_grid%configuration%vertical%cza(k-1)
           end do          
         end do        
       end do           
     end if

Here is the caller graph for this function:

Variable Documentation

◆ cos_x

real(kind=default_precision), dimension(:), allocatable fftsolver_mod::cos_x

private

Definition at line 21 of file fftsolver.F90.

21 real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: cos_x, cos_y

◆ cos_y

real(kind=default_precision), dimension(:), allocatable fftsolver_mod::cos_y

private

Definition at line 21 of file fftsolver.F90.

◆ fourier_space_sizes

integer, dimension(3) fftsolver_mod::fourier_space_sizes

private

Definition at line 24 of file fftsolver.F90.

24 integer :: fourier_space_sizes(3)

◆ halo_swap_state

type(halo_communication_type), save fftsolver_mod::halo_swap_state

private

Definition at line 25 of file fftsolver.F90.

25 type(halo_communication_type), save :: halo_swap_state

◆ pi

real(kind=default_precision) fftsolver_mod::pi

private

Definition at line 22 of file fftsolver.F90.

22 real(kind=DEFAULT_PRECISION) :: pi

◆ pt

real(kind=default_precision), dimension(:,:,:), allocatable fftsolver_mod::pt

private

Definition at line 23 of file fftsolver.F90.

23 real(kind=DEFAULT_PRECISION), dimension(:,:,:), allocatable :: pt

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ complete_psrce_calculation()

◆ copy_halo_buffer_to_p()

◆ copy_p_to_halo_buffer()

◆ fftsolver_get_descriptor()

◆ finalisation_callback()

◆ initialisation_callback()

◆ perform_local_data_copy_for_p()

◆ timestep_callback()

◆ tridiagonal_solver()

Variable Documentation

◆ cos_x

◆ cos_y

◆ fourier_space_sizes

◆ halo_swap_state

◆ pi

◆ pt