extract code from print_cfl into dss_hvtensor routine

HOMME fortran init code refactor: important code hidden in print_cfl()
is extracted into a new routine, dss_hvtensor

dss_hvtensor improves HV results by making the HV coefficients
smooth accross elements.  It also does a bilinear projection
of the coefficients within each element to minimize oscillations
in the tensor coefficieints.

components/homme/src/preqx/viscosity_mod.F90

@@ -5,7 +5,7 @@
 module viscosity_mod
   use viscosity_base, only: compute_zeta_C0, compute_div_C0, compute_zeta_C0_contra, compute_div_C0_contra, make_c0, &
        make_c0_vector, neighbor_minmax
-  use viscosity_base, only: biharmonic_wk_scalar, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis
+  use viscosity_base, only: biharmonic_wk_scalar, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis, dss_hvtensor
   use viscosity_preqx_base, only: biharmonic_wk_dp3d
   implicit none
 end module viscosity_mod

components/homme/src/preqx_acc/viscosity_mod.F90

@@ -5,7 +5,7 @@
 
 module viscosity_mod
   use viscosity_base, only: compute_zeta_C0, compute_div_C0, compute_zeta_C0_contra, compute_div_C0_contra, make_c0, make_c0_vector
-  use viscosity_base, only: biharmonic_wk_scalar,neighbor_minmax, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis
+  use viscosity_base, only: biharmonic_wk_scalar,neighbor_minmax, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis, dss_hvtensor
   use viscosity_preqx_base, only: biharmonic_wk_dp3d
   use thread_mod, only : omp_get_num_threads
   use kinds, only : real_kind, iulog
@@ -24,7 +24,7 @@ module viscosity_mod
   public :: biharmonic_wk_scalar, neighbor_minmax, neighbor_minmax_start,neighbor_minmax_finish, biharmonic_wk_dp3d
   public :: biharmonic_wk_scalar_openacc
   public :: neighbor_minmax_openacc
-  public :: smooth_phis
+  public :: smooth_phis, dss_hvtensor
 
 
 contains

components/homme/src/share/global_norms_mod.F90

@@ -5,7 +5,6 @@
 module global_norms_mod
 
   use kinds, only : iulog
-  use edgetype_mod, only : EdgeBuffer_t
 
   implicit none
   private
@@ -25,7 +24,6 @@ module global_norms_mod
   public :: wrap_repro_sum
 
   private :: global_maximum
-  type (EdgeBuffer_t), private :: edgebuf
 
 contains
 
@@ -111,9 +109,6 @@ subroutine test_global_integral(elem,hybrid,nets,nete,mindxout)
     use reduction_mod, only : ParallelMin,ParallelMax
     use physical_constants, only : scale_factor,dd_pi
     use parallel_mod, only : abortmp, global_shared_buf, global_shared_sum
-    use edgetype_mod, only : EdgeBuffer_t
-    use edge_mod, only :  initedgebuffer, FreeEdgeBuffer, edgeVpack, edgeVunpack
-    use bndry_mod, only : bndry_exchangeV
     use control_mod, only : geometry
 
     type(element_t)      , intent(inout) :: elem(:)
@@ -131,7 +126,7 @@ subroutine test_global_integral(elem,hybrid,nets,nete,mindxout)
     real (kind=real_kind) :: min_min_dx, max_min_dx, avg_min_dx
     real (kind=real_kind) :: min_normDinv, max_normDinv
     real (kind=real_kind) :: min_len
-    integer :: ie,corner, i, j,nlon
+    integer :: ie, i, j
 
 
     h(:,:,nets:nete)=1.0D0
@@ -261,9 +256,6 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)
                             hypervis_scaling, dcmip16_mu,dcmip16_mu_s,dcmip16_mu_q
     use control_mod, only : tstep_type
     use parallel_mod, only : abortmp, global_shared_buf, global_shared_sum
-    use edgetype_mod, only : EdgeBuffer_t 
-    use edge_mod, only : initedgebuffer, FreeEdgeBuffer, edgeVpack, edgeVunpack
-    use bndry_mod, only : bndry_exchangeV
     use time_mod, only : tstep
 
     type(element_t)      , intent(inout) :: elem(:)
@@ -276,10 +268,8 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)
     real (kind=real_kind) :: max_normDinv  ! used for CFL
     real (kind=real_kind) :: min_hypervis, max_hypervis, avg_hypervis, stable_hv
     real (kind=real_kind) :: normDinv_hypervis
-    real (kind=real_kind) :: x, y, noreast, nw, se, sw
-    real (kind=real_kind), dimension(np,np,nets:nete) :: zeta
     real (kind=real_kind) :: lambda_max, lambda_vis, min_gw, lambda, nu_div_actual, nu_top_actual
-    integer :: ie,corner, i, j, rowind, colind
+    integer :: ie, i, j
     type (quadrature_t)    :: gp
 
 
@@ -326,6 +316,8 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)
 
     gp=gausslobatto(np)
     min_gw = minval(gp%weights)
+    deallocate(gp%points)
+    deallocate(gp%weights)
 
     max_normDinv=0
     min_max_dx=1d99
@@ -350,58 +342,6 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)
     endif
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  TENSOR, RESOLUTION-AWARE HYPERVISCOSITY
-!  The tensorVisc() array is computed in cube_mod.F90
-!  this block of code will DSS it so the tensor if C0
-!  and also make it bilinear in each element.
-!  Oksana Guba
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   if (hypervis_scaling /= 0) then
-    call initEdgeBuffer(hybrid%par,edgebuf,elem,1)
-    do rowind=1,2
-      do colind=1,2
-	do ie=nets,nete
-	  zeta(:,:,ie) = elem(ie)%tensorVisc(:,:,rowind,colind)*elem(ie)%spheremp(:,:)
-	  call edgeVpack(edgebuf,zeta(1,1,ie),1,0,ie)
-	end do
-
-	call bndry_exchangeV(hybrid,edgebuf)
-	do ie=nets,nete
-	  call edgeVunpack(edgebuf,zeta(1,1,ie),1,0,ie)
-          elem(ie)%tensorVisc(:,:,rowind,colind) = zeta(:,:,ie)*elem(ie)%rspheremp(:,:)
-	end do
-      enddo !rowind
-    enddo !colind
-    call FreeEdgeBuffer(edgebuf)
-
-!IF BILINEAR MAP OF V NEEDED
-    do rowind=1,2
-      do colind=1,2
-    ! replace hypervis w/ bilinear based on continuous corner values
-	do ie=nets,nete
-	  noreast = elem(ie)%tensorVisc(np,np,rowind,colind)
-	  nw = elem(ie)%tensorVisc(1,np,rowind,colind)
-	  se = elem(ie)%tensorVisc(np,1,rowind,colind)
-	  sw = elem(ie)%tensorVisc(1,1,rowind,colind)
-	  do i=1,np
-	    x = gp%points(i)
-	    do j=1,np
-		y = gp%points(j)
-		elem(ie)%tensorVisc(i,j,rowind,colind) = 0.25d0*( &
-					(1.0d0-x)*(1.0d0-y)*sw + &
-					(1.0d0-x)*(y+1.0d0)*nw + &
-					(x+1.0d0)*(1.0d0-y)*se + &
-					(x+1.0d0)*(y+1.0d0)*noreast)
-	    end do
-	  end do
-	end do
-      enddo !rowind
-    enddo !colind
-    endif
-    deallocate(gp%points)
-    deallocate(gp%weights)
-
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !    
@@ -463,6 +403,8 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)
 
   end subroutine print_cfl
 
+
+
   ! ================================
   ! global_maximum:
   !

components/homme/src/share/prim_driver_base.F90

@@ -749,6 +749,7 @@ subroutine prim_init2(elem, hybrid, nets, nete, tl, hvcoord)
     use model_init_mod,       only: model_init2
     use time_mod,             only: timelevel_t, tstep, timelevel_init, nendstep, smooth, nsplit, TimeLevel_Qdp
     use control_mod,          only: smooth_phis_numcycle
+    use viscosity_mod,        only: dss_hvtensor
 
 #ifdef TRILINOS
     use prim_derived_type_mod ,only : derived_type, initialize
@@ -980,9 +981,10 @@ end subroutine noxinit
     endif
 
     call model_init2(elem(:), hybrid,deriv1,hvcoord,tl,nets,nete)
+    ! apply dss and bilinear projection to tensor coefficients 
+    call dss_hvtensor(elem,hybrid,nets,nete)
 
     ! advective and viscious CFL estimates
-    ! may also adjust tensor coefficients based on CFL
     call print_cfl(elem,hybrid,nets,nete,dtnu)
 
     ! Use the flexible time stepper if dt_remap_factor == 0 (vertically Eulerian

components/homme/src/share/viscosity_base.F90

@@ -37,7 +37,7 @@ module viscosity_base
 public :: neighbor_minmax_start,neighbor_minmax_finish
 public :: smooth_phis
 #endif
-
+public :: dss_hvtensor
 
 !
 ! compute vorticity/divergence and then project to make continious
@@ -64,7 +64,6 @@ module viscosity_base
 public :: compute_div_C0_contra     ! velocity around in contra-coordinates
 public :: compute_eta_C0_contra
 
-type (EdgeBuffer_t)          :: edge1
 
 contains
 
@@ -856,8 +855,97 @@ subroutine neighbor_minmax(elem,hybrid,edge3,nets,nete,nt,min_neigh,max_neigh,mi
 end subroutine
 
 
+#endif
+
+
+  subroutine dss_hvtensor(elem,hybrid,nets,nete)
+!
+!   estimate various CFL limits
+!   also, for variable resolution viscosity coefficient, make sure
+!   worse viscosity CFL (given by dtnu) is not violated by reducing 
+!   viscosity coefficient in regions where CFL is violated
+!
+    use kinds,       only : real_kind
+    use hybrid_mod,  only : hybrid_t
+    use element_mod, only : element_t
+    use dimensions_mod, only : np,ne
+    use quadrature_mod, only : gausslobatto, quadrature_t
+
+    use control_mod, only : hypervis_scaling
+    use edgetype_mod, only : EdgeBuffer_t 
+    use edge_mod, only : initedgebuffer, FreeEdgeBuffer, edgeVpack, edgeVunpack
+    use bndry_mod, only : bndry_exchangeV
+
+
+    type(element_t)      , intent(inout) :: elem(:)
+    integer              , intent(in) :: nets,nete
+    type (hybrid_t)      , intent(in) :: hybrid
+
+
+    real (kind=real_kind) :: x, y, noreast, nw, se, sw
+    real (kind=real_kind), dimension(np,np,nets:nete) :: zeta
+    integer :: ie,corner, i, j, rowind, colind
+    type (quadrature_t)    :: gp
+    type (EdgeBuffer_t)          :: edgebuf
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  TENSOR, RESOLUTION-AWARE HYPERVISCOSITY
+!  The tensorVisc() array is computed in cube_mod.F90
+!  this block of code will DSS it so the tensor if C0
+!  and also make it bilinear in each element.
+!  Oksana Guba
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   if (hypervis_scaling /= 0) then
+    gp=gausslobatto(np)
+    call initEdgeBuffer(hybrid%par,edgebuf,elem,1)
+    do rowind=1,2
+      do colind=1,2
+	do ie=nets,nete
+	  zeta(:,:,ie) = elem(ie)%tensorVisc(:,:,rowind,colind)*elem(ie)%spheremp(:,:)
+	  call edgeVpack(edgebuf,zeta(1,1,ie),1,0,ie)
+	end do
+
+	call bndry_exchangeV(hybrid,edgebuf)
+	do ie=nets,nete
+	  call edgeVunpack(edgebuf,zeta(1,1,ie),1,0,ie)
+          elem(ie)%tensorVisc(:,:,rowind,colind) = zeta(:,:,ie)*elem(ie)%rspheremp(:,:)
+	end do
+      enddo !rowind
+    enddo !colind
+    call FreeEdgeBuffer(edgebuf)
+
+!IF BILINEAR MAP OF V NEEDED
+    do rowind=1,2
+      do colind=1,2
+    ! replace hypervis w/ bilinear based on continuous corner values
+	do ie=nets,nete
+	  noreast = elem(ie)%tensorVisc(np,np,rowind,colind)
+	  nw = elem(ie)%tensorVisc(1,np,rowind,colind)
+	  se = elem(ie)%tensorVisc(np,1,rowind,colind)
+	  sw = elem(ie)%tensorVisc(1,1,rowind,colind)
+	  do i=1,np
+	    x = gp%points(i)
+	    do j=1,np
+		y = gp%points(j)
+		elem(ie)%tensorVisc(i,j,rowind,colind) = 0.25d0*( &
+					(1.0d0-x)*(1.0d0-y)*sw + &
+					(1.0d0-x)*(y+1.0d0)*nw + &
+					(x+1.0d0)*(1.0d0-y)*se + &
+					(x+1.0d0)*(y+1.0d0)*noreast)
+	    end do
+	  end do
+	end do
+      enddo !rowind
+    enddo !colind
+    deallocate(gp%points)
+    deallocate(gp%weights)
+    endif
+  end subroutine dss_hvtensor
+
+
+
+
 
 
 
-#endif
 end module viscosity_base

components/homme/src/sweqx/sweq_mod.F90

@@ -18,6 +18,8 @@ subroutine sweq(elem,edge1,edge2,edge3,red,par,ithr,nets,nete)
     !-----------------
     use derivative_mod, only : derivative_t, derivinit, allocate_subcell_integration_matrix
     !-----------------
+    use viscosity_mod, only : dss_hvtensor
+    !-----------------
     use dimensions_mod, only : np, nlev, npsq, nelemd
     !-----------------
     use shallow_water_mod, only : tc1_init_state, tc2_init_state, tc5_init_state, tc6_init_state, tc5_invariants, &
@@ -192,7 +194,7 @@ end subroutine noxfinish
     hybrid = hybrid_create(par,ithr,hthreads)
     simday=0
        call test_global_integral(elem,hybrid,nets,nete)
-
+       call dss_hvtensor(elem,hybrid,nets,nete)
        dtnu = 2.0d0*tstep*max(nu,nu_s)/hypervis_subcycle
        call print_cfl(elem,hybrid,nets,nete,dtnu)
 
@@ -829,6 +831,8 @@ subroutine sweq_rk(elem, edge1,edge2,edge3,red,par,ithr,nets,nete)
     !-----------------
     use derivative_mod, only : derivative_t, derivinit
     !-----------------
+    use viscosity_mod, only : dss_hvtensor
+    !-----------------
     use dimensions_mod, only :  np, nlev, npsq
     !-----------------
     use shallow_water_mod, only : tc1_init_state, tc2_init_state, tc5_init_state, &
@@ -932,6 +936,7 @@ subroutine sweq_rk(elem, edge1,edge2,edge3,red,par,ithr,nets,nete)
     hybrid = hybrid_create(par,ithr,hthreads)
 
        call test_global_integral(elem,hybrid,nets,nete,mindx)
+       call dss_hvtensor(elem,hybrid,nets,nete)
        dtnu = (tstep/rk_stage_user)*max(nu,nu_s)/hypervis_subcycle
        call print_cfl(elem,hybrid,nets,nete,dtnu)
 

components/homme/src/sweqx/viscosity_mod.F90

@@ -5,7 +5,7 @@
 module viscosity_mod
   use viscosity_base, only: compute_zeta_C0, compute_div_C0, &
     compute_zeta_C0_contra, compute_eta_C0_contra, compute_div_C0_contra, make_c0, neighbor_minmax, &
-    make_c0_vector
+    make_c0_vector, dss_hvtensor
 
 use dimensions_mod, only : np, nlev,qsize,nelemd
 use hybrid_mod, only : hybrid_t
@@ -23,6 +23,7 @@ module viscosity_mod
 implicit none
 
 public :: compute_pv_C0_contra
+public :: dss_hvtensor
 
 contains
 

components/homme/src/theta-l/share/viscosity_mod.F90

@@ -4,6 +4,6 @@
 
 module viscosity_mod
   use viscosity_base, only: compute_zeta_C0, compute_div_C0, compute_zeta_C0_contra, compute_div_C0_contra, make_c0, make_c0_vector, neighbor_minmax
-  use viscosity_base, only: biharmonic_wk_scalar, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis
+  use viscosity_base, only: biharmonic_wk_scalar, neighbor_minmax_start,neighbor_minmax_finish, smooth_phis, dss_hvtensor
   implicit none
 end module viscosity_mod