scale_multigrid_solver_2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Multigrid / Solver 2D
!!
!! @par Description
!!      Manage a multigrid solver for 2D domain
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_multigrid_solver_2d
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_sparsemat, only: sparsemat
 
  use scale_element_base, only: elementbase2d
  use scale_element_quadrilateral, only: quadrilateralelement
 
  use scale_localmesh_2d, only: localmesh2d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_rectdom2d, only: meshrectdom2d
 
  use scale_meshfield_base, only: meshfield2d
  use scale_meshfieldcomm_rectdom2d, only: meshfieldcommrectdom2d
 
  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_2d, only: &
    meshhierarchy2d,      &
    meshhierarchylevel2d, &
    meshhierarchylocalmgdata2d
  
  use scale_multigrid_solver_base, only: multigridsolverbase
  use scale_multigrid_smoother_base, only: &
    mgsmootherbase2d, &
    mgsmoother_pre_id, mgsmoother_post_id
  use scale_multigrid_fieldset_base, only: mgfieldset2d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type for multigrid solver in 2D domain
  type, extends(multigridsolverbase), public :: multigridsolver2d
    type(meshhierarchy2d), pointer :: mesh_hierarchy_ptr  !< Pointer to an object to manage 2D mesh hierarchy
    class(mgsmootherbase2d), pointer :: mg_smoother_ptr   !< Pointer to an object to manage 2D MG smoother
 
    type(mgfieldset2d), allocatable :: fields_h(:)         !< Array of objects to manage MG field sets for h-MG
    type(mgfieldset2d), allocatable :: fields_p(:)         !< Array of objects to manage MG field sets for p-MG
  contains
    procedure :: init => multigridsolver2d_init
    procedure :: final => multigridsolver2d_final
    procedure :: solve => multigridsolver2d_solve
    procedure :: do_vcycle => multigridsolver2d_do_vcycle
    !-
    procedure :: do_hmg_vcycle => multigridsolver2d_do_hmg_vcycle
    procedure :: operate_pmg_restriction => multigridsolver2d_operate_pmg_restriction
    procedure :: operate_pmg_correction => multigridsolver2d_operate_pmg_correction
    procedure :: operate_hmg_restriction => multigridsolver2d_operate_hmg_restriction
    procedure :: operate_hmg_correction => multigridsolver2d_operate_hmg_correction
  end type multigridsolver2d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
contains
  !> Initialize an object for multigrid solver in 2D domain
!OCL SERIAL
  subroutine multigridsolver2d_init( this, &
    mesh_hierarchy, mg_smoother,    &
    aux_var_num, aux_vec_num        )
 
    use scale_multigrid_solver_base, only: multigridsolverbase_init
    implicit none
    class(multigridsolver2d), intent(inout) :: this
    class(meshhierarchy2d), intent(in), target :: mesh_hierarchy
    class(mgsmootherbase2d), 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 multigridsolver2d_init
 
  !> Finalize an object for multigrid solver in 2D domain
!OCL SERIAL
  subroutine multigridsolver2d_final(this)
    use scale_multigrid_solver_base, only: multigridsolverbase_final
    implicit none
    class(multigridsolver2d), intent(inout) :: this
 
    integer :: lev_p
    integer :: lev_h
    !-------------------------------------------------------------
 
    call multigridsolverbase_final( this )
 
    if ( allocated(this%fields_p) ) then
      do lev_p=1, size(this%fields_p)
        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, size(this%fields_h)
        call this%fields_h(lev_h)%Final()
      end do
    end if
    deallocate( this%fields_h )
 
    return
  end subroutine multigridsolver2d_final
 
  !> Solve a linear system using multigrid method in 2D domain
!OCL SERIAL
  subroutine multigridsolver2d_solve(this, q, &
    f )
    implicit none
    class(multigridsolver2d), intent(inout) :: this
 
    class(meshfield2d), intent(inout), target :: q
    class(meshfield2d), intent(inout) :: f
 
    integer :: vcyc_itr
 
    class(localmesh2d), 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("MultiGridSolver2D_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("MultiGridSolver2D_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 multigridsolver2d_solve
 
  !> Do a V-cycle in 2D domain
!OCL SERIAL
  recursive subroutine multigridsolver2d_do_vcycle(this, mg_level, f_in, vcyc_itr)
    implicit none
    class(multigridsolver2d), intent(inout), target :: this
    integer, intent(in) :: mg_level
    type(meshfield2d), intent(inout) :: f_in
    integer, intent(in) :: vcyc_itr
    logical :: invoke_hMG
 
    class(mgfieldset2d), pointer :: fs_p
    class(meshhierarchy2d), 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("MultiGridSolver2D_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%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 ( mg_level == mesh_hierarchy%NUM_pMG_LEVEL .and. mesh_hierarchy%NUM_hMG_LEVEL == 0 ) then
      log_info("MultiGridSolver2D_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%Lift, mesh_hierarchy%p_mesh_list(mg_level)%ptr, &
      this%current_p_lev, this%current_h_lev,  mgsmoother_post_id )
 
    call this%mg_smoother_ptr%Output_residual_history()
    log_info("MultiGridSolver2D_do_Vcycle",*) "End: p_level=", this%current_p_lev
 
    this%current_p_lev = this%current_p_lev - 1 
    return
  end subroutine multigridsolver2d_do_vcycle
 
!OCL SERIAL
  recursive subroutine multigridsolver2d_do_hmg_vcycle( this, mg_level, f_in )
    implicit none
    class(multigridsolver2d), intent(inout), target :: this
    integer, intent(in) :: mg_level
    type(meshfield2d), intent(inout) :: f_in
 
    class(meshhierarchy2d), pointer :: mesh_hierarchy
    class(mgfieldset2d), 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("MultiGridSolver2D_do_hMG_Vcycle",*) "Start: h_level=", this%current_h_lev
 
    if ( mg_level == mesh_hierarchy%NUM_hMG_LEVEL ) then
      ! It should be replaced with a direct solver
      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%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("MultiGridSolver2D_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%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%Lift, mesh_hierarchy%h_mesh_list(mg_level)%ptr, &
        this%current_p_lev, this%current_h_lev, mgsmoother_post_id             )
 
    log_info("MultiGridSolver2D_do_hMG_Vcycle",*) "End: h_level=", this%current_h_lev
    this%current_h_lev = this%current_h_lev - 1
    return
  end subroutine multigridsolver2d_do_hmg_vcycle
 
!--
!OCL SERIAL
  subroutine multigridsolver2d_operate_pmg_restriction( this, res_c, &
    res, p_lev )
    implicit none
    class(multigridsolver2d), intent(in), target :: this
    class(meshfield2d), intent(inout) :: res_c
    class(meshfield2d), intent(in) :: res
    integer, intent(in) :: p_lev
 
    class(meshhierarchy2d), pointer :: hierarchy
    integer :: ldomID
    class(meshbase2d), pointer :: mesh2D
    !------------------------------------------------
 
    if ( this%mesh_hierarchy_ptr%NUM_pMG_LEVEL <  p_lev+1 ) then
      call prc_abort()
    end if
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh2d => hierarchy%p_mesh_list(p_lev)%ptr
    do ldomid=1, mesh2d%LOCAL_MESH_NUM
      call multigridsolver2d_pmg_operation( res_c%local(ldomid)%val,     &
        res%local(ldomid)%val,                                           &
        hierarchy%elem2D_list(p_lev), hierarchy%elem2D_list(p_lev+1),    &
        mesh2d%lcmesh_list(ldomid), hierarchy%p_level(p_lev)%pMat1D_f2c, &
        .false.  )
    end do
    return
  end subroutine multigridsolver2d_operate_pmg_restriction
 
!OCL SERIAL
  subroutine multigridsolver2d_operate_pmg_correction( this, dq, &
    cor_c, p_lev )
    implicit none
    class(multigridsolver2d), intent(in), target :: this
    class(meshfield2d), intent(inout) :: dq
    class(meshfield2d), intent(in) :: cor_c
    integer, intent(in) :: p_lev
 
    class(meshhierarchy2d), pointer :: hierarchy
 
    integer :: ldomID
    class(meshbase2d), pointer :: mesh2D
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh2d => hierarchy%p_mesh_list(p_lev)%ptr
 
    do ldomid=1, mesh2d%LOCAL_MESH_NUM
      call multigridsolver2d_pmg_operation( dq%local(ldomid)%val,        &
        cor_c%local(ldomid)%val,                                         &
        hierarchy%elem2D_list(p_lev+1), hierarchy%elem2D_list(p_lev),    &
        mesh2d%lcmesh_list(ldomid), hierarchy%p_level(p_lev)%pMat1D_c2f, &
        .true. )
    end do
    return
  end subroutine multigridsolver2d_operate_pmg_correction
 
!OCL SERIAL
  subroutine multigridsolver2d_operate_hmg_restriction( this, res_c, &
    res, h_lev )
    use scale_meshfield_base, only: meshfield2d
    implicit none
    class(multigridsolver2d), intent(in), target :: this
    class(meshfield2d), intent(inout) :: res_c
    class(meshfield2d), intent(in) :: res
    integer, intent(in) :: h_lev
 
    class(meshhierarchy2d), pointer :: hierarchy
 
    integer :: ldomID
    class(meshbase2d), pointer :: mesh2D
    class(meshbase2d), pointer :: mesh2D_c
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh2d => hierarchy%h_mesh_list(h_lev)%ptr
    mesh2d_c => hierarchy%h_mesh_list(h_lev+1)%ptr
 
    do ldomid=1, mesh2d%LOCAL_MESH_NUM
      call multigridsolver2d_hmg_restriction_core( res_c%local(ldomid)%val,          &
        res%local(ldomid)%val,                                                       &
        hierarchy%h_level(h_lev)%mg_local(ldomid), mesh2d%lcmesh_list(ldomid),       &
        hierarchy%elem2D_list(hierarchy%NUM_pMG_LEVEL), mesh2d_c%lcmesh_list(ldomid) )
    end do
    return
  end subroutine multigridsolver2d_operate_hmg_restriction
 
!OCL SERIAL
  subroutine multigridsolver2d_operate_hmg_correction( this, dq, &
    cor_c, h_lev )
    use scale_meshfield_base, only: meshfield2d
    implicit none
    class(multigridsolver2d), intent(in), target :: this
    class(meshfield2d), intent(inout) :: dq
    class(meshfield2d), intent(in) :: cor_c
    integer, intent(in) :: h_lev
 
    integer :: ldomID
    class(meshbase2d), pointer :: mesh2D
    class(meshhierarchy2d), pointer :: hierarchy
    !------------------------------------------------
 
    hierarchy => this%mesh_hierarchy_ptr
    mesh2d => hierarchy%h_mesh_list(h_lev)%ptr
 
    do ldomid=1, mesh2d%LOCAL_MESH_NUM
      call multigridsolver2d_hmg_correction_core( dq%local(ldomid)%val,     &
        cor_c%local(ldomid)%val, hierarchy%h_level(h_lev)%mg_local(ldomid), &
        mesh2d%lcmesh_list(ldomid), mesh2d%refElem2D                        )
    end do
    return
  end subroutine multigridsolver2d_operate_hmg_correction
 
!-- private --------------------------------------------------------------
 
 
!OCL SERIAL
  subroutine multigridsolver2d_hmg_restriction_core( res_c_lc,      &
      res_lc, mg_local, lmesh, elem, lmesh_c )
    implicit none
    class(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), intent(in) :: elem
    class(localmesh2d), 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(meshhierarchylocalmgdata2d), intent(in) :: mg_local
 
    integer :: k, ke, ke_c
    real(RP) :: Ic2fT_lc(4,4)
    real(RP) :: tmp_c(elem%Np)
    real(RP) :: tmp2(elem%Np)
 
    real(RP) :: int_tmp(lmesh_c%Ne)
    real(RP) :: res_cor
    !-----------------------------------------
 
    !$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) )
        ! tmp2(:) = 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%refElem2D%invM, tmp_c(:)) * 0.25_rp        
      ! res_c_lc(:,ke_c) = tmp_c(:) * 0.25_RP
    end do
    !$omp end parallel
    return
  end subroutine multigridsolver2d_hmg_restriction_core
 
!OCL SERIAL
  subroutine multigridsolver2d_hmg_correction_core( dq_lc,  &
    dq_c_lc, mg_local, lmesh, elem )
    implicit none
    class(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), 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(meshhierarchylocalmgdata2d), 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 multigridsolver2d_hmg_correction_core
 
!OCL SERIAL
  subroutine multigridsolver2d_pmg_operation( q_o, &
    q_i, elem2D_i, elem2D_o, lcmesh, pMat1D, is_added )
    implicit none
    class(elementbase2d), intent(in) :: elem2D_i
    class(elementbase2d), intent(in) :: elem2D_o
    class(localmesh2d), intent(in) :: lcmesh
    real(RP), intent(inout) :: q_o(elem2D_o%Nfp,elem2D_o%Nfp,lcmesh%NeA)
    real(RP), intent(in) :: q_i(elem2D_i%Nfp,elem2D_i%Nfp,lcmesh%NeA)
    real(RP), intent(in) :: pMat1D(elem2D_o%Nfp,elem2D_i%Nfp)
    logical, intent(in) :: is_added
 
    integer :: ke
 
    integer :: px, py
    integer :: pxx, pyy
    real(RP) :: tmp1
    real(RP) :: tmp2(elem2D_o%Nfp,elem2D_i%Nfp)
    real(RP) :: tmp3(elem2D_o%Nfp)
 
    real(RP) :: mat_tr(elem2D_i%Nfp,elem2D_o%Nfp)
    !-------------------------------------------
 
    mat_tr(:,:) = transpose(pmat1d)
 
    !$omp parallel do private(tmp1, tmp2, tmp3)
    do ke=lcmesh%NeS, lcmesh%NeE
      do py=1, elem2d_i%Nfp
      do pxx=1, elem2d_o%Nfp
        tmp1 = 0.0_rp
        do px=1, elem2d_i%Nfp
          tmp1 = tmp1 + mat_tr(px,pxx) * q_i(px,py,ke)
        end do
        tmp2(pxx,py) = tmp1
      end do
      end do
 
      if ( is_added ) then
        do pyy=1, elem2d_o%Nfp
          tmp3(:) = 0.0_rp
          do py=1, elem2d_i%Nfp
            do px=1, elem2d_o%Nfp
              tmp3(px) = tmp3(px) + mat_tr(py,pyy) * tmp2(px,py)
            end do
          end do
          q_o(:,pyy,ke) = q_o(:,pyy,ke)  + tmp3(:)
        end do
      else
        do pyy=1, elem2d_o%Nfp
          tmp3(:) = 0.0_rp
          do py=1, elem2d_i%Nfp
            do px=1, elem2d_o%Nfp
              tmp3(px) = tmp3(px) + mat_tr(py,pyy) * tmp2(px,py)
            end do
          end do
          q_o(:,pyy,ke) = tmp3(:)
        end do
      end if
    end do
    return
  end subroutine multigridsolver2d_pmg_operation
end module scale_multigrid_solver_2d