iterativesolver.F90

Go to the documentation of this file.

 
module iterativesolver_mod
  use monc_component_mod, only : component_descriptor_type
  use collections_mod, only : map_type
  use optionsdatabase_mod, only : options_get_real, options_get_integer, options_get_logical
  use grids_mod, only : z_index, y_index, x_index, local_grid_type, grid_configuration_type
  use state_mod, only : model_state_type
  use datadefn_mod, only : default_precision, precision_type
  use logging_mod, only : log_warn, log_debug, log_log, log_get_logging_level
  use conversions_mod, only : conv_to_string
  use communication_types_mod, only : halo_communication_type, neighbour_description_type, field_data_wrapper_type
  use halo_communication_mod, only : copy_buffer_to_field, copy_field_to_buffer, perform_local_data_copy_for_field, &
       init_halo_communication, finalise_halo_communication, initiate_nonblocking_halo_swap, complete_nonblocking_halo_swap, &
       blocking_halo_swap, get_single_field_per_halo_cell
  use registry_mod, only : is_component_enabled
  use logging_mod, only : log_error, log_master_log
  use mpi, only : mpi_max, mpi_sum, mpi_comm_world, mpi_request_null, mpi_statuses_ignore
  implicit none
 
#ifndef TEST_MODE
  private
#endif
 
  type matrix_type
     real(kind=default_precision) :: n, s, e, w
     real(kind=default_precision), dimension(:), allocatable :: u, d, p, lu_d, lu_u, vol
  end type matrix_type
 
  real(kind=default_precision) :: tolerance, relaxation     
  integer :: max_iterations, &     !< Maximum number of BiCGStab iterations
       preconditioner_iterations
  logical :: symm_prob
  
  real(kind=default_precision), parameter :: tiny = 1.0e-16 
 
  type(halo_communication_type), save :: halo_swap_state
  real(kind=default_precision), dimension(:,:,:), allocatable :: psource, prev_p      
  logical :: first_run=.true.
  type(matrix_type) :: a
 
  public iterativesolver_get_descriptor
contains
 
  type(component_descriptor_type) function iterativesolver_get_descriptor()
    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
  end function iterativesolver_get_descriptor
 
  subroutine initialisation_callback(current_state)
    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))
  end subroutine initialisation_callback
 
  subroutine timestep_callback(current_state)
    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
  end subroutine timestep_callback
 
  subroutine finalisation_callback(current_state)
    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)
  end subroutine finalisation_callback  
 
  subroutine bicgstab(current_state, A, x, b, i_strt, i_end, j_strt, j_end, k_end)
    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
  end subroutine bicgstab
 
  subroutine cg_solver(current_state, A, x, b, i_strt, i_end, j_strt, j_end, k_end)
    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
   end subroutine cg_solver
 
  subroutine precond(current_state, A, s, r, preits)
    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
  end subroutine precond
 
  subroutine calc_ax(current_state, A, x, Ax)
    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
  end subroutine calc_ax
 
  real(kind=default_precision) function inner_prod(current_state, x, y, i_strt, i_end, j_strt, j_end, k_end)
    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
  end function inner_prod
 
  function inner_prod_three_way(current_state, t, s, i_strt, i_end, j_strt, j_end, k_end)
    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
  end function inner_prod_three_way 
 
  subroutine set_matrix_for_poisson(grid_configuration, A, z_size)
    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
  end subroutine set_matrix_for_poisson
 
  subroutine deduce_global_divmax(current_state)
    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)
  end subroutine deduce_global_divmax   
 
  subroutine copy_p_to_halo_buffer(current_state, neighbour_description, dim, source_index, &
       pid_location, current_page, source_data)
    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
  end subroutine copy_p_to_halo_buffer
 
  subroutine copy_calc_ax_to_halo_buffer(current_state, neighbour_description, dim, source_index, &
       pid_location, current_page, source_data)
    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
  end subroutine copy_calc_ax_to_halo_buffer
 
  subroutine copy_halo_buffer_to_p(current_state, neighbour_description, dim, target_index, &
       neighbour_location, current_page, source_data)
    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
  end subroutine copy_halo_buffer_to_p
 
  subroutine copy_halo_buffer_to_calc_ax(current_state, neighbour_description, dim, target_index, &
       neighbour_location, current_page, source_data)
    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
  end subroutine copy_halo_buffer_to_calc_ax
 
  subroutine perform_local_data_copy_for_p(current_state, halo_depth, involve_corners, source_data)
    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)
  end subroutine perform_local_data_copy_for_p
 
  subroutine perform_local_data_copy_for_calc_ax(current_state, halo_depth, involve_corners, source_data)
    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)
  end subroutine perform_local_data_copy_for_calc_ax
 
  function create_problem_matrix(z_size)
    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))
  end function create_problem_matrix
 
  subroutine complete_psrce_calculation(current_state, y_halo_size, x_halo_size)
    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)
  end subroutine complete_psrce_calculation 
end module iterativesolver_mod