added more versatile inflow BCs

and the Kovaznay flow case

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

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

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]

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

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)

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