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

---

Generated by   1.8.14