2008-12-08 19:07:42 -06:00
|
|
|
!!*****************************************************************************
|
|
|
|
!!
|
|
|
|
!! module: scheme - handling the actual solver of the set of equations
|
|
|
|
!!
|
|
|
|
!! Copyright (C) 2008 Grzegorz Kowal <kowal@astro.wisc.edu>
|
|
|
|
!!
|
|
|
|
!!*****************************************************************************
|
|
|
|
!!
|
|
|
|
!! This file is part of Godunov-AMR.
|
|
|
|
!!
|
|
|
|
!! Godunov-AMR is free software; you can redistribute it and/or modify
|
|
|
|
!! it under the terms of the GNU General Public License as published by
|
|
|
|
!! the Free Software Foundation; either version 3 of the License, or
|
|
|
|
!! (at your option) any later version.
|
|
|
|
!!
|
|
|
|
!! Godunov-AMR is distributed in the hope that it will be useful,
|
|
|
|
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
!! GNU General Public License for more details.
|
|
|
|
!!
|
|
|
|
!! You should have received a copy of the GNU General Public License
|
|
|
|
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
!!
|
|
|
|
!!*****************************************************************************
|
|
|
|
!!
|
|
|
|
!
|
|
|
|
module scheme
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
contains
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!==============================================================================
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
! update: subroutine sweeps over all directions and integrates the directional
|
|
|
|
! derivatives of the flux in order to get the increment of solution
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!==============================================================================
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
2008-12-08 20:03:01 -06:00
|
|
|
subroutine update(u, du, dxi, dyi, dzi)
|
2008-12-08 19:07:42 -06:00
|
|
|
|
2008-12-08 20:03:01 -06:00
|
|
|
use blocks, only : nvars, idn, imx, imy, imz, ien
|
|
|
|
use config, only : igrids, jgrids, kgrids, ngrids
|
2008-12-08 19:07:42 -06:00
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input arguments
|
|
|
|
!
|
|
|
|
real, dimension(nvars,igrids,jgrids,kgrids), intent(in) :: u
|
|
|
|
real, dimension(nvars,igrids,jgrids,kgrids), intent(out) :: du
|
2008-12-08 20:03:01 -06:00
|
|
|
real , intent(in) :: dxi, dyi, dzi
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i, j, k
|
|
|
|
|
|
|
|
! local arrays
|
|
|
|
!
|
|
|
|
real, dimension(nvars,ngrids) :: ul, fl
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
|
|
|
du(:,:,:,:) = 0.0
|
|
|
|
|
2008-12-08 20:03:01 -06:00
|
|
|
! update along X-direction
|
|
|
|
!
|
|
|
|
do k = 1, kgrids
|
|
|
|
do j = 1, jgrids
|
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
|
|
|
do i = 1, igrids
|
|
|
|
ul(1,i) = u(idn,i,j,k)
|
|
|
|
ul(2,i) = u(imx,i,j,k)
|
|
|
|
ul(3,i) = u(imy,i,j,k)
|
|
|
|
ul(4,i) = u(imz,i,j,k)
|
|
|
|
#ifdef ADI
|
|
|
|
ul(5,i) = u(ien,i,j,k)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2008-12-08 20:29:13 -06:00
|
|
|
call hll(nvars, igrids, ul, fl)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
|
|
|
do i = 1, igrids
|
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dxi * fl(1,i)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dxi * fl(2,i)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dxi * fl(3,i)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dxi * fl(4,i)
|
|
|
|
#ifdef ADI
|
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dxi * fl(5,i)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
! update along Y-direction
|
|
|
|
!
|
|
|
|
do k = 1, kgrids
|
|
|
|
do i = 1, igrids
|
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
|
|
|
do j = 1, jgrids
|
|
|
|
ul(1,j) = u(idn,i,j,k)
|
|
|
|
ul(2,i) = u(imy,i,j,k)
|
|
|
|
ul(3,j) = u(imz,i,j,k)
|
|
|
|
ul(4,j) = u(imx,i,j,k)
|
|
|
|
#ifdef ADI
|
|
|
|
ul(5,j) = u(ien,i,j,k)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2008-12-08 20:29:13 -06:00
|
|
|
call hll(nvars, jgrids, ul, fl)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
|
|
|
do j = 1, jgrids
|
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dyi * fl(1,j)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dyi * fl(4,j)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dyi * fl(2,j)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dyi * fl(3,j)
|
|
|
|
#ifdef ADI
|
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dyi * fl(5,j)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
#if NDIMS == 3
|
|
|
|
! update along Z-direction
|
|
|
|
!
|
|
|
|
do j = 1, jgrids
|
|
|
|
do i = 1, igrids
|
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
|
|
|
do k = 1, kgrids
|
|
|
|
ul(1,k) = u(idn,i,j,k)
|
|
|
|
ul(2,k) = u(imz,i,j,k)
|
|
|
|
ul(3,k) = u(imx,i,j,k)
|
|
|
|
ul(4,k) = u(imy,i,j,k)
|
|
|
|
#ifdef ADI
|
|
|
|
ul(5,k) = u(ien,i,j,k)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2008-12-08 20:29:13 -06:00
|
|
|
call hll(nvars, kgrids, ul, fl)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
|
|
|
do k = 1, kgrids
|
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dzi * fl(1,k)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dzi * fl(3,k)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dzi * fl(4,k)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dzi * fl(2,k)
|
|
|
|
#ifdef ADI
|
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dzi * fl(5,k)
|
|
|
|
#endif /* ADI */
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
#endif /* NDIMS == 3 */
|
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2008-12-08 19:07:42 -06:00
|
|
|
end subroutine update
|
2008-12-08 20:16:37 -06:00
|
|
|
#ifdef HLL
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!===============================================================================
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
! hll: subroutine computes the approximated flux using HLL method
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!===============================================================================
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
subroutine hll(m, n, uc, f)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
2008-12-08 20:53:29 -06:00
|
|
|
use interpolation, only : reconstruct
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
integer , intent(in) :: m, n
|
|
|
|
real, dimension(m,n), intent(in) :: uc
|
|
|
|
real, dimension(m,n), intent(out) :: f
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
integer :: p, i
|
|
|
|
real, dimension(m,n) :: ul, ur, ql, qr, qc
|
|
|
|
real, dimension(m,n) :: fl, fr, fx
|
|
|
|
real, dimension(n) :: cl, cr
|
|
|
|
real :: al, ar, ap, div
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
|
|
|
#ifdef CONREC
|
|
|
|
! reconstruct left and right states of conserved variables
|
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
do p = 1, m
|
2008-12-08 20:53:29 -06:00
|
|
|
call reconstruct(n,uc(p,:),ul(p,:),ur(p,:))
|
2008-12-08 20:16:37 -06:00
|
|
|
enddo
|
|
|
|
|
|
|
|
! calculate primitive variables
|
|
|
|
!
|
2008-12-08 20:38:46 -06:00
|
|
|
call cons2prim(m,n,ul,ql)
|
|
|
|
call cons2prim(m,n,ur,qr)
|
2008-12-08 20:16:37 -06:00
|
|
|
#else /* CONREC */
|
|
|
|
|
|
|
|
! calculate primitive variables
|
|
|
|
!
|
2008-12-08 20:38:46 -06:00
|
|
|
call cons2prim(m,n,uc,qc)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! reconstruct left and right states of primitive variables
|
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
do p = 1, m
|
2008-12-08 20:53:29 -06:00
|
|
|
call reconstruct(n,qc(p,:),ql(p,:),qr(p,:))
|
2008-12-08 20:16:37 -06:00
|
|
|
enddo
|
|
|
|
|
|
|
|
! calculate conservative variables at states
|
|
|
|
!
|
2008-12-08 20:38:46 -06:00
|
|
|
call prim2cons(m,n,ql,ul)
|
|
|
|
call prim2cons(m,n,qr,ur)
|
2008-12-08 20:16:37 -06:00
|
|
|
#endif /* CONREC */
|
|
|
|
|
|
|
|
! calculate fluxes and speeds
|
|
|
|
!
|
2008-12-08 20:29:13 -06:00
|
|
|
call fluxspeed(m,n,ql,ul,fl,cl)
|
|
|
|
call fluxspeed(m,n,qr,ur,fr,cr)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! iterate over all points
|
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
|
|
|
|
! calculate min and max and intermediate speeds: eq. (67)
|
|
|
|
!
|
|
|
|
al = min(ql(2,i) - cl(i),qr(2,i) - cr(i))
|
|
|
|
ar = max(ql(2,i) + cl(i),qr(2,i) + cr(i))
|
|
|
|
|
|
|
|
! calculate HLL flux
|
|
|
|
!
|
|
|
|
if (al .ge. 0.0) then
|
|
|
|
fx(:,i) = fl(:,i)
|
|
|
|
else if (ar .le. 0.0) then
|
|
|
|
fx(:,i) = fr(:,i)
|
|
|
|
else
|
|
|
|
ap = ar * al
|
|
|
|
div = 1.0/(ar - al)
|
|
|
|
|
|
|
|
fx(:,i) = div*(ar*fl(:,i) - al*fr(:,i) + ap*(ur(:,i) - ul(:,i)))
|
|
|
|
endif
|
|
|
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
! calculate numerical flux
|
|
|
|
!
|
|
|
|
f(:,2:n) = - fx(:,2:n) + fx(:,1:n-1)
|
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
|
|
|
end subroutine hll
|
|
|
|
#endif /* HLL */
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
! fluxspeed: subroutine computes fluxes and speeds for a given set of equations
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
subroutine fluxspeed(m, n, q, u, f, c)
|
|
|
|
|
|
|
|
use config, only : gamma, csnd, csnd2
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
|
|
|
integer , intent(in) :: m, n
|
|
|
|
real, dimension(m,n), intent(in) :: q, u
|
|
|
|
real, dimension(m,n), intent(out) :: f
|
|
|
|
real, dimension(n) , intent(out) :: c
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
|
|
|
real :: cs
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
|
|
|
! sweep over all points
|
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
|
|
|
|
! compute fluxes
|
|
|
|
!
|
|
|
|
f(1,i) = u(2,i)
|
|
|
|
#ifdef ADI
|
|
|
|
f(2,i) = q(2,i)*u(2,i) + q(5,i)
|
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
f(2,i) = q(2,i)*u(2,i) + q(1,i)*csnd2
|
|
|
|
#endif /* ISO */
|
|
|
|
f(3,i) = q(2,i)*u(3,i)
|
|
|
|
f(4,i) = q(2,i)*u(4,i)
|
|
|
|
#ifdef ADI
|
|
|
|
f(5,i) = q(2,i)*(u(5,i) + q(5,i))
|
|
|
|
#endif /* ADI */
|
|
|
|
|
|
|
|
! compute speeds
|
|
|
|
!
|
|
|
|
#ifdef ADI
|
|
|
|
c(i) = sqrt(gamma*q(5,i)/q(1,i))
|
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
c(i) = csnd
|
|
|
|
#endif /* ISO */
|
|
|
|
enddo
|
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
|
|
|
end subroutine fluxspeed
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
! cons2prim: subroutine converts primitive variables to conservative
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
subroutine cons2prim(m, n, u, q)
|
|
|
|
|
|
|
|
use config , only : gammam1
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
|
|
|
integer , intent(in) :: m, n
|
|
|
|
real, dimension(m,n), intent(in) :: u
|
|
|
|
real, dimension(m,n), intent(out) :: q
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
|
|
|
real :: dni
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
dni = 1.0 / u(1,i)
|
|
|
|
|
|
|
|
q(1,i) = u(1,i)
|
|
|
|
q(2,i) = dni*u(2,i)
|
|
|
|
q(3,i) = dni*u(3,i)
|
|
|
|
q(4,i) = dni*u(4,i)
|
|
|
|
#ifdef ADI
|
|
|
|
q(5,i) = u(5,i) &
|
|
|
|
- 0.5*(u(2,i)*q(2,i) + u(3,i)*q(3,i) + u(4,i)*q(4,i))
|
|
|
|
q(5,i) = gammam1*q(5,i)
|
|
|
|
#endif /* ADI */
|
|
|
|
enddo
|
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
end subroutine cons2prim
|
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
! prim2cons: subroutine converts primitive variables to conservative
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
subroutine prim2cons(m, n, q, u)
|
|
|
|
|
|
|
|
use config, only : gammam1i
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
|
|
|
integer , intent(in) :: m, n
|
|
|
|
real, dimension(m,n), intent(in) :: q
|
|
|
|
real, dimension(m,n), intent(out) :: u
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
u(1,i) = q(1,i)
|
|
|
|
u(2,i) = q(1,i)*q(2,i)
|
|
|
|
u(3,i) = q(1,i)*q(3,i)
|
|
|
|
u(4,i) = q(1,i)*q(4,i)
|
|
|
|
#ifdef ADI
|
|
|
|
u(5,i) = gammam1i*q(5,i) &
|
|
|
|
+ 0.5*(u(2,i)*q(2,i) + u(3,i)*q(3,i) + u(4,i)*q(4,i))
|
|
|
|
#endif /* ADI */
|
|
|
|
enddo
|
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
end subroutine prim2cons
|
2008-12-08 19:07:42 -06:00
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!===============================================================================
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
|
|
|
end module
|