Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

extract code from print_cfl into dss_hvtensor routine #6723

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

extract code from print_cfl into dss_hvtensor routine #6723

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion components/homme/src/preqx/viscosity_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions components/homme/src/preqx_acc/viscosity_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
70 changes: 6 additions & 64 deletions components/homme/src/share/global_norms_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
module global_norms_mod

use kinds, only : iulog
use edgetype_mod, only : EdgeBuffer_t

implicit none
private
Expand All @@ -25,7 +24,6 @@ module global_norms_mod
public :: wrap_repro_sum

private :: global_maximum
type (EdgeBuffer_t), private :: edgebuf

contains

Expand Down Expand Up @@ -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(:)
Expand All @@ -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
Expand Down Expand Up @@ -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(:)
Expand All @@ -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


Expand Down Expand Up @@ -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
Expand All @@ -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)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
Expand Down Expand Up @@ -463,6 +403,8 @@ subroutine print_cfl(elem,hybrid,nets,nete,dtnu)

end subroutine print_cfl



! ================================
! global_maximum:
!
Expand Down
4 changes: 3 additions & 1 deletion components/homme/src/share/prim_driver_base.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
94 changes: 91 additions & 3 deletions components/homme/src/share/viscosity_base.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
7 changes: 6 additions & 1 deletion components/homme/src/sweqx/sweq_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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)

Expand Down
3 changes: 2 additions & 1 deletion components/homme/src/sweqx/viscosity_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -23,6 +23,7 @@ module viscosity_mod
implicit none

public :: compute_pv_C0_contra
public :: dss_hvtensor

contains

Expand Down
2 changes: 1 addition & 1 deletion components/homme/src/theta-l/share/viscosity_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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