This is the iterative pressure solver and uses a Jacobi preconditioned BiCGStab which we implement here. More...

Data Types
type	matrix_type
	A helper type to abstract the concrete details of the matrix. More...

Functions/Subroutines
type(component_descriptor_type) function, public	iterativesolver_get_descriptor ()
	Descriptor of the iterative solver component used by the registry. More...

subroutine	initialisation_callback (current_state)
	Initialisation callback hook which will set up the halo swapping state and allocate some data. More...

subroutine	timestep_callback (current_state)
	Timestep callback, this ignores all but the last column where it calls the solver. More...

subroutine	finalisation_callback (current_state)
	Called as MONC is shutting down and frees the halo swap state and deallocates local data. More...

subroutine	bicgstab (current_state, A, x, b, i_strt, i_end, j_strt, j_end, k_end)
	Performs the BiCGStab KS method. More...

subroutine	cg_solver (current_state, A, x, b, i_strt, i_end, j_strt, j_end, k_end)
	Performs the preconditioned conjugate gradient method. More...

subroutine	precond (current_state, A, s, r, preits)
	Jacobi preconditioner. More...

subroutine	calc_ax (current_state, A, x, Ax)
	Calculates A * x. More...

real(kind=default_precision) function	inner_prod (current_state, x, y, i_strt, i_end, j_strt, j_end, k_end)
	Returns the global inner product of two vectors, ignoring the halo cells. More...

real(kind=default_precision) function, dimension(3)	inner_prod_three_way (current_state, t, s, i_strt, i_end, j_strt, j_end, k_end)
	Returns the global inner product of a pair of vectors, ignoring the halo cells for three separate pairs. This call is for optimisation to bunch up the comms in a BiCGStab solver per iteration. More...

subroutine	set_matrix_for_poisson (grid_configuration, A, z_size)
	Sets the values of the provided matrix to solve the poisson equation. More...

subroutine	deduce_global_divmax (current_state)
	Determines the global divmax which is written into the current state. 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_calc_ax_to_halo_buffer (current_state, neighbour_description, dim, source_index, pid_location, current_page, source_data)
	Copies the source field data to halo buffers for a specific process in a dimension and halo cell - for the calc_Ax halo swaps. More...

subroutine	copy_halo_buffer_to_p (current_state, neighbour_description, dim, target_index, neighbour_location, current_page, source_data)
	Copies the halo buffer to halo location for the p field. More...

subroutine	copy_halo_buffer_to_calc_ax (current_state, neighbour_description, dim, target_index, neighbour_location, current_page, source_data)
	Copies the halo buffer to halo location for the source field as required in the calc_Ax procedure. More...

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...

subroutine	perform_local_data_copy_for_calc_ax (current_state, halo_depth, involve_corners, source_data)
	Does a local data copy for halo swapping cells with wrap around (to maintain periodic boundary condition) More...

type(matrix_type) function	create_problem_matrix (z_size)
	Creates a problem matrix, allocating the required data based upon the column size. 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...

Variables
real(kind=default_precision)	tolerance

real(kind=default_precision)	relaxation
	Solving tollerance. More...

integer	max_iterations
	Maximum number of BiCGStab iterations. More...

integer	preconditioner_iterations
	Number of preconditioner iterations to perform per call. More...

logical	symm_prob

real(kind=default_precision), parameter	tiny = 1.0e-16
	Minimum residual - if we go below this then something has gone wrong. More...

type(halo_communication_type), save	halo_swap_state
	The halo swap state as initialised by that module. More...

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

real(kind=default_precision), dimension(:,:,:), allocatable	prev_p
	Passed to BiCGStab as the RHS. More...

logical	first_run =.true.

type(matrix_type)	a

Detailed Description

This is the iterative pressure solver and uses a Jacobi preconditioned BiCGStab which we implement here.

Function/Subroutine Documentation

◆ bicgstab()

subroutine iterativesolver_mod::bicgstab	(	type(model_state_type), intent(inout), target	current_state,
		type(matrix_type), intent(inout)	A,
		real(kind=default_precision), dimension(:,:,:), intent(inout)	x,
		real(kind=default_precision), dimension(:,:,:), intent(in)	b,
		integer, intent(in)	i_strt,
		integer, intent(in)	i_end,
		integer, intent(in)	j_strt,
		integer, intent(in)	j_end,
		integer, intent(in)	k_end
	)

private

Performs the BiCGStab KS method.

Parameters

current_state	The current model state
A	The matrix
x	The solution
b	The RHS

Definition at line 135 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     type(matrix_type), intent(inout) :: A    
     real(kind=default_precision), dimension(:,:,:), intent(inout) :: x
     real(kind=default_precision), dimension(:,:,:), intent(in) :: b
     integer, intent(in) :: i_strt, i_end, j_strt, j_end, k_end
  
     integer :: it, i, j, k
     real(kind=default_precision) :: sc_err, alf, omg, nrm, my_rho, bet, tt, ts, ss, err, init_err, inner_prod_results(3)
     real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX) + &
         current_state%local_grid%halo_size(Z_INDEX) * 2, &
         current_state%local_grid%size(Y_INDEX) + current_state%local_grid%halo_size(Y_INDEX) * 2, &
         current_state%local_grid%size(X_INDEX) + current_state%local_grid%halo_size(X_INDEX) * 2) :: ax, r, cr, pp, v, t, s, cs
  
     ! Calculate scale factor for error
     sc_err = sqrt(inner_prod(current_state, b, b, i_strt, i_end, j_strt, j_end, k_end))
     sc_err = max(sc_err, 0.0001_default_precision)
  
     ! Calculate initial residual
     call calc_ax(current_state, a, x, ax)
  
     do i = i_strt, i_end
       do j = j_strt, j_end
         do k = 2, k_end
           r(k,j,i) = b(k,j,i) - ax(k,j,i)
           cr(k,j,i) = r(k,j,i)
         end do
       end do
     end do
  
     my_rho = inner_prod(current_state, r, r, i_strt, i_end, j_strt, j_end, k_end)
     err = sqrt(my_rho)/sc_err
     init_err = err
  
     alf = 1.0_default_precision
     omg = 1.0_default_precision
     nrm = 1.0_default_precision
  
     if (err .ge. tolerance) then
       do it=1, max_iterations
         if (it > 1) my_rho = inner_prod(current_state, r, cr, i_strt, i_end, j_strt, j_end, k_end)
         bet = (my_rho/nrm) * (alf/omg)
         if (it == 1) then
           call precond(current_state, a, pp, r, preconditioner_iterations)
         else
           do i = i_strt, i_end
             do j = j_strt, j_end
               do k = 2, k_end
                 t(k,j,i) = r(k,j,i) - bet*omg*v(k,j,i)
               end do
             end do
           end do
           call precond(current_state, a, s, t, preconditioner_iterations)
           do i = i_strt, i_end
             do j = j_strt, j_end
               do k = 2, k_end
                 pp(k,j,i) = s(k,j,i) + bet*pp(k,j,i)
               end do
             end do
           end do
         end if       
         call calc_ax(current_state, a, pp, v)
         nrm = inner_prod(current_state, cr, v, i_strt, i_end, j_strt, j_end, k_end)
         alf = my_rho / nrm
  
         do i = i_strt, i_end
           do j = j_strt, j_end
             do k = 2, k_end
               s(k,j,i) = r(k,j,i) - alf*v(k,j,i)
             end do
           end do
         end do
  
         call precond(current_state, a, cs, s, preconditioner_iterations)
         call calc_ax(current_state, a, cs, t)
  
         inner_prod_results=inner_prod_three_way(current_state, t, s, i_strt, i_end, j_strt, j_end, k_end)
         tt = inner_prod_results(1)
         ts = inner_prod_results(2)
         ss = inner_prod_results(3)
         omg = ts/tt
         x = x + alf*pp + omg*cs
         do i = i_strt, i_end
           do j = j_strt, j_end
             do k = 2, k_end
               r(k,j,i) = s(k,j,i) - omg*t(k,j,i)
             end do
           end do
         end do
         nrm = my_rho
  
         if (abs(omg) < tiny) then
           call log_log(log_warn, "Convergence problem, omega="//conv_to_string(omg))
         endif
  
         err = sqrt(ss - 2*omg*ts + omg**2 *tt)/sc_err        
         if (err < tolerance) exit        
       end do
     end if
     
     if (err > tolerance) then
       call log_log(log_warn, "Convergence failed, RNorm="//conv_to_string(err, exponent=.true.))
     else if (current_state%parallel%my_rank==0 .and. log_get_logging_level() .eq. log_debug) then
       call log_log(log_debug, "Converged in "//trim(conv_to_string(it))//" iterations with RNorm="//&
            trim(conv_to_string(err, 5, .true.))//" initial norm="//trim(conv_to_string(init_err, 5, .true.)))
     end if
  

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calc_ax()

subroutine iterativesolver_mod::calc_ax	(	type(model_state_type), intent(inout), target	current_state,
		type(matrix_type), intent(in)	A,
		real(kind=default_precision), dimension(:,:,:), intent(inout), target	x,
		real(kind=default_precision), dimension(:,:,:), intent(inout), target	Ax
	)

private

Calculates A * x.

Parameters

current_state	The current model state
A	The matrix
x	Vector to multiply with
Ax	Result of A*x

Definition at line 378 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     type(matrix_type), intent(in) :: A
     real(kind=default_precision), dimension(:,:,:), target, intent(inout) :: x, ax
  
     integer :: i, k, j, n, istart, iend, jstart, jend
     type(field_data_wrapper_type) :: source_data   
  
     source_data%data=>x
  
     call initiate_nonblocking_halo_swap(current_state, halo_swap_state, &
          copy_calc_ax_to_halo_buffer, source_data=(/source_data/))
  
     ax(1,:,:) = 0.0_default_precision
     if (symm_prob) then
       do n=1, 5
         if (n==1) then
           istart=current_state%local_grid%local_domain_start_index(x_index)+1
           iend=current_state%local_grid%local_domain_end_index(x_index)-1
           jstart=current_state%local_grid%local_domain_start_index(y_index)+1
           jend=current_state%local_grid%local_domain_end_index(y_index)-1
         else if (n==2) then
           istart=current_state%local_grid%local_domain_start_index(x_index)
           iend=current_state%local_grid%local_domain_start_index(x_index)
         else if (n==3) then
           istart=current_state%local_grid%local_domain_end_index(x_index)
           iend=current_state%local_grid%local_domain_end_index(x_index)
         else if (n==4) then
           jstart=current_state%local_grid%local_domain_start_index(y_index)
           jend=current_state%local_grid%local_domain_start_index(y_index)
           istart=current_state%local_grid%local_domain_start_index(x_index)
           iend=current_state%local_grid%local_domain_end_index(x_index)
         else if (n==5) then
           jstart=current_state%local_grid%local_domain_end_index(y_index)
           jend=current_state%local_grid%local_domain_end_index(y_index)
         end if
         do i=istart, iend
           do j=jstart, jend
             k=2
             ax(k,j,i)=a%vol(k)*(a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i)))+ a%u(k)*x(k+1,j,i)+a%p(k)*x(k,j,i)
             do k=3,current_state%local_grid%size(z_index)-1
               ax(k,j,i)=a%vol(k)*(a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i)))+&
                    a%u(k)*x(k+1,j,i)+a%d(k)*x(k-1,j,i)+a%p(k)*x(k,j,i)
             end do
             k=current_state%local_grid%size(z_index)
             ax(k,j,i) = a%vol(k)*(a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i)))+ a%d(k)*x(k-1,j,i)+a%p(k)*x(k,j,i)
           end do
         end do
         if (n==1) then
           call complete_nonblocking_halo_swap(current_state, halo_swap_state, perform_local_data_copy_for_calc_ax, &
                copy_halo_buffer_to_calc_ax, source_data=(/source_data/))
         end if
       end do
     else
       do n=1, 5
         if (n==1) then
           istart=current_state%local_grid%local_domain_start_index(x_index)+1
           iend=current_state%local_grid%local_domain_end_index(x_index)-1
           jstart=current_state%local_grid%local_domain_start_index(y_index)+1
           jend=current_state%local_grid%local_domain_end_index(y_index)-1
         else if (n==2) then
           istart=current_state%local_grid%local_domain_start_index(x_index)
           iend=current_state%local_grid%local_domain_start_index(x_index)
         else if (n==3) then
           istart=current_state%local_grid%local_domain_end_index(x_index)
           iend=current_state%local_grid%local_domain_end_index(x_index)
         else if (n==4) then
           jstart=current_state%local_grid%local_domain_start_index(y_index)
           jend=current_state%local_grid%local_domain_start_index(y_index)
           istart=current_state%local_grid%local_domain_start_index(x_index)
           iend=current_state%local_grid%local_domain_end_index(x_index)
         else if (n==5) then
           jstart=current_state%local_grid%local_domain_end_index(y_index)
           jend=current_state%local_grid%local_domain_end_index(y_index)
         end if
         do i=istart, iend
           do j=jstart, jend
             k=2
             ax(k,j,i)=a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i))+ a%u(k)*x(k+1,j,i)+a%p(k)*x(k,j,i)
             do k=3,current_state%local_grid%size(z_index)-1
               ax(k,j,i)=a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i))+&
                    a%u(k)*x(k+1,j,i)+a%d(k)*x(k-1,j,i)+a%p(k)*x(k,j,i)
             end do
             k=current_state%local_grid%size(z_index)
             ax(k,j,i) = a%n*(x(k,j,i+1)+x(k,j,i-1))+a%e*(x(k,j+1,i)+x(k,j-1,i))+ a%d(k)*x(k-1,j,i)+a%p(k)*x(k,j,i)
           end do
         end do
         if (n==1) then
           call complete_nonblocking_halo_swap(current_state, halo_swap_state, perform_local_data_copy_for_calc_ax, &
                copy_halo_buffer_to_calc_ax, source_data=(/source_data/))
         end if
       end do
     endif

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cg_solver()

subroutine iterativesolver_mod::cg_solver	(	type(model_state_type), intent(inout), target	current_state,
		type(matrix_type), intent(inout)	A,
		real(kind=default_precision), dimension(:,:,:), intent(inout)	x,
		real(kind=default_precision), dimension(:,:,:), intent(in)	b,
		integer, intent(in)	i_strt,
		integer, intent(in)	i_end,
		integer, intent(in)	j_strt,
		integer, intent(in)	j_end,
		integer, intent(in)	k_end
	)

private

Performs the preconditioned conjugate gradient method.

Parameters

current_state	The current model state
A	The matrix
x	The solution
b	The RHS

Definition at line 248 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     type(matrix_type), intent(inout) :: A
     real(kind=default_precision), dimension(:,:,:), intent(inout) :: x
     real(kind=default_precision), dimension(:,:,:), intent(in) :: b
     integer, intent(in) :: i_strt, i_end, j_strt, j_end, k_end
  
     integer :: it, k, i, j
     real(kind=default_precision) :: sc_err, alf, bet, err, init_err, rho
     real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX) + &
         current_state%local_grid%halo_size(Z_INDEX) * 2, &
         current_state%local_grid%size(Y_INDEX) + current_state%local_grid%halo_size(Y_INDEX) * 2, &
         current_state%local_grid%size(X_INDEX) + current_state%local_grid%halo_size(X_INDEX) * 2) :: ax, r, z, p
  
     ! first rescale RHS for symmetry (this could be done when p_source is calculated
     do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
       do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
  
         r(1,j,i) = 0.0_default_precision
         do k=2,current_state%local_grid%size(z_index)
            r(k,j,i) = b(k,j,i) * a%vol(k)
         end do
       end do
     end do
  
     ! Calculate scale factor for error
  
     call calc_ax(current_state, a, x, ax)
  
     sc_err = sqrt(inner_prod(current_state, r, r, i_strt, i_end, j_strt, j_end, k_end))
     sc_err = max(sc_err, 0.0001_default_precision)
     r = r - ax
     init_err = sqrt(inner_prod(current_state, r, r, i_strt, i_end, j_strt, j_end, k_end))/sc_err
  
     do it=1, max_iterations
        if( it == 1 ) then
           call precond(current_state, a, p, r, preconditioner_iterations)
           rho = inner_prod(current_state, p, r, i_strt, i_end, j_strt, j_end, k_end)
           alf = rho
        else
           call precond(current_state, a, z, r, preconditioner_iterations)
           alf = inner_prod(current_state, z, r, i_strt, i_end, j_strt, j_end, k_end)
           bet = alf/rho
           rho = alf
           p   = z + bet*p
        end if
  
        call calc_ax(current_state, a, p, ax)
        alf = alf/inner_prod(current_state, p, ax, i_strt, i_end, j_strt, j_end, k_end)
        x   = x + alf*p
        r   = r - alf*ax
  
        err = sqrt(inner_prod(current_state, r, r, i_strt, i_end, j_strt, j_end, k_end))/sc_err
        if (err < tolerance) exit
     end do
  
     if( current_state%parallel%my_rank == 0 ) print*,it, err, init_err
  
     if (err > tolerance) then
       call log_log(log_warn, "Convergence failed, RNorm="//conv_to_string(err, exponent=.true.))
     else if (current_state%parallel%my_rank==0 .and. log_get_logging_level() .eq. log_debug) then
       call log_log(log_debug, "Converged in "//trim(conv_to_string(it))//" iterations with RNorm="//&
            trim(conv_to_string(err, 5, .true.))//" initial norm="//trim(conv_to_string(init_err, 5, .true.)))
     end if
  

Here is the call graph for this function:

Here is the caller graph for this function:

◆ complete_psrce_calculation()

subroutine iterativesolver_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 735 of file iterativesolver.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_calc_ax_to_halo_buffer()

subroutine iterativesolver_mod::copy_calc_ax_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 source field data to halo buffers for a specific process in a dimension and halo cell - for the calc_Ax halo swaps.

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 623 of file iterativesolver.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
  
     type(field_data_wrapper_type) :: selected_source
  
     selected_source=source_data(1)
  
     call copy_field_to_buffer(current_state%local_grid, neighbour_description%send_halo_buffer, selected_source%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:

◆ copy_halo_buffer_to_calc_ax()

subroutine iterativesolver_mod::copy_halo_buffer_to_calc_ax	(	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

Copies the halo buffer to halo location for the source field as required in the calc_Ax procedure.

Parameters

current_state	The current model state
neighbour_description	The halo swapping description of the neighbour we are accessing the buffer of
dim	The dimension we receive for
target_index	The target index for the dimension we are receiving for
neighbour_location	The location in the local neighbour data stores of this neighbour
current_page	The current, next, halo swap page to read from (all previous have been read and copied already)
source_data	Optional source data which is written into

Definition at line 671 of file iterativesolver.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
  
     type(field_data_wrapper_type) :: selected_source
  
     selected_source=source_data(1)
  
     call copy_buffer_to_field(current_state%local_grid, neighbour_description%recv_halo_buffer, selected_source%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_halo_buffer_to_p()

subroutine iterativesolver_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

Copies the halo buffer to halo location for the p field.

Parameters

current_state	The current model state
neighbour_description	The halo swapping description of the neighbour we are accessing the buffer of
dim	The dimension we receive for
target_index	The target index for the dimension we are receiving for
neighbour_location	The location in the local neighbour data stores of this neighbour
current_page	The current, next, halo swap page to read from (all previous have been read and copied already)
source_data	Optional source data which is written into

Definition at line 649 of file iterativesolver.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 iterativesolver_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 601 of file iterativesolver.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:

◆ create_problem_matrix()

type(matrix_type) function iterativesolver_mod::create_problem_matrix ( integer, intent(in) z_size )

private

Creates a problem matrix, allocating the required data based upon the column size.

Parameters

z_size Number of elements in the vertical column

Returns: The allocated matrix ready to be used

Definition at line 722 of file iterativesolver.F90.

     integer, intent(in) :: z_size
     type(matrix_type) :: create_problem_matrix
  
     allocate(create_problem_matrix%u(z_size), create_problem_matrix%d(z_size), create_problem_matrix%p(z_size), &
          create_problem_matrix%lu_u(z_size), create_problem_matrix%lu_d(z_size), create_problem_matrix%vol(z_size))

Here is the caller graph for this function:

◆ deduce_global_divmax()

subroutine iterativesolver_mod::deduce_global_divmax ( type(model_state_type), intent(inout), target current_state )

private

Determines the global divmax which is written into the current state.

Parameters

current_state The current model state

Definition at line 584 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     
     integer :: ierr
     
     call mpi_allreduce(current_state%local_divmax, current_state%global_divmax, 1, precision_type, mpi_max, &
          current_state%parallel%monc_communicator, ierr)

Here is the caller graph for this function:

◆ finalisation_callback()

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

private

Called as MONC is shutting down and frees the halo swap state and deallocates local data.

Parameters

current_state The current model state

Definition at line 123 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
  
     call finalise_halo_communication(halo_swap_state)
     deallocate(psource, prev_p, a%u, a%d, a%p, a%lu_u, a%lu_d, a%vol)

Here is the caller graph for this function:

◆ initialisation_callback()

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

private

Initialisation callback hook which will set up the halo swapping state and allocate some data.

Parameters

current_state The current model state

Definition at line 57 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
  
     if (.not. is_component_enabled(current_state%options_database, "diverr")) then
       call log_master_log(log_error, "The iterative solver component requires the diverr component to be enabled")
     end if
  
     tolerance=options_get_real(current_state%options_database, "tolerance")
     max_iterations=options_get_integer(current_state%options_database, "max_iterations")
     preconditioner_iterations=options_get_integer(current_state%options_database, "preconditioner_iterations")
     symm_prob=options_get_logical(current_state%options_database, "symm_prob")
  
     call init_halo_communication(current_state, get_single_field_per_halo_cell, halo_swap_state, 1, .false.)
  
     allocate(psource(current_state%local_grid%size(z_index) + current_state%local_grid%halo_size(z_index) * 2, &
          current_state%local_grid%size(y_index) + current_state%local_grid%halo_size(y_index) * 2, &
          current_state%local_grid%size(x_index) + current_state%local_grid%halo_size(x_index) * 2),&
          prev_p(current_state%local_grid%size(z_index) + current_state%local_grid%halo_size(z_index) * 2, &
          current_state%local_grid%size(y_index) + current_state%local_grid%halo_size(y_index) * 2, &
          current_state%local_grid%size(x_index) + current_state%local_grid%halo_size(x_index) * 2))
  
     a=create_problem_matrix(current_state%local_grid%size(z_index))
     call set_matrix_for_poisson(current_state%global_grid%configuration, a, current_state%local_grid%size(z_index))

Here is the call graph for this function:

Here is the caller graph for this function:

◆ inner_prod()

real(kind=default_precision) function iterativesolver_mod::inner_prod	(	type(model_state_type), intent(inout), target	current_state,
		real(kind=default_precision), dimension(:,:,:), intent(in)	x,
		real(kind=default_precision), dimension(:,:,:), intent(in)	y,
		integer, intent(in)	i_strt,
		integer, intent(in)	i_end,
		integer, intent(in)	j_strt,
		integer, intent(in)	j_end,
		integer, intent(in)	k_end
	)

private

Returns the global inner product of two vectors, ignoring the halo cells.

Parameters

current_state	The current model state
x	First vector @praam y Second vector

Returns: Global inner product of the two input vectors

Definition at line 478 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     real(kind=default_precision), dimension(:,:,:), intent(in) :: x, y
     integer, intent(in) :: i_strt, i_end, j_strt, j_end, k_end
  
     real(kind=default_precision) :: local_sum, global_sum
     integer :: ierr, i, j, k
  
     local_sum=0.0_default_precision
  
      do i=i_strt, i_end
       do j=j_strt, j_end
         do k=2, k_end
           local_sum=local_sum+x(k,j,i)*y(k,j,i)
         end do        
       end do      
     end do
  
     call mpi_allreduce(local_sum, global_sum, 1, precision_type, mpi_sum, current_state%parallel%monc_communicator, ierr)
     inner_prod=global_sum

Here is the call graph for this function:

Here is the caller graph for this function:

◆ inner_prod_three_way()

real(kind=default_precision) function, dimension(3) iterativesolver_mod::inner_prod_three_way	(	type(model_state_type), intent(inout), target	current_state,
		real(kind=default_precision), dimension(:,:,:), intent(in)	t,
		real(kind=default_precision), dimension(:,:,:), intent(in)	s,
		integer, intent(in)	i_strt,
		integer, intent(in)	i_end,
		integer, intent(in)	j_strt,
		integer, intent(in)	j_end,
		integer, intent(in)	k_end
	)

private

Returns the global inner product of a pair of vectors, ignoring the halo cells for three separate pairs. This call is for optimisation to bunch up the comms in a BiCGStab solver per iteration.

Parameters

current_state	The current model state
x	First vector @praam y Second vector

Returns: Global inner product of the two input vectors

Definition at line 506 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     integer, intent(in) :: i_strt, i_end, j_strt, j_end, k_end
     real(kind=default_precision), dimension(:,:,:), intent(in) :: t, s
     real(kind=default_precision), dimension(3) :: inner_prod_three_way
  
     real(kind=default_precision), dimension(3) :: local_sum, global_sum
     integer :: ierr, i, j, k
  
     local_sum(1)=0.0_default_precision
     local_sum(2)=0.0_default_precision
     local_sum(3)=0.0_default_precision
  
     do i=i_strt, i_end
       do j=j_strt, j_end
         do k=2, k_end
           local_sum(1)=local_sum(1)+t(k,j,i)*t(k,j,i)
           local_sum(2)=local_sum(2)+t(k,j,i)*s(k,j,i)
           local_sum(3)=local_sum(3)+s(k,j,i)*s(k,j,i)
         end do
       end do
     end do
  
     call mpi_allreduce(local_sum, global_sum, 3, precision_type, mpi_sum, current_state%parallel%monc_communicator, ierr)
     inner_prod_three_way=global_sum

Here is the caller graph for this function:

◆ iterativesolver_get_descriptor()

type(component_descriptor_type) function, public iterativesolver_mod::iterativesolver_get_descriptor

Descriptor of the iterative solver component used by the registry.

Returns: The iterative solver component descriptor

Definition at line 47 of file iterativesolver.F90.

     iterativesolver_get_descriptor%name="iterativesolver"
     iterativesolver_get_descriptor%version=0.1
     iterativesolver_get_descriptor%initialisation=>initialisation_callback
     iterativesolver_get_descriptor%timestep=>timestep_callback
     iterativesolver_get_descriptor%finalisation=>finalisation_callback

Here is the call graph for this function:

◆ perform_local_data_copy_for_calc_ax()

subroutine iterativesolver_mod::perform_local_data_copy_for_calc_ax	(	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 a local data copy for halo swapping cells with wrap around (to maintain periodic boundary condition)

Parameters

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

Definition at line 705 of file iterativesolver.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
  
     type(field_data_wrapper_type) :: selected_source
  
     selected_source=source_data(1)
  
     call perform_local_data_copy_for_field(selected_source%data, current_state%local_grid, &
          current_state%parallel%my_rank, halo_depth, involve_corners)

Here is the caller graph for this function:

◆ perform_local_data_copy_for_p()

subroutine iterativesolver_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 692 of file iterativesolver.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:

◆ precond()

subroutine iterativesolver_mod::precond	(	type(model_state_type), intent(inout), target	current_state,
		type(matrix_type), intent(inout)	A,
		real(kind=default_precision), dimension(:,:,:), intent(inout)	s,
		real(kind=default_precision), dimension(:,:,:), intent(in)	r,
		integer, intent(in)	preits
	)

private

Jacobi preconditioner.

Parameters

current_state	The current model state
A	The matrix
s	Written into as result of preconditioning
r	Input values to preconditioner
preits	Number of iterations of the preconditioner to perform per call

Definition at line 320 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
     real(kind=default_precision), dimension(:,:,:), intent(in) :: r
     real(kind=default_precision), dimension(:,:,:), intent(inout) :: s
     integer, intent(in) :: preits
     type(matrix_type), intent(inout) :: A
  
     real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX) + &
         current_state%local_grid%halo_size(Z_INDEX) * 2, current_state%local_grid%size(Y_INDEX) + &
         current_state%local_grid%halo_size(Y_INDEX) * 2, current_state%local_grid%size(X_INDEX) + &
         current_state%local_grid%halo_size(X_INDEX) * 2) :: t
     real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX)) :: s_k
     integer :: it, i, j, k
  
     if (preits .lt. 0) then
       s=r
       return
     end if    
  
     do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
       do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
         s(1,j,i) = 0.0_default_precision
         k=2 
         s(k,j,i)=r(k,j,i)*a%LU_d(k)
         do k=3,current_state%local_grid%size(z_index)
           s(k,j,i)=(r(k,j,i) - a%d(k)*s(k-1,j,i))*a%lu_d(k)
         end do
         do k=current_state%local_grid%size(z_index)-1, 2, -1
           s(k,j,i)=s(k,j,i) - a%lu_u(k)*s(k+1,j,i)
         end do        
       end do
     end do
  
     do it=1, preits
       call calc_ax(current_state, a, s, t)
       do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
       do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
           k=2 
           s_k(k)=(r(k,j,i) - t(k,j,i))*a%lu_d(k) 
           do k=3,current_state%local_grid%size(z_index) 
             s_k(k)=(r(k,j,i) - t(k,j,i) - a%d(k)*s_k(k-1))*a%lu_d(k) 
           end do
           k=current_state%local_grid%size(z_index)
           s(k,j,i)=s(k,j,i)+s_k(k) 
           do k=current_state%local_grid%size(z_index)-1, 2, -1 
             s_k(k)=s_k(k) - a%lu_u(k)*s_k(k+1) 
             s(k,j,i)=s(k,j,i) + relaxation*s_k(k) 
           end do
         end do
       end do
     end do
  

Here is the call graph for this function:

Here is the caller graph for this function:

◆ set_matrix_for_poisson()

subroutine iterativesolver_mod::set_matrix_for_poisson	(	type(grid_configuration_type), intent(inout)	grid_configuration,
		type(matrix_type), intent(inout)	A,
		integer, intent(in)	z_size
	)

private

Sets the values of the provided matrix to solve the poisson equation.

Parameters

grid_configuration	Configuration of the vertical and horizontal grids
A	The matrix that the values are written into
z_size	Number of elements in a column

Definition at line 537 of file iterativesolver.F90.

     type(grid_configuration_type), intent(inout) :: grid_configuration
     type(matrix_type), intent(inout) :: A
     integer, intent(in) :: z_size    
  
     integer :: k
     real(kind=default_precision) :: d_sc, concat_scalars    
  
     a%n=grid_configuration%horizontal%cx*grid_configuration%horizontal%cx
     a%s=a%n
     a%e=grid_configuration%horizontal%cy*grid_configuration%horizontal%cy
     a%w=a%e
     concat_scalars=a%n+a%s+a%e+a%w
     do k=2, z_size      
       if (symm_prob) then
          a%vol(k)=grid_configuration%vertical%dz(k)
          d_sc=1.0/grid_configuration%vertical%rhon(k)
       else
          d_sc=grid_configuration%vertical%rdz(k) / grid_configuration%vertical%rhon(k)
          a%vol(k)=1.0
       endif
  
       if (k==z_size) then
         a%u(k)=0.0_default_precision
       else
         a%u(k)=grid_configuration%vertical%rho(k)*grid_configuration%vertical%rdzn(k+1)
       end if
       if (k==2) then
         a%d(k)=0.0_default_precision
       else
         a%d(k)=grid_configuration%vertical%rho(k-1)*grid_configuration%vertical%rdzn(k)
       end if
       a%p(k) = d_sc * (-(a%u(k) + a%d(k))) - concat_scalars * a%vol(k)
       a%u(k)=d_sc * a%u(k)
       a%d(k)=d_sc * a%d(k)
     end do
      k=2 
      a%lu_d(k)=1.0_default_precision/a%p(k) 
      a%lu_u(k)=a%lu_d(k)*a%u(k) 
      do k=3, z_size 
        a%lu_d(k)=1.0_default_precision/(a%p(k) - a%d(k)*a%lu_u(k-1))
        a%lu_u(k)=a%u(k)*a%lu_d(k) 
      end do

Here is the caller graph for this function:

◆ timestep_callback()

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

private

Timestep callback, this ignores all but the last column where it calls the solver.

Parameters

current_state The current model state

Definition at line 84 of file iterativesolver.F90.

     type(model_state_type), target, intent(inout) :: current_state
  
     integer :: i_strt, i_end, j_strt, j_end, k_end
  
     i_strt = current_state%local_grid%local_domain_start_index(x_index)
     i_end  = current_state%local_grid%local_domain_end_index(x_index)
     j_strt = current_state%local_grid%local_domain_start_index(y_index)
     j_end  = current_state%local_grid%local_domain_end_index(y_index)
     k_end  = current_state%local_grid%size(z_index)
  
     call complete_psrce_calculation(current_state, current_state%local_grid%halo_size(y_index), &
          current_state%local_grid%halo_size(x_index))
      
     call initiate_nonblocking_halo_swap(current_state, halo_swap_state, copy_p_to_halo_buffer)
     call deduce_global_divmax(current_state)    
     call complete_nonblocking_halo_swap(current_state, halo_swap_state, perform_local_data_copy_for_p, copy_halo_buffer_to_p)
  
     psource=current_state%p%data
     if (first_run) then
       ! If first timestep then initial guess is zero
       current_state%p%data=0.0_default_precision
       first_run=.false.
     else
       ! Initial guess is set to previous timesteps p
       current_state%p%data=prev_p
     end if
     
     if (symm_prob) then
       call cg_solver(current_state, a, current_state%p%data, psource, i_strt, i_end, j_strt, j_end, k_end)
     else
       call bicgstab(current_state, a, current_state%p%data, psource, i_strt, i_end, j_strt, j_end, k_end)
     end if
  
     prev_p=current_state%p%data

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

◆ a

type(matrix_type) iterativesolver_mod::a

private

Definition at line 40 of file iterativesolver.F90.

40 type(matrix_type) :: A

◆ first_run

logical iterativesolver_mod::first_run =.true.

private

Definition at line 39 of file iterativesolver.F90.

39 logical :: first_run=.true.

◆ halo_swap_state

type(halo_communication_type), save iterativesolver_mod::halo_swap_state

private

The halo swap state as initialised by that module.

Definition at line 37 of file iterativesolver.F90.

37 type(halo_communication_type), save :: halo_swap_state

◆ max_iterations

integer iterativesolver_mod::max_iterations

private

Maximum number of BiCGStab iterations.

Definition at line 31 of file iterativesolver.F90.

31 integer :: max_iterations, & !< Maximum number of BiCGStab iterations

32 preconditioner_iterations

◆ preconditioner_iterations

integer iterativesolver_mod::preconditioner_iterations

private

Number of preconditioner iterations to perform per call.

Definition at line 31 of file iterativesolver.F90.

◆ prev_p

real(kind=default_precision), dimension(:,:,:), allocatable iterativesolver_mod::prev_p

private

Passed to BiCGStab as the RHS.

Definition at line 38 of file iterativesolver.F90.

◆ psource

real(kind=default_precision), dimension(:,:,:), allocatable iterativesolver_mod::psource

private

Definition at line 38 of file iterativesolver.F90.

38 real(kind=default_precision), dimension(:,:,:), allocatable :: psource, prev_p

◆ relaxation

real(kind=default_precision) iterativesolver_mod::relaxation

private

Solving tollerance.

Definition at line 30 of file iterativesolver.F90.

◆ symm_prob

logical iterativesolver_mod::symm_prob

private

Definition at line 33 of file iterativesolver.F90.

33 logical :: symm_prob

◆ tiny

real(kind=default_precision), parameter iterativesolver_mod::tiny = 1.0e-16

private

Minimum residual - if we go below this then something has gone wrong.

Definition at line 35 of file iterativesolver.F90.

35 real(kind=default_precision), parameter :: tiny = 1.0e-16

◆ tolerance

real(kind=default_precision) iterativesolver_mod::tolerance

private

Definition at line 30 of file iterativesolver.F90.

30 real(kind=default_precision) :: tolerance, relaxation

Data Types

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ bicgstab()

◆ calc_ax()

◆ cg_solver()

◆ complete_psrce_calculation()

◆ copy_calc_ax_to_halo_buffer()

◆ copy_halo_buffer_to_calc_ax()

◆ copy_halo_buffer_to_p()

◆ copy_p_to_halo_buffer()

◆ create_problem_matrix()

◆ deduce_global_divmax()

◆ finalisation_callback()

◆ initialisation_callback()

◆ inner_prod()

◆ inner_prod_three_way()

◆ iterativesolver_get_descriptor()

◆ perform_local_data_copy_for_calc_ax()

◆ perform_local_data_copy_for_p()

◆ precond()

◆ set_matrix_for_poisson()

◆ timestep_callback()

Variable Documentation

◆ a

◆ first_run

◆ halo_swap_state

◆ max_iterations

◆ preconditioner_iterations

◆ prev_p

◆ psource

◆ relaxation

◆ symm_prob

◆ tiny

◆ tolerance