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

Fix initial forward Euler step under rk_*. #19

Merged
merged 4 commits into from
Jan 17, 2024
Merged

Fix initial forward Euler step under rk_*. #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4050797
Fix initial forward Euler step under `rk_*`.
p-costa Jan 12, 2024
6890121
Buoyancy term missed diffision by `rho`.
p-costa Jan 12, 2024
e768c9f
Removed redundant line.
p-costa Jan 12, 2024
d81352f
Forgot to uncomment some code.
p-costa Jan 16, 2024
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
3 changes: 1 addition & 2 deletions src/main.f90
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
Loading