Skip to content

Commit

Permalink
Merge pull request #19 from p-costa/fix-euler-forward
Browse files Browse the repository at this point in the history
Fix initial forward Euler step under `rk_*`.
  • Loading branch information
gianlupo authored Jan 17, 2024
2 parents 5da3f00 + d81352f commit 2373663
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 21 deletions.
3 changes: 1 addition & 2 deletions src/main.f90
Original file line number Diff line number Diff line change
Expand Up @@ -375,7 +375,6 @@ program cans
call updatep(pp,p)
call boundp(cbcpre,n,bcpre,nb,is_bound,dl,dzc,p)
end if
dto = dt
!
! check simulation stopping criteria
!
Expand All @@ -389,6 +388,7 @@ program cans
tw = (MPI_WTIME()-twi)/3600.
if(tw >= tw_max ) is_done = is_done.or..true.
end if
dto = dt
if(mod(istep,icheck) == 0) then
if(myid == 0) print*, 'Calculating maximum velocity to set ACDI gamma parameter...'
call acdi_set_gamma(n,acdi_gam_factor,u,v,w,gam)
Expand All @@ -403,7 +403,6 @@ program cans
is_done = .true.
kill = .true.
end if
dto = dt
dti = 1./dt
call chkdiv(lo,hi,dli,dzfi,u,v,w,divtot,divmax)
if(myid == 0) print*, 'Total divergence = ', divtot, '| Maximum divergence = ', divmax
Expand Down
12 changes: 6 additions & 6 deletions src/mom.f90
Original file line number Diff line number Diff line change
Expand Up @@ -826,9 +826,9 @@ subroutine mom_xyz_oth(nx,ny,nz,dxi,dyi,dzci,rho12,beta12,bforce,gacc,sigma,rho0
factorxp = rhobeta + drhobeta*psixp
factoryp = rhobeta + drhobeta*psiyp
factorzp = rhobeta + drhobeta*psizp
dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)
dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)
dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)
dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)/rhoxp
dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)/rhoyp
dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)/rhozp
#endif
dudt(i,j,k) = dudt_aux
dvdt(i,j,k) = dvdt_aux
Expand Down Expand Up @@ -1101,9 +1101,9 @@ subroutine mom_xyz_all(nx,ny,nz,dxi,dyi,dzci,dzfi,rho12,mu12,beta12,bforce,gacc,
factorxp = rhobeta + drhobeta*psixp
factoryp = rhobeta + drhobeta*psiyp
factorzp = rhobeta + drhobeta*psizp
dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)
dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)
dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)
dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)/rhoxp
dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)/rhoyp
dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)/rhozp
#endif
dudt(i,j,k) = dudt_aux
dvdt(i,j,k) = dvdt_aux
Expand Down
25 changes: 12 additions & 13 deletions src/rk.f90
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +51,6 @@ subroutine rk(rkpar,n,dli,dzci,dzfi,dt, &
allocate(dudtrko_t(n(1),n(2),n(3)),dvdtrko_t(n(1),n(2),n(3)),dwdtrko_t(n(1),n(2),n(3)))
!$acc enter data create(dudtrk_t ,dvdtrk_t ,dwdtrk_t ) async(1)
!$acc enter data create(dudtrko_t,dvdtrko_t,dwdtrko_t) async(1)
!$acc kernels default(present) async(1) ! not really necessary
dudtrko_t(:,:,:) = dudtrk_t(:,:,:)
dvdtrko_t(:,:,:) = dvdtrk_t(:,:,:)
dwdtrko_t(:,:,:) = dwdtrk_t(:,:,:)
!$acc end kernels
dudtrk => dudtrk_t
dvdtrk => dvdtrk_t
dwdtrk => dwdtrk_t
Expand Down Expand Up @@ -146,16 +141,18 @@ subroutine rk_scal(rkpar,n,dli,dzci,dzfi,dt, &
factor12 = factor1 + factor2
rhocp12 = rho12(:)*cp12(:)
if(is_first) then ! leverage save attribute to allocate these arrays on the device only once
is_first = .false.
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)
!$acc kernels default(present) async(1) ! not really necessary
dsdtrko_t(:,:,:) = 0._rp
!$acc end kernels
dsdtrk => dsdtrk_t
dsdtrko => dsdtrko_t
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 ! use Euler forward
!$acc kernels
dsdtrko(:,:,:) = dsdtrk(:,:,:)
!$acc end kernels
is_first = .false.
end if
!$acc parallel loop collapse(3) default(present) async(1)
do k=1,n(3)
do j=1,n(2)
Expand Down Expand Up @@ -195,16 +192,18 @@ subroutine rk_2fl(rkpar,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi)
factor2 = rkpar(2)*dt
factor12 = factor1 + factor2
if(is_first) then ! leverage save attribute to allocate these arrays on the device only once
is_first = .false.
allocate(dpsidtrk_t(n(1),n(2),n(3)),dpsidtrko_t(n(1),n(2),n(3)))
!$acc enter data create(dpsidtrk_t,dpsidtrko_t) async(1)
!$acc kernels default(present) async(1) ! not really necessary
dpsidtrko_t(:,:,:) = 0._rp
!$acc end kernels
dpsidtrk => dpsidtrk_t
dpsidtrko => dpsidtrko_t
end if
call acdi_transport_pf(n(1),n(2),n(3),dli(1),dli(2),dli(3),dzci,dzfi,gam,seps,u,v,w,psi,dpsidtrk)
if(is_first) then ! use Euler forward
!$acc kernels
dpsidtrko(:,:,:) = dpsidtrk(:,:,:)
!$acc end kernels
is_first = .false.
end if
!$acc parallel loop collapse(3) default(present) async(1)
!$OMP PARALLEL DO COLLAPSE(3) DEFAULT(shared)
do k=1,n(3)
Expand Down

0 comments on commit 2373663

Please sign in to comment.