Skip to content

Commit

Permalink
polished outflow-related branch w/ several fixes along the way
Browse files Browse the repository at this point in the history 
  * such as fixing integer overflow under load.f90
  * some fixes for the traction outflow BCs
  * remove printing of wall forces by default
  * add more general inflow routines
  * add Kovaznay benchmark
  • Loading branch information
@p-costa
p-costa committed Feb 10, 2022
1 parent 9e69a02 commit b9d0d78
Show file tree
Hide file tree
Showing 4 changed files with 77 additions and 100 deletions.
100 changes: 34 additions & 66 deletions src/bound.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,7 +265,7 @@ subroutine set_bc(ctype,ibound,lo,hi,idir,centered,rvalue,dr,p)
end select
end subroutine set_bc
!
subroutine set_open_bc_uvw(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr_x,tr_y,tr_z)
subroutine set_open_bc_uvw(ibound,lo,hi,idir,visc,dr,u,v,w,p,up,vp,wp,tr_x,tr_y,tr_z)
!
! a zero or estimated-traction open BC (Bozonnet et al., JCP 2021)
! the latter serves well as a robust outflow BC;
Expand All @@ -277,14 +277,13 @@ subroutine set_open_bc_uvw(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr_x,tr_y,tr_z)
integer , intent(in ), dimension(3) :: lo,hi
integer , intent(in ) :: idir
real(rp), intent(in ) :: visc,dr
real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: p
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u,v,w
real(rp), intent(out ), dimension(lo(2)-1:,lo(3)-1:,0:), optional :: tr_x
real(rp), intent(out ), dimension(lo(1)-1:,lo(3)-1:,0:), optional :: tr_y
real(rp), intent(out ), dimension(lo(1)-1:,lo(2)-1:,0:), optional :: tr_z
real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u ,v ,w ,p
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: up,vp,wp
real(rp), intent(inout), dimension(lo(2)-1:,lo(3)-1:,0:), optional :: tr_x
real(rp), intent(inout), dimension(lo(1)-1:,lo(3)-1:,0:), optional :: tr_y
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,0:), optional :: tr_z
real(rp) :: factor,norm
integer :: q
logical :: is_estimated_traction
!
if( ibound == 0) then
norm = 1._rp
Expand All @@ -294,90 +293,59 @@ subroutine set_open_bc_uvw(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr_x,tr_y,tr_z)
q = hi(idir)
end if
factor = -norm*dr/(2*visc)
is_estimated_traction = .false.
!
select case(idir)
case(1)
if(present(tr_x)) is_estimated_traction = .true.
if (ibound == 0) then
!$OMP WORKSHARE
u(q,:,:) = u(q+1,:,:) + factor*(.5_rp*min(0._rp,norm*u(q+1,:,:))**2 + p(q+1,:,:))
tr_x(:,:,ibound) = tr_x(:,:,ibound) + .5_rp*min(0._rp,norm*u(q+1,:,:))**2
up(q,:,:) = up(q+1,:,:) + factor*(tr_x(:,:,ibound) + p(q+1,:,:))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
u(q,:,:) = u(q,:,:) + factor*tr_x(:,:,ibound)
!$OMP END WORKSHARE
end if
else if(ibound == 1) then
!$OMP WORKSHARE
u(q,:,:) = u(q-1,:,:) + factor*(.5_rp*min(0._rp,norm*u(q-1,:,:))**2 + p(q ,:,:))
tr_x(:,:,ibound) = tr_x(:,:,ibound) + .5_rp*min(0._rp,norm*u(q-1,:,:))**2
up(q,:,:) = up(q-1,:,:) + factor*(tr_x(:,:,ibound) + p(q ,:,:))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
u(q,:,:) = u(q,:,:) + factor*tr_x(:,:,ibound)
!$OMP END WORKSHARE
end if
!$OMP WORKSHARE
u(hi(idir)+1,:,:) = u(hi(idir) ,:,:) ! not needed
up(hi(idir)+1,:,:) = up(hi(idir) ,:,:) ! not needed
!$OMP END WORKSHARE
end if
case(2)
if(present(tr_y)) is_estimated_traction = .true.
if (ibound == 0) then
!$OMP WORKSHARE
v(:,q,:) = v(:,q+1,:) + factor*(.5_rp*min(0._rp,norm*v(:,q+1,:))**2 + p(:,q+1,:))
tr_y(:,:,ibound) = tr_y(:,:,ibound) + .5_rp*min(0._rp,norm*v(:,q+1,:))**2
vp(:,q,:) = vp(:,q+1,:) + factor*(tr_y(:,:,ibound) + p(:,q+1,:))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
v(:,q,:) = v(:,q,:) + factor*tr_y(:,:,ibound)
!$OMP END WORKSHARE
end if
else if(ibound == 1) then
!$OMP WORKSHARE
v(:,q,:) = v(:,q-1,:) + factor*(.5_rp*min(0._rp,norm*v(:,q-1,:))**2 + p(:,q ,:))
tr_y(:,:,ibound) = tr_y(:,:,ibound) + .5_rp*min(0._rp,norm*v(:,q-1,:))**2
vp(:,q,:) = vp(:,q-1,:) + factor*(tr_y(:,:,ibound) + p(:,q ,:))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
v(:,q,:) = v(:,q,:) + factor*tr_y(:,:,ibound)
!$OMP END WORKSHARE
end if
!$OMP WORKSHARE
v(:,hi(idir)+1,:) = v(:,hi(idir) ,:) ! not needed
vp(:,hi(idir)+1,:) = vp(:,hi(idir) ,:) ! not needed
!$OMP END WORKSHARE
end if
case(3)
if(present(tr_z)) is_estimated_traction = .true.
if (ibound == 0) then
!$OMP WORKSHARE
w(:,:,q) = w(:,:,q+1) + factor*(.5_rp*min(0._rp,norm*w(:,:,q+1))**2 + p(:,:,q+1))
tr_z(:,:,ibound) = tr_z(:,:,ibound) + .5_rp*min(0._rp,norm*w(:,:,q+1))**2
wp(:,:,q) = wp(:,:,q+1) + factor*(tr_z(:,:,ibound) + p(:,:,q+1))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
w(:,:,q) = w(:,:,q) + factor*tr_z(:,:,ibound)
!$OMP END WORKSHARE
end if
else if(ibound == 1) then
!$OMP WORKSHARE
w(:,:,q) = w(:,:,q-1) + factor*(.5_rp*min(0._rp,norm*w(:,:,q-1))**2 + p(:,:,q ))
tr_z(:,:,ibound) = tr_z(:,:,ibound) + .5_rp*min(0._rp,norm*w(:,:,q-1))**2
wp(:,:,q) = wp(:,:,q-1) + factor*(tr_z(:,:,ibound) + p(:,:,q ))
!$OMP END WORKSHARE
if(is_estimated_traction) then
!$OMP WORKSHARE
w(:,:,q) = w(:,:,q) + factor*tr_z(:,:,ibound)
!$OMP END WORKSHARE
end if
!$OMP WORKSHARE
w(:,:,hi(idir)+1) = w(:,:,hi(idir) ) ! not needed
wp(:,:,hi(idir)+1) = wp(:,:,hi(idir) ) ! not needed
!$OMP END WORKSHARE
end if
end select
end subroutine set_open_bc_uvw
!
subroutine set_open_bc_p(ibound,lo,hi,idir,drc,drf,alpha,p)
!
! a zero or estimated-traction open BC (Bozonnet et al., JCP 2021)
! the latter serves well as a robust outflow BC;
! the estimated traction can be computed in subroutine
! `cmpt_estimated_traction`
! traction BC for the correction pressure (Bozonnet et al., JCP 2021)
!
implicit none
integer , intent(in ) :: ibound
Expand Down Expand Up @@ -427,7 +395,7 @@ subroutine set_open_bc_p(ibound,lo,hi,idir,drc,drf,alpha,p)
end select
end subroutine set_open_bc_p
!
subroutine cmpt_estimated_traction(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr_x,tr_y,tr_z)
subroutine cmpt_estimated_traction(ibound,lo,hi,idir,visc,dr,u,v,w,p,tr_x,tr_y,tr_z)
!
! computes the estimated traction from an interior grid cell
! adjacent to the boundary grid cell, to be used for computing the
Expand All @@ -438,8 +406,7 @@ subroutine cmpt_estimated_traction(ibound,lo,hi,idir,visc,dr,p,u,v,w,tr_x,tr_y,t
integer , intent(in ), dimension(3) :: lo,hi
integer , intent(in ) :: idir
real(rp), intent(in ) :: visc,dr ! dr -- grid spacing of the interior grid cell
real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: p
real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u,v,w
real(rp), intent(in ), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u,v,w,p
real(rp), intent(inout), dimension(lo(2)-1:,lo(3)-1:,0:), optional :: tr_x
real(rp), intent(inout), dimension(lo(1)-1:,lo(3)-1:,0:), optional :: tr_y
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,0:), optional :: tr_z
Expand Down Expand Up @@ -527,17 +494,18 @@ subroutine inflow(is_inflow_bound,is_correc,lo,hi,uin_x,vin_x,win_x, &
end do
end subroutine inflow
!
subroutine outflow(is_outflow_bound,is_estimated_traction,lo,hi,dl,visc,tr_x,tr_y,tr_z,u,v,w,p)
subroutine outflow(is_outflow_bound,is_estimated_traction,lo,hi,dl,dl_h,visc,u,v,w,p,tr_x,tr_y,tr_z,up,vp,wp)
!
implicit none
logical , intent(in ), dimension(0:1,1:3) :: is_outflow_bound,is_estimated_traction
integer , intent(in ), dimension(3 ) :: lo,hi
real(rp), intent(in ), dimension(0:1,1:3) :: dl
real(rp), intent(in ), dimension(0:1,1:3) :: dl,dl_h
real(rp), intent(in ) :: visc
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u ,v ,w ,p
real(rp), intent(inout), dimension(lo(2)-1:,lo(3)-1:,0:) :: tr_x
real(rp), intent(inout), dimension(lo(1)-1:,lo(3)-1:,0:) :: tr_y
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,0:) :: tr_z
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: u,v,w,p
real(rp), intent(inout), dimension(lo(1)-1:,lo(2)-1:,lo(3)-1:) :: up,vp,wp
integer :: ib,idir
!
do idir=1,3
Expand All @@ -547,18 +515,18 @@ subroutine outflow(is_outflow_bound,is_estimated_traction,lo,hi,dl,visc,tr_x,tr_
case(1)
tr_x(:,:,ib) = 0._rp
if( is_estimated_traction(ib,idir) ) &
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_x=tr_x)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_x=tr_x)
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl_h(ib,idir),u,v,w,p,tr_x=tr_x)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),u,v,w,p,up,vp,wp,tr_x=tr_x)
case(2)
tr_y(:,:,ib) = 0._rp
if( is_estimated_traction(ib,idir) ) &
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_y=tr_y)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_y=tr_y)
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl_h(ib,idir),u,v,w,p,tr_y=tr_y)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),u,v,w,p,up,vp,wp,tr_y=tr_y)
case(3)
tr_z(:,:,ib) = 0._rp
if( is_estimated_traction(ib,idir) ) &
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_z=tr_z)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),p,u,v,w,tr_z=tr_z)
call cmpt_estimated_traction(ib,lo,hi,idir,visc,dl_h(ib,idir),u,v,w,p,tr_z=tr_z)
call set_open_bc_uvw(ib,lo,hi,idir,visc,dl(ib,idir),u,v,w,p,up,vp,wp,tr_z=tr_z)
end select
end if
end do
Expand Down
6 changes: 3 additions & 3 deletions src/initflow.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,9 +120,9 @@ subroutine initflow(inivel,is_wallturb,lo,hi,lo_g,hi_g,l,uref,lref,visc,bforce,
end do
end if
if(is_noise) then
call add_noise(lo,hi,lo_g,hi_g,123,.5_rp,u)
call add_noise(lo,hi,lo_g,hi_g,456,.5_rp,v)
call add_noise(lo,hi,lo_g,hi_g,789,.5_rp,w)
call add_noise(lo,hi,lo_g,hi_g,123,.05_rp,u)
call add_noise(lo,hi,lo_g,hi_g,456,.05_rp,v)
call add_noise(lo,hi,lo_g,hi_g,789,.05_rp,w)
end if
if(is_mean) then
call set_mean(lo,hi,l,dxc,dyf,dzf,uref,u)
Expand Down
43 changes: 26 additions & 17 deletions src/main.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ program snac
type(MPI_DATATYPE) , dimension( 3) :: halos
integer , dimension( 3) :: ng,lo_g,hi_g,lo_1,hi_1
real(rp), allocatable, dimension(:,:,:) :: u,v,w,p,up,vp,wp,pp,po
real(rp), allocatable, dimension(:,:,:) :: velin_x,velin_y,velin_z,tr_x,tr_y,tr_z
real(rp), allocatable, dimension(:,:,:) :: tr_x,tr_y,tr_z
#ifdef _IMPDIFF
real(rp), allocatable, dimension(:,:,:) :: uo,vo,wo
#endif
Expand All @@ -92,17 +92,21 @@ program snac
#ifdef _IMPDIFF
type(rhs_bound) :: rhsu,rhsv,rhsw,bcu,bcv,bcw
#endif
real(rp), dimension( 0:1,3) :: dl,dl_outflow
real(rp), dimension( 0:1,3) :: dl,dl_outflow,dl_outflow_h
real(rp), dimension(0:1,0:1,3) :: dlc_outflow
logical :: outflow_exists
#ifdef _IMPDIFF
real(rp), dimension(0:1,3) :: dlu,dlv,dlw
integer , dimension( 3) :: hiu,hiv,hiw
#endif
type(hypre_solver) :: psolver
logical :: is_symm_matrix_p
#ifdef _IMPDIFF
type(hypre_solver) :: usolver,vsolver,wsolver
real(rp) :: alphai,alphaoi
logical :: is_symm_matrix_u, &
is_symm_matrix_v, &
is_symm_matrix_w
#endif
!
real(rp) :: dt,dtmax,time,dtrk,divtot,divmax
Expand Down Expand Up @@ -605,15 +609,18 @@ program snac
!
! outflow BCs
!
allocate(tr_x,mold=velin_x)
allocate(tr_y,mold=velin_y)
allocate(tr_z,mold=velin_z)
allocate(tr_x,mold=u_in%x)
allocate(tr_y,mold=v_in%y)
allocate(tr_z,mold=w_in%z)
tr_x(:,:,:) = 0._rp
tr_y(:,:,:) = 0._rp
tr_z(:,:,:) = 0._rp
dl_outflow = reshape([dxf_g(lo_g(1)),dxf_g(hi_g(1)), &
dyf_g(lo_g(2)),dyf_g(hi_g(2)), &
dzf_g(lo_g(3)),dzf_g(hi_g(3))],shape(dl_outflow))
dl_outflow_h = reshape([dxf_g(lo_g(1)+1),dxf_g(hi_g(1)-1), &
dyf_g(lo_g(2)+1),dyf_g(hi_g(2)-1), &
dzf_g(lo_g(3)+1),dzf_g(hi_g(3)-1)],shape(dl_outflow))
dlc_outflow(0,:,:) = reshape([dxc_g(lo_g(1)-1),dxc_g(hi_g(1)-1), &
dyc_g(lo_g(2)-1),dyc_g(hi_g(2)-1), &
dzc_g(lo_g(3)-1),dzc_g(hi_g(3)-1)],shape(dl_outflow))
Expand All @@ -628,7 +635,6 @@ program snac
end do
end do
call MPI_ALLREDUCE(any(is_bound_outflow(:,:)),outflow_exists,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD)
call outflow(is_bound_outflow,is_estimated_traction,lo,hi,dl_outflow,visc,tr_x,tr_y,tr_z,u,v,w,p)
alpha_bc( :,:) = 0._rp
alpha_bc_o(:,:) = alpha_bc(:,:)
!
Expand Down Expand Up @@ -739,6 +745,7 @@ program snac
allocate(comms_fft(hi_a(idir)-lo_a(idir)+1))
allocate(lambda_p_a(hi_a(idir)-lo_a(idir)+1))
is_bound_a(:,:) = is_bound(:,:)
is_symm_matrix_p = is_uniform_grid.and..not.outflow_exists
#ifndef _FFT_USE_SLABS
comms_fft(:) = MPI_COMM_WORLD
lambda_p_a(:) = lambda_p
Expand All @@ -753,17 +760,17 @@ program snac
#endif
#ifndef _FFT_USE_SLICED_PENCILS
call init_n_2d_matrices(cbcpre(:,il:iu:iskip),bcpre(:,il:iu:iskip),dl(:,il:iu:iskip), &
is_uniform_grid,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
is_symm_matrix_p,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
lo_a(idir),hi_a(idir),lo_a(il:iu:iskip),hi_a(il:iu:iskip),periods(il:iu:iskip), &
dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,lambda_p_a,comms_fft,psolver_fft)
#else
call init_n_3d_matrices(idir,nslices,cbcpre,bcpre,dl,is_uniform_grid,is_bound,is_centered,lo,periods, &
call init_n_3d_matrices(idir,nslices,cbcpre,bcpre,dl,is_symm_matrix_p,is_bound,is_centered,lo,periods, &
lo_sp,hi_sp,dxc,dxf,dyc,dyf,dzc,dzf,alpha,alpha_bc,lambda_p_a,psolver_fft)
#endif
call create_n_solvers(npsolvers,hypre_maxiter,hypre_tol,hypre_solver_i,psolver_fft)
call setup_n_solvers(npsolvers,psolver_fft)
#else
call init_matrix_3d(cbcpre,bcpre,dl,is_uniform_grid,is_bound,is_centered,lo,hi,periods, &
call init_matrix_3d(cbcpre,bcpre,dl,is_symm_matrix_p,is_bound,is_centered,lo,hi,periods, &
dxc,dxf,dyc,dyf,dzc,dzf,alpha,alpha_bc,psolver,is_bound_outflow=is_bound_outflow)
call create_solver(hypre_maxiter,hypre_tol,hypre_solver_i,psolver)
call setup_solver(psolver)
Expand All @@ -775,6 +782,7 @@ program snac
hiu(:) = hi(:)
if(is_bound(1,1)) hiu(:) = hiu(:)-[1,0,0]
is_centered(:) = [.false.,.true.,.true.]
is_symm_matrix_u = is_uniform_grid
call init_bc_rhs(cbcvel(:,:,1),bcvel(:,:,1),dlu,is_bound,is_centered,lo,hiu,periods, &
dxf,dxc,dyc,dyf,dzc,dzf,rhsu%x,rhsu%y,rhsu%z,bcu%x,bcu%y,bcu%z)
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
Expand All @@ -792,11 +800,11 @@ program snac
if(periods(1) == 0) hiu_a(:) = hiu_a(:)-[1,0,0]
#endif
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,1),bcvel(:,il:iu:iskip,1),dlu(:,il:iu:iskip), &
is_uniform_grid,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
is_symm_matrix_u,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
lo_a(idir),hiu_a(idir),lo_a(il:iu:iskip),hiu_a(il:iu:iskip),periods(il:iu:iskip), &
dlu1_1,dlu1_2,dlu2_1,dlu2_2,alpha,alpha_bc,lambda_u_a,comms_fft,usolver_fft)
#else
call init_matrix_3d(cbcvel(:,:,1),bcvel(:,:,1),dlu,is_uniform_grid,is_bound,is_centered,lo,hiu,periods, &
call init_matrix_3d(cbcvel(:,:,1),bcvel(:,:,1),dlu,is_symm_matrix_u,is_bound,is_centered,lo,hiu,periods, &
dxf,dxc,dyc,dyf,dzc,dzf,alpha,alpha_bc,usolver)
#endif
dlv = reshape([dxc_g(lo_g(1)-1),dxc_g(hi_g(1)), &
Expand All @@ -820,13 +828,14 @@ program snac
lambda_v_a(:) = lambda_v(lo_s(idir)-lo(idir)+1:hi_s(idir)-lo(idir)+1)
hiv_a(:) = hi_s(:)
if(periods(2) == 0) hiv_a(:) = hiv_a(:)-[0,1,0]
is_symm_matrix_v = is_uniform_grid
#endif
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,2),bcvel(:,il:iu:iskip,2),dlv(:,il:iu:iskip), &
is_uniform_grid,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
is_symm_matrix_v,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
lo_a(idir),hiv_a(idir),lo_a(il:iu:iskip),hiv_a(il:iu:iskip),periods(il:iu:iskip), &
dlv1_1,dlv1_2,dlv2_1,dlv2_2,alpha,alpha_bc,lambda_v_a,comms_fft,vsolver_fft)
#else
call init_matrix_3d(cbcvel(:,:,2),bcvel(:,:,2),dlv,is_uniform_grid,is_bound,is_centered,lo,hiv,periods, &
call init_matrix_3d(cbcvel(:,:,2),bcvel(:,:,2),dlv,is_symm_matrix_v,is_bound,is_centered,lo,hiv,periods, &
dxc,dxf,dyf,dyc,dzc,dzf,alpha,alpha_bc,vsolver)
#endif
dlw = reshape([dxc_g(lo_g(1)-1),dxc_g(hi_g(1)), &
Expand All @@ -850,13 +859,14 @@ program snac
lambda_w_a(:) = lambda_w(lo_s(idir)-lo(idir)+1:hi_s(idir)-lo(idir)+1)
hiw_a(:) = hi_s(:)
if(periods(3) == 0) hiw_a(:) = hiw_a(:)-[0,0,1]
is_symm_matrix_w = is_uniform_grid
#endif
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,3),bcvel(:,il:iu:iskip,3),dlw(:,il:iu:iskip), &
is_uniform_grid,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
is_symm_matrix_w,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
lo_a(idir),hiw_a(idir),lo_a(il:iu:iskip),hiw_a(il:iu:iskip),periods(il:iu:iskip), &
dlw1_1,dlw1_2,dlw2_1,dlw2_2,alpha,alpha_bc,lambda_w_a,comms_fft,wsolver_fft)
#else
call init_matrix_3d(cbcvel(:,:,3),bcvel(:,:,3),dlw,is_uniform_grid,is_bound,is_centered,lo,hiw,periods, &
call init_matrix_3d(cbcvel(:,:,3),bcvel(:,:,3),dlw,is_symm_matrix_w,is_bound,is_centered,lo,hiw,periods, &
dxc,dxf,dyc,dyf,dzf,dzc,alpha,alpha_bc,wsolver)
#endif
#endif
Expand Down Expand Up @@ -983,7 +993,7 @@ program snac
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)
call outflow(is_bound_outflow,is_estimated_traction,lo,hi,dl_outflow,visc,tr_x,tr_y,tr_z,u,v,w,p)
call outflow(is_bound_outflow,is_estimated_traction,lo,hi,dl_outflow,dl_outflow_h,visc,u,v,w,p,tr_x,tr_y,tr_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
Expand Down Expand Up @@ -1057,7 +1067,6 @@ program snac
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)
call chkdiv(lo,hi,dxf,dyf,dzf,u,v,w,vol_all,MPI_COMM_WORLD,divtot,divmax)
end do
!
! check simulation stopping criteria
Expand Down
Loading

0 comments on commit b9d0d78

Please sign in to comment.