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
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : nvr, nqt
|
|
|
|
use blocks , only : idn, imx, imy, imz, ien
|
2010-02-11 23:30:46 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
2010-02-28 18:35:57 -03:00
|
|
|
use config , only : im, jm, km
|
|
|
|
#ifdef FLUXCT
|
|
|
|
use interpolation, only : magtocen
|
|
|
|
#endif /* FLUXCT */
|
2008-12-08 19:07:42 -06:00
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
real, dimension(nqt,im,jm,km), intent(in) :: u
|
|
|
|
real, dimension(nqt,im,jm,km), intent(out) :: du
|
|
|
|
real , intent(in) :: dxi, dyi, dzi
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
2009-10-28 00:12:18 -02:00
|
|
|
integer :: i, j, k, im1, jm1, km1, ip1, jp1, kp1
|
2008-12-08 20:03:01 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
! local temporary arrays
|
2008-12-08 20:03:01 -06:00
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
real, dimension(nvr,im) :: ux
|
|
|
|
real, dimension(nqt,im) :: fx
|
|
|
|
real, dimension(nvr,jm) :: uy
|
|
|
|
real, dimension(nqt,jm) :: fy
|
|
|
|
#if NDIMS == 3
|
|
|
|
real, dimension(nvr,km) :: uz
|
|
|
|
real, dimension(nqt,km) :: fz
|
|
|
|
#endif /* NDIMS == 3 */
|
|
|
|
real, dimension(nvr,im,jm,km) :: w
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2010-02-28 18:35:57 -03:00
|
|
|
!
|
|
|
|
! copy variables to a new array
|
|
|
|
!
|
|
|
|
w(1:nqt,:,:,:) = u(1:nqt,:,:,:)
|
|
|
|
#ifdef FLUXCT
|
|
|
|
! interpolate cell centered components of the magnetic field
|
|
|
|
!
|
|
|
|
call magtocen(w)
|
|
|
|
#endif /* FLUXCT */
|
|
|
|
|
|
|
|
! reset the increment array
|
2008-12-08 19:07:42 -06:00
|
|
|
!
|
|
|
|
du(:,:,:,:) = 0.0
|
|
|
|
|
2008-12-08 20:03:01 -06:00
|
|
|
! update along X-direction
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do k = 1, km
|
2009-10-28 00:12:18 -02:00
|
|
|
#if NDIMS == 3
|
|
|
|
#ifdef MHD
|
|
|
|
km1 = max(k - 1, 1)
|
|
|
|
kp1 = min(k + 1,km)
|
|
|
|
#endif /* MHD */
|
|
|
|
#endif /* NDIMS == 3 */
|
2008-12-09 22:28:29 -06:00
|
|
|
do j = 1, jm
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
jm1 = max(j - 1, 1)
|
|
|
|
jp1 = min(j + 1,jm)
|
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do i = 1, im
|
2010-02-28 18:35:57 -03:00
|
|
|
ux(idn,i) = w(idn,i,j,k)
|
|
|
|
ux(imx,i) = w(imx,i,j,k)
|
|
|
|
ux(imy,i) = w(imy,i,j,k)
|
|
|
|
ux(imz,i) = w(imz,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
ux(ien,i) = w(ien,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
ux(ibx,i) = w(ibx,i,j,k)
|
|
|
|
ux(iby,i) = w(iby,i,j,k)
|
|
|
|
ux(ibz,i) = w(ibz,i,j,k)
|
|
|
|
#ifdef FLUXCT
|
|
|
|
ux(icx,i) = w(icx,i,j,k)
|
|
|
|
ux(icy,i) = w(icy,i,j,k)
|
|
|
|
ux(icz,i) = w(icz,i,j,k)
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2010-02-28 18:35:57 -03:00
|
|
|
call hll (im, ux, fx)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
2008-12-16 14:38:04 -06:00
|
|
|
#ifdef HLLC
|
2010-02-28 18:35:57 -03:00
|
|
|
call hllc(im, ux, fx)
|
2008-12-16 14:38:04 -06:00
|
|
|
#endif /* HLLC */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do i = 1, im
|
2010-02-28 18:35:57 -03:00
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dxi * fx(idn,i)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dxi * fx(imx,i)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dxi * fx(imy,i)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dxi * fx(imz,i)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dxi * fx(ien,i)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
! update magnetic variables
|
|
|
|
!
|
|
|
|
im1 = max(i - 1, 1)
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FIELDCD
|
2010-02-22 16:21:20 -03:00
|
|
|
ip1 = min(i + 1,im)
|
2009-10-28 00:12:18 -02:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
du(iby,i,j,k) = du(iby,i,j,k) - 0.5 * dxi * (fx(iby,ip1) - fx(iby,im1))
|
2009-10-28 00:12:18 -02:00
|
|
|
#if NDIMS == 3
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ibz,i,j,k) = du(ibz,i,j,k) - 0.5 * dxi * (fx(ibz,ip1) - fx(ibz,im1))
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* NDIMS == 3 */
|
2010-02-22 16:21:20 -03:00
|
|
|
#endif /* FIELDCD */
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
|
|
|
du(ibx,i ,jm1,k ) = du(ibx,i ,jm1,k ) + dyi * fx(iby,i)
|
|
|
|
du(ibx,i ,jp1,k ) = du(ibx,i ,jp1,k ) - dyi * fx(iby,i)
|
|
|
|
du(iby,i ,j ,k ) = du(iby,i ,j ,k ) - dxi * (fx(iby,i) - fx(iby,im1))
|
|
|
|
du(iby,i ,jm1,k ) = du(iby,i ,jm1,k ) - dxi * (fx(iby,i) - fx(iby,im1))
|
|
|
|
#if NDIMS == 3
|
|
|
|
du(ibx,i ,j ,km1) = du(ibx,i ,j ,km1) + dzi * fx(ibz,i)
|
|
|
|
du(ibx,i ,j ,kp1) = du(ibx,i ,j ,kp1) - dzi * fx(ibz,i)
|
|
|
|
du(ibz,i ,j ,k ) = du(ibz,i ,j ,k ) - dxi * (fx(ibz,i) - fx(ibz,im1))
|
|
|
|
du(ibz,i ,j ,km1) = du(ibz,i ,j ,km1) - dxi * (fx(ibz,i) - fx(ibz,im1))
|
|
|
|
#endif /* NDIMS == 3 */
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
! update along Y-direction
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do k = 1, km
|
2009-10-28 00:12:18 -02:00
|
|
|
#if NDIMS == 3
|
|
|
|
#ifdef MHD
|
|
|
|
km1 = max(k - 1, 1)
|
|
|
|
kp1 = min(k + 1,km)
|
|
|
|
#endif /* MHD */
|
|
|
|
#endif /* NDIMS == 3 */
|
2008-12-09 22:28:29 -06:00
|
|
|
do i = 1, im
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
im1 = max(i - 1, 1)
|
|
|
|
ip1 = min(i + 1,im)
|
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do j = 1, jm
|
2010-02-28 18:35:57 -03:00
|
|
|
uy(idn,j) = w(idn,i,j,k)
|
|
|
|
uy(imx,j) = w(imy,i,j,k)
|
|
|
|
uy(imy,j) = w(imz,i,j,k)
|
|
|
|
uy(imz,j) = w(imx,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
uy(ien,j) = w(ien,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
uy(ibx,j) = w(iby,i,j,k)
|
|
|
|
uy(iby,j) = w(ibz,i,j,k)
|
|
|
|
uy(ibz,j) = w(ibx,i,j,k)
|
|
|
|
#ifdef FLUXCT
|
|
|
|
uy(icx,j) = w(icy,i,j,k)
|
|
|
|
uy(icy,j) = w(icz,i,j,k)
|
|
|
|
uy(icz,j) = w(icx,i,j,k)
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2010-02-28 18:35:57 -03:00
|
|
|
call hll (jm, uy, fy)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
2008-12-16 14:38:04 -06:00
|
|
|
#ifdef HLLC
|
2010-02-28 18:35:57 -03:00
|
|
|
call hllc(jm, uy, fy)
|
2008-12-16 14:38:04 -06:00
|
|
|
#endif /* HLLC */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do j = 1, jm
|
2010-02-28 18:35:57 -03:00
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dyi * fy(idn,j)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dyi * fy(imz,j)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dyi * fy(imx,j)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dyi * fy(imy,j)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dyi * fy(ien,j)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
! update magnetic variables
|
|
|
|
!
|
|
|
|
jm1 = max(j - 1, 1)
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FIELDCD
|
2010-02-22 16:21:20 -03:00
|
|
|
jp1 = min(j + 1,jm)
|
|
|
|
|
2009-10-28 00:12:18 -02:00
|
|
|
#if NDIMS == 3
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ibz,i,j,k) = du(ibz,i,j,k) - 0.5 * dyi * (fy(iby,jp1) - fy(iby,jm1))
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* NDIMS == 3 */
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ibx,i,j,k) = du(ibx,i,j,k) - 0.5 * dyi * (fy(ibz,jp1) - fy(ibz,jm1))
|
2010-02-22 16:21:20 -03:00
|
|
|
#endif /* FIELDCD */
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
|
|
|
#if NDIMS == 3
|
|
|
|
du(iby,i ,j ,km1) = du(iby,i ,j ,km1) + dzi * fy(iby,j)
|
|
|
|
du(iby,i ,j ,kp1) = du(iby,i ,j ,kp1) - dzi * fy(iby,j)
|
|
|
|
du(ibz,i ,j ,k ) = du(ibz,i ,j ,k ) - dyi * (fy(iby,j) - fy(iby,jm1))
|
|
|
|
du(ibz,i ,j ,km1) = du(ibz,i ,j ,km1) - dyi * (fy(iby,j) - fy(iby,jm1))
|
|
|
|
#endif /* NDIMS == 3 */
|
|
|
|
du(iby,im1,j ,k ) = du(iby,im1,j ,k ) + dxi * fy(ibz,j)
|
|
|
|
du(iby,ip1,j ,k ) = du(iby,ip1,j ,k ) - dxi * fy(ibz,j)
|
|
|
|
du(ibx,i ,j ,k ) = du(ibx,i ,j ,k ) - dyi * (fy(ibz,j) - fy(ibz,jm1))
|
|
|
|
du(ibx,im1,j ,k ) = du(ibx,im1,j ,k ) - dyi * (fy(ibz,j) - fy(ibz,jm1))
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
#if NDIMS == 3
|
|
|
|
! update along Z-direction
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do j = 1, jm
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
jm1 = max(j - 1, 1)
|
|
|
|
jp1 = min(j + 1,jm)
|
|
|
|
#endif /* MHD */
|
2008-12-09 22:28:29 -06:00
|
|
|
do i = 1, im
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
im1 = max(i - 1, 1)
|
|
|
|
ip1 = min(i + 1,im)
|
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! copy directional vectors of variables for the one dimensional solver
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do k = 1, km
|
2010-02-28 18:35:57 -03:00
|
|
|
uz(idn,k) = w(idn,i,j,k)
|
|
|
|
uz(imx,k) = w(imz,i,j,k)
|
|
|
|
uz(imy,k) = w(imx,i,j,k)
|
|
|
|
uz(imz,k) = w(imy,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
uz(ien,k) = w(ien,i,j,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
uz(ibx,k) = w(ibz,i,j,k)
|
|
|
|
uz(iby,k) = w(ibx,i,j,k)
|
|
|
|
uz(ibz,k) = w(iby,i,j,k)
|
|
|
|
#ifdef FLUXCT
|
|
|
|
uz(icx,k) = w(icz,i,j,k)
|
|
|
|
uz(icy,k) = w(icx,i,j,k)
|
|
|
|
uz(icz,k) = w(icy,i,j,k)
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
end do
|
|
|
|
|
|
|
|
! execute solver (returns fluxes for the update)
|
|
|
|
!
|
|
|
|
#ifdef HLL
|
2010-02-28 18:35:57 -03:00
|
|
|
call hll (km, uz, fz)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* HLL */
|
2008-12-16 14:38:04 -06:00
|
|
|
#ifdef HLLC
|
2010-02-28 18:35:57 -03:00
|
|
|
call hllc(km, uz, fz)
|
2008-12-16 14:38:04 -06:00
|
|
|
#endif /* HLLC */
|
2008-12-08 20:03:01 -06:00
|
|
|
|
|
|
|
! update the arrays of increments
|
|
|
|
!
|
2008-12-09 22:28:29 -06:00
|
|
|
do k = 1, km
|
2010-02-28 18:35:57 -03:00
|
|
|
du(idn,i,j,k) = du(idn,i,j,k) + dzi * fz(idn,k)
|
|
|
|
du(imx,i,j,k) = du(imx,i,j,k) + dzi * fz(imy,k)
|
|
|
|
du(imy,i,j,k) = du(imy,i,j,k) + dzi * fz(imz,k)
|
|
|
|
du(imz,i,j,k) = du(imz,i,j,k) + dzi * fz(imx,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ien,i,j,k) = du(ien,i,j,k) + dzi * fz(ien,k)
|
2008-12-08 20:03:01 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
! update magnetic variables
|
|
|
|
!
|
|
|
|
km1 = max(k - 1, 1)
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FIELDCD
|
2010-02-22 16:21:20 -03:00
|
|
|
kp1 = min(k + 1,km)
|
2009-10-28 00:12:18 -02:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
du(ibx,i,j,k) = du(ibx,i,j,k) - 0.5 * dzi * (fz(iby,kp1) - fz(iby,km1))
|
|
|
|
du(iby,i,j,k) = du(iby,i,j,k) - 0.5 * dzi * (fz(ibz,kp1) - fz(ibz,km1))
|
2010-02-22 16:21:20 -03:00
|
|
|
#endif /* FIELDCD */
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
|
|
|
du(ibz,im1,j ,k ) = du(ibz,im1,j ,k ) + dxi * fz(iby,k)
|
|
|
|
du(ibz,ip1,j ,k ) = du(ibz,ip1,j ,k ) - dxi * fz(iby,k)
|
|
|
|
du(ibx,i ,j ,k ) = du(ibx,i ,j ,k ) - dzi * (fz(iby,k) - fz(iby,km1))
|
|
|
|
du(ibx,im1,j ,k ) = du(ibx,im1,j ,k ) - dzi * (fz(iby,k) - fz(iby,km1))
|
|
|
|
|
|
|
|
du(ibz,i ,jm1,k ) = du(ibz,i ,jm1,k ) + dyi * fz(ibz,k)
|
|
|
|
du(ibz,i ,jp1,k ) = du(ibz,i ,jp1,k ) - dyi * fz(ibz,k)
|
|
|
|
du(iby,i ,j ,k ) = du(iby,i ,j ,k ) - dzi * (fz(ibz,k) - fz(ibz,km1))
|
|
|
|
du(iby,i ,jm1,k ) = du(iby,i ,jm1,k ) - dzi * (fz(ibz,k) - fz(ibz,km1))
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:03:01 -06:00
|
|
|
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
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
subroutine hll(n, u, f)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : nvr, nfl, nqt
|
|
|
|
use blocks , only : idn, imx, imy, imz, ivx, ivy, ivz, ipr, ien
|
2010-02-11 23:30:46 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
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
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(nvr,n), intent(in) :: u
|
|
|
|
real, dimension(nqt,n), intent(out) :: f
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer :: p, i
|
|
|
|
real, dimension(nvr,n) :: q, ql, qr, ul, ur
|
|
|
|
real, dimension(nqt,n) :: fl, fr, fn
|
|
|
|
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
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
! calculate the primitive variables
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call cons2prim(n, u, q)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
! reconstruct the left and right states of the primitive variables
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
do p = 1, nfl
|
2009-10-28 00:12:18 -02:00
|
|
|
call reconstruct(n, q(p,:), ql(p,:), qr(p,:))
|
|
|
|
end do
|
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef MHD
|
|
|
|
#ifdef FLUXCT
|
|
|
|
! reconstruct the left and right states of the magnetic field components
|
|
|
|
!
|
|
|
|
ql(ibx,:) = q(ibx,:)
|
|
|
|
qr(ibx,:) = q(ibx,:)
|
|
|
|
call reconstruct(n, q(icy,:), ql(iby,:), qr(iby,:))
|
|
|
|
call reconstruct(n, q(icz,:), ql(ibz,:), qr(ibz,:))
|
|
|
|
ql(icx:icz,:) = ql(ibx:ibz,:)
|
|
|
|
qr(icx:icz,:) = qr(ibx:ibz,:)
|
|
|
|
#endif /* FLUXCT */
|
|
|
|
#endif /* MHD */
|
|
|
|
|
2008-12-08 20:16:37 -06:00
|
|
|
! calculate conservative variables at states
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call prim2cons(n, ql, ul)
|
|
|
|
call prim2cons(n, qr, ur)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! calculate fluxes and speeds
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call fluxspeed(n, ql, ul, fl, cl)
|
|
|
|
call fluxspeed(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)
|
|
|
|
!
|
2009-10-28 00:12:18 -02:00
|
|
|
al = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i))
|
|
|
|
ar = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i))
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! calculate HLL flux
|
|
|
|
!
|
|
|
|
if (al .ge. 0.0) then
|
2010-02-28 18:35:57 -03:00
|
|
|
fn(:,i) = fl(:,i)
|
2008-12-08 20:16:37 -06:00
|
|
|
else if (ar .le. 0.0) then
|
2010-02-28 18:35:57 -03:00
|
|
|
fn(:,i) = fr(:,i)
|
2008-12-08 20:16:37 -06:00
|
|
|
else
|
|
|
|
ap = ar * al
|
2009-10-28 00:12:18 -02:00
|
|
|
div = 1.0 / (ar - al)
|
2008-12-08 20:16:37 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
fn(1:nqt,i) = div * (ar * fl(1:nqt,i) - al * fr(1:nqt,i) &
|
|
|
|
+ ap * (ur(1:nqt,i) - ul(1:nqt,i)))
|
2009-10-28 00:12:18 -02:00
|
|
|
end if
|
|
|
|
end do
|
2008-12-08 20:16:37 -06:00
|
|
|
|
|
|
|
! calculate numerical flux
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1)
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2010-02-22 16:21:20 -03:00
|
|
|
#ifdef FIELDCD
|
|
|
|
call emf(n, q(ivx:ivz,:), q(ibx:ibz,:), f(ibx:ibz,:))
|
|
|
|
#endif /* FIELDCD */
|
|
|
|
#ifdef FLUXCT
|
2010-02-28 18:35:57 -03:00
|
|
|
f(ibx:ibz,1:n) = 0.25 * fn(ibx:ibz,1:n)
|
2010-02-22 16:21:20 -03:00
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:16:37 -06:00
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:16:37 -06:00
|
|
|
!
|
|
|
|
end subroutine hll
|
|
|
|
#endif /* HLL */
|
2008-12-16 14:38:04 -06:00
|
|
|
#ifdef HLLC
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
! hllc: subroutine to compute flux approximated by HLLC method (HYDRO only)
|
|
|
|
! ([1] Batten et al., 1997, JSC, 18, 6, 1553)
|
|
|
|
! ([2] Miyoshi & Kusano, 2005, JCP, 208, 315)
|
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
subroutine hllc(n, u, f)
|
2008-12-16 14:38:04 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : nvr, nfl, nqt
|
2008-12-16 14:38:04 -06:00
|
|
|
use interpolation, only : reconstruct
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(nvr,n), intent(in) :: u
|
|
|
|
real, dimension(nqt,n), intent(out) :: f
|
2008-12-16 14:38:04 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer :: p, i
|
|
|
|
real, dimension(nvr,n) :: ql, qr, qc, ul, ur
|
|
|
|
real, dimension(nfl,n) :: fl, fr, fx
|
|
|
|
real, dimension(n) :: cl, cr, cm
|
|
|
|
real :: sl, sr, sm, sml, smr, srmv, slmv, srmm, slmm &
|
|
|
|
, smvl, smvr, div, pt
|
|
|
|
real, dimension(nvr) :: q1l, q1r, u1l, u1r
|
2008-12-16 14:38:04 -06:00
|
|
|
!
|
|
|
|
!----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
f (:,:) = 0.0
|
|
|
|
fx(:,:) = 0.0
|
|
|
|
|
|
|
|
! calculate primitive variables
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call cons2prim(nvr, n, uc, qc)
|
2008-12-16 14:38:04 -06:00
|
|
|
|
|
|
|
! reconstruct left and right states of primitive variables
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
do p = 1, nfl
|
|
|
|
call reconstruct(n, qc(p,:), ql(p,:), qr(p,:))
|
2008-12-16 14:38:04 -06:00
|
|
|
enddo
|
|
|
|
|
|
|
|
! calculate conservative variables at states
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call prim2cons(nvr, n, ql, ul)
|
|
|
|
call prim2cons(nvr, n, qr, ur)
|
2008-12-16 14:38:04 -06:00
|
|
|
|
|
|
|
! calculate fluxes and speeds
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
call fluxspeed(n, ql, ul, fl, cl)
|
|
|
|
call fluxspeed(n, qr, ur, fr, cr)
|
2008-12-16 14:38:04 -06:00
|
|
|
|
|
|
|
! calculate fluxes
|
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
|
|
|
|
! calculate min and max and intermediate speeds: eq. (67)
|
|
|
|
!
|
|
|
|
sl = min(ql(2,i) - cl(i),qr(2,i) - cr(i))
|
|
|
|
sr = max(ql(2,i) + cl(i),qr(2,i) + cr(i))
|
|
|
|
|
|
|
|
! all speeds >= 0, left side flux
|
|
|
|
!
|
|
|
|
if (sl .ge. 0.0) then
|
|
|
|
|
|
|
|
fx(:,i) = fl(:,i)
|
|
|
|
|
|
|
|
! all speeds <= 0, right side flux
|
|
|
|
!
|
|
|
|
else if (sr .le. 0.0) then
|
|
|
|
|
|
|
|
fx(:,i) = fr(:,i)
|
|
|
|
|
|
|
|
! intermediate states
|
|
|
|
!
|
|
|
|
else ! sl < 0 & sr > 0
|
|
|
|
|
|
|
|
! useful differences
|
|
|
|
!
|
|
|
|
slmv = sl - ql(2,i)
|
|
|
|
srmv = sr - qr(2,i)
|
|
|
|
|
|
|
|
! speed of contact discontinuity (eq. 34 [1], 14 [2])
|
|
|
|
!
|
|
|
|
div = srmv*qr(1,i) - slmv*ql(1,i)
|
|
|
|
sml = (srmv*ur(2,i) - slmv*ul(2,i) - qr(5,i) + ql(5,i))/div
|
|
|
|
div = slmv*ql(1,i) - srmv*qr(1,i)
|
|
|
|
smr = (slmv*ul(2,i) - srmv*ur(2,i) - ql(5,i) + qr(5,i))/div
|
|
|
|
sm = 0.5d0 * (sml + smr)
|
|
|
|
|
|
|
|
if (sm .eq. 0.0d0) then
|
|
|
|
|
|
|
|
! calculate left intermediate state
|
|
|
|
!
|
|
|
|
pt = ql(5,i) - ul(2,i)*slmv
|
|
|
|
u1l(1) = ql(1,i)*slmv/sl
|
|
|
|
u1l(2) = 0.0d0
|
|
|
|
u1l(3) = u1l(1)*ql(3,i)
|
|
|
|
u1l(4) = u1l(1)*ql(4,i)
|
|
|
|
if (sl .eq. 0.0) then
|
|
|
|
u1l(5) = ul(5,i)
|
|
|
|
else
|
|
|
|
u1l(5) = (slmv*ul(5,i) - ql(5,i)*ql(2,i))/sl
|
|
|
|
endif
|
|
|
|
|
|
|
|
! calculate right intermediate state
|
|
|
|
!
|
|
|
|
pt = qr(5,i) - ur(2,i)*srmv
|
|
|
|
u1r(1) = qr(1,i)*srmv/sr
|
|
|
|
u1r(2) = 0.0d0
|
|
|
|
u1r(3) = u1r(1)*qr(3,i)
|
|
|
|
u1r(4) = u1r(1)*qr(4,i)
|
|
|
|
if (sr .eq. 0.0) then
|
|
|
|
u1r(5) = ur(5,i)
|
|
|
|
else
|
|
|
|
u1r(5) = (srmv*ur(5,i) - qr(5,i)*qr(2,i))/sr
|
|
|
|
endif
|
|
|
|
|
|
|
|
! calculate intermediate flux
|
|
|
|
!
|
|
|
|
fx(:,i) = 0.5*(fl(:,i) + sl*(u1l(:) - ul(:,i)) &
|
|
|
|
+ fr(:,i) + sr*(u1r(:) - ur(:,i)))
|
|
|
|
|
|
|
|
else
|
|
|
|
|
|
|
|
! useful differences
|
|
|
|
!
|
|
|
|
slmm = sl - sm
|
|
|
|
srmm = sr - sm
|
|
|
|
smvl = sm - ql(2,i)
|
|
|
|
smvr = sm - qr(2,i)
|
|
|
|
|
|
|
|
! intermediate discontinuities
|
|
|
|
!
|
|
|
|
if (sm .gt. 0.0) then
|
|
|
|
|
|
|
|
! pressure of intermediate states (eq. 36 [1], 16 [2])
|
|
|
|
!
|
|
|
|
pt = ql(5,i) + ql(1,i)*slmv*smvl
|
|
|
|
|
|
|
|
! calculate left intermediate state
|
|
|
|
!
|
|
|
|
u1l(1) = ql(1,i)*slmv/slmm ! eq. (35 [1], 17 [2])
|
|
|
|
u1l(2) = u1l(1)*sm ! eq. (37 [1])
|
|
|
|
u1l(3) = u1l(1)*ql(3,i) ! eq. (38 [1], 18 [2])
|
|
|
|
u1l(4) = u1l(1)*ql(4,i) ! eq. (39 [1], 19 [2])
|
|
|
|
if (slmm .eq. 0.0) then
|
|
|
|
u1l(5) = ul(5,i)
|
|
|
|
else
|
|
|
|
u1l(5) = (slmv*ul(5,i) - ql(5,i)*ql(2,i) + pt*sm)/slmm ! eq. (40 [1], 20 [2])
|
|
|
|
endif
|
|
|
|
|
|
|
|
! calculate left intermediate flux
|
|
|
|
!
|
|
|
|
fx(:,i) = fl(:,i) + sl*(u1l(:) - ul(:,i)) ! eq. (29 [1], 15 [2])
|
|
|
|
|
|
|
|
else if (sm .lt. 0.0) then
|
|
|
|
|
|
|
|
! pressure of intermediate states (eq. 36 [1], 16 [2])
|
|
|
|
!
|
|
|
|
pt = qr(5,i) + qr(1,i)*srmv*smvr
|
|
|
|
|
|
|
|
! calculate right intermediate state
|
|
|
|
!
|
|
|
|
u1r(1) = qr(1,i)*srmv/srmm ! eq. (35 [1], 17 [2])
|
|
|
|
u1r(2) = u1r(1)*sm ! eq. (37 [1])
|
|
|
|
u1r(3) = u1r(1)*qr(3,i) ! eq. (38 [1], 18 [2])
|
|
|
|
u1r(4) = u1r(1)*qr(4,i) ! eq. (39 [1], 19 [2])
|
|
|
|
if (srmm .eq. 0.0) then
|
|
|
|
u1r(5) = ur(5,i)
|
|
|
|
else
|
|
|
|
u1r(5) = (srmv*ur(5,i) - qr(5,i)*qr(2,i) + pt*sm)/srmm ! eq. (40 [1], 20 [2])
|
|
|
|
endif
|
|
|
|
|
|
|
|
! calculate right intermediate flux
|
|
|
|
!
|
|
|
|
fx(:,i) = fr(:,i) + sr*(u1r(:) - ur(:,i)) ! eq. (30 [1], 15 [2])
|
|
|
|
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
! calculate numerical flux
|
|
|
|
!
|
|
|
|
f(:,2:n) = - fx(:,2:n) + fx(:,1:n-1)
|
|
|
|
|
|
|
|
!-------------------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
end subroutine hllc
|
|
|
|
#endif /* HLLC */
|
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
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
subroutine fluxspeed(n, q, u, f, c)
|
2008-12-08 20:29:13 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : nvr, nqt
|
2010-02-11 23:30:46 -02:00
|
|
|
use blocks, only : idn, imx, imy, imz, ivx, ivy, ivz, ipr, ien
|
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:29:13 -06:00
|
|
|
use config, only : gamma, csnd, csnd2
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(nvr,n), intent(in) :: q, u
|
|
|
|
real, dimension(nqt,n), intent(out) :: f
|
|
|
|
real, dimension(n) , intent(out) :: c
|
2008-12-08 20:29:13 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
2009-10-28 00:12:18 -02:00
|
|
|
real :: bb, pm, vb, cs, cb, ca
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
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
|
|
|
|
!
|
2009-10-28 00:12:18 -02:00
|
|
|
f(idn,i) = u(imx,i)
|
2008-12-08 20:29:13 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
f(imx,i) = q(ivx,i) * u(imx,i) + q(ipr,i)
|
2008-12-08 20:29:13 -06:00
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
2009-10-28 00:12:18 -02:00
|
|
|
f(imx,i) = q(ivx,i) * u(imx,i) + q(idn,i) * csnd2
|
2008-12-08 20:29:13 -06:00
|
|
|
#endif /* ISO */
|
2009-10-28 00:12:18 -02:00
|
|
|
f(imy,i) = q(ivx,i) * u(imy,i)
|
|
|
|
f(imz,i) = q(ivx,i) * u(imz,i)
|
|
|
|
#ifdef ADI
|
|
|
|
f(ien,i) = q(ivx,i) * (u(ien,i) + q(ipr,i))
|
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef MHD
|
|
|
|
bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
|
|
|
pm = 0.5 * bb
|
|
|
|
vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
|
|
|
|
f(imx,i) = f(imx,i) - q(ibx,i) * q(ibx,i) + pm
|
|
|
|
f(imy,i) = f(imy,i) - q(ibx,i) * q(iby,i)
|
|
|
|
f(imz,i) = f(imz,i) - q(ibx,i) * q(ibz,i)
|
2008-12-08 20:29:13 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
f(ien,i) = f(ien,i) + q(ivx,i) * pm - q(ibx,i) * vb
|
2008-12-08 20:29:13 -06:00
|
|
|
#endif /* ADI */
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
2009-10-28 00:12:18 -02:00
|
|
|
f(ibx,i) = 0.0
|
|
|
|
f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
|
|
|
|
f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
|
2010-02-28 18:35:57 -03:00
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:29:13 -06:00
|
|
|
|
|
|
|
! compute speeds
|
|
|
|
!
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
#ifdef ADI
|
|
|
|
cs = gamma * q(ipr,i)
|
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
cs = csnd2 * q(idn,i)
|
|
|
|
#endif /* ISO */
|
|
|
|
cb = cs + bb
|
|
|
|
ca = q(ibx,i) * q(ibx,i)
|
|
|
|
c(i) = sqrt(0.5 * (cb + sqrt(max(0.0, cb * cb - 4.0 * cs * ca))) / q(idn,i))
|
|
|
|
#else /* MHD */
|
2008-12-08 20:29:13 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
c(i) = sqrt(gamma * q(ipr,i) / q(idn,i))
|
2008-12-08 20:29:13 -06:00
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
c(i) = csnd
|
|
|
|
#endif /* ISO */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
|
|
|
end do
|
2008-12-08 20:29:13 -06:00
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:29:13 -06:00
|
|
|
!
|
|
|
|
end subroutine fluxspeed
|
2010-02-22 16:21:20 -03:00
|
|
|
#ifdef MHD
|
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
! emf: subroutine computes magnetic fluxes (electromotive force)
|
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
subroutine emf(n, v, b, f)
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(3,n), intent(in) :: v, b
|
|
|
|
real, dimension(3,n), intent(out) :: f
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
|
|
|
!
|
|
|
|
!-------------------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
! sweep over all points
|
|
|
|
!
|
|
|
|
do i = 1, n
|
|
|
|
f(1,i) = 0.0
|
|
|
|
f(2,i) = v(1,i) * b(2,i) - b(1,i) * v(2,i)
|
|
|
|
f(3,i) = v(1,i) * b(3,i) - b(1,i) * v(3,i)
|
|
|
|
end do
|
|
|
|
!
|
|
|
|
!-------------------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
end subroutine emf
|
|
|
|
#endif /* MHD */
|
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
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
subroutine cons2prim(n, u, q)
|
2008-12-08 20:38:46 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : nvr
|
2010-02-11 23:30:46 -02:00
|
|
|
use blocks, only : idn, imx, imy, imz, ivx, ivy, ivz, ipr, ien
|
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
2009-10-28 00:12:18 -02:00
|
|
|
use config, only : gammam1
|
2008-12-08 20:38:46 -06:00
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(nvr,n), intent(in) :: u
|
|
|
|
real, dimension(nvr,n), intent(out) :: q
|
2008-12-08 20:38:46 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
2009-10-28 00:12:18 -02:00
|
|
|
real :: dni, ei, ek, em
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
do i = 1, n
|
2009-10-28 00:12:18 -02:00
|
|
|
dni = 1.0 / u(idn,i)
|
2008-12-08 20:38:46 -06:00
|
|
|
|
2009-10-28 00:12:18 -02:00
|
|
|
q(idn,i) = u(idn,i)
|
|
|
|
q(ivx,i) = dni * u(imx,i)
|
|
|
|
q(ivy,i) = dni * u(imy,i)
|
|
|
|
q(ivz,i) = dni * u(imz,i)
|
2008-12-08 20:38:46 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
ek = 0.5 * sum(u(imx:imz,i) * q(ivx:ivz,i))
|
|
|
|
ei = u(ien,i) - ek
|
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
em = 0.5 * sum(u(icx:icz,i) * u(icx:icz,i))
|
2009-10-28 00:12:18 -02:00
|
|
|
ei = ei - em
|
|
|
|
#endif /* MHD */
|
|
|
|
q(ipr,i) = gammam1 * ei
|
2008-12-08 20:38:46 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
q(ibx,i) = u(ibx,i)
|
|
|
|
q(iby,i) = u(iby,i)
|
|
|
|
q(ibz,i) = u(ibz,i)
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
|
|
|
q(icx,i) = u(icx,i)
|
|
|
|
q(icy,i) = u(icy,i)
|
|
|
|
q(icz,i) = u(icz,i)
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
|
|
|
end do
|
2008-12-08 20:38:46 -06:00
|
|
|
|
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
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
subroutine prim2cons(n, q, u)
|
2008-12-08 20:38:46 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : nvr
|
2010-02-11 23:30:46 -02:00
|
|
|
use blocks, only : idn, imx, imy, imz, ivx, ivy, ivz, ipr, ien
|
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks, only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
2008-12-08 20:38:46 -06:00
|
|
|
use config, only : gammam1i
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input/output arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
integer , intent(in) :: n
|
|
|
|
real, dimension(nvr,n), intent(in) :: q
|
|
|
|
real, dimension(nvr,n), intent(out) :: u
|
2008-12-08 20:38:46 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i
|
2009-10-28 00:12:18 -02:00
|
|
|
real :: ei, ek, em
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
do i = 1, n
|
2009-10-28 00:12:18 -02:00
|
|
|
u(idn,i) = q(idn,i)
|
|
|
|
u(imx,i) = q(idn,i) * q(ivx,i)
|
|
|
|
u(imy,i) = q(idn,i) * q(ivy,i)
|
|
|
|
u(imz,i) = q(idn,i) * q(ivz,i)
|
2008-12-08 20:38:46 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
ei = gammam1i * q(ipr,i)
|
|
|
|
ek = 0.5 * sum(u(imx:imz,i) * q(ivx:ivz,i))
|
|
|
|
u(ien,i) = ei + ek
|
2008-12-08 20:38:46 -06:00
|
|
|
#endif /* ADI */
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
#ifdef ADI
|
2010-02-28 18:35:57 -03:00
|
|
|
em = 0.5 * sum(q(icx:icz,i) * q(icx:icz,i))
|
2009-10-28 00:12:18 -02:00
|
|
|
u(ien,i) = u(ien,i) + em
|
|
|
|
#endif /* ADI */
|
|
|
|
u(ibx,i) = q(ibx,i)
|
|
|
|
u(iby,i) = q(iby,i)
|
|
|
|
u(ibz,i) = q(ibz,i)
|
2010-02-28 18:35:57 -03:00
|
|
|
#ifdef FLUXCT
|
|
|
|
u(icx,i) = q(icx,i)
|
|
|
|
u(icy,i) = q(icy,i)
|
|
|
|
u(icz,i) = q(icz,i)
|
|
|
|
#endif /* FLUXCT */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
|
|
|
end do
|
2008-12-08 20:38:46 -06:00
|
|
|
|
2008-12-08 21:04:20 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
2008-12-08 20:38:46 -06:00
|
|
|
!
|
|
|
|
end subroutine prim2cons
|
2008-12-09 14:51:33 -06:00
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
! maxspeed: function to calculate maximum speed in the system
|
|
|
|
!
|
|
|
|
!===============================================================================
|
|
|
|
!
|
|
|
|
function maxspeed(u)
|
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : nvr, nqt
|
|
|
|
use blocks , only : idn, ivx, ivz, ipr
|
2010-02-11 23:30:46 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
use blocks , only : ibx, iby, ibz, icx, icy, icz
|
2010-02-11 23:30:46 -02:00
|
|
|
#endif /* MHD */
|
2010-02-28 18:35:57 -03:00
|
|
|
use config , only : im, jm, km, ib, ie, jb, je, kb, ke, gamma
|
|
|
|
#if defined MHD && defined FLUXCT
|
|
|
|
use interpolation, only : magtocen
|
|
|
|
#endif /* MHD && FLUXCT */
|
2008-12-09 14:51:33 -06:00
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
! input arguments
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
real, dimension(nqt,im,jm,km), intent(in) :: u
|
2008-12-09 14:51:33 -06:00
|
|
|
|
|
|
|
! local variables
|
|
|
|
!
|
|
|
|
integer :: i, j, k
|
|
|
|
real :: vv, v, c
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
|
|
|
real :: bb
|
|
|
|
#endif /* MHD */
|
2008-12-09 14:51:33 -06:00
|
|
|
real :: maxspeed
|
|
|
|
|
|
|
|
! local arrays
|
|
|
|
!
|
2010-02-28 18:35:57 -03:00
|
|
|
real, dimension(nvr,im,jm,km) :: w
|
|
|
|
real, dimension(nvr,im) :: q
|
2008-12-09 14:51:33 -06:00
|
|
|
!
|
|
|
|
!----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
maxspeed = 0.0
|
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
! copy variables to a new array
|
|
|
|
!
|
|
|
|
w(1:nqt,1:im,1:jm,1:km) = u(1:nqt,1:im,1:jm,1:km)
|
|
|
|
#ifdef FLUXCT
|
|
|
|
! interpolate cell centered components of the magnetic field
|
|
|
|
!
|
|
|
|
call magtocen(w)
|
|
|
|
#endif /* FLUXCT */
|
|
|
|
|
2008-12-09 14:51:33 -06:00
|
|
|
! iterate over all points and calculate maximum speed
|
|
|
|
!
|
2008-12-13 21:05:51 -06:00
|
|
|
do k = kb, ke
|
|
|
|
do j = jb, je
|
2008-12-09 14:51:33 -06:00
|
|
|
|
2010-02-28 18:35:57 -03:00
|
|
|
call cons2prim(im, w(1:nvr,1:im,j,k), q(1:nvr,1:im))
|
2008-12-09 14:51:33 -06:00
|
|
|
|
2008-12-13 21:05:51 -06:00
|
|
|
do i = ib, ie
|
2008-12-09 14:51:33 -06:00
|
|
|
|
|
|
|
! calculate the velocity
|
|
|
|
!
|
2008-12-13 21:05:51 -06:00
|
|
|
vv = sum(q(ivx:ivz,i)**2)
|
2008-12-09 14:51:33 -06:00
|
|
|
v = sqrt(vv)
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2010-02-28 18:35:57 -03:00
|
|
|
bb = sum(q(icx:icz,i)**2)
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-09 14:51:33 -06:00
|
|
|
|
|
|
|
! calculate the maximum characteristic speed
|
|
|
|
!
|
2009-10-28 00:12:18 -02:00
|
|
|
#ifdef MHD
|
2008-12-09 14:51:33 -06:00
|
|
|
#ifdef ADI
|
2009-10-28 00:12:18 -02:00
|
|
|
c = sqrt((gamma * q(ipr,i) + bb) / q(idn,i))
|
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
c = sqrt(csnd2 + bb / q(idn,i))
|
|
|
|
#endif /* ISO */
|
|
|
|
#else /* MHD */
|
|
|
|
#ifdef ADI
|
|
|
|
c = sqrt(gamma * q(ipr,i) / q(idn,i))
|
2008-12-09 14:51:33 -06:00
|
|
|
#endif /* ADI */
|
|
|
|
#ifdef ISO
|
|
|
|
c = csnd
|
|
|
|
#endif /* ISO */
|
2009-10-28 00:12:18 -02:00
|
|
|
#endif /* MHD */
|
2008-12-09 14:51:33 -06:00
|
|
|
|
|
|
|
! calculate maximum of the speed
|
|
|
|
!
|
|
|
|
maxspeed = max(maxspeed, v + c)
|
2009-10-28 00:12:18 -02:00
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
2010-02-28 18:35:57 -03:00
|
|
|
!
|
2008-12-09 14:51:33 -06:00
|
|
|
!-------------------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
end function maxspeed
|
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
|