scale_multigrid_smoother_base.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Multigrid / Smoother base
!!
!! @par Description
!!      A module to provide derived type to manage multigrid smoother
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_multigrid_smoother_base
   !-----------------------------------------------------------------------------
   !
   !++ used modules
   !
   !
   use scale_precision
   use scale_io
   use scale_prc, only: prc_abort
 
   use scale_sparsemat, only: sparsemat
 
   use scale_element_quadrilateral, only: quadrilateralelement
   use scale_element_hexahedral, only: hexahedralelement
 
   use scale_mesh_base2d, only: meshbase2d
   use scale_localmesh_2d, only: localmesh2d
   use scale_mesh_base3d, only: meshbase3d
   use scale_localmesh_3d, only: localmesh3d
   use scale_mesh_rectdom2d, only: meshrectdom2d
   use scale_mesh_cubedom3d, only: meshcubedom3d
 
   use scale_meshfield_base, only: &
     meshfield2d, meshfield3d
   use scale_meshfieldcomm_base, only: meshfieldcommbase
   use scale_meshfieldcomm_rectdom2d, only: meshfieldcommrectdom2d
   use scale_meshfieldcomm_cubedom3d, only: meshfieldcommcubedom3d
   !-----------------------------------------------------------------------------
   implicit none
   private
 
   !-----------------------------------------------------------------------------
   !
   !++ Public type & procedure
   ! 
   integer, public, parameter :: mgsmoother_pre_id  = 1  !< ID to represent pre-smoothing
   integer, public, parameter :: mgsmoother_post_id = 2  !< ID to represent post-smoothing
 
   !> Base type for multigrid smoother
   type, public :: mgsmootherbase
      integer :: itration           !< Current smoothing iteration
      integer :: num_smooth_itr_max !< Maximum number of smoothing iterations
 
      real(rp) :: threshold_ratio_residual_l2  !< Threshold of decreasing ratio of residual for smoothing (L2 norm)
      real(rp) :: threshold_residual_l2  !< Residual threshold for smoothing (L2 norm)
      real(rp) :: threshold_residual_max !< Residual threshold for smoothing (max norm)
 
      integer :: history_residual_step                 !< Step recording residual history
      integer :: history_residual_rstep                !< Remaining step to record residual history
      integer :: history_residual_count                !< Count of recorded residuals
      real(rp), allocatable :: history_residual_l2(:)  !< Residual history of L2 norm
      real(rp) :: history_residual_l2_initial          !< Initial residual L2 norm
      real(rp), allocatable :: history_residual_max(:) !< Residual history of max norm
      real(rp) :: history_residual_max_initial         !< Initial residual max norm
   contains
      procedure :: set_smoothing_param => mgsmootherbase_set_smoothing_param
      procedure :: output_residual_history => mgsmootherbase_output_residual_history
      procedure :: is_converged => mgsmootherbase_is_converged
      procedure :: get_initial_residual_statistics => mgsmootherbase_get_initial_residual_statistics
      procedure :: get_current_residual_statistics => mgsmootherbase_get_current_residual_statistics
      procedure, private :: is_residual_check_step => mgsmootherbase_is_residual_check_step
      procedure, private :: eval_global_statistic => mgsmootherbase_eval_global_statistic
      
   end type mgsmootherbase
   public :: mgsmootherbase_init
   public :: mgsmootherbase_final
   public :: mgsmootherbase_put_residual_history
   public :: mgsmootherbase_get_statistic_lc
 
   integer, public, parameter :: mgsmoother_default_num_smooth_ite_max = 10
   real(rp), public, parameter :: mgsmoother_default_threshold_ratio_residual_l2 = 1.0e-3_rp
   real(rp), public, parameter :: mgsmoother_default_threshold_residual_l2  = 1.0e-3_rp
   real(rp), public, parameter :: mgsmoother_default_threshold_residual_max = 1.0e-3_rp
 
 
   !> Derived type for 2D multigrid smoother
   type, abstract, extends(mgsmootherbase), public :: mgsmootherbase2d
   contains
     procedure :: do_smoothing => mgsmootherbase2d_do_smoothing
     procedure(mgsmootherbase2d_advance_itr_1step), deferred :: advance_itr_1step
     procedure :: eval_residual_statistic => mgsmootherbase2d_eval_residual_statistic
     procedure :: set_initial_residual => mgsmootherbase2d_set_initial_residual
     procedure, private :: put_residual => mgsmootherbase2d_put_residual
   end type mgsmootherbase2d
   public :: mgsmootherbase2d_init
   public :: mgsmootherbase2d_final
 
   abstract interface
    subroutine mgsmootherbase2d_advance_itr_1step( this, q, res, &
      f, aux_var, itr, var_comm, aux_comm, Dx, Dy, Lift, mesh2D, &
      cal_res_flag, zero_initial_guess, mg_p_level, mg_h_level,  &
      pre_or_post_smooth )
      import meshbase2d
      import :: meshfield2d
      import meshfieldcommbase
      import sparsemat
      import :: mgsmootherbase2d
      class(mgsmootherbase2d), intent(inout) :: this
      type(meshfield2d), intent(inout), target :: q
      type(meshfield2d), intent(inout) :: res
      type(meshfield2d), intent(in) :: f
      type(meshfield2d), intent(inout), target :: aux_var(:)
      integer, intent(in) :: itr
      class(meshfieldcommbase), intent(inout) :: var_comm
      class(meshfieldcommbase), intent(inout) :: aux_comm
      type(sparsemat), intent(in) :: Dx
      type(sparsemat), intent(in) :: Dy
      type(sparsemat), intent(in) :: Lift
      class(meshbase2d), intent(in), target :: mesh2D
      logical, intent(in) :: cal_res_flag
      logical, intent(in) :: zero_initial_guess
      integer, intent(in) :: mg_p_level
      integer, intent(in) :: mg_h_level
      integer, intent(in) :: pre_or_post_smooth
    end subroutine mgsmootherbase2d_advance_itr_1step
   end interface  
 
   !> Derived type for 3D multigrid smoother
   type, abstract, extends(mgsmootherbase), public :: mgsmootherbase3d
   contains
    procedure, public :: do_smoothing => mgsmootherbase3d_do_smoothing
    procedure(mgsmootherbase3d_advance_itr_1step), deferred :: advance_itr_1step
    procedure :: eval_residual_statistic => mgsmootherbase3d_eval_residual_statistic
    procedure :: set_initial_residual => mgsmootherbase3d_set_initial_residual
    procedure, private :: put_residual => mgsmootherbase3d_put_residual
   end type mgsmootherbase3d
   public :: mgsmootherbase3d_init
   public :: mgsmootherbase3d_final
 
  abstract interface
    subroutine mgsmootherbase3d_advance_itr_1step( this, q, res,      &
      f, aux_var, itr, var_comm, aux_comm, Dx, Dy, Dz, Lift, mesh3D,  &
      cal_res_flag, zero_initial_guess, mg_p_level, mg_h_level,       &
      pre_or_post_smooth )
      import meshbase3d
      import :: meshfield3d
      import meshfieldcommbase
      import sparsemat
      import :: mgsmootherbase3d
      class(mgsmootherbase3d), intent(inout) :: this
      type(meshfield3d), intent(inout), target :: q
      type(meshfield3d), intent(inout) :: res
      type(meshfield3d), intent(in) :: f
      type(meshfield3d), intent(inout), target :: aux_var(:)
      integer, intent(in) :: itr
      class(meshfieldcommbase), intent(inout) :: var_comm
      class(meshfieldcommbase), intent(inout) :: aux_comm
      type(sparsemat), intent(in) :: Dx
      type(sparsemat), intent(in) :: Dy
      type(sparsemat), intent(in) :: Dz
      type(sparsemat), intent(in) :: Lift
      class(meshbase3d), intent(in), target :: mesh3D
      logical, intent(in) :: cal_res_flag
      logical, intent(in) :: zero_initial_guess
      integer, intent(in) :: mg_p_level
      integer, intent(in) :: mg_h_level
      integer, intent(in) :: pre_or_post_smooth
    end subroutine mgsmootherbase3d_advance_itr_1step
  end interface  
 
contains
!-- Base
!OCL SERIAL
  subroutine mgsmootherbase_init(this)
    implicit none
    class(mgsmootherbase), intent(inout) :: this
    !---------------------------------
    this%num_smooth_itr_max  = mgsmoother_default_num_smooth_ite_max
    this%threshold_ratio_residual_l2 = mgsmoother_default_threshold_ratio_residual_l2
    this%threshold_residual_l2  = mgsmoother_default_threshold_residual_l2
    this%threshold_residual_max = mgsmoother_default_threshold_residual_max
    return
  end subroutine mgsmootherbase_init
 
!OCL SERIAL
  subroutine mgsmootherbase_final(this)
    implicit none
    class(mgsmootherbase), intent(inout) :: this
    !---------------------------------
    return
  end subroutine mgsmootherbase_final
 
  !> Set smoothing parameters
!OCL SERIAL
  subroutine mgsmootherbase_set_smoothing_param( this, &
    num_smooth_itr_max,       &
    threshold_ratio_residual_l2, &
    threshold_residual_l2,    &
    threshold_residual_max,   &
    history_residual_step     )
    implicit none
    class(mgsmootherbase), intent(inout) :: this
    integer, intent(in) :: num_smooth_itr_max      !< Maximum number of smoothing iterations
    real(rp), intent(in) :: threshold_ratio_residual_l2  !< Threshold of decreasing ratio of residual for smoothing (L2 norm)
    real(rp), intent(in) :: threshold_residual_l2  !< Residual threshold for smoothing (L2 norm)
    real(rp), intent(in) :: threshold_residual_max !< Residual threshold for smoothing (max norm)
    integer, intent(in) :: history_residual_step   !< Step recording residual history
 
    integer :: num_history_residual
    !---------------------------------
 
    this%num_smooth_itr_max = num_smooth_itr_max
    this%threshold_ratio_residual_l2 = threshold_ratio_residual_l2
    this%threshold_residual_l2  = threshold_residual_l2
    this%threshold_residual_max = threshold_residual_max 
 
    !-
    this%history_residual_step = history_residual_step
    num_history_residual = this%num_smooth_itr_max / this%history_residual_step + 1
 
    allocate( this%history_residual_l2(num_history_residual) )
    allocate( this%history_residual_max(num_history_residual) )
 
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "  Smoothing parameters set:"
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    num_smooth_itr_max          =", this%num_smooth_itr_max
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    threshold_ratio_residual_l2 =", this%threshold_ratio_residual_l2
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    threshold_residual_l2       =", this%threshold_residual_l2
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    threshold_residual_max      =", this%threshold_residual_max
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    history_residual_step       =", this%history_residual_step
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "    num of history residual     =", num_history_residual
    log_info( 'MGSmootherBase_set_smoothing_param',* ) "----------------------------------------"
    return
  end subroutine mgsmootherbase_set_smoothing_param
 
!OCL SERIAL
  function mgsmootherbase_is_residual_check_step( this ) result( is_check_step )
    implicit none
    class(mgsmootherbase), intent(in) :: this
    logical :: is_check_step
    !---------------------------------
 
    if ( this%history_residual_rstep == 1 ) then
      is_check_step = .true.
    else
      is_check_step = .false.
    end if
    return
  end function mgsmootherbase_is_residual_check_step
 
!OCL SERIAL
  subroutine mgsmootherbase_get_initial_residual_statistics( this, res_l2_init, res_max_init )
    implicit none
    class(mgsmootherbase), intent(in) :: this
    real(rp), intent(out) :: res_l2_init
    real(rp), intent(out) :: res_max_init
    !---------------------------------
 
    res_l2_init = this%history_residual_l2_initial
    res_max_init = this%history_residual_max_initial
    return
  end subroutine mgsmootherbase_get_initial_residual_statistics
 
!OCL SERIAL
  subroutine mgsmootherbase_get_current_residual_statistics( this, res_l2, res_max )
    implicit none
    class(mgsmootherbase), intent(in) :: this
    real(rp), intent(out) :: res_l2
    real(rp), intent(out) :: res_max
    !---------------------------------
 
    res_l2 = this%history_residual_l2( this%history_residual_count )
    res_max = this%history_residual_max( this%history_residual_count )
    return
  end subroutine mgsmootherbase_get_current_residual_statistics
 
!OCL SERIAL
  function mgsmootherbase_is_converged( this, is_Vcycle_end ) result( conv_flag )
    implicit none
    class(mgsmootherbase), intent(in) :: this
    logical, intent(in) :: is_vcycle_end
    logical :: conv_flag
 
    real(rp) :: res_l2, res_max
    !---------------------------------
 
    call this%Get_current_residual_statistics( res_l2, res_max )
 
    conv_flag = .false.
    if ( is_vcycle_end ) then
      if ( (      res_l2 < this%threshold_residual_l2         &
                .and. res_max < this%threshold_residual_max ) ) then
        conv_flag = .true.
      end if
    else
      if ( ( res_l2 / this%history_residual_l2_initial < this%threshold_ratio_residual_l2 ) &
          .or. (      res_l2 < this%threshold_residual_l2                                   &
                .and. res_max < this%threshold_residual_max ) ) then
        conv_flag = .true.
      end if
    end if
    return
  end function mgsmootherbase_is_converged
 
!OCL SERIAL
  subroutine mgsmootherbase_output_residual_history( this )
    implicit none
    class(mgsmootherbase), intent(in) :: this
    
    integer :: i
    !---------------------------------
 
    log_info( 'MGSmootherBase_output_residual_history','(A,2ES12.4)') "Smoothing residual history (ini, L2 norm, Max norm) : ", this%history_residual_l2_initial, this%history_residual_max_initial
    log_info( 'MGSmootherBase_output_residual_history',*) "Smoothing residual history (itr, L2 norm, Max norm) : "
    do i=1, this%history_residual_count
      log_info( 'MGSmootherBase_output_residual_history', '(I8, 2ES12.4)' ) i*this%history_residual_step, &
        this%history_residual_l2(i), this%history_residual_max(i) 
    end do
    return
  end subroutine mgsmootherbase_output_residual_history
 
!-- 2D -----------------------------------------
 
!OCL SERIAL
  subroutine mgsmootherbase2d_init(this)
    implicit none
    class(mgsmootherbase2d), intent(inout) :: this
    !---------------------------------
    call mgsmootherbase_init( this )
    return
  end subroutine mgsmootherbase2d_init
 
!OCL SERIAL
  subroutine mgsmootherbase2d_final(this)
    implicit none
    class(mgsmootherbase2d), intent(inout) :: this
    !---------------------------------
    call mgsmootherbase_final( this )
    return
  end subroutine mgsmootherbase2d_final
 
!OCL SERIAL
  subroutine mgsmootherbase2d_do_smoothing( this, q, res,  &
    f, aux_var, var_comm, aux_comm, Dx, Dy, Lift, mesh2D,  &
    mg_p_level, mg_h_level,                                &
    pre_or_post_smooth )
    implicit none
    class(mgsmootherbase2d), intent(inout) :: this
    type(meshfield2d), intent(inout), target :: q
    type(meshfield2d), intent(inout) :: res
    type(meshfield2d), intent(in) :: f
    type(meshfield2d), intent(inout), target :: aux_var(:)
    class(meshfieldcommbase), intent(inout) :: var_comm
    class(meshfieldcommbase), intent(inout) :: aux_comm
    type(sparsemat), intent(in) :: dx
    type(sparsemat), intent(in) :: dy
    type(sparsemat), intent(in) :: lift
    class(meshbase2d), intent(in), target :: mesh2d
    integer, intent(in) :: mg_p_level
    integer, intent(in) :: mg_h_level
    integer, intent(in) :: pre_or_post_smooth
 
    integer :: itr
    logical :: zero_initial_guess
    logical :: cal_res_flag
    logical :: conv_flag
    !--------------------------------------
 
    this%history_residual_count = 0
    this%history_residual_rstep = this%history_residual_step
    
    do itr=1, this%num_smooth_itr_max
      zero_initial_guess = ( pre_or_post_smooth == mgsmoother_pre_id .and. mg_p_level /= 1 .and. itr == 1 )
      cal_res_flag = this%Is_residual_check_step()
 
      call this%Advance_itr_1step( q, res,                          &
        f, aux_var, itr, var_comm, aux_comm, dx, dy, lift, mesh2d,  &
        cal_res_flag, zero_initial_guess, mg_p_level, mg_h_level,   &
        pre_or_post_smooth )
 
      call this%put_residual( res, conv_flag )
      if ( cal_res_flag ) then
        if ( conv_flag ) then
          exit
        else if ( itr == this%num_smooth_itr_max ) then
          log_info( 'MGSmootherBase2D_Do_smoothing',*) "Smoothing did not converge in the maximum number of iterations.."
          log_info( 'MGSmootherBase2D_Do_smoothing',*) "Residual L2, max norm   : ", this%history_residual_l2(this%history_residual_count), this%history_residual_max(this%history_residual_count)
          exit
        end if
      end if
    end do
    return
  end subroutine mgsmootherbase2d_do_smoothing
 
!OCL SERIAL
  subroutine mgsmootherbase2d_eval_residual_statistic( this, res, &
    res_l2, res_max )
    implicit none
    class(mgsmootherbase2d), intent(in) :: this
    type(meshfield2d), intent(in), target :: res
    real(rp), intent(out) :: res_l2
    real(rp), intent(out) :: res_max
 
    integer :: ldom
    class(meshbase2d), pointer :: mesh2d
    class(localmesh2d), pointer :: lmesh2d
 
    real(rp) :: res_l2_lc, res_max_lc
    !---------------------------------
 
    mesh2d => res%mesh
    res_l2_lc = 0.0_rp; res_max_lc = 0.0_rp
 
    do ldom=1, mesh2d%LOCAL_MESH_NUM
      lmesh2d => mesh2d%lcmesh_list(ldom)
      call mgsmootherbase_get_statistic_lc( res_l2_lc, res_max_lc, &
        res%local(ldom)%val, lmesh2d, mesh2d%refElem2D )
    end do
    call this%Eval_global_statistic( res_l2, res_max, &
      res_l2_lc, res_max_lc )
    return
  end subroutine mgsmootherbase2d_eval_residual_statistic
 
!OCL SERIAL
  subroutine mgsmootherbase2d_set_initial_residual( this, res )
    implicit none
    class(mgsmootherbase2d), intent(inout) :: this
    type(meshfield2d), intent(in) :: res
    !----------------------------------------
    call this%Eval_residual_statistic( res, & ! (in)
      this%history_residual_l2_initial, & ! (out)
      this%history_residual_max_initial ) ! (out)
    return
  end subroutine mgsmootherbase2d_set_initial_residual
  
  !> Put residual history for 2D smoother.
  !! To update internal counters, this subroutine should be called after each smoothing iteration 
  !! even when the evaluation of residual is not performed. 
!OCL SERIAL
  subroutine mgsmootherbase2d_put_residual( this, res, conv_flag)
    implicit none
    class(mgsmootherbase2d), intent(inout) :: this
    type(meshfield2d), intent(in), target :: res !< Residual field
    logical, intent(out) :: conv_flag
 
    real(rp) :: res_l2, res_max
    logical :: is_residual_check_step
    !---------------------------------
 
    is_residual_check_step = this%Is_residual_check_step()
    if ( is_residual_check_step ) then
      call this%Eval_residual_statistic( res, & ! (in)
        res_l2, res_max ) ! (out)
    end if
    call mgsmootherbase_put_residual_history( this, &
      conv_flag,                              & ! (out)
      is_residual_check_step, res_l2, res_max ) ! (in)
    return
  end subroutine mgsmootherbase2d_put_residual
 
!-- 3D -----------------------------------------
 
!OCL SERIAL
  subroutine mgsmootherbase3d_init(this)
    implicit none
    class(mgsmootherbase3d), intent(inout) :: this
    !---------------------------------
    call mgsmootherbase_init( this )
    return
  end subroutine mgsmootherbase3d_init
 
!OCL SERIAL
  subroutine mgsmootherbase3d_final(this)
    implicit none
    class(mgsmootherbase3d), intent(inout) :: this
    !---------------------------------
    call mgsmootherbase_final( this )
    return
  end subroutine mgsmootherbase3d_final
 
!OCL SERIAL
  subroutine mgsmootherbase3d_do_smoothing( this, q, res,      &
    f, aux_var, var_comm, aux_comm, Dx, Dy, Dz, Lift, mesh3D,  &
    mg_p_level, mg_h_level,                                    &
    pre_or_post_smooth )
    implicit none
    class(mgsmootherbase3d), intent(inout) :: this
    type(meshfield3d), intent(inout), target :: q
    type(meshfield3d), intent(inout) :: res
    type(meshfield3d), intent(in) :: f
    type(meshfield3d), intent(inout), target :: aux_var(:)
    class(meshfieldcommbase), intent(inout) :: var_comm
    class(meshfieldcommbase), intent(inout) :: aux_comm
    type(sparsemat), intent(in) :: dx
    type(sparsemat), intent(in) :: dy
    type(sparsemat), intent(in) :: dz
    type(sparsemat), intent(in) :: lift
    class(meshbase3d), intent(in), target :: mesh3d
    integer, intent(in) :: mg_p_level
    integer, intent(in) :: mg_h_level
    integer, intent(in) :: pre_or_post_smooth
 
    integer :: itr
    logical :: zero_initial_guess
    logical :: cal_res_flag
    logical :: conv_flag
    !--------------------------------------
 
    this%history_residual_count = 0
    this%history_residual_rstep = this%history_residual_step
    
    do itr=1, this%num_smooth_itr_max
      cal_res_flag = this%Is_residual_check_step()
      zero_initial_guess = ( pre_or_post_smooth == mgsmoother_pre_id .and. mg_p_level /= 1 .and. itr == 1 )
 
      call this%Advance_itr_1step( q, res,                              &
        f, aux_var, itr, var_comm, aux_comm, dx, dy, dz, lift, mesh3d,  &
        cal_res_flag, zero_initial_guess, mg_p_level, mg_h_level,       &
        pre_or_post_smooth )
 
      call this%put_residual( res, conv_flag )
      if ( cal_res_flag ) then
        if ( conv_flag ) then
          exit
        else if ( itr == this%num_smooth_itr_max ) then
          log_info( 'MGSmootherBase3D_Do_smoothing',*) "Smoothing did not converge in the maximum number of iterations.."
          log_info( 'MGSmootherBase3D_Do_smoothing',*) "Residual L2, max norm   : ", this%history_residual_l2(this%history_residual_count), this%history_residual_max(this%history_residual_count)
          exit
        end if
      end if
    end do
    return
  end subroutine mgsmootherbase3d_do_smoothing
 
!OCL SERIAL
  subroutine mgsmootherbase3d_eval_residual_statistic( this, res, &
    res_l2, res_max )
    implicit none
    class(mgsmootherbase3d), intent(in) :: this
    type(meshfield3d), intent(in) :: res
    real(rp), intent(out) :: res_l2
    real(rp), intent(out) :: res_max
 
    integer :: ldom
    class(meshbase3d), pointer :: mesh3d
    class(localmesh3d), pointer :: lmesh3d
 
    real(rp) :: res_l2_lc, res_max_lc
    !---------------------------------
 
    mesh3d => res%mesh
    res_l2_lc = 0.0_rp; res_max_lc = 0.0_rp
 
    do ldom=1, mesh3d%LOCAL_MESH_NUM
      lmesh3d => mesh3d%lcmesh_list(ldom)
      call mgsmootherbase_get_statistic_lc( res_l2_lc, res_max_lc, &
        res%local(ldom)%val, lmesh3d, mesh3d%refElem3D )
    end do
    call this%Eval_global_statistic( res_l2, res_max, &
      res_l2_lc, res_max_lc )
    return
  end subroutine mgsmootherbase3d_eval_residual_statistic
 
!OCL SERIAL
  subroutine mgsmootherbase3d_set_initial_residual( this, res )
    implicit none
    class(mgsmootherbase3d), intent(inout) :: this
    type(meshfield3d), intent(in) :: res
    !----------------------------------------
    call this%Eval_residual_statistic( res, & ! (in)
      this%history_residual_l2_initial, & ! (out)
      this%history_residual_max_initial ) ! (out)
    return
  end subroutine mgsmootherbase3d_set_initial_residual
 
!OCL SERIAL
  subroutine mgsmootherbase3d_put_residual( this, res, conv_flag )
    implicit none
    class(mgsmootherbase3d), intent(inout) :: this
    type(meshfield3d), intent(in), target :: res
    logical, intent(out) :: conv_flag
 
    real(rp) :: res_l2, res_max
    logical :: is_residual_check_step
    !---------------------------------
 
    is_residual_check_step = this%Is_residual_check_step()
    if ( is_residual_check_step ) then
      call this%Eval_residual_statistic( res, & ! (in)
        res_l2, res_max ) ! (out)
    end if
    call mgsmootherbase_put_residual_history( this, &
      conv_flag,                              & ! (out)
      is_residual_check_step, res_l2, res_max ) ! (in)
    return
  end subroutine mgsmootherbase3d_put_residual
 
!-- private subroutines ----------------------------------------------------
 
  !> Put residual history and evaluate whether convergence is achieved.
  !! This subroutine also updates internal counters.
!OCL SERIAL
  subroutine mgsmootherbase_put_residual_history( this, conv_flag, &
    is_residual_check_step, res_l2, res_max )
    implicit none
    class(mgsmootherbase), intent(inout) :: this
    logical, intent(out) :: conv_flag
    logical, intent(in) :: is_residual_check_step
    real(rp), intent(in) :: res_l2
    real(rp), intent(in) :: res_max
    !---------------------------------
 
    if ( is_residual_check_step ) then
      !- Store the statistics of residual to history
      this%history_residual_count = this%history_residual_count + 1
      this%history_residual_l2 (this%history_residual_count) = res_l2
      this%history_residual_max(this%history_residual_count) = res_max
 
      !- Reset remaining step
      this%history_residual_rstep = this%history_residual_step
 
      !- Evaluate convergence
      conv_flag = this%Is_converged( .false. )
    else
      this%history_residual_rstep = this%history_residual_rstep - 1
      conv_flag = .false.
    end if
    return
  end subroutine mgsmootherbase_put_residual_history
 
  !> Get residual statistics on local mesh
!OCL SERIAL
  subroutine mgsmootherbase_get_statistic_lc( res_l2_lc, res_max_lc, &
    res_lc, lmesh, elem )
    use scale_localmesh_base, only: localmeshbase
    use scale_element_base, only: elementbase
    implicit none
    class(localmeshbase), intent(in), target :: lmesh
    class(elementbase), intent(in), target :: elem
    real(rp), intent(inout) :: res_l2_lc
    real(rp), intent(inout) :: res_max_lc
    real(rp), intent(in) :: res_lc(elem%np,lmesh%nea)
 
    integer :: ke
    real(rp) :: res_(elem%np)
    !---------------------------------
 
    !$omp parallel do private(ke,res_) reduction(+:res_l2_lc) reduction(max:res_max_lc)
    do ke=lmesh%NeS, lmesh%NeE
      res_(:) = res_lc(:,ke)
      res_l2_lc = res_l2_lc + sum( res_(:)**2 )
      res_max_lc = max(res_max_lc, maxval( abs(res_(:)) ))
    end do
    return
  end subroutine mgsmootherbase_get_statistic_lc
 
!OCL SERIAL
  subroutine mgsmootherbase_eval_global_statistic( this, res_l2_global, res_max_global, &
    res_l2_lc, res_max_lc)
    use mpi
    use scale_prc, only: &
      prc_local_comm_world
    implicit none
    class(mgsmootherbase), intent(in) :: this
    real(rp), intent(out) :: res_l2_global
    real(rp), intent(out) :: res_max_global
    real(rp), intent(in) :: res_l2_lc
    real(rp), intent(in) :: res_max_lc
 
    integer :: ierr
    !---------------------------------
 
    !- Gather L2 residual from all processes
    call mpi_allreduce( res_l2_lc, res_l2_global, &
      1, mpi_double_precision, mpi_sum, prc_local_comm_world, ierr )
      
    res_l2_global = sqrt( res_l2_global )
 
    !- Gather max residual from all processes
    call mpi_allreduce( res_max_lc, res_max_global, &
      1, mpi_double_precision, mpi_max, prc_local_comm_world, ierr )
    return
  end subroutine mgsmootherbase_eval_global_statistic
 
end module scale_multigrid_smoother_base