Feature/speedup trcadv

FElib/src/element/scale_element_operation_base.F90

@@ -95,19 +95,14 @@ subroutine ElementOperationBase_DxDyDzLift( this, vec_in, vec_in_lift, vec_out_d
       real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
     end subroutine ElementOperationBase_DxDyDzLift
 
-    subroutine ElementOperationBase_Div( this, vec_in_x, vec_in_y, vec_in_z, vec_in_lift, &
-      vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
+    subroutine ElementOperationBase_Div( this, vec_in, vec_in_lift, &
+      vec_out )
       import ElementOperationBase3D      
       import RP
       class(ElementOperationBase3D), intent(in) :: this
-      real(RP), intent(in) :: vec_in_x(this%elem3D%Np)
-      real(RP), intent(in) :: vec_in_y(this%elem3D%Np)
-      real(RP), intent(in) :: vec_in_z(this%elem3D%Np)
+      real(RP), intent(in) :: vec_in(this%elem3D%Np,3)
       real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
-      real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
-      real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
-      real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
-      real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
+      real(RP), intent(out) :: vec_out(this%elem3D%Np,4)
     end subroutine ElementOperationBase_Div
 
     subroutine ElementOperationBase_Div_var5( this, vec_in, vec_in_lift, &

FElib/src/element/scale_element_operation_general.F90

@@ -264,24 +264,19 @@ end subroutine element_operation_general_DxDyDzLift
 !> Calculate the 3D gradient
 !!
 !OCL SERIAL
-  subroutine element_operation_general_Div( this, vec_in_x, vec_in_y, vec_in_z, vec_in_lift, &
-    vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
+  subroutine element_operation_general_Div( this, vec_in, vec_in_lift, &
+    vec_out )
     implicit none
     class(ElementOperationGeneral), intent(in) :: this
-    real(RP), intent(in) :: vec_in_x(this%elem3D%Np)
-    real(RP), intent(in) :: vec_in_y(this%elem3D%Np)
-    real(RP), intent(in) :: vec_in_z(this%elem3D%Np)
+    real(RP), intent(in) :: vec_in(this%elem3D%Np,3)
     real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
-    real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
+    real(RP), intent(out) :: vec_out(this%elem3D%Np,4)
     !---------------------------------------------------------------
 
-    call sparsemat_matmul( this%Dx_sm, vec_in_x, vec_out_dx )
-    call sparsemat_matmul( this%Dy_sm, vec_in_y, vec_out_dy )
-    call sparsemat_matmul( this%Dz_sm, vec_in_z, vec_out_dz )
-    call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out_lift )
+    call sparsemat_matmul( this%Dx_sm, vec_in(:,1), vec_out(:,1) )
+    call sparsemat_matmul( this%Dy_sm, vec_in(:,2), vec_out(:,2) )
+    call sparsemat_matmul( this%Dz_sm, vec_in(:,3), vec_out(:,3) )
+    call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out(:,4) )
     return
   end subroutine element_operation_general_Div    
 

FElib/src/element/scale_element_operation_tensorprod3D.F90.erb

@@ -286,31 +286,30 @@ contains
     return
   end subroutine element_operation_tensorprod3D_DxDyDzLift_P<%=p%>
 
-!> Calculate the 3D gradient
+!> Calculate the 3D divergence
+!!
+!! @param vec_in Array storing flux data (Fx,Fy,Fz) at nodes that we apply the differential matrices Dx, Dy, and Dz, respectively
+!! @param vec_in_lift Array storing surface flux data (Fs) at nodes that we apply the lifting operator
+!! @param vec_out Array storing the resulting data (Dx Fx, Dy Fy, Dz Fz, Lift Fs)
 !!
 !OCL SERIAL
-  subroutine element_operation_tensorprod3D_Div_P<%=p%>(  this, vec_in_x, vec_in_y, vec_in_z, vec_in_lift, &
-    vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
+  subroutine element_operation_tensorprod3D_Div_P<%=p%>(  this, vec_in, vec_in_lift, &
+    vec_out )
     use scale_element_operation_tensorprod3D_kernel, only: &
       element_operation_kernel_matvec_divlike_dirXYZ_P<%=p%>, &
       element_operation_kernel_matvec_Lift_hexahedral_P<%=p%>
     implicit none
     class(ElementOperationTensorProd3D_P<%=p%>), intent(in) :: this
-    real(RP), intent(in) :: vec_in_x(this%elem3D%Np)
-    real(RP), intent(in) :: vec_in_y(this%elem3D%Np)
-    real(RP), intent(in) :: vec_in_z(this%elem3D%Np)
+    real(RP), intent(in) :: vec_in(this%elem3D%Np,3)
     real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
-    real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
-    real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
+    real(RP), intent(out) :: vec_out(this%elem3D%Np,4)
     !----------------------------------------------------------
 
     call element_operation_kernel_matvec_Lift_hexahedral_P<%=p%>( this%Lift_mat, vec_in_lift, &
-      vec_out_lift )      
+      vec_out(:,4) )      
 
-    call element_operation_kernel_matvec_divlike_dirXYZ_P<%=p%>( this%D1D, this%D1D_tr, vec_in_x, vec_in_y, vec_in_z, &
-      vec_out_dx, vec_out_dy, vec_out_dz )
+    call element_operation_kernel_matvec_divlike_dirXYZ_P<%=p%>( this%D1D, this%D1D_tr, vec_in(:,1), vec_in(:,2), vec_in(:,3), &
+      vec_out(:,1), vec_out(:,2), vec_out(:,3) )
     return
   end subroutine element_operation_tensorprod3D_Div_P<%=p%>
 

FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90

@@ -454,7 +454,7 @@ subroutine AtmDynDGMDriver_trcadv3d_update( this, &
               this%alphaDensM%local(n)%face_val, this%alphaDensP%local(n)%face_val,         & ! (in)
               DENS_hyd%local(n)%val, DDENS_TRC%local(n)%val, DDENS0_TRC%local(n)%val,       & ! (in)
               this%tint(n)%coef_c_ex(rkstage), dttmp_trc,                                   & ! (in) 
-              Dx, Dy, Dz, Sx, Sy, Sz, Lift, this%FaceIntMat,                                & ! (in)
+              this%FaceIntMat,                                                              & ! (in)
               lcmesh3D, lcmesh3D%refElem3D, lcmesh3D%lcmesh2D, lcmesh3D%lcmesh2D%refElem2D, & ! (in)
               this%disable_limiter                                                          ) ! (in)
           end do
@@ -478,7 +478,7 @@ subroutine AtmDynDGMDriver_trcadv3d_update( this, &
               this%alphaDensM%local(n)%face_val, this%alphaDensP%local(n)%face_val,         & ! (in)
               this%AUX_TRCVARS3D(AUXTRCVARS3D_FCTCOEF_ID)%local(n)%val,                     & ! (out)
               RHOQ_tp%local(n)%val,                                                         & ! (in)
-              Dx, Dy, Dz, Sx, Sy, Sz, Lift, this%FaceIntMat,                                & ! (in)
+              element_operation, this%FaceIntMat,                                           & ! (in)
               lcmesh3D, lcmesh3D%refElem3D, lcmesh3D%lcmesh2D, lcmesh3D%lcmesh2D%refElem2D  ) ! (in)
           else
             !$omp parallel do

FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_nonhydro3d_common.F90

@@ -552,9 +552,9 @@ subroutine atm_dyn_dgm_nonhydro3d_common_calc_phyd_hgrad_lc( DPhydDx, DPhydDy, &
 
     integer :: ke, ke2D, p  
 
-    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np,2), DFlux(elem%Np,4,2)
+    real(RP) :: Flux(elem%Np,3), Fz(elem%Np), DFlux(elem%Np,4,2)
     real(RP) :: del_flux_hyd(elem%NfpTot,2,lmesh%Ne)
-    real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np)
+    real(RP) :: GsqrtV, RGsqrtV(elem%Np)
 
     real(RP) :: E33
     real(RP) :: GradPhyd_x, GradPhyd_y
@@ -568,26 +568,27 @@ subroutine atm_dyn_dgm_nonhydro3d_common_calc_phyd_hgrad_lc( DPhydDx, DPhydDy, &
       lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D                     ) ! (in)
 
     !$omp parallel do private( ke, ke2D, p, &
-    !$omp Fx, Fy, Fz, DFlux, GsqrtV, RGsqrtV, E33, &
+    !$omp Flux, Fz, DFlux, GsqrtV, RGsqrtV, E33, &
     !$omp GradPhyd_x, GradPhyd_y )      
     do ke = lmesh%NeS, lmesh%NeE
       ke2d = lmesh%EMap3Dto2D(ke)
 
       do p=1, elem%Np
-        GsqrtV(p)  = lmesh%Gsqrt(p,ke) / ( lmesh%gam(p,ke)**2 * lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) )
-        RGsqrtV(p) = 1.0_RP / GsqrtV(p)
+        GsqrtV  = lmesh%Gsqrt(p,ke) / ( lmesh%gam(p,ke)**2 * lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) )
+
+        RGsqrtV(p) = 1.0_RP / GsqrtV
+        Flux(p,1) = GsqrtV * ( PRES_hyd(p,ke) - PRES_hyd_ref(p,ke) )
       end do
 
-      do p=1, elem%Np
-        Fx(p) = GsqrtV(p) * ( PRES_hyd(p,ke) - PRES_hyd_ref(p,ke) )
-        Fy(p) = Fx(p)
-        Fz(p,1) = lmesh%GI3(p,ke,1) * Fx(p)
-        Fz(p,2) = lmesh%GI3(p,ke,2) * Fx(p)
+      do p=1, elem%Np        
+        Flux(p,2) = Flux(p,1)
+        Flux(p,3) = lmesh%GI3(p,ke,1) * Flux(p,1)
+        Fz  (p)   = lmesh%GI3(p,ke,2) * Flux(p,1)
       end do   
 
-      call element3D_operation%Div( Fx, Fy, Fz(:,1), del_flux_hyd(:,1,ke), &
-        DFlux(:,1,1), DFlux(:,2,1), DFlux(:,3,1), DFlux(:,4,1) )
-      call element3D_operation%Dz( Fz(:,2), DFlux(:,3,2) )
+      call element3D_operation%Div( Flux, del_flux_hyd(:,1,ke), &
+        DFlux(:,:,1) )
+      call element3D_operation%Dz( Fz, DFlux(:,3,2) )
       call element3D_operation%Lift( del_flux_hyd(:,2,ke), DFlux(:,4,2) )
 
       do p=1, elem%Np

FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_trcadvect3d_heve.F90

@@ -31,6 +31,7 @@ module scale_atm_dyn_dgm_trcadvect3d_heve
   use scale_element_base, only: &
     ElementBase2D, ElementBase3D
   use scale_element_hexahedral, only: HexahedralElement
+  use scale_element_operation_base, only: ElementOperationBase3D
   use scale_localmesh_2d, only: LocalMesh2D  
   use scale_localmesh_3d, only: LocalMesh3D
   use scale_mesh_base2d, only: MeshBase2D
@@ -150,14 +151,13 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_cal_tend( &
     QTRC_, MOMX_, MOMY_, MOMZ_,                     & ! (in)
     alphDENS_M, alphDENS_P, fct_coef,               & ! (in)
     RHOQ_tp,                                        & ! (in)
-    Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat,       & ! (in) 
+    element3D_operation, FaceIntMat,                & ! (in) 
     lmesh, elem, lmesh2D, elem2D                    ) ! (in)
 
     class(LocalMesh3D), intent(in) :: lmesh
     class(ElementBase3D), intent(in) :: elem
     class(LocalMesh2D), intent(in) :: lmesh2D
     class(ElementBase2D), intent(in) :: elem2D
-    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat
 
     real(RP), intent(out) :: QTRC_dt(elem%Np,lmesh%NeA)
     real(RP), intent(in) :: QTRC_(elem%Np,lmesh%NeA)
@@ -168,11 +168,15 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_cal_tend( &
     real(RP), intent(in) :: alphDENS_P(elem%NfpTot,lmesh%Ne)
     real(RP), intent(in) :: fct_coef(elem%Np,lmesh%NeA)
     real(RP), intent(in) :: RHOQ_tp(elem%Np,lmesh%NeA)
+    class(ElementOperationBase3D), intent(in) :: element3D_operation
+    type(SparseMat), intent(in) :: FaceIntMat
 
-    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)    
+    real(RP) :: Flux(elem%Np,3), DFlux(elem%Np,4)
     real(RP) :: del_flux(elem%NfpTot,lmesh%Ne)
-    real(RP) :: momwt_(elem%Np)    
-
+    real(RP) :: momwt_
+    real(RP) :: Gsqrt_
+    real(RP) :: RGsqrt(elem%Np), RGsqrtV(elem%Np)
+
     integer :: ke, ke2d 
 
     real(RP) :: Q0, Q1, vol
@@ -191,26 +195,34 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_cal_tend( &
     call PROF_rapend('cal_trcadv_tend_bndflux', 3)
 
     call PROF_rapstart('cal_trcadv_tend_interior', 3)
-    !$omp parallel do private( &
-    !$omp momwt_, Fx, Fy, Fz, LiftDelFlx,  &    
-    !$omp ke, ke2d                         )  
+    !$omp parallel do private( ke, ke2D, &
+    !$omp momwt_, Flux, DFlux, Gsqrt_, RGsqrt, RGsqrtV )  
     do ke=lmesh%NeS, lmesh%NeE
       ke2d = lmesh%EMap3Dto2D(ke)
 
-      momwt_(:) = MOMZ_(:,ke) * lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d) / lmesh%Gsqrt(:,ke) &
-                + lmesh%GI3(:,ke,1) * MOMX_(:,ke) + lmesh%GI3(:,ke,2) * MOMY_(:,ke)
-
-      call sparsemat_matmul(Dx, lmesh%Gsqrt(:,ke) * MOMX_(:,ke) * QTRC_(:,ke), Fx)
-      call sparsemat_matmul(Dy, lmesh%Gsqrt(:,ke) * MOMY_(:,ke) * QTRC_(:,ke), Fy)
-      call sparsemat_matmul(Dz, lmesh%Gsqrt(:,ke) * momwt_(:)   * QTRC_(:,ke), Fz)
-      call sparsemat_matmul(Lift, lmesh%Fscale(:,ke) * del_flux(:,ke), LiftDelFlx)
-
-      QTRC_dt(:,ke) = ( &
-          - lmesh%Escale(:,ke,1,1) * Fx(:)     &
-          - lmesh%Escale(:,ke,2,2) * Fy(:)     &
-          - lmesh%Escale(:,ke,3,3) * Fz(:)     &
-          - LiftDelFlx(:)                   )  &
-          / lmesh%Gsqrt(:,ke)                  &
+      do p=1, elem%Np
+        RGsqrt(p) = 1.0_RP / lmesh%Gsqrt(p,ke)
+        RGsqrtV(p) = lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) / RGsqrt(p)
+      end do
+
+      do p=1, elem%Np
+        Gsqrt_ = lmesh%Gsqrt(p,ke)
+        Flux(p,1) = Gsqrt_ * MOMX_(p,ke) * QTRC_(p,ke)
+        Flux(p,2) = Gsqrt_ * MOMY_(p,ke) * QTRC_(p,ke)
+        Flux(p,3) = Gsqrt_ * ( &
+            MOMZ_(p,ke) * RGsqrtV(p)        &
+          + lmesh%GI3(p,ke,1) * MOMX_(p,ke) &
+          + lmesh%GI3(p,ke,2) * MOMY_(p,ke) ) * QTRC_(p,ke)
+      end do
+      call element3D_operation%Div( &
+        Flux, del_flux(:,ke), & ! (in)
+        DFlux                 ) ! (out)
+
+      QTRC_dt(:,ke) = - ( &
+            lmesh%Escale(:,ke,1,1) * DFlux(:,1) &
+          + lmesh%Escale(:,ke,2,2) * DFlux(:,2) &
+          + lmesh%Escale(:,ke,3,3) * DFlux(:,3) &
+          + DFlux(:,4) ) / lmesh%Gsqrt(:,ke)    &
           + RHOQ_tp(:,ke)
     end do
     call PROF_rapend('cal_trcadv_tend_interior', 3)
@@ -223,8 +235,7 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef( &
     fct_coef,                                                     & ! (out)
     QTRC_, MOMX_, MOMY_, MOMZ_, RHOQ_tp_, AlphDens_M, AlphDens_P, & ! (in)
     DENS_hyd, DDENS_, DDENS0_, rk_c_ssm1, dt,                     & ! (in)
-    Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat,                     & ! (in)
-    lmesh, elem, lmesh2D, elem2D,                                 & ! (in)
+    FaceIntMat, lmesh, elem, lmesh2D, elem2D,                     & ! (in)
     disable_limiter                                               ) ! (in)
 
     implicit none
@@ -245,7 +256,7 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef( &
     real(RP), intent(in) :: DDENS0_(elem%Np,lmesh%NeA)
     real(RP), intent(in) :: rk_c_ssm1
     real(RP), intent(in) :: dt
-    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat
+    type(SparseMat), intent(in) :: FaceIntMat
     logical, intent(in), optional :: disable_limiter
 
     real(RP) :: netOutwardFlux(lmesh%Ne)
@@ -644,15 +655,17 @@ subroutine atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
       do f=1, elem%Nfaces_h
         do p=1, elem%Nfp_h
           fp = p + (f-1)*elem%Nfp_h
-          del_flux(fp,ke) = numflux(fp) * 0.5_RP * ( R_P(fp) + R_M(fp) - ( R_P(fp) - R_M(fp) ) * sign( 1.0_RP, outward_flux_tmp(f) ) ) &
-                          - QTRC_M(fp) * MomFlxM(fp)    
+          del_flux(fp,ke) = lmesh%Fscale(fp,ke) * &
+                        (   numflux(fp) * 0.5_RP * ( R_P(fp) + R_M(fp) - ( R_P(fp) - R_M(fp) ) * sign( 1.0_RP, outward_flux_tmp(f) ) ) &
+                          - QTRC_M(fp) * MomFlxM(fp) )
         end do
       end do 
       do f=1, elem%Nfaces_v
         do p=1, elem%Nfp_v
           fp = p + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
-          del_flux(fp,ke) = numflux(fp) * 0.5_RP * ( R_P(fp) + R_M(fp) - ( R_P(fp) - R_M(fp) ) * sign( 1.0_RP, outward_flux_tmp(elem%Nfaces_h+f) ) ) &
-                          - QTRC_M(fp) * MomFlxM(fp) 
+          del_flux(fp,ke) = lmesh%Fscale(fp,ke) * & 
+                        (   numflux(fp) * 0.5_RP * ( R_P(fp) + R_M(fp) - ( R_P(fp) - R_M(fp) ) * sign( 1.0_RP, outward_flux_tmp(elem%Nfaces_h+f) ) ) &
+                          - QTRC_M(fp) * MomFlxM(fp) )
         end do
       end do
     end do