diff --git a/components/homme/src/preqx/viscosity_mod.F90 b/components/homme/src/preqx/viscosity_mod.F90 index 7d826b8a9be6..1068637e42f3 100644 --- a/components/homme/src/preqx/viscosity_mod.F90 +++ b/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 diff --git a/components/homme/src/preqx_acc/viscosity_mod.F90 b/components/homme/src/preqx_acc/viscosity_mod.F90 index ead8ce2e770c..9abc920fc58c 100644 --- a/components/homme/src/preqx_acc/viscosity_mod.F90 +++ b/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 diff --git a/components/homme/src/share/global_norms_mod.F90 b/components/homme/src/share/global_norms_mod.F90 index c3084e32803f..f9ca6b9479b9 100644 --- a/components/homme/src/share/global_norms_mod.F90 +++ b/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: ! diff --git a/components/homme/src/share/prim_driver_base.F90 b/components/homme/src/share/prim_driver_base.F90 index df83e1cc386b..cba9ef642072 100644 --- a/components/homme/src/share/prim_driver_base.F90 +++ b/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 diff --git a/components/homme/src/share/viscosity_base.F90 b/components/homme/src/share/viscosity_base.F90 index 9c42f57158f7..a0cb2056d1bf 100644 --- a/components/homme/src/share/viscosity_base.F90 +++ b/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 diff --git a/components/homme/src/sweqx/sweq_mod.F90 b/components/homme/src/sweqx/sweq_mod.F90 index 50fbc72ba41a..dc01dab92769 100644 --- a/components/homme/src/sweqx/sweq_mod.F90 +++ b/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) diff --git a/components/homme/src/sweqx/viscosity_mod.F90 b/components/homme/src/sweqx/viscosity_mod.F90 index 56fab03084af..ca6027487373 100644 --- a/components/homme/src/sweqx/viscosity_mod.F90 +++ b/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 diff --git a/components/homme/src/theta-l/share/viscosity_mod.F90 b/components/homme/src/theta-l/share/viscosity_mod.F90 index 0627ea2c9d9d..630c71d5475f 100644 --- a/components/homme/src/theta-l/share/viscosity_mod.F90 +++ b/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