From c260994a1d80836933a5b680dfa016a1d881d87e Mon Sep 17 00:00:00 2001 From: gianlupo Date: Wed, 6 Nov 2024 17:45:16 +0100 Subject: [PATCH] Enforce BC for the normal at walls --- src/bound.f90 | 2 +- src/main.f90 | 13 +++++++------ src/param.f90 | 44 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 7 deletions(-) diff --git a/src/bound.f90 b/src/bound.f90 index de8d638..ea04b69 100644 --- a/src/bound.f90 +++ b/src/bound.f90 @@ -80,7 +80,7 @@ end subroutine bounduvw ! subroutine boundp(cbc,n,bc,nb,is_bound,dl,dzc,p) ! - ! imposes pressure boundary conditions + ! imposes pressure and scalar boundary conditions ! implicit none character(len=1), intent(in), dimension(0:1,3) :: cbc diff --git a/src/main.f90 b/src/main.f90 index 66d30ae..75e932c 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -49,6 +49,7 @@ program cans use mod_output , only: out0d,gen_alias,out1d,out1d_chan,out2d,out3d,write_log_output,write_visu_2d,write_visu_3d use mod_param , only: rkcoeff,small, & nb,is_bound,cbcvel,bcvel,cbcpre,bcpre,cbcsca,bcsca,cbcpsi,bcpsi, & + cbcnox,bcnox,cbcnoy,bcnoy,cbcnoz,bcnoz, & icheck,iout0d,iout1d,iout2d,iout3d,isave, & nstep,time_max,tw_max,stop_type,restart,is_overwrite_save,nsaves_max, & datadir, & @@ -327,9 +328,9 @@ program cans call acdi_cmpt_phi(n,seps,psi,phi) call cmpt_norm_curv(n,dli,dzci,dzfi,phi,normx,normy,normz,kappa) call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,kappa) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normx) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normy) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normz) + call boundp(cbcnox,n,bcnox,nb,is_bound,dl,dzc,normx) + call boundp(cbcnoy,n,bcnoy,nb,is_bound,dl,dzc,normy) + call boundp(cbcnoz,n,bcnoz,nb,is_bound,dl,dzc,normz) ! #if defined(_CONSERVATIVE_MOMENTUM) !$acc kernels default(present) async(1) @@ -391,9 +392,9 @@ program cans call acdi_cmpt_phi(n,seps,psi,phi) call cmpt_norm_curv(n,dli,dzci,dzfi,phi,normx,normy,normz,kappa) call boundp(cbcpsi,n,bcpre,nb,is_bound,dl,dzc,kappa) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normx) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normy) - call boundp(cbcpsi,n,bcpsi,nb,is_bound,dl,dzc,normz) + call boundp(cbcnox,n,bcnox,nb,is_bound,dl,dzc,normx) + call boundp(cbcnoy,n,bcnoy,nb,is_bound,dl,dzc,normy) + call boundp(cbcnoz,n,bcnoz,nb,is_bound,dl,dzc,normz) end if #if defined(_SCALAR) call tm_scal(tm_coeff,n,dli,dzci,dzfi,dt,ssource,rho12,ka12,cp12,psi,u,v,w,s) diff --git a/src/param.f90 b/src/param.f90 index df4f201..6de016c 100644 --- a/src/param.f90 +++ b/src/param.f90 @@ -55,6 +55,12 @@ module mod_param real(rp) , protected, dimension(0:1,3) :: bcpsi character(len=1), protected, dimension(0:1,3) :: cbcsca real(rp) , protected, dimension(0:1,3) :: bcsca +character(len=1), protected, dimension(0:1,3) :: cbcnox +real(rp) , protected, dimension(0:1,3) :: bcnox +character(len=1), protected, dimension(0:1,3) :: cbcnoy +real(rp) , protected, dimension(0:1,3) :: bcnoy +character(len=1), protected, dimension(0:1,3) :: cbcnoz +real(rp) , protected, dimension(0:1,3) :: bcnoz ! real(rp), protected, dimension(3) :: bforce,gacc real(rp), protected :: ssource @@ -167,6 +173,44 @@ subroutine read_input(myid) #if defined(_CONSTANT_COEFFS_POISSON) rho0 = minval(rho12(:)) #endif + ! + ! set normals BC based on cbcpsi + ! + block + integer :: b,c + logical :: is_wall + cbcnox = cbcpsi + cbcnoy = cbcpsi + cbcnoz = cbcpsi + bcnox = bcpsi + bcnoy = bcpsi + bcnoz = bcpsi + do b = 0,1 + do c = 1,3 + is_wall = cbcpsi(b,c)=='N' + select case(c) + case(1) + if (is_wall) then + cbcnox(b,1) = 'D' + cbcnoy(b,1) = 'I' + cbcnoz(b,1) = 'I' + end if + case(2) + if (is_wall) then + cbcnox(b,2) = 'I' + cbcnoy(b,2) = 'D' + cbcnoz(b,2) = 'I' + end if + case(3) + if (is_wall) then + cbcnox(b,3) = 'I' + cbcnoy(b,3) = 'I' + cbcnoz(b,3) = 'D' + end if + end select + end do + end do + end block #if defined(_OPENACC) ! ! read cuDecomp parameter file cudecomp.in, if it exists