Skip to content

Commit

Permalink
Bugfix in scalar transport and minor polishing.
Browse files Browse the repository at this point in the history
  • Loading branch information
@p-costa
p-costa committed Nov 13, 2024
1 parent 3f92b3f commit 4063c50
Show file tree
Hide file tree
Showing 4 changed files with 27 additions and 24 deletions.
26 changes: 14 additions & 12 deletions src/main.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,8 +188,8 @@ program cans
rhsbp%z(n(1),n(2),0:1))
allocate(psi,phi,kappa,normx,normy,normz,mold=pp)
#if defined(_CONSERVATIVE_MOMENTUM)
allocate(psio,mold=pp)
allocate(acdi_rgx,acdi_rgy,acdi_rgz,mold=pp)
if(.not. allocated(psio)) allocate(psio(0:n(1)+1,0:n(2)+1,0:n(3)+1))
#endif
#if defined(_DEBUG)
if(myid == 0) print*, 'This executable of CaNS was built with compiler: ', compiler_version()
Expand Down Expand Up @@ -318,13 +318,15 @@ program cans
!$acc end kernels
#endif
#if defined(_SCALAR)
!$acc enter data copyin(s)
!$acc enter data copyin(s) async(1)
call boundp(cbcsca,n,bcsca,nb,is_bound,dl,dzc,s)
!$acc wait
#endif
!$acc enter data copyin(psi) create(phi,kappa,normx,normy,normz)
!$acc enter data create(acdi_rgx,acdi_rgy,acdi_rgz)
!$acc enter data create(psio)
!$acc enter data copyin(psi) create(phi,kappa,normx,normy,normz) async(1)
!$acc enter data create(acdi_rgx,acdi_rgy,acdi_rgz) async(1)
!$acc enter data create(psio) async(1)
call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi)
!$acc wait
!
call acdi_cmpt_phi(n,seps,psi,phi)
call cmpt_norm_curv(n,dli,dzci,dzfi,phi,normx,normy,normz,kappa)
Expand All @@ -347,7 +349,7 @@ program cans
!
write(fldnum,'(i7.7)') istep
!$acc wait
!$acc update self(u,v,w,p,psi,kappa)
!$acc update self(u,v,w,p,psi,kappa,s)
#include "out1d.h90"
#include "out2d.h90"
#include "out3d.h90"
Expand Down Expand Up @@ -387,7 +389,7 @@ program cans
if(is_track_interface) then
call tm_2fl(tm_coeff,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,normx,normy,normz,phi,psi,acdi_rgx,acdi_rgy,acdi_rgz)
#if defined(_CONSERVATIVE_MOMENTUM)
call bounduvw(cbcvel,n,bcvel,nb,is_bound,.false.,dl,dzc,dzf,acdi_rgx,acdi_rgy,acdi_rgz)
call bounduvw(cbcnor,n,bcnor,nb,is_bound,.false.,dl,dzc,dzf,acdi_rgx,acdi_rgy,acdi_rgz)
#endif
call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi)
call acdi_cmpt_phi(n,seps,psi,phi)
Expand All @@ -412,7 +414,7 @@ program cans
call bulk_mean_12(n,grid_vol_ratio_c,psi,rho12,rho_av)
end if
call tm(tm_coeff,n,dli,dzci,dzfi,dt,dt_r, &
bforce,gacc,sigma,rho_av,rho12,mu12,beta12,rho0,psi,kappa,s,p,pn,po,psio, &
bforce,gacc,sigma,rho_av,rho12,mu12,beta12,rho0,psi,kappa,p,pn,po,psio,s, &
acdi_rgx,acdi_rgy,acdi_rgz,u,v,w)
if(is_forced_hit) then
call lscale_forcing(2,lo,hi,0.5_rp,dtrk,l,dl,zc,zf,u,v,w)
Expand Down Expand Up @@ -501,17 +503,17 @@ program cans
write(fldnum,'(i7.7)') istep
if(mod(istep,iout1d) == 0) then
!$acc wait
!$acc update self(u,v,w,p,psi,kappa)
!$acc update self(u,v,w,p,psi,kappa,s)
#include "out1d.h90"
end if
if(mod(istep,iout2d) == 0) then
!$acc wait
!$acc update self(u,v,w,p,psi,kappa)
!$acc update self(u,v,w,p,psi,kappa,s)
#include "out2d.h90"
end if
if(mod(istep,iout3d) == 0) then
!$acc wait
!$acc update self(u,v,w,p,psi,kappa)
!$acc update self(u,v,w,p,psi,kappa,s)
#include "out3d.h90"
end if
if(mod(istep,isave ) == 0.or.(is_done.and..not.kill)) then
Expand All @@ -531,7 +533,7 @@ program cans
end if
end if
!$acc wait
!$acc update self(u,v,w,p,psi)
!$acc update self(u,v,w,p,psi,s)
call load_one('w',trim(datadir)//trim(filename)//'_'//trim(fexts(1))//'.bin',MPI_COMM_WORLD,ng,[1,1,1],lo,hi,u,time,istep)
call load_one('w',trim(datadir)//trim(filename)//'_'//trim(fexts(2))//'.bin',MPI_COMM_WORLD,ng,[1,1,1],lo,hi,v,time,istep)
call load_one('w',trim(datadir)//trim(filename)//'_'//trim(fexts(3))//'.bin',MPI_COMM_WORLD,ng,[1,1,1],lo,hi,w,time,istep)
Expand Down
4 changes: 2 additions & 2 deletions src/mom.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -992,8 +992,8 @@ subroutine mom_xyz_oth(n,dli,dzci,dzfi,dt_r,rho12,beta12,bforce,gacc,sigma,rho0,
real(rp), intent(in ), dimension(3) :: bforce,gacc
real(rp), intent(in ) :: sigma
real(rp), intent(in ) :: rho0,rho_av
real(rp), intent(in ), dimension(0:,0:,0:) :: p,psi,kappa,s
real(rp), intent(in ), dimension(0:,0:,0:), optional :: pn,po
real(rp), intent(in ), dimension(0:,0:,0:) :: p,psi,kappa
real(rp), intent(in ), dimension(0:,0:,0:), optional :: s,pn,po
real(rp), intent(inout), dimension( :, :, :) :: dudt,dvdt,dwdt
integer :: i,j,k
real(rp) :: rho,drho,rhobeta,drhobeta
Expand Down
13 changes: 7 additions & 6 deletions src/rk.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module mod_rk
public rk,rk_scal,rk_2fl
contains
subroutine rk(rkpar,n,dli,dzci,dzfi,dt,dt_r, &
bforce,gacc,sigma,rho_av,rho12,mu12,beta12,rho0,psi,kappa,s,p,pn,po,psio, &
bforce,gacc,sigma,rho_av,rho12,mu12,beta12,rho0,psi,kappa,p,pn,po,psio,s, &
acdi_rgx,acdi_rgy,acdi_rgz,u,v,w)
!
! RK3 scheme for time integration of the momentum equations
Expand All @@ -29,9 +29,9 @@ subroutine rk(rkpar,n,dli,dzci,dzfi,dt,dt_r, &
real(rp), intent(in ) :: sigma,rho_av
real(rp), intent(in ), dimension(2) :: rho12,mu12,beta12
real(rp), intent(in ) :: rho0
real(rp), intent(in ), dimension(0:,0:,0:) :: psi,kappa,s
real(rp), intent(in ), dimension(0:,0:,0:) :: psi,kappa
real(rp), intent(in ), dimension(0:,0:,0:) :: p
real(rp), intent(in ), dimension(0:,0:,0:), optional :: pn,po,psio
real(rp), intent(in ), dimension(0:,0:,0:), optional :: pn,po,psio,s
real(rp), intent(in ), dimension(0:,0:,0:), optional :: acdi_rgx,acdi_rgy,acdi_rgz
real(rp), intent(inout), dimension(0:,0:,0:) :: u,v,w
real(rp), target , allocatable, dimension(:,:,:), save :: dudtrk_t ,dvdtrk_t ,dwdtrk_t , &
Expand Down Expand Up @@ -63,7 +63,7 @@ subroutine rk(rkpar,n,dli,dzci,dzfi,dt,dt_r, &
!
! advection, diffusion and regularization terms
!
call mom_xyz_ad(n,dli,dzci,dzfi,rho12,mu12,acdi_rgx,acdi_rgy,acdi_rgz,u,v,w,psi,psio(:,:,:),dudtrk,dvdtrk,dwdtrk)
call mom_xyz_ad(n,dli,dzci,dzfi,rho12,mu12,acdi_rgx,acdi_rgy,acdi_rgz,u,v,w,psi,psio,dudtrk,dvdtrk,dwdtrk)
!
if(is_first) then
!$acc kernels default(present) async(1)
Expand Down Expand Up @@ -156,17 +156,18 @@ subroutine rk_scal(rkpar,n,dli,dzci,dzfi,dt, &
logical, save :: is_first = .true.
real(rp) :: factor1,factor2,factor12
integer :: i,j,k
real(rp), dimension(2) :: rhocp12
real(rp), dimension(2), save :: rhocp12
!
factor1 = rkpar(1)*dt
factor2 = rkpar(2)*dt
factor12 = factor1 + factor2
rhocp12 = rho12(:)*cp12(:)
if(is_first) then ! leverage save attribute to allocate these arrays on the device only once
allocate(dsdtrk_t(n(1),n(2),n(3)),dsdtrko_t(n(1),n(2),n(3)))
!$acc enter data create(dsdtrk_t,dsdtrko_t) async(1)
dsdtrk => dsdtrk_t
dsdtrko => dsdtrko_t
rhocp12(:) = rho12(:)*cp12(:)
!$acc wait(1)
end if
call scal_ad(n(1),n(2),n(3),dli(1),dli(2),dzci,dzfi,ssource,ka12,rhocp12,psi,u,v,w,s,dsdtrk)
if(is_first) then
Expand Down
8 changes: 4 additions & 4 deletions src/scal.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ subroutine scal_ad(nx,ny,nz,dxi,dyi,dzci,dzfi,ssource,ka12,rhocp12,psi,u,v,w,s,d
integer :: i,j,k
real(rp) :: usip,usim,vsjp,vsjm,wskp,wskm
real(rp) :: dsdxp,dsdxm,dsdyp,dsdym,dsdzp,dsdzm
real(rp) :: kaxp,kaxm,kayp,kaym,kazp,kazm
real(rp) :: kaxp,kaxm,kayp,kaym,kazp,kazm,rhocp_c
real(rp) :: ka,dka,rhocp,drhocp
!
ka = ka12(2) ; dka = ka12(1)-ka12(2)
Expand All @@ -32,7 +32,7 @@ subroutine scal_ad(nx,ny,nz,dxi,dyi,dzci,dzfi,ssource,ka12,rhocp12,psi,u,v,w,s,d
!$acc parallel loop collapse(3) default(present) &
!$acc private(usip,usim,vsjp,vsjm,wskp,wskm) &
!$acc private(dsdxp,dsdxm,dsdyp,dsdym,dsdzp,dsdzm) &
!$acc private(kaxp,kaxm,kayp,kaym,kazp,kazm,rhocp) async(1)
!$acc private(kaxp,kaxm,kayp,kaym,kazp,kazm,rhocp_c) async(1)
do k=1,nz
do j=1,ny
do i=1,nx
Expand All @@ -57,14 +57,14 @@ subroutine scal_ad(nx,ny,nz,dxi,dyi,dzci,dzfi,ssource,ka12,rhocp12,psi,u,v,w,s,d
dsdzp = (s(i,j,k+1)-s(i,j,k ))*dzci(k )
dsdzm = (s(i,j,k )-s(i,j,k-1))*dzci(k-1)
!
rhocp = rhocp+drhocp*psi(i,j,k)
rhocp_c = rhocp+drhocp*psi(i,j,k)
!
dsdt(i,j,k) = dxi*( -usip + usim ) + &
dyi*( -vsjp + vsjm ) + &
dzfi(k)*( -wskp + wskm ) + &
( (kaxp*dsdxp-kaxm*dsdxm)*dxi + &
(kayp*dsdyp-kaym*dsdym)*dyi + &
(kazp*dsdzp-kazm*dsdzm)*dzfi(k) + ssource )/rhocp
(kazp*dsdzp-kazm*dsdzm)*dzfi(k) + ssource )/rhocp_c
end do
end do
end do
Expand Down

0 comments on commit 4063c50

Please sign in to comment.