From dccbfddc4cc480cf387444e51aa2382f3ec62005 Mon Sep 17 00:00:00 2001 From: Pedro Costa Date: Mon, 6 Dec 2021 13:02:15 +0000 Subject: [PATCH] added more versatile inflow BCs and the Kovaznay flow case --- .../kovasznay_flow/dns.in | 8 ++ .../kovasznay_flow/geo/block.001 | 18 ++++ .../kovasznay_flow/hypre.in | 1 + src/bound.f90 | 32 ++++++-- src/initflow.f90 | 82 +++++++++++++++---- src/main.f90 | 63 ++++++++++---- 6 files changed, 161 insertions(+), 43 deletions(-) create mode 100644 examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/dns.in create mode 100644 examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/geo/block.001 create mode 100644 examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/hypre.in diff --git a/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/dns.in b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/dns.in new file mode 100644 index 0000000..c21ebcd --- /dev/null +++ b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/dns.in @@ -0,0 +1,8 @@ +.95 1.0e5 ! cfl, dtmin +1. 1. .025 ! uref, lref, rey +50000 100. 0.1 ! nstep,time_max,tw_max +T F F ! stop_type(1:3) +F T ! restart,is_overwrite_save +10 10 200000 50000000 100 200 ! icheck, iout0d, iout1d, iout2d, iout3d, isave +0. 0. 0. ! bforce(1:3) +4 ! nthreadsmax diff --git a/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/geo/block.001 b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/geo/block.001 new file mode 100644 index 0000000..b5bfc0e --- /dev/null +++ b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/geo/block.001 @@ -0,0 +1,18 @@ +1 1 1 ! dims(1:3) +1 1 1 ! lo(1:3) +320 2 64 ! hi(1:3) +-0.5 0. -0.5 ! lmin(1:3) +4.5 .03125 0.5 ! lmax(1:3) +0 0 0 ! gt(1:3) +0. 0. 0. ! gr(1:3) +D N F F F F ! cbcvel(0:1,1:3,1) [u BC type] +D N F F F F ! cbcvel(0:1,1:3,2) [v BC type] +D N F F F F ! cbcvel(0:1,1:3,3) [w BC type] +N D F F F F ! cbcpre(0:1,1:3 ) [p BC type] +1. 0. 1. 1. 1. 1. ! bcvel(0:1,1:3,1) [u BC value] +0. 0. 1. 1. 1. 1. ! bcvel(0:1,1:3,2) [v BC value] +0. 0. 1. 1. 1. 1. ! bcvel(0:1,1:3,3) [w BC value] +0. 0. 1. 1. 1. 1. ! bcpre(0:1,1:3 ) [p BC value] +2 0 0 0 0 0 ! is_inflow(0:1,1:3) +kov ! inivel +1 ! id diff --git a/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/hypre.in b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/hypre.in new file mode 100644 index 0000000..c6aa43d --- /dev/null +++ b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/hypre.in @@ -0,0 +1 @@ +1 1.e-4 50 ! hypre_solver_i,hypre_tol,hypre_maxiter | note: hypre_solver_i = [1->SMG;2->PFMG;3->GMRES;4->BiCGSTAB] diff --git a/src/bound.f90 b/src/bound.f90 index 2634a6c..6a2c86b 100644 --- a/src/bound.f90 +++ b/src/bound.f90 @@ -423,27 +423,41 @@ subroutine cmpt_estimated_traction(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr) end select end subroutine cmpt_estimated_traction ! - subroutine inflow(is_inflow_bound,lo,hi,velx,vely,velz,u,v,w) + subroutine inflow(is_inflow_bound,is_correc,lo,hi,uin_x,vin_x,win_x, & + uin_y,vin_y,win_y, & + uin_z,vin_z,win_z,u,v,w) ! implicit none logical , intent(in ), dimension(0:1,1:3) :: is_inflow_bound + logical , intent(in ) :: is_correc integer , intent(in ), dimension(3 ) :: lo,hi - real(rp), intent(in ), dimension(lo(2)-1:,lo(3)-1:,0:) :: velx - real(rp), intent(in ), dimension(lo(1)-1:,lo(3)-1:,0:) :: vely - real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,0:) :: velz + real(rp), intent(in ), dimension(lo(2)-1:,lo(3)-1:,0:) :: uin_x,vin_x,win_x + real(rp), intent(in ), dimension(lo(1)-1:,lo(3)-1:,0:) :: uin_y,vin_y,win_y + real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,0:) :: uin_z,vin_z,win_z real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u,v,w - integer :: ibound,idir,q + integer :: ibound,idir,qn,qt,dq do idir=1,3 do ibound=0,1 if(is_inflow_bound(ibound,idir)) then - q = (1-ibound)*(lo(idir)-1)+ibound*hi(idir) + qn = (1-ibound)*(lo(idir)-1)+ibound*(hi(idir) ) + qt = (1-ibound)*(lo(idir)-1)+ibound*(hi(idir)+1) + if(ibound == 0) dq = 1; if(ibound == 1) dq = -1 select case(idir) case(1) - u(q,:,:) = velx(:,:,ibound) + if(.not.is_correc) & + u(qn,:,:) = uin_x(:,:,ibound) + v(qt,:,:) = -v(qt+dq,:,:) + 2._rp*vin_x(:,:,ibound) + w(qt,:,:) = -w(qt+dq,:,:) + 2._rp*win_x(:,:,ibound) case(2) - v(:,q,:) = vely(:,:,ibound) + u(:,qt,:) = -u(:,qt+dq,:) + 2._rp*uin_y(:,:,ibound) + if(.not.is_correc) & + v(:,qn,:) = vin_y(:,:,ibound) + w(:,qt,:) = -w(:,qt+dq,:) + 2._rp*win_y(:,:,ibound) case(3) - w(:,:,q) = velz(:,:,ibound) + u(:,:,qt) = -u(:,:,qt+dq) + 2._rp*uin_z(:,:,ibound) + v(:,:,qt) = -v(:,:,qt+dq) + 2._rp*vin_z(:,:,ibound) + if(.not.is_correc) & + w(:,:,qn) = win_z(:,:,ibound) end select end if end do diff --git a/src/initflow.f90 b/src/initflow.f90 index 6aea33a..4144c20 100644 --- a/src/initflow.f90 +++ b/src/initflow.f90 @@ -26,7 +26,7 @@ subroutine initflow(inivel,is_wallturb,lo,hi,lo_g,hi_g,l,uref,lref,visc,bforce, real(rp), allocatable, dimension(:) :: u1d integer :: i,j,k logical :: is_noise,is_mean,is_pair - real(rp) :: reb,retau + real(rp) :: reb,retau,lambda real(rp) :: xcl,xfl,ycl,yfl,zcl,zfl ! allocate(u1d(lo(3):hi(3))) @@ -70,6 +70,25 @@ subroutine initflow(inivel,is_wallturb,lo,hi,lo_g,hi_g,l,uref,lref,visc,bforce, end do end do end do + case('kov') + reb = uref*lref/visc + lambda = reb/2-sqrt(reb**2/4+4*pi**2) + do k=lo(3),hi(3) + zcl = zc(k) + zfl = zf(k) + do j=lo(2),hi(2) + ycl = yc(j) + yfl = yf(j) + do i=lo(1),hi(1) + xcl = xc(i) + xfl = xf(i) + u(i,j,k) = 1._rp -exp(lambda*xfl)*cos(2*pi*zcl) + v(i,j,k) = 0._rp + w(i,j,k) = lambda/(2._rp*pi)*exp(lambda*xcl)*sin(2*pi*zfl) + p(i,j,k) = 0._rp !.5_rp*(1._rp-exp(2._rp*lambda*xcl)) + end do + end do + end do case('pdc') if(is_wallturb) then ! turbulent flow retau = sqrt(bforce*lref)*uref/visc @@ -88,7 +107,7 @@ subroutine initflow(inivel,is_wallturb,lo,hi,lo_g,hi_g,l,uref,lref,visc,bforce, call MPI_FINALIZE() error stop end select - if(inivel /= 'tgv') then + if(all(inivel /= ['tgv','kov'])) then do k=lo(3),hi(3) do j=lo(2),hi(2) do i=lo(1),hi(1) @@ -268,7 +287,8 @@ function poiseuille_square(r,lmin,lmax) result(vel) vel = 6._rp*(rr*(1._rp-rr)) end function poiseuille_square ! - subroutine init_inflow(periods,lo,hi,lmin,lmax,x1c,x2c,uref,vel) + subroutine init_inflow(inflow_type,periods,lo,hi,lmin,lmax,x1c,x1f,x2c,x2f, & + uref,lref,visc,veln,velt1,velt2) ! ! below an inflow is calculated for a 2d plane ! later this subroutine can be generalized with other shapes of @@ -278,27 +298,55 @@ subroutine init_inflow(periods,lo,hi,lmin,lmax,x1c,x2c,uref,vel) ! there is no danger of experiencing a singularity '0._rp**0'. ! implicit none + integer , parameter :: INFLOW_POISEUILLE = 1, & + INFLOW_KOVASZNAY = 2 + ! + integer , intent(in ) :: inflow_type integer , intent(in ), dimension(2 ) :: periods integer , intent(in ), dimension(2 ) :: lo,hi real(rp), intent(in ), dimension(2 ) :: lmin,lmax - real(rp), intent(in ), dimension(lo(1)-1:) :: x1c - real(rp), intent(in ), dimension(lo(2)-1:) :: x2c - real(rp), intent(in ) :: uref ! target bulk velocity - real(rp), intent(out), dimension(lo(1)-1:,lo(2)-1:) :: vel + real(rp), intent(in ), dimension(lo(1)-1:) :: x1c,x1f + real(rp), intent(in ), dimension(lo(2)-1:) :: x2c,x2f + real(rp), intent(in ) :: uref,lref,visc ! target bulk velocity + real(rp), intent(out), dimension(lo(1)-1:,lo(2)-1:) :: veln,velt1,velt2 integer, dimension(2) :: q integer :: i1,i2 - ! real(rp), external :: poiseuille_square - ! procedure (), pointer :: inflow_function => null() - ! if(inflow_type == POISEUILLE) inflow_function => poiseuille_square + real(rp) :: lambda,reb,x0 ! - q(:) = 1 - where(periods(:) /= 0) q(:) = 0 - do i2 = lo(2)-1,hi(2)+1 - do i1 = lo(1)-1,hi(1)+1 - vel(i1,i2) = uref*poiseuille_square(x1c(i1),lmin(1),lmax(1))**q(1) * & - poiseuille_square(x2c(i2),lmin(2),lmax(2))**q(2) + select case(inflow_type) + case(INFLOW_POISEUILLE) + q(:) = 1 + where(periods(:) /= 0) q(:) = 0 + do i2 = lo(2)-1,hi(2)+1 + do i1 = lo(1)-1,hi(1)+1 + veln(i1,i2) = uref*poiseuille_square(x1c(i1),lmin(1),lmax(1))**q(1) * & + poiseuille_square(x2c(i2),lmin(2),lmax(2))**q(2) + velt1(i1,i2) = 0._rp + velt2(i1,i2) = 0._rp + end do end do - end do + case(INFLOW_KOVASZNAY) + reb = uref*lref/visc + lambda = reb/2-sqrt(reb**2/4+4*pi**2) + x0 = -0.5 + do i2 = lo(2)-1,hi(2)+1 + do i1 = lo(1)-1,hi(1)+1 + veln( i1,i2) = 1._rp- exp(lambda*x0)*cos(2*pi*x2c(i2)) + velt1(i1,i2) = 0._rp + velt2(i1,i2) = lambda/(2._rp*pi)*exp(lambda*x0)*sin(2*pi*x2f(i2)) + end do + end do + case default + q(:) = 1 + where(periods(:) /= 0) q(:) = 0 + do i2 = lo(2)-1,hi(2)+1 + do i1 = lo(1)-1,hi(1)+1 + veln(i1,i2) = 0._rp + velt1(i1,i2) = 0._rp + velt2(i1,i2) = 0._rp + end do + end do + end select end subroutine init_inflow ! subroutine log_profile(lo,hi,zc,l,uref,lref,visc,p) diff --git a/src/main.f90 b/src/main.f90 index a7fbadd..36e7873 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -85,6 +85,7 @@ program snac real(rp), allocatable, dimension(:,:,:) :: z end type rhs_bound type(rhs_bound) :: rhsp + type(rhs_bound) :: u_in,v_in,w_in real(rp) :: alpha,alpha_bc(0:1,1:3) #ifdef _IMPDIFF type(rhs_bound) :: rhsu,rhsv,rhsw,bcu,bcv,bcw @@ -220,9 +221,15 @@ program snac wp(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1), & pp(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1), & po(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1)) - allocate(velin_x( lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1,0:1), & - velin_y( lo(1)-1:hi(1)+1,lo(3)-1:hi(3)+1,0:1), & - velin_z( lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,0:1)) + allocate(u_in%x(lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1,0:1), & + u_in%y(lo(1)-1:hi(1)+1,lo(3)-1:hi(3)+1,0:1), & + u_in%z(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,0:1), & + v_in%x(lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1,0:1), & + v_in%y(lo(1)-1:hi(1)+1,lo(3)-1:hi(3)+1,0:1), & + v_in%z(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,0:1), & + w_in%x(lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1,0:1), & + w_in%y(lo(1)-1:hi(1)+1,lo(3)-1:hi(3)+1,0:1), & + w_in%z(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,0:1)) #ifdef _IMPDIFF allocate(uo(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1), & vo(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1), & @@ -521,9 +528,15 @@ program snac call bounduvw(cbcvel,lo,hi,bcvel,.false.,halos,is_bound,nb, & dxc,dxf,dyc,dyf,dzc,dzf,u,v,w) call boundp( cbcpre,lo,hi,bcpre,halos,is_bound,nb,dxc,dyc,dzc,p) - velin_x(:,:,:) = 0._rp - velin_y(:,:,:) = 0._rp - velin_z(:,:,:) = 0._rp + u_in%x(:,:,:) = 0._rp + v_in%x(:,:,:) = 0._rp + w_in%x(:,:,:) = 0._rp + u_in%y(:,:,:) = 0._rp + v_in%y(:,:,:) = 0._rp + w_in%y(:,:,:) = 0._rp + u_in%z(:,:,:) = 0._rp + v_in%z(:,:,:) = 0._rp + w_in%z(:,:,:) = 0._rp #ifdef _IMPDIFF do ib=0,1 bcu%x(:,:,ib) = bcvel(ib,1,1) @@ -544,30 +557,41 @@ program snac select case(idir) case(1) il = 2;iu = 3;iskip = 1 - call init_inflow(periods(il:iu:iskip),lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & - yc,zc,bcvel(ib,idir,idir),velin_x(:,:,ib)) + call init_inflow(inflow_type(ib,idir),periods(il:iu:iskip), & + lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & + yc,yf,zc,zf,bcvel(ib,idir,idir),lref,visc,u_in%x(:,:,ib),v_in%x(:,:,ib),w_in%x(:,:,ib)) #ifdef _IMPDIFF - bcu%x(:,:,ib) = velin_x(lo(2):hi(2),lo(3):hi(3),ib) + bcu%x(:,:,ib) = u_in%x(lo(2):hi(2),lo(3):hi(3),ib) + bcv%x(:,:,ib) = v_in%x(lo(2):hi(2),lo(3):hi(3),ib) + bcw%x(:,:,ib) = w_in%x(lo(2):hi(2),lo(3):hi(3),ib) #endif case(2) il = 1;iu = 3;iskip = 2 - call init_inflow(periods(il:iu:iskip),lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & - xc,zc,bcvel(ib,idir,idir),velin_y(:,:,ib)) + call init_inflow(inflow_type(ib,idir),periods(il:iu:iskip), & + lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & + xc,xf,zc,zf,bcvel(ib,idir,idir),lref,visc,v_in%y(:,:,ib),u_in%y(:,:,ib),w_in%y(:,:,ib)) #ifdef _IMPDIFF - bcv%y(:,:,ib) = velin_y(lo(1):hi(1),lo(3):hi(3),ib) + bcu%y(:,:,ib) = u_in%y(lo(1):hi(1),lo(3):hi(3),ib) + bcv%y(:,:,ib) = v_in%y(lo(1):hi(1),lo(3):hi(3),ib) + bcw%y(:,:,ib) = w_in%y(lo(1):hi(1),lo(3):hi(3),ib) #endif case(3) il = 1;iu = 2;iskip = 1 - call init_inflow(periods(il:iu:iskip),lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & - xc,yc,bcvel(ib,idir,idir),velin_z(:,:,ib)) + call init_inflow(inflow_type(ib,idir),periods(il:iu:iskip), & + lo(il:iu:iskip),hi(il:iu:iskip),lmin(il:iu:iskip),lmax(il:iu:iskip), & + xc,xf,yc,yf,bcvel(ib,idir,idir),lref,visc,w_in%z(:,:,ib),u_in%z(:,:,ib),v_in%z(:,:,ib)) #ifdef _IMPDIFF - bcw%z(:,:,ib) = velin_z(lo(1):hi(1),lo(2):hi(2),ib) + bcu%z(:,:,ib) = u_in%z(lo(1):hi(1),lo(2):hi(2),ib) + bcv%z(:,:,ib) = v_in%z(lo(1):hi(1),lo(2):hi(2),ib) + bcw%z(:,:,ib) = w_in%z(lo(1):hi(1),lo(2):hi(2),ib) #endif end select end if end do end do - call inflow(is_bound_inflow,lo,hi,velin_x,velin_y,velin_z,u,v,w) + call inflow(is_bound_inflow,.false.,lo,hi,u_in%x,v_in%x,w_in%x, & + u_in%y,v_in%y,w_in%y, & + u_in%z,v_in%z,w_in%z,u,v,w) up(:,:,:) = 0._rp vp(:,:,:) = 0._rp wp(:,:,:) = 0._rp @@ -909,7 +933,9 @@ program snac #endif call bounduvw(cbcvel,lo,hi,bcvel,.false.,halos,is_bound,nb, & dxc,dxf,dyc,dyf,dzc,dzf,up,vp,wp) - call inflow(is_bound_inflow,lo,hi,velin_x,velin_y,velin_z,up,vp,wp) + call inflow(is_bound_inflow,.false.,lo,hi,u_in%x,v_in%x,w_in%x, & + u_in%y,v_in%y,w_in%y, & + u_in%z,v_in%z,w_in%z,up,vp,wp) #if !defined(_IMPDIFF) && defined(_ONE_PRESS_CORR) dtrk = dt if(irk < 3) then ! pressure correction only at the last RK step @@ -947,6 +973,9 @@ program snac call correc(lo,hi,dxc,dyc,dzc,dtrk,pp,up,vp,wp,u,v,w) call bounduvw(cbcvel,lo,hi,bcvel,.true.,halos,is_bound,nb, & dxc,dxf,dyc,dyf,dzc,dzf,u,v,w) + call inflow(is_bound_inflow,.true.,lo,hi,u_in%x,v_in%x,w_in%x, & + u_in%y,v_in%y,w_in%y, & + u_in%z,v_in%z,w_in%z,u,v,w) call updt_pressure(lo,hi,dxc,dxf,dyc,dyf,dzc,dzf,alpha,pp,p) call boundp( cbcpre,lo,hi,bcpre,halos,is_bound,nb,dxc,dyc,dzc,p) end do