momentum.f90 4.47 KiB
!####################################################################
!> @author Holger Grosshans
!> @brief one sweep to relax momentum eq. using Jacobi method
!> @param err: L2 error norm
subroutine momentum(err)
use var
real(kind=pr),dimension(ii,jj,ll) :: u1,v1,w1
real(kind=pr) :: err,syncSum
real(kind=pr) :: dum,dvm,dwm
real(kind=pr) :: du,dv,dw,dua,dva,dwa
real(kind=pr) :: tu,tv,tw,qu,qv,qw
integer :: i,j,l
err=0._pr
u1=0._pr; v1=0._pr; w1=0._pr
dum=0._pr; dvm=0._pr; dwm=0._pr
do 1 i=imin,imax,1
do 1 j=jmin,jmax,1
do 1 l=lmin,lmax,1
if (celltype(i,j,l).ne.active) cycle
dua=0._pr; dva=0._pr; dwa=0._pr
if (.not.((bcx.eq.'w').and.(myid.eq.nrprocs-1).and.(i.eq.imax))) then
call momx1(tu,qu,i,j,l)
du=(defect_u(i,j,l)+tu)/qu*urfu
u1(i,j,l)=u(i,j,l)+du
dua=abs(du)
if (dum.lt.dua) dum=dua
endif
if (.not.((bcy.eq.'w').and.(j.eq.jmax))) then
call momy1(tv,qv,i,j,l)
dv=(defect_v(i,j,l)+tv)/qv*urfv
v1(i,j,l)=v(i,j,l)+dv
dva=abs(dv)
if (dvm.lt.dva) dvm=dva
endif
if (.not.((bcz.eq.'w').and.(l.eq.lmax))) then
call momz1(tw,qw,i,j,l)
dw=(defect_w(i,j,l)+tw)/qw*urfw
w1(i,j,l)=w(i,j,l)+dw
dwa=abs(dw)
if (dwm.lt.dwa) dwm=dwa
endif
12 err=err+dua*dua+dva*dva+dwa*dwa
1 enddo
u=u1; v=v1; w=w1
call sync(u); call sync(v); call sync(w)
call bcUVW
err=(syncSum(err)/dimgptot/3._pr)**(0.5_pr)
! write(*,'(x,a,es9.2e2,a,3(es9.2e2))') &
! 'res. inner it. mom =',err,', max du,dv,dw =',dum,dvm,dwm
return
end
!####################################################################!> @author Holger Grosshans
!> @brief second-order time integration
subroutine ddt
use var
integer :: i,j,l
dudt=-(1._pr+tau01+tau02)*u+tau02*u01
dvdt=-(1._pr+tau01+tau02)*v+tau02*v01
dwdt=-(1._pr+tau01+tau02)*w+tau02*w01
u01=u
v01=v
w01=w
return
end
!####################################################################
!> @author Holger Grosshans
!> @brief HO-LO terms for deferred correction method
subroutine deferredCorrection
use var
real(kind=pr) :: &
tu,tv,tw,qu,qv,qw,ra,raHO,tuHO,tvHO,twHO,quHO,qvHO,qwHO
integer :: i,j,l
! velocity components are defined on ranges:
! u(imin:imax, jmin:jmax, lmin:lmax)
! v(imin:imax, jmin:jmax-1, lmin:lmax)
! w(imin:imax, jmin:jmax, lmin:lmax-1)
do 1 i=imin,imax
do 1 j=jmin+1,jmax-1
do 1 l=lmin+1,lmax-1
if (celltype(i,j,l).ne.active) cycle
if (((bcx.eq.'p'.or.bcx.eq.'i').or.((bcx.eq.'w').and.(i-2.ge.imin-1).and.(i+1.le.imax))).and. &
((bcy.eq.'p').or.((bcy.eq.'w').and.(j-2.ge.jmin-1).and.(j+1.le.jmax))).and. &
((bcz.eq.'p').or.((bcz.eq.'w').and.(l-2.ge.lmin-1).and.(l+1.le.lmax)))) then
call mass2(ra,i,j,l)
call mass4(raHO,i,j,l)
defect_c(i,j,l)=raHO-ra
else
defect_c(i,j,l)=0._pr
endif
if (((bcx.eq.'p'.or.bcx.eq.'i').or.((bcx.eq.'w').and.(i-3.ge.imin).and.(i+3.le.imax-1))).and. &
((bcy.eq.'p').or.((bcy.eq.'w').and.(j-3.ge.jmin).and.(j+3.le.jmax))).and. &
((bcz.eq.'p').or.((bcz.eq.'w').and.(l-3.ge.lmin).and.(l+3.le.lmax)))) then
call momx1(tu,qu,i,j,l)
! call momx3(tuHO,quHO,i,j,l)
call momx5(tuHO,quHO,i,j,l)
defect_u(i,j,l)=tuHO-tu
else
defect_u(i,j,l)=0._pr
endif
if (((bcx.eq.'p'.or.bcx.eq.'i').or.((bcx.eq.'w').and.(i-3.ge.imin).and.(i+3.le.imax))).and. &
((bcy.eq.'p').or.((bcy.eq.'w').and.(j-3.ge.jmin).and.(j+3.le.jmax-1))).and. &
((bcz.eq.'p').or.((bcz.eq.'w').and.(l-3.ge.lmin).and.(l+3.le.lmax)))) then
call momy1(tv,qv,i,j,l)
! call momy3(tvHO,qvHO,i,j,l)
call momy5(tvHO,qvHO,i,j,l)
defect_v(i,j,l)=tvHO-tv
else
defect_v(i,j,l)=0._pr
endif
if (((bcx.eq.'p'.or.bcx.eq.'i').or.((bcx.eq.'w').and.(i-3.ge.imin).and.(i+3.le.imax))).and. &
((bcy.eq.'p').or.((bcy.eq.'w').and.(j-3.ge.jmin).and.(j+3.le.jmax))).and. &
((bcz.eq.'p').or.((bcz.eq.'w').and.(l-3.ge.lmin).and.(l+3.le.lmax-1)))) then
call momz1(tw,qw,i,j,l)
! call momz3(twHO,qwHO,i,j,l)
call momz5(twHO,qwHO,i,j,l)
defect_w(i,j,l)=twHO-tw
else
defect_w(i,j,l)=0._pr
endif
1 enddo
return
end