scale_multigrid_solver_3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Multigrid / Solver 3D
!!
!! @par Description
!!      Manage a multigrid solver for 3D domain
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_multigrid_solver_3d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_sparsemat, only: sparsemat
 
  use scale_element_base, only: elementbase3d
  use scale_element_quadrilateral, only: quadrilateralelement
  use scale_element_hexahedral, only: hexahedralelement
 
  use scale_localmesh_3d, only: localmesh3d
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_mesh_rectdom2d, only: meshrectdom2d
  use scale_mesh_cubedom3d, only: meshcubedom3d
 
  use scale_meshfield_base, only: &
    meshfield2d, meshfield3d
  use scale_meshfieldcomm_rectdom2d, only: meshfieldcommrectdom2d
  use scale_meshfieldcomm_cubedom3d, only: meshfieldcommcubedom3d
 
  use scale_mesh_hierarchy_base, only: &
    pmg_finest_level => mesh_hierarchy_pmg_finest_level, &
    hmg_finest_level => mesh_hierarchy_hmg_finest_level, &
    mesh_hierarchy_type_pmg,                             &
    mesh_hierarchy_type_hmg
  use scale_mesh_hierarchy_3d, only: &
    meshhierarchy3d,      &
    meshhierarchylevel3d, &
    meshhierarchylocalmgdata3d
  
  use scale_multigrid_solver_base, only: multigridsolverbase
  use scale_multigrid_smoother_base, only: &
    mgsmootherbase3d, mgsmoother_pre_id, mgsmoother_post_id
  use scale_multigrid_fieldset_base, only: mgfieldset3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type for multigrid solver in 3D domain
  type, extends(multigridsolverbase), public :: multigridsolver3d
    type(meshhierarchy3d), pointer :: mesh_hierarchy_ptr  !< Pointer to an object to manage 3D mesh hierarchy
    class(mgsmootherbase3d), pointer :: mg_smoother_ptr   !< Pointer to an object to manage 3D MG smoother
 
    type(mgfieldset3d), allocatable :: fields_h(:)        !< Array of objects to manage MG field sets for h-MG
    type(mgfieldset3d), allocatable :: fields_p(:)        !< Array of objects to manage MG field sets for p-MG
  contains
    procedure :: init => multigridsolver3d_init
    procedure :: final => multigridsolver3d_final
    procedure :: solve => multigridsolver3d_solve
    procedure :: do_vcycle => multigridsolver3d_do_vcycle
    !-
    procedure :: do_hmg_vcycle => multigridsolver3d_do_hmg_vcycle
    procedure :: operate_pmg_restriction => multigridsolver3d_operate_pmg_restriction
    procedure :: operate_pmg_correction => multigridsolver3d_operate_pmg_correction
    procedure :: operate_hmg_restriction => multigridsolver3d_operate_hmg_restriction
    procedure :: operate_hmg_correction => multigridsolver3d_operate_hmg_correction
  end type multigridsolver3d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
contains
  !> Initialize an object for multigrid solver in 3D domain
!OCL SERIAL
  subroutine multigridsolver3d_init( this, &
    mesh_hierarchy, mg_smoother,    &
    aux_var_num, aux_vec_num        )
    use scale_multigrid_solver_base, only: multigridsolverbase_init
    implicit none
    class(multigridsolver3d), intent(inout) :: this
    class(meshhierarchy3d), intent(in), target :: mesh_hierarchy
    class(mgsmootherbase3d), intent(in), target :: mg_smoother
    integer, intent(in) :: aux_var_num
    integer, intent(in) :: aux_vec_num
 
    integer :: lev_p
    integer :: lev_h
    !-------------------------------------------------------------
 
    call multigridsolverbase_init( this, mesh_hierarchy )
 
    this%mesh_hierarchy_ptr => mesh_hierarchy
    this%mg_smoother_ptr => mg_smoother
 
    !- Prepare p-MG field sets
    allocate( this%fields_p( mesh_hierarchy%NUM_pMG_LEVEL ) )
    do lev_p=1, mesh_hierarchy%NUM_pMG_LEVEL
      call this%fields_p(lev_p)%Init( mesh_hierarchy%p_mesh_list(lev_p)%ptr, &
        aux_var_num, aux_vec_num, lev_p )
    end do
 
    !- Prepare h-MG field sets
    allocate( this%fields_h( mesh_hierarchy%NUM_hMG_LEVEL ) )
    do lev_h=1, mesh_hierarchy%NUM_hMG_LEVEL
      call this%fields_h(lev_h)%Init( mesh_hierarchy%h_mesh_list(lev_h)%ptr, &
        aux_var_num, aux_vec_num, lev_h )
    end do
    return
  end subroutine multigridsolver3d_init
 
  !> Finalize an object for multigrid solver in 3D domain
!OCL SERIAL
  subroutine multigridsolver3d_final(this)
    use scale_multigrid_solver_base, only: multigridsolverbase_final
    implicit none
    class(multigridsolver3d), intent(inout) :: this
 
    integer :: lev_p
    integer :: lev_h
    !-------------------------------------------------------------
 
    call multigridsolverbase_final( this )
 
    if (allocated(this%fields_p) ) then
      do lev_p=1, this%mesh_hierarchy_ptr%NUM_pMG_LEVEL
        call this%fields_p(lev_p)%Final()
      end do
    end if
    deallocate( this%fields_p )
 
    if (allocated(this%fields_h) ) then
      do lev_h=1, this%mesh_hierarchy_ptr%NUM_hMG_LEVEL
        call this%fields_h(lev_h)%Final()
      end do
    end if
    deallocate( this%fields_h )
 
    return
  end subroutine multigridsolver3d_final
 
  !> Solve a linear system using multigrid method in 3D domain
!OCL SERIAL
  subroutine multigridsolver3d_solve(this, q, &
    f )
    implicit none
    class(multigridsolver3d), intent(inout) :: this
 
    class(meshfield3d), intent(inout), target :: q
    class(meshfield3d), intent(inout) :: f
 
    integer :: vcyc_itr
 
    class(localmesh3d), pointer :: lmesh
    integer :: ldomID
    integer :: ke
    !-------------------------------------------------------------
 
    do ldomid=1, q%mesh%LOCAL_MESH_NUM
      lmesh => q%mesh%lcmesh_list(ldomid)
      do ke=lmesh%NeS, lmesh%NeE
        this%fields_p(pmg_finest_level)%dq%local(ldomid)%val(:,ke) = q%local(ldomid)%val(:,ke)
      end do
    end do
 
    this%current_p_lev = pmg_finest_level-1
    this%current_h_lev = hmg_finest_level-1
 
    !-
    do vcyc_itr=1, this%vcyc_num_max
      log_info("MultiGridSolver3D_solve",*) "V-cycle iteration:", vcyc_itr
 
      call this%do_Vcycle( pmg_finest_level, f, vcyc_itr )
      if ( this%Is_converged(this%mg_smoother_ptr) ) then
        log_info("MultiGridSolver3D_solve",*) "V-cycle converged: vcyc_itr=", vcyc_itr
        exit
      end if
    end do
 
    !-
    do ldomid=1, q%mesh%LOCAL_MESH_NUM
      lmesh => q%mesh%lcmesh_list(ldomid)
      do ke=lmesh%NeS, lmesh%NeE
        q%local(ldomid)%val(:,ke) = this%fields_p(pmg_finest_level)%dq%local(ldomid)%val(:,ke)
      end do
    end do
    return
  end subroutine multigridsolver3d_solve
 
  !> Do a V-cycle in 3D domain
!OCL SERIAL
  recursive subroutine multigridsolver3d_do_vcycle( this, &
    mg_level, f_in, vcyc_itr )
    implicit none
    class(multigridsolver3d), intent(inout), target :: this
    integer, intent(in) :: mg_level
    type(meshfield3d), intent(inout) :: f_in
    integer, intent(in) :: vcyc_itr
 
    logical :: invoke_hMG
 
    class(mgfieldset3d), pointer :: fs_p
    class(meshhierarchy3d), pointer :: mesh_hierarchy
    !---------------------------------------------------------------------
 
    this%current_p_lev = this%current_p_lev + 1    
 
    mesh_hierarchy => this%mesh_hierarchy_ptr
    fs_p => this%fields_p(mg_level)
 
    log_info("MultiGridSolver3D_do_Vcycle",*) "Start: p_level=", this%current_p_lev
 
    !- Pre-relaxation
    call this%mg_smoother_ptr%Do_smoothing( fs_p%dq, fs_p%res,                        &
      f_in, fs_p%aux_var, fs_p%var_comm_ptr, fs_p%aux_comm_ptr,                       &
      fs_p%Dx, fs_p%Dy, fs_p%Dz, fs_p%Lift, mesh_hierarchy%p_mesh_list(mg_level)%ptr, & 
      this%current_p_lev, this%current_h_lev, mgsmoother_pre_id )
 
    if ( vcyc_itr == 1 .and. this%current_p_lev == pmg_finest_level ) then
      call this%mg_smoother_ptr%Get_initial_residual_statistics( &
        this%history_residual_l2_initial, this%history_residual_max_initial )
    end if
    call this%mg_smoother_ptr%Output_residual_history()
 
    ! if ( vcyc_itr == 1 ) call Output_tmp_data(this, fs_p, f_in, vcyc_itr, "_pre")
 
    if ( mg_level == mesh_hierarchy%NUM_pMG_LEVEL .and. mesh_hierarchy%NUM_hMG_LEVEL == 0 ) then
      ! if ( vcyc_itr == 1 ) call Output_tmp_data(this, fs_p, f_in, vcyc_itr, "_post")
      log_info("MultiGridSolver3D_do_Vcycle",*) "End: p_level=", this%current_p_lev
      this%current_p_lev = this%current_p_lev - 1 
      return
    end if
 
    !-
    invoke_hmg = ( mg_level+1 >= mesh_hierarchy%NUM_pMG_LEVEL .and. mesh_hierarchy%NUM_hMG_LEVEL > 0 )
 
    if ( invoke_hmg ) then
      !- Restriction
      call this%Operate_pMG_restriction( this%fields_h(hmg_finest_level)%f, &
        this%fields_p(mg_level)%res, mg_level )
 
      !- Advance node in the V-cycle
      call this%do_hMG_Vcycle( hmg_finest_level, this%fields_h(hmg_finest_level)%f )
      
      !- Correction
      call this%Operate_pMG_correction( this%fields_p(mg_level)%dq, &
        this%fields_h(hmg_finest_level)%dq, mg_level )
    else
      !- Restriction
      call this%Operate_pMG_restriction( this%fields_p(mg_level+1)%f, &
        this%fields_p(mg_level)%res, mg_level )
      
      !- Advance node in the V-cycle
      call this%do_Vcycle( mg_level+1, this%fields_p(mg_level+1)%f, vcyc_itr )
 
      !- Correction
      call this%Operate_pMG_correction( this%fields_p(mg_level)%dq, &
        this%fields_p(mg_level+1)%dq, mg_level )
    end if
 
    !- Post-relaxation
    call this%mg_smoother_ptr%Do_smoothing( fs_p%dq, fs_p%res,                        &
      f_in, fs_p%aux_var, fs_p%var_comm_ptr, fs_p%aux_comm_ptr,                       &
      fs_p%Dx, fs_p%Dy, fs_p%Dz, fs_p%Lift, mesh_hierarchy%p_mesh_list(mg_level)%ptr, &
      this%current_p_lev, this%current_h_lev, mgsmoother_post_id )
 
    ! if ( vcyc_itr == 1 ) call Output_tmp_data(this, fs_p, f_in, vcyc_itr, "_post")
 
    call this%mg_smoother_ptr%Output_residual_history()
    log_info("MultiGridSolver3D_do_Vcycle",*) "End: p_level=", this%current_p_lev
 
    this%current_p_lev = this%current_p_lev - 1
    return
  end subroutine multigridsolver3d_do_vcycle
 
!OCL SERIAL
  recursive subroutine multigridsolver3d_do_hmg_vcycle( this, mg_level, f_in )
    implicit none
    class(multigridsolver3d), intent(inout), target :: this
    integer, intent(in) :: mg_level
    type(meshfield3d), intent(inout) :: f_in
 
    class(meshhierarchy3d), pointer :: mesh_hierarchy
    class(mgfieldset3d), pointer :: fs_h
    !----------------------------------------------
 
    mesh_hierarchy => this%mesh_hierarchy_ptr
    fs_h => this%fields_h(mg_level)
 
    this%current_h_lev = this%current_h_lev + 1    
 
    log_info("MultiGridSolver3D_do_hMG_Vcycle",*) "Start: h_level=", this%current_h_lev
 
    if ( mg_level == mesh_hierarchy%NUM_hMG_LEVEL ) then
      ! It should be replaced by a direct solver in the future
      call this%mg_smoother_ptr%Do_smoothing( fs_h%dq, fs_h%res,                        &
        f_in, fs_h%aux_var, fs_h%var_comm_ptr, fs_h%aux_comm_ptr,                       &
        fs_h%Dx, fs_h%Dy, fs_h%Dz, fs_h%Lift, mesh_hierarchy%h_mesh_list(mg_level)%ptr, &
        this%current_p_lev, this%current_h_lev, mgsmoother_pre_id                       )
 
      call this%mg_smoother_ptr%Output_residual_history()
 
      log_info("MultiGridSolver3D_do_hMG_Vcycle",*) "End: h_level=", this%current_h_lev
      this%current_h_lev = this%current_h_lev - 1
      return
    end if
 
    !- Pre-relaxation
    call this%mg_smoother_ptr%Do_smoothing( fs_h%dq, fs_h%res,                        &
      f_in, fs_h%aux_var, fs_h%var_comm_ptr, fs_h%aux_comm_ptr,                       &
      fs_h%Dx, fs_h%Dy, fs_h%Dz, fs_h%Lift, mesh_hierarchy%h_mesh_list(mg_level)%ptr, &
      this%current_p_lev, this%current_h_lev, mgsmoother_pre_id                       )
 
    !- Restriction
    call this%Operate_hMG_restriction( this%fields_h(mg_level+1)%f, &
      this%fields_h(mg_level)%res, mg_level )
 
    !- Advance node in the V-cycle
    call this%do_hMG_Vcycle( mg_level+1, this%fields_h(mg_level+1)%f )
 
    !- Correction
    call this%Operate_hMG_correction( this%fields_h(mg_level)%dq, &
      this%fields_h(mg_level+1)%dq, mg_level )
 
    !- Post-relaxation
    call this%mg_smoother_ptr%Do_smoothing( fs_h%dq, fs_h%res,                        &
      f_in, fs_h%aux_var, fs_h%var_comm_ptr, fs_h%aux_comm_ptr,                       &
      fs_h%Dx, fs_h%Dy, fs_h%Dz, fs_h%Lift, mesh_hierarchy%h_mesh_list(mg_level)%ptr, &
      this%current_p_lev, this%current_h_lev, mgsmoother_post_id                      )
 
    log_info("MultiGridSolver3D_do_hMG_Vcycle",*) "End: h_level=", this%current_h_lev
    this%current_h_lev = this%current_h_lev - 1
    return
  end subroutine multigridsolver3d_do_hmg_vcycle
 
!--
!OCL SERIAL
  subroutine multigridsolver3d_operate_pmg_restriction( this, res_c, &
    res, p_lev )
    implicit none
    class(multigridsolver3d), intent(in), target :: this
    class(meshfield3d), intent(inout) :: res_c
    class(meshfield3d), intent(in) :: res
    integer, intent(in) :: p_lev
 
    class(meshhierarchy3d), pointer :: hierarchy
    integer :: ldomID
    class(meshbase3d), pointer :: mesh3D
    !------------------------------------------------
 
    if ( this%mesh_hierarchy_ptr%NUM_pMG_LEVEL <  p_lev+1 ) then
      call prc_abort()
    end if
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh3d => hierarchy%p_mesh_list(p_lev)%ptr
    do ldomid=1, mesh3d%LOCAL_MESH_NUM
      call multigridsolver3d_pmg_operation( res_c%local(ldomid)%val,     &
        res%local(ldomid)%val,                                           &
        hierarchy%elem3D_list(p_lev), hierarchy%elem3D_list(p_lev+1),    &
        mesh3d%lcmesh_list(ldomid), hierarchy%p_level(p_lev)%pMat1D_f2c, &
        .false.  )
    end do
    return
  end subroutine multigridsolver3d_operate_pmg_restriction
 
!OCL SERIAL
  subroutine multigridsolver3d_operate_pmg_correction( this, dq, &
    cor_c, p_lev )
    implicit none
    class(multigridsolver3d), intent(in), target :: this
    class(meshfield3d), intent(inout) :: dq
    class(meshfield3d), intent(in) :: cor_c
    integer, intent(in) :: p_lev
 
    class(meshhierarchy3d), pointer :: hierarchy
 
    integer :: ldomID
    class(meshbase3d), pointer :: mesh3D
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh3d => hierarchy%p_mesh_list(p_lev)%ptr
 
    do ldomid=1, mesh3d%LOCAL_MESH_NUM
      call multigridsolver3d_pmg_operation( dq%local(ldomid)%val,        &
        cor_c%local(ldomid)%val,                                         &
        hierarchy%elem3D_list(p_lev+1), hierarchy%elem3D_list(p_lev),    &
        mesh3d%lcmesh_list(ldomid), hierarchy%p_level(p_lev)%pMat1D_c2f, &
        .true. )
    end do
    return
  end subroutine multigridsolver3d_operate_pmg_correction
 
!OCL SERIAL
  subroutine multigridsolver3d_operate_hmg_restriction( this, res_c, &
    res, h_lev )
    use scale_meshfield_base, only: meshfield3d
    implicit none
    class(multigridsolver3d), intent(in), target :: this
    class(meshfield3d), intent(inout) :: res_c
    class(meshfield3d), intent(in) :: res
    integer, intent(in) :: h_lev
 
    class(meshhierarchy3d), pointer :: hierarchy
 
    integer :: ldomID
    class(meshbase3d), pointer :: mesh3D
    class(meshbase3d), pointer :: mesh3D_c
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh3d => hierarchy%h_mesh_list(h_lev)%ptr
    mesh3d_c => hierarchy%h_mesh_list(h_lev+1)%ptr
 
    do ldomid=1, mesh3d%LOCAL_MESH_NUM
      call multigridsolver3d_hmg_restriction_core( res_c%local(ldomid)%val,          &
        res%local(ldomid)%val,                                                       &
        hierarchy%h_level(h_lev)%mg_local(ldomid), mesh3d%lcmesh_list(ldomid),       &
        hierarchy%elem3D_list(hierarchy%NUM_pMG_LEVEL), mesh3d_c%lcmesh_list(ldomid) )
    end do
    return
  end subroutine multigridsolver3d_operate_hmg_restriction
 
!OCL SERIAL
  subroutine multigridsolver3d_operate_hmg_correction( this, dq, &
    cor_c, h_lev )
    use scale_meshfield_base, only: meshfield3d
    implicit none
    class(multigridsolver3d), intent(in), target :: this
    class(meshfield3d), intent(inout) :: dq
    class(meshfield3d), intent(in) :: cor_c
    integer, intent(in) :: h_lev
 
    integer :: ldomID
    class(meshbase3d), pointer :: mesh3D
    class(meshhierarchy3d), pointer :: hierarchy
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh3d => hierarchy%h_mesh_list(h_lev)%ptr
 
    do ldomid=1, mesh3d%LOCAL_MESH_NUM
      call multigridsolver3d_hmg_correction_core( dq%local(ldomid)%val,     &
        cor_c%local(ldomid)%val, hierarchy%h_level(h_lev)%mg_local(ldomid), &
        mesh3d%lcmesh_list(ldomid), mesh3d%refElem3D                        )
    end do
    return
  end subroutine multigridsolver3d_operate_hmg_correction
 
!-- private --------------------------------------------------------------
 
!OCL SERIAL
  subroutine output_tmp_data( this, fs, fin, vcyc_itr, postfix )
    use scale_prc, only: prc_myrank
    use scale_file_base_meshfield, only: file_base_meshfield
    use scale_mesh_base3d, only: meshbase3d_dimtypeid_xyz
    implicit none
    class(multigridsolver3d), intent(in) :: this
    class(mgfieldset3d), intent(in), target :: fs
    class(meshfield3d), intent(in) :: fin
    integer, intent(in) :: vcyc_itr
    character(len=*), intent(in) :: postfix
 
    type(file_base_meshfield) :: file
    class(meshbase3d), pointer :: mesh
    character(len=H_MID) :: fname
    logical :: fileexisted
 
    integer, parameter :: DQ_VID = 1
    integer, parameter :: RES_VID = 2
    integer, parameter :: FIN_VID = 3
    !-----------------------------------------------------
 
    select type (mesh => fs%dq%mesh)
    type is (meshcubedom3d)
      call file%Init(2, mesh3d=mesh)
    end select
 
    write(fname,'(a,i2.2,a,i2.2,a)') "tmp_plev", this%current_p_lev, "_vcyc", vcyc_itr, trim(postfix)
    call file%Create( fname, "MG", "REAL8", fileexisted, myrank=prc_myrank )
    call file%Def_Var( "dq", "1", "dq", dq_vid, meshbase3d_dimtypeid_xyz, "REAL8")
    call file%Def_Var( "res", "1", "res", res_vid, meshbase3d_dimtypeid_xyz, "REAL8")
    call file%Def_Var( "fin", "1", "fin", fin_vid, meshbase3d_dimtypeid_xyz, "REAL8")
    call file%End_def()
 
    call file%Write_var3D( dq_vid, fs%dq, 0.0_rp, 1.0_rp)
    call file%Write_var3D( res_vid, fs%res, 0.0_rp, 1.0_rp)
    call file%Write_var3D( fin_vid, fin, 0.0_rp, 1.0_rp)
    call file%Close()
    call file%Final()
    return
  end subroutine output_tmp_data
 
!OCL SERIAL
  subroutine multigridsolver3d_hmg_restriction_core( res_c_lc,      &
      res_lc, mg_local, lmesh, elem, lmesh_c )
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh3d), intent(in) :: lmesh_c
    real(RP), intent(out) :: res_c_lc(elem%Np,lmesh_c%NeA)
    real(RP), intent(in) :: res_lc(elem%Np,lmesh%NeA)
    class(meshhierarchylocalmgdata3d), intent(in) :: mg_local
 
    integer :: k, ke, ke_c
    real(RP) :: Ic2fT_lc(8,8)
    real(RP) :: tmp_c(elem%Np)
    real(RP) :: tmp2(elem%Np)
    !-----------------------------------------
 
    !$omp parallel private(ke, ke_c, tmp2, tmp_c, Ic2fT_lc)
    !$omp do
    do ke_c=lmesh_c%NeS, lmesh_c%NeE
      
      tmp_c(:) = 0.0_rp
      do k=1, 4
        ke = mg_local%If2c_emap(k,ke_c)
 
        tmp2(:) = matmul( elem%M, res_lc(:,ke) )
        ic2ft_lc(:,:) = transpose( mg_local%Ic2f(:,:,ke) )
        tmp_c(:) = tmp_c(:) + matmul( ic2ft_lc, tmp2(:) )
      end do
 
      res_c_lc(:,ke_c) = matmul(lmesh_c%refElem3D%invM, tmp_c(:))  * 0.25_rp        
    end do
    !$omp end parallel
    return
  end subroutine multigridsolver3d_hmg_restriction_core
 
!OCL SERIAL
  subroutine multigridsolver3d_hmg_correction_core( dq_lc,  &
    dq_c_lc, mg_local, lmesh, elem )
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    real(RP), intent(out) :: dq_lc(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: dq_c_lc(elem%Np,lmesh%NeA)
    class(meshhierarchylocalmgdata3d), intent(in) :: mg_local
 
    integer :: k, ke, ke_c
    !-----------------------------------------
 
    !$omp parallel do private(ke, ke_c)
    do ke=lmesh%NeS, lmesh%NeE
      ke_c = mg_local%Ic2f_emap(ke)
      dq_lc(:,ke) = dq_lc(:,ke) + &
            matmul( mg_local%Ic2f(:,:,ke), dq_c_lc(:,ke_c) )
    end do
    return
  end subroutine multigridsolver3d_hmg_correction_core
 
!OCL SERIAL
  subroutine multigridsolver3d_pmg_operation( q_o, &
    q_i, elem3D_i, elem3D_o, lcmesh, pMat1D, is_added )
    implicit none
    class(elementbase3d), intent(in) :: elem3D_i
    class(elementbase3d), intent(in) :: elem3D_o
    class(localmesh3d), intent(in) :: lcmesh
    real(RP), intent(inout) :: q_o(elem3D_o%Nnode_h1D,elem3D_o%Nnode_h1D,elem3D_o%Nnode_v,lcmesh%NeA)
    real(RP), intent(in) :: q_i(elem3D_i%Nnode_h1D,elem3D_i%Nnode_h1D,elem3D_i%Nnode_v,lcmesh%NeA)
    real(RP), intent(in) :: pMat1D(elem3D_o%Nnode_h1D,elem3D_i%Nnode_h1D)
    logical, intent(in) :: is_added
 
    integer :: ke
 
    integer :: px, py, pz
    integer :: pxx, pyy, pzz
    real(RP) :: tmp1
    real(RP) :: tmp2(elem3D_o%Nnode_h1D,elem3D_i%Nnode_h1D)
    real(RP) :: tmp3(elem3D_o%Nnode_h1D,elem3D_o%Nnode_h1D,elem3D_i%Nnode_v)
    real(RP) :: tmp4(elem3D_o%Nnode_h1D)
    real(RP) :: tmp_h(elem3D_o%Nnode_h1D,elem3D_o%Nnode_h1D)
 
    real(RP) :: mat_tr(elem3D_i%Nnode_h1D,elem3D_o%Nnode_h1D)
    !-------------------------------------------
 
    mat_tr(:,:) = transpose(pmat1d)
 
    !$omp parallel do private(tmp1, tmp2, tmp3, tmp4, tmp_h)
    do ke=lcmesh%NeS, lcmesh%NeE
 
      do pz=1, elem3d_i%Nnode_v
        do py=1, elem3d_i%Nnode_h1D
        do pxx=1, elem3d_o%Nnode_h1D
          tmp1 = 0.0_rp
          do px=1, elem3d_i%Nnode_h1D
            tmp1 = tmp1 + mat_tr(px,pxx) * q_i(px,py,pz,ke)
          end do
          tmp2(pxx,py) = tmp1
        end do
        end do
 
        do pyy=1, elem3d_o%Nnode_h1D
          tmp4(:) = 0.0_rp
          do py=1, elem3d_i%Nnode_h1D
            do px=1, elem3d_o%Nnode_h1D
              tmp4(px) = tmp4(px) + mat_tr(py,pyy) * tmp2(px,py)
            end do
          end do
          tmp3(:,pyy,pz) = tmp4(:)
        end do
      end do
 
      if ( is_added ) then
        do pzz=1, elem3d_o%Nnode_v
          tmp_h(:,:) = 0.0_rp
          do pz=1, elem3d_i%Nnode_v
            do py=1, elem3d_o%Nnode_h1D
            do px=1, elem3d_o%Nnode_h1D
              tmp_h(px,py) = tmp_h(px,py) + mat_tr(pz,pzz) * tmp3(px,py,pz)
            end do
            end do
          end do
          q_o(:,:,pzz,ke) = q_o(:,:,pzz,ke) + tmp_h(:,:)
        end do
      else
        do pzz=1, elem3d_o%Nnode_v
          tmp_h(:,:) = 0.0_rp
          do pz=1, elem3d_i%Nnode_v
            do py=1, elem3d_o%Nnode_h1D
            do px=1, elem3d_o%Nnode_h1D
              tmp_h(px,py) = tmp_h(px,py) + mat_tr(pz,pzz) * tmp3(px,py,pz)
            end do
            end do
          end do
          q_o(:,:,pzz,ke) = tmp_h(:,:)
        end do
      end if
    end do
    return
  end subroutine multigridsolver3d_pmg_operation
end module scale_multigrid_solver_3d