From 4050797091c67b4374c53d2009b330dee773a9b7 Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Fri, 12 Jan 2024 23:05:25 +0100 Subject: [PATCH 1/4] Fix initial forward Euler step under `rk_*`. --- src/rk.f90 | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/rk.f90 b/src/rk.f90 index 9a2d1b8..baa2aba 100644 --- a/src/rk.f90 +++ b/src/rk.f90 @@ -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 @@ -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) @@ -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) From 6890121da2913062adc020d8689cfe9f3d31dcdb Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Fri, 12 Jan 2024 23:44:20 +0100 Subject: [PATCH 2/4] Buoyancy term missed diffision by `rho`. --- src/mom.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mom.f90 b/src/mom.f90 index 2542861..59f3a62 100644 --- a/src/mom.f90 +++ b/src/mom.f90 @@ -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 @@ -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 From e768c9f606f45dbc6a1ab322a82f08aa9c203ced Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Sat, 13 Jan 2024 00:00:10 +0100 Subject: [PATCH 3/4] Removed redundant line. --- src/main.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/main.f90 b/src/main.f90 index 4f36c1e..eaa41c1 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -260,7 +260,7 @@ program cans nb,is_bound,cbcvel,cbcpre,bcvel,bcpre) #endif ! - call acdi_set_epsilon(dl,dzfi,acdi_eps_factor,seps) + !call acdi_set_epsilon(dl,dzfi,acdi_eps_factor,seps) ! fexts(1) = 'u' fexts(2) = 'v' @@ -335,10 +335,10 @@ program cans ! ! VoF update comes here! (discuss) ! - call tm_2fl(tm_coeff,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi) - call cmpt_norm_curv(n,dl,dli,dzc,dzf,dzci,dzfi,psi,kappa) - call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa) +! call tm_2fl(tm_coeff,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi) +! call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi) +! call cmpt_norm_curv(n,dl,dli,dzc,dzf,dzci,dzfi,psi,kappa) +! call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa) #if defined(_SCALAR) call tm_scal(tm_coeff,n,dli,dzci,dzfi,dt,0._rp,rho12,ka12,cp12,psi,u,v,w,s) call boundp(cbcsca,n,bcsca,nb,is_bound,dl,dzc,s) @@ -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 ! @@ -389,9 +388,10 @@ 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) + !call acdi_set_gamma(n,acdi_gam_factor,u,v,w,gam) if(myid == 0) print*, 'Gamma = ', gam, 'Epsilon = ', seps if(myid == 0) print*, 'Checking stability and divergence...' call chkdt(n,dl,dzci,dzfi,is_solve_ns,mu12,rho12,sigma,gacc,u,v,w,dtmax,gam,seps) !add the scalar time step check @@ -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 From d81352ff8b9baa2c566dbe08592fa56f00fc1cda Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Tue, 16 Jan 2024 14:06:06 +0100 Subject: [PATCH 4/4] Forgot to uncomment some code. --- src/main.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/main.f90 b/src/main.f90 index eaa41c1..000c1d6 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -260,7 +260,7 @@ program cans nb,is_bound,cbcvel,cbcpre,bcvel,bcpre) #endif ! - !call acdi_set_epsilon(dl,dzfi,acdi_eps_factor,seps) + call acdi_set_epsilon(dl,dzfi,acdi_eps_factor,seps) ! fexts(1) = 'u' fexts(2) = 'v' @@ -335,10 +335,10 @@ program cans ! ! VoF update comes here! (discuss) ! -! call tm_2fl(tm_coeff,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi) -! call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi) -! call cmpt_norm_curv(n,dl,dli,dzc,dzf,dzci,dzfi,psi,kappa) -! call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa) + call tm_2fl(tm_coeff,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi) + call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,psi) + call cmpt_norm_curv(n,dl,dli,dzc,dzf,dzci,dzfi,psi,kappa) + call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa) #if defined(_SCALAR) call tm_scal(tm_coeff,n,dli,dzci,dzfi,dt,0._rp,rho12,ka12,cp12,psi,u,v,w,s) call boundp(cbcsca,n,bcsca,nb,is_bound,dl,dzc,s) @@ -391,7 +391,7 @@ program cans 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) + call acdi_set_gamma(n,acdi_gam_factor,u,v,w,gam) if(myid == 0) print*, 'Gamma = ', gam, 'Epsilon = ', seps if(myid == 0) print*, 'Checking stability and divergence...' call chkdt(n,dl,dzci,dzfi,is_solve_ns,mu12,rho12,sigma,gacc,u,v,w,dtmax,gam,seps) !add the scalar time step check