From f1378c6e4ea9a09df11ee1a3cf12b1fce4aa4abe Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 5 Mar 2014 19:04:23 -0300 Subject: [PATCH 1/7] SCHEMES: Group subroutines according to the equations. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 1124 ++++++++++++++++++++++++----------------------- 1 file changed, 567 insertions(+), 557 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index 47904f3..7589e10 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -451,6 +451,8 @@ module schemes !! !=============================================================================== ! +!***** ISOTHERMAL HYDRODYNAMICS ***** +! !=============================================================================== ! ! subroutine UPDATE_FLUX_HD_ISO: @@ -467,7 +469,6 @@ module schemes ! q - the array of primitive variables; ! f - the array of numerical fluxes; ! -! !=============================================================================== ! subroutine update_flux_hd_iso(idir, dx, q, f) @@ -491,14 +492,14 @@ module schemes ! local variables ! - integer :: i, j, k + integer :: i, j, k ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, fx + real(kind=8), dimension(nv,jm) :: qy, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -616,557 +617,6 @@ module schemes ! !=============================================================================== ! -! subroutine UPDATE_FLUX_HD_ADI: -! ----------------------------- -! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. -! -! Arguments: -! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! -! -!=============================================================================== -! - subroutine update_flux_hd_adi(idir, dx, q, f) - -! include external variables -! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f - -! local variables -! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for flux update -! - call start_timer(imf) -#endif /* PROFILE */ - -! reset the flux array -! - f(:,:,:,:) = 0.0d+00 - -! select the directional flux to compute -! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ipr,1:im) = q(ipr,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ien,1:im,j,k) = fx(ien,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ipr,1:jm) = q(ipr,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ien,i,1:jm,k) = fy(ien,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ipr,1:km) = q(ipr,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ien,i,j,1:km) = fz(ien,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select - -#ifdef PROFILE -! stop accounting time for flux update -! - call stop_timer(imf) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_flux_hd_adi -! -!=============================================================================== -! -! subroutine UPDATE_FLUX_MHD_ADI: -! ------------------------------ -! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. -! -! Arguments: -! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! -! -!=============================================================================== -! - subroutine update_flux_mhd_iso(idir, dx, q, f) - -! include external variables -! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz - use equations , only : ibx, iby, ibz, ibp - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f - -! local variables -! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for flux update -! - call start_timer(imf) -#endif /* PROFILE */ - -! reset the flux array -! - f(:,:,:,:) = 0.0d+00 - -! select the directional flux to compute -! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ibx,1:im) = q(ibx,1:im,j,k) - qx(iby,1:im) = q(iby,1:im,j,k) - qx(ibz,1:im) = q(ibz,1:im,j,k) - qx(ibp,1:im) = q(ibp,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ibx,1:im,j,k) = fx(ibx,1:im) - f(iby,1:im,j,k) = fx(iby,1:im) - f(ibz,1:im,j,k) = fx(ibz,1:im) - f(ibp,1:im,j,k) = fx(ibp,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ibx,1:jm) = q(iby,i,1:jm,k) - qy(iby,1:jm) = q(ibz,i,1:jm,k) - qy(ibz,1:jm) = q(ibx,i,1:jm,k) - qy(ibp,1:jm) = q(ibp,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ibx,i,1:jm,k) = fy(ibz,1:jm) - f(iby,i,1:jm,k) = fy(ibx,1:jm) - f(ibz,i,1:jm,k) = fy(iby,1:jm) - f(ibp,i,1:jm,k) = fy(ibp,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ibx,1:km) = q(ibz,i,j,1:km) - qz(iby,1:km) = q(ibx,i,j,1:km) - qz(ibz,1:km) = q(iby,i,j,1:km) - qz(ibp,1:km) = q(ibp,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ibx,i,j,1:km) = fz(iby,1:km) - f(iby,i,j,1:km) = fz(ibz,1:km) - f(ibz,i,j,1:km) = fz(ibx,1:km) - f(ibp,i,j,1:km) = fz(ibp,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select - -#ifdef PROFILE -! stop accounting time for flux update -! - call stop_timer(imf) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_flux_mhd_iso -! -!=============================================================================== -! -! subroutine UPDATE_FLUX_MHD_ADI: -! ------------------------------ -! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. -! -! Arguments: -! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! -! -!=============================================================================== -! - subroutine update_flux_mhd_adi(idir, dx, q, f) - -! include external variables -! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien - use equations , only : ibx, iby, ibz, ibp - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f - -! local variables -! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for flux update -! - call start_timer(imf) -#endif /* PROFILE */ - -! reset the flux array -! - f(:,:,:,:) = 0.0d+00 - -! select the directional flux to compute -! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ibx,1:im) = q(ibx,1:im,j,k) - qx(iby,1:im) = q(iby,1:im,j,k) - qx(ibz,1:im) = q(ibz,1:im,j,k) - qx(ibp,1:im) = q(ibp,1:im,j,k) - qx(ipr,1:im) = q(ipr,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ibx,1:im,j,k) = fx(ibx,1:im) - f(iby,1:im,j,k) = fx(iby,1:im) - f(ibz,1:im,j,k) = fx(ibz,1:im) - f(ibp,1:im,j,k) = fx(ibp,1:im) - f(ien,1:im,j,k) = fx(ien,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ibx,1:jm) = q(iby,i,1:jm,k) - qy(iby,1:jm) = q(ibz,i,1:jm,k) - qy(ibz,1:jm) = q(ibx,i,1:jm,k) - qy(ibp,1:jm) = q(ibp,i,1:jm,k) - qy(ipr,1:jm) = q(ipr,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ibx,i,1:jm,k) = fy(ibz,1:jm) - f(iby,i,1:jm,k) = fy(ibx,1:jm) - f(ibz,i,1:jm,k) = fy(iby,1:jm) - f(ibp,i,1:jm,k) = fy(ibp,1:jm) - f(ien,i,1:jm,k) = fy(ien,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ibx,1:km) = q(ibz,i,j,1:km) - qz(iby,1:km) = q(ibx,i,j,1:km) - qz(ibz,1:km) = q(iby,i,j,1:km) - qz(ibp,1:km) = q(ibp,i,j,1:km) - qz(ipr,1:km) = q(ipr,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ibx,i,j,1:km) = fz(iby,1:km) - f(iby,i,j,1:km) = fz(ibz,1:km) - f(ibz,i,j,1:km) = fz(ibx,1:km) - f(ibp,i,j,1:km) = fz(ibp,1:km) - f(ien,i,j,1:km) = fz(ien,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select - -#ifdef PROFILE -! stop accounting time for flux update -! - call stop_timer(imf) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_flux_mhd_adi -! -!=============================================================================== -! ! subroutine RIEMANN_HD_ISO_HLL: ! ----------------------------- ! @@ -1299,6 +749,178 @@ module schemes ! !=============================================================================== ! +!***** ADIABATIC HYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_HD_ADI: +! ----------------------------- +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_hd_adi(idir, dx, q, f) + +! include external variables +! + use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, fx + real(kind=8), dimension(nv,jm) :: qy, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ipr,1:im) = q(ipr,1:im,j,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ien,1:im,j,k) = fx(ien,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ipr,1:jm) = q(ipr,i,1:jm,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ien,i,1:jm,k) = fy(ien,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ipr,1:km) = q(ipr,i,j,1:km) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ien,i,j,1:km) = fz(ien,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_hd_adi +! +!=============================================================================== +! ! subroutine RIEMANN_HD_ADI_HLL: ! ----------------------------- ! @@ -1623,8 +1245,199 @@ module schemes ! !=============================================================================== ! +!***** ISOTHERMAL MAGNETOHYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_MHD_ISO: +! ------------------------------ +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_mhd_iso(idir, dx, q, f) + +! include external variables +! + use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz + use equations , only : ibx, iby, ibz, ibp + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, fx + real(kind=8), dimension(nv,jm) :: qy, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ibx,1:im) = q(ibx,1:im,j,k) + qx(iby,1:im) = q(iby,1:im,j,k) + qx(ibz,1:im) = q(ibz,1:im,j,k) + qx(ibp,1:im) = q(ibp,1:im,j,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ibx,1:im,j,k) = fx(ibx,1:im) + f(iby,1:im,j,k) = fx(iby,1:im) + f(ibz,1:im,j,k) = fx(ibz,1:im) + f(ibp,1:im,j,k) = fx(ibp,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ibx,1:jm) = q(iby,i,1:jm,k) + qy(iby,1:jm) = q(ibz,i,1:jm,k) + qy(ibz,1:jm) = q(ibx,i,1:jm,k) + qy(ibp,1:jm) = q(ibp,i,1:jm,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ibx,i,1:jm,k) = fy(ibz,1:jm) + f(iby,i,1:jm,k) = fy(ibx,1:jm) + f(ibz,i,1:jm,k) = fy(iby,1:jm) + f(ibp,i,1:jm,k) = fy(ibp,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ibx,1:km) = q(ibz,i,j,1:km) + qz(iby,1:km) = q(ibx,i,j,1:km) + qz(ibz,1:km) = q(iby,i,j,1:km) + qz(ibp,1:km) = q(ibp,i,j,1:km) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ibx,i,j,1:km) = fz(iby,1:km) + f(iby,i,j,1:km) = fz(ibz,1:km) + f(ibz,i,j,1:km) = fz(ibx,1:km) + f(ibp,i,j,1:km) = fz(ibp,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_mhd_iso +! +!=============================================================================== +! ! subroutine RIEMANN_MHD_ISO_HLL: -! ----------------------------- +! ------------------------------ ! ! Subroutine solves one dimensional Riemann problem using ! the Harten-Lax-van Leer (HLL) method. @@ -2548,6 +2361,203 @@ module schemes ! !=============================================================================== ! +!***** ADIABATIC MAGNETOHYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_MHD_ADI: +! ------------------------------ +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_mhd_adi(idir, dx, q, f) + +! include external variables +! + use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + use equations , only : ibx, iby, ibz, ibp + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, fx + real(kind=8), dimension(nv,jm) :: qy, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ibx,1:im) = q(ibx,1:im,j,k) + qx(iby,1:im) = q(iby,1:im,j,k) + qx(ibz,1:im) = q(ibz,1:im,j,k) + qx(ibp,1:im) = q(ibp,1:im,j,k) + qx(ipr,1:im) = q(ipr,1:im,j,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ibx,1:im,j,k) = fx(ibx,1:im) + f(iby,1:im,j,k) = fx(iby,1:im) + f(ibz,1:im,j,k) = fx(ibz,1:im) + f(ibp,1:im,j,k) = fx(ibp,1:im) + f(ien,1:im,j,k) = fx(ien,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ibx,1:jm) = q(iby,i,1:jm,k) + qy(iby,1:jm) = q(ibz,i,1:jm,k) + qy(ibz,1:jm) = q(ibx,i,1:jm,k) + qy(ibp,1:jm) = q(ibp,i,1:jm,k) + qy(ipr,1:jm) = q(ipr,i,1:jm,k) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ibx,i,1:jm,k) = fy(ibz,1:jm) + f(iby,i,1:jm,k) = fy(ibx,1:jm) + f(ibz,i,1:jm,k) = fy(iby,1:jm) + f(ibp,i,1:jm,k) = fy(ibp,1:jm) + f(ien,i,1:jm,k) = fy(ien,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ibx,1:km) = q(ibz,i,j,1:km) + qz(iby,1:km) = q(ibx,i,j,1:km) + qz(ibz,1:km) = q(iby,i,j,1:km) + qz(ibp,1:km) = q(ibp,i,j,1:km) + qz(ipr,1:km) = q(ipr,i,j,1:km) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ibx,i,j,1:km) = fz(iby,1:km) + f(iby,i,j,1:km) = fz(ibz,1:km) + f(ibz,i,j,1:km) = fz(ibx,1:km) + f(ibp,i,j,1:km) = fz(ibp,1:km) + f(ien,i,j,1:km) = fz(ien,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_mhd_adi +! +!=============================================================================== +! ! subroutine RIEMANN_MHD_ADI_HLL: ! ------------------------------ ! From a49dad5ad041f306c35dec56665b1384a82cc9ec Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Mar 2014 12:15:44 -0300 Subject: [PATCH 2/7] SCHEMES: Separate Riemann state reconstruction from Riemann solvers. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 893 +++++++++++++++++++++++++++++------------------- 1 file changed, 535 insertions(+), 358 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index 7589e10..e7f600a 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -46,13 +46,17 @@ module schemes #ifdef PROFILE ! timer indices ! - integer , save :: imi, imu, imf, imr + integer , save :: imi, imu, imf, ims, imr #endif /* PROFILE */ ! pointer to the flux update procedure ! procedure(update_flux_hd_iso), pointer, save :: update_flux => null() +! pointer to the state reconstruction +! + procedure(states_hd_iso) , pointer, save :: states => null() + ! pointer to the Riemann solver ! procedure(riemann_hd_iso_hll), pointer, save :: riemann => null() @@ -93,8 +97,8 @@ module schemes ! include external procedures and variables ! - use equations , only : eqsys, eos - use parameters, only : get_parameter_string + use equations , only : eqsys, eos + use parameters , only : get_parameter_string ! local variables are not implicit by default ! @@ -118,6 +122,7 @@ module schemes call set_timer('schemes:: initialization' , imi) call set_timer('schemes:: increment update', imu) call set_timer('schemes:: flux update' , imf) + call set_timer('schemes:: Riemann states' , ims) call set_timer('schemes:: Riemann solver' , imr) ! start accounting time for module initialization/finalization @@ -146,6 +151,7 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_hd_iso + states => states_hd_iso ! select the Riemann solver ! @@ -171,6 +177,7 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_hd_adi + states => states_hd_adi ! select the Riemann solver ! @@ -215,6 +222,7 @@ module schemes ! set pointers to the subroutines ! update_flux => update_flux_mhd_iso + states => states_mhd_iso ! select the Riemann solver ! @@ -259,6 +267,7 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_mhd_adi + states => states_mhd_adi ! select the Riemann solver ! @@ -354,6 +363,7 @@ module schemes ! nullify procedure pointers ! nullify(update_flux) + nullify(states) nullify(riemann) #ifdef PROFILE @@ -386,8 +396,8 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv ! local variables are not implicit by default ! @@ -475,9 +485,9 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz ! local variables are not implicit by default ! @@ -496,10 +506,10 @@ module schemes ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -531,9 +541,13 @@ module schemes qx(ivy,1:im) = q(ivy,1:im,j,k) qx(ivz,1:im) = q(ivz,1:im,j,k) +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) ! update the array of fluxes ! @@ -559,9 +573,13 @@ module schemes qy(ivy,1:jm) = q(ivz,i,1:jm,k) qy(ivz,1:jm) = q(ivx,i,1:jm,k) +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) ! update the array of fluxes ! @@ -588,9 +606,13 @@ module schemes qz(ivy,1:km) = q(ivx,i,j,1:km) qz(ivz,1:km) = q(ivy,i,j,1:km) +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) ! update the array of fluxes ! @@ -617,36 +639,27 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_HD_ISO_HLL: -! ----------------------------- +! subroutine STATES_HD_ISO: +! ------------------------ ! -! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! Subroutine reconstructs the Riemann states. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; -! -! References: -! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; ! !=============================================================================== ! - subroutine riemann_hd_iso_hll(n, h, q, f) + subroutine states_hd_iso(n, h, q, ql, qr) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -657,16 +670,90 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr + +! local variables +! + integer :: p +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the state reconstruction +! + call start_timer(ims) +#endif /* PROFILE */ + +! reconstruct the left and right states of primitive variables +! + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do + +! check if the reconstruction gives negative values of density, +! if so, correct the states +! + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) + +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_hd_iso +! +!=============================================================================== +! +! subroutine RIEMANN_HD_ISO_HLL: +! ----------------------------- +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_hd_iso_hll(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -678,17 +765,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -773,9 +849,9 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien ! local variables are not implicit by default ! @@ -794,10 +870,10 @@ module schemes ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -830,9 +906,13 @@ module schemes qx(ivz,1:im) = q(ivz,1:im,j,k) qx(ipr,1:im) = q(ipr,1:im,j,k) +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) ! update the array of fluxes ! @@ -860,9 +940,13 @@ module schemes qy(ivz,1:jm) = q(ivx,i,1:jm,k) qy(ipr,1:jm) = q(ipr,i,1:jm,k) +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) ! update the array of fluxes ! @@ -891,9 +975,13 @@ module schemes qz(ivz,1:km) = q(ivy,i,j,1:km) qz(ipr,1:km) = q(ipr,i,j,1:km) +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) ! update the array of fluxes ! @@ -921,36 +1009,27 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_HD_ADI_HLL: -! ----------------------------- +! subroutine STATES_HD_ADI: +! ------------------------ ! -! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! Subroutine reconstructs the Riemann states. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; -! -! References: -! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; ! !=============================================================================== ! - subroutine riemann_hd_adi_hll(n, h, q, f) + subroutine states_hd_adi(n, h, q, ql, qr) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ipr, ivx - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ipr + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -961,16 +1040,91 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr + +! local variables +! + integer :: p +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the state reconstruction +! + call start_timer(ims) +#endif /* PROFILE */ + +! reconstruct the left and right states of primitive variables +! + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do + +! check if the reconstruction gives negative values of density or pressure, +! if so, correct the states +! + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) + call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) + +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_hd_adi +! +!=============================================================================== +! +! subroutine RIEMANN_HD_ADI_HLL: +! ----------------------------- +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_hd_adi_hll(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -982,18 +1136,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1063,10 +1205,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1076,14 +1217,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_hd_adi_hllc(n, h, q, f) + subroutine riemann_hd_adi_hllc(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -1092,20 +1232,19 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm real(kind=8) :: srml, slmm, srmm real(kind=8) :: dn, pr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -1117,18 +1256,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1269,10 +1396,10 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz - use equations , only : ibx, iby, ibz, ibp + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz + use equations , only : ibx, iby, ibz, ibp ! local variables are not implicit by default ! @@ -1291,10 +1418,10 @@ module schemes ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -1330,9 +1457,13 @@ module schemes qx(ibz,1:im) = q(ibz,1:im,j,k) qx(ibp,1:im) = q(ibp,1:im,j,k) +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) ! update the array of fluxes ! @@ -1366,9 +1497,13 @@ module schemes qy(ibz,1:jm) = q(ibx,i,1:jm,k) qy(ibp,1:jm) = q(ibp,i,1:jm,k) +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) ! update the array of fluxes ! @@ -1403,9 +1538,13 @@ module schemes qz(ibz,1:km) = q(iby,i,j,1:km) qz(ibp,1:km) = q(ibp,i,j,1:km) +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) ! update the array of fluxes ! @@ -1436,37 +1575,28 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_MHD_ISO_HLL: -! ------------------------------ +! subroutine STATES_MHD_ISO: +! ------------------------- ! -! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! Subroutine reconstructs the Riemann states. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; -! -! References: -! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; ! !=============================================================================== ! - subroutine riemann_mhd_iso_hll(n, h, q, f) + subroutine states_mhd_iso(n, h, q, ql, qr) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ibx, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ibx, ibp + use equations , only : cmax + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -1477,25 +1607,19 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q - real(kind=8), dimension(nv,n), intent(out) :: f + real(kind=8), dimension(nv,n), intent(out) :: ql, qr ! local variables ! - integer :: p, i - real(kind=8) :: sl, sr, srml - -! local arrays to store the states -! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr - real(kind=8), dimension(nv) :: wl, wr - real(kind=8), dimension(n) :: cl, cr + integer :: i, p + real(kind=8) :: bx, bp ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for Riemann solver +! start accounting time for the state reconstruction ! - call start_timer(imr) + call start_timer(ims) #endif /* PROFILE */ ! reconstruct the left and right states of primitive variables @@ -1506,18 +1630,95 @@ module schemes ! obtain the state values for Bx and Psi for the GLM-MHD equations ! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) + do i = 1, n -! check if the reconstruction doesn't give the negative density or pressure, + bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) & + - (qr(ibp,i) - ql(ibp,i)) / cmax) + bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) & + - (qr(ibx,i) - ql(ibx,i)) * cmax) + + ql(ibx,i) = bx + qr(ibx,i) = bx + ql(ibp,i) = bp + qr(ibp,i) = bp + + end do ! i = 1, n + +! check if the reconstruction gives negative values of density, ! if so, correct the states ! call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_mhd_iso +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ISO_HLL: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_mhd_iso_hll(n, ql, qr, f) + +! include external variables and procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: i + real(kind=8) :: sl, sr, srml + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(nv) :: wl, wr + real(kind=8), dimension(n) :: cl, cr +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1586,10 +1787,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1600,15 +1800,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_iso_hlld(n, h, q, f) + subroutine riemann_mhd_iso_hlld(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -1617,19 +1815,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -1641,26 +1838,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1979,10 +2156,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1992,15 +2168,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_iso_hlldm(n, h, q, f) + subroutine riemann_mhd_iso_hlldm(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -2009,19 +2183,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr, ca ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -2033,26 +2206,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -2385,10 +2538,10 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien - use equations , only : ibx, iby, ibz, ibp + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + use equations , only : ibx, iby, ibz, ibp ! local variables are not implicit by default ! @@ -2407,10 +2560,10 @@ module schemes ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -2447,9 +2600,13 @@ module schemes qx(ibp,1:im) = q(ibp,1:im,j,k) qx(ipr,1:im) = q(ipr,1:im,j,k) +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) ! update the array of fluxes ! @@ -2485,9 +2642,13 @@ module schemes qy(ibp,1:jm) = q(ibp,i,1:jm,k) qy(ipr,1:jm) = q(ipr,i,1:jm,k) +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) ! update the array of fluxes ! @@ -2524,9 +2685,13 @@ module schemes qz(ibp,1:km) = q(ibp,i,j,1:km) qz(ipr,1:km) = q(ipr,i,j,1:km) +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) ! update the array of fluxes ! @@ -2558,37 +2723,28 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_MHD_ADI_HLL: -! ------------------------------ +! subroutine STATES_MHD_ADI: +! ------------------------- ! -! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! Subroutine reconstructs the Riemann states. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; -! -! References: -! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; ! !=============================================================================== ! - subroutine riemann_mhd_adi_hll(n, h, q, f) + subroutine states_mhd_adi(n, h, q, ql, qr) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ipr, ivx, ibx, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ipr, ibx, ibp + use equations , only : cmax + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -2599,16 +2755,108 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr + +! local variables +! + integer :: i, p + real(kind=8) :: bx, bp +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the state reconstruction +! + call start_timer(ims) +#endif /* PROFILE */ + +! reconstruct the left and right states of primitive variables +! + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do ! p = 1, nv + +! obtain the state values for Bx and Psi for the GLM-MHD equations +! + do i = 1, n + + bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) & + - (qr(ibp,i) - ql(ibp,i)) / cmax) + bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) & + - (qr(ibx,i) - ql(ibx,i)) * cmax) + + ql(ibx,i) = bx + qr(ibx,i) = bx + ql(ibp,i) = bp + qr(ibp,i) = bp + + end do ! i = 1, n + +! check if the reconstruction gives negative values of density or density, +! if so, correct the states +! + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) + call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) + +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_mhd_adi +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ADI_HLL: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_mhd_adi_hll(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -2620,27 +2868,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -2710,10 +2937,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -2729,16 +2955,14 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_adi_hllc(n, h, q, f) + subroutine riemann_mhd_adi_hllc(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr - use equations , only : imx, imy, imz, ien - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -2747,19 +2971,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -2771,27 +2994,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -3014,10 +3216,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -3028,16 +3229,14 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_adi_hlld(n, h, q, f) + subroutine riemann_mhd_adi_hlld(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr - use equations , only : imx, imy, imz, ien - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -3046,20 +3245,19 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb, dv real(kind=8) :: dnl, dnr, cal, car, sml, smr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -3071,27 +3269,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) From a8d575e471d137b1ad0530c9a6698aa7147d20dd Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Mar 2014 12:59:51 -0300 Subject: [PATCH 3/7] EQUATIONS, SCHEMES: Implement Roe's solver for isothermal HD. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 148 ++++++++++++++++++++++++++++++++++++++++--- src/schemes.F90 | 156 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 295 insertions(+), 9 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 108b75b..0cb82be 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -63,16 +63,20 @@ module equations ! pointers to the conversion procedures ! - procedure(prim2cons_hd_iso), pointer, save :: prim2cons => null() - procedure(cons2prim_hd_iso), pointer, save :: cons2prim => null() + procedure(prim2cons_hd_iso) , pointer, save :: prim2cons => null() + procedure(cons2prim_hd_iso) , pointer, save :: cons2prim => null() ! pointer to the flux procedure ! - procedure(fluxspeed_hd_iso), pointer, save :: fluxspeed => null() + procedure(fluxspeed_hd_iso) , pointer, save :: fluxspeed => null() ! pointer to the maxspeed procedure ! - procedure(maxspeed_hd_iso) , pointer, save :: maxspeed => null() + procedure(maxspeed_hd_iso) , pointer, save :: maxspeed => null() + +! pointer to the Roe eigensystem procedure +! + procedure(esystem_roe_hd_iso), pointer, save :: eigensystem_roe => null() ! the system of equations and the equation of state @@ -97,6 +101,10 @@ module equations ! character(len=4), dimension(:), allocatable, save :: pvars, cvars +! eigenvectors +! + real(kind=8), dimension(:,:,:), allocatable, save :: evroe + ! adiabatic heat ratio ! real(kind=8) , save :: gamma = 5.0d+00 / 3.0d+00 @@ -123,6 +131,7 @@ module equations public :: prim2cons, cons2prim public :: fluxspeed public :: maxspeed, reset_maxspeed, get_maxspeed + public :: eigensystem_roe public :: update_primitive_variables public :: gamma public :: csnd @@ -237,10 +246,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_iso - cons2prim => cons2prim_hd_iso - fluxspeed => fluxspeed_hd_iso - maxspeed => maxspeed_hd_iso + prim2cons => prim2cons_hd_iso + cons2prim => cons2prim_hd_iso + fluxspeed => fluxspeed_hd_iso + maxspeed => maxspeed_hd_iso + eigensystem_roe => esystem_roe_hd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -435,6 +445,10 @@ module equations ! csnd2 = csnd * csnd +! allocate space for Roe eigenvectors +! + allocate(evroe(2,nv,nv)) + ! print information about the equation module ! if (verbose) then @@ -490,12 +504,17 @@ module equations if (allocated(pvars)) deallocate(pvars) if (allocated(cvars)) deallocate(cvars) +! deallocate Roe eigenvectors +! + if (allocated(evroe)) deallocate(evroe) + ! release the procedure pointers ! nullify(prim2cons) nullify(cons2prim) nullify(fluxspeed) - nullify(maxspeed ) + nullify(maxspeed) + nullify(eigensystem_roe) #ifdef PROFILE ! stop accounting time for module initialization/finalization @@ -911,6 +930,117 @@ module equations ! end function maxspeed_hd_iso ! +!=============================================================================== +! +! subroutine ESYSTEM_ROE_HD_ISO: +! ----------------------------- +! +! Subroutine computes eigenvalues and eigenvectors for a given set of +! equations and input variables. +! +! Arguments: +! +! q - the intermediate Roe state vector; +! c - the vector of eigenvalues; +! r - the matrix of right eigenvectors; +! l - the matrix of left eigenvectors; +! +! References: +! +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! +!=============================================================================== +! + subroutine esystem_roe_hd_iso(q, c, r, l) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), dimension(nv) , intent(in) :: q + real(kind=8), dimension(nv) , intent(inout) :: c + real(kind=8), dimension(nv,nv), intent(inout) :: l, r + +! local variables +! + logical , save :: first = .true. + real(kind=8), save :: ch +! +!------------------------------------------------------------------------------- +! +! prepare the internal arrays at the first run +! + if (first) then + +! prepare constants +! + ch = 0.5d+00 / csnd + +! reset all elements +! + evroe(:, : ,:) = 0.0d+00 + +! initiate the matrix of left eigenvectors +! + evroe(1,ivx,1) = - ch + evroe(1,ivy,2) = 1.0d+00 + evroe(1,ivz,3) = 1.0d+00 + evroe(1,ivx,4) = ch + +! initiate the matrix of right eigenvectors +! + evroe(2,1,idn) = 1.0d+00 + evroe(2,2,ivy) = 1.0d+00 + evroe(2,3,ivz) = 1.0d+00 + evroe(2,4,idn) = 1.0d+00 + +! unset the first execution flag +! + first = .false. + + end if ! first execution + +! prepare eigenvalues +! + c(1) = q(ivx) - csnd + c(2) = q(ivx) + c(3) = q(ivx) + c(4) = q(ivx) + csnd + +! update the varying elements of the matrix of left eigenvectors +! + evroe(1,idn,1) = ch * c(4) + evroe(1,idn,2) = - q(ivy) + evroe(1,idn,3) = - q(ivz) + evroe(1,idn,4) = - ch * c(1) + +! update the varying elements of the matrix of right eigenvectors +! + evroe(2,1,ivx) = c(1) + evroe(2,1,ivy) = q(ivy) + evroe(2,1,ivz) = q(ivz) + + evroe(2,4,ivx) = c(4) + evroe(2,4,ivy) = q(ivy) + evroe(2,4,ivz) = q(ivz) + +! copy matrices of eigenvectors +! + l(1:nv,1:nv) = evroe(1,1:nv,1:nv) + r(1:nv,1:nv) = evroe(2,1:nv,1:nv) + +!------------------------------------------------------------------------------- +! + end subroutine esystem_roe_hd_iso +! !******************************************************************************* ! ! ADIABATIC HYDRODYNAMIC EQUATIONS diff --git a/src/schemes.F90 b/src/schemes.F90 index e7f600a..0e82819 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -157,6 +157,15 @@ module schemes ! select case(trim(solver)) + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_hd_iso_roe ! in the case of unknown Riemann solver, revert to HLL ! @@ -825,6 +834,153 @@ module schemes ! !=============================================================================== ! +! subroutine RIEMANN_HD_ISO_ROE: +! ----------------------------- +! +! Subroutine solves one dimensional Riemann problem using +! the Roe's method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 +! +!=============================================================================== +! + subroutine riemann_hd_iso_roe(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + else + do p = nv, 1, - 1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_hd_iso_roe +! +!=============================================================================== +! !***** ADIABATIC HYDRODYNAMICS ***** ! !=============================================================================== From 78b15d5596d0256c718b52d34258a43c5bfa576a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Mar 2014 13:16:57 -0300 Subject: [PATCH 4/7] EQUATIONS, SCHEMES: Implement Roe's solver for adiabatic HD. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 158 +++++++++++++++++++++++++++++++++++++++++++-- src/schemes.F90 | 159 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 313 insertions(+), 4 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 0cb82be..7927385 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -269,10 +269,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_adi - cons2prim => cons2prim_hd_adi - fluxspeed => fluxspeed_hd_adi - maxspeed => maxspeed_hd_adi + prim2cons => prim2cons_hd_adi + cons2prim => cons2prim_hd_adi + fluxspeed => fluxspeed_hd_adi + maxspeed => maxspeed_hd_adi + eigensystem_roe => esystem_roe_hd_adi ! warn about the unimplemented equation of state ! @@ -1335,6 +1336,155 @@ module equations ! end function maxspeed_hd_adi ! +!=============================================================================== +! +! subroutine ESYSTEM_ROE_HD_ADI: +! ----------------------------- +! +! Subroutine computes eigenvalues and eigenvectors for a given set of +! equations and input variables. +! +! Arguments: +! +! q - the intermediate Roe state vector; +! c - the vector of eigenvalues; +! r - the matrix of right eigenvectors; +! l - the matrix of left eigenvectors; +! +! References: +! +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! +!=============================================================================== +! + subroutine esystem_roe_hd_adi(q, c, r, l) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), dimension(nv) , intent(in) :: q + real(kind=8), dimension(nv) , intent(inout) :: c + real(kind=8), dimension(nv,nv), intent(inout) :: l, r + +! local variables +! + logical, save :: first = .true. + real(kind=8) :: vv, vh, c2, na, cc, vc, ng, nd, nw, nh, nc +! +!------------------------------------------------------------------------------- +! +! prepare the internal arrays at the first run +! + if (first) then + +! reset all elements +! + evroe(:, : ,:) = 0.0d+00 + +! initiate the matrix of left eigenvectors +! + evroe(1,ivy,2) = 1.0d+00 + evroe(1,ivz,3) = 1.0d+00 + +! initiate the matrix of right eigenvectors +! + evroe(2,1,idn) = 1.0d+00 + evroe(2,2,ivy) = 1.0d+00 + evroe(2,3,ivz) = 1.0d+00 + evroe(2,4,idn) = 1.0d+00 + evroe(2,5,idn) = 1.0d+00 + +! unset the first execution flag +! + first = .false. + + end if ! first execution + +! calculate characteristic speeds and useful variables +! + vv = sum(q(ivx:ivz)**2) + vh = 0.5d+00 * vv + c2 = gammam1 * (q(ien) - vh) + na = 0.5d+00 / c2 + cc = sqrt(c2) + vc = q(ivx) * cc + ng = na * gammam1 + nd = 2.0d+00 * ng + nw = na * vc + nh = na * gammam1 * vh + nc = na * cc + +! prepare eigenvalues +! + c(1) = q(ivx) - cc + c(2) = q(ivx) + c(3) = q(ivx) + c(4) = q(ivx) + c(5) = q(ivx) + cc + +! update the varying elements of the matrix of left eigenvectors +! + evroe(1,idn,1) = nh + nw + evroe(1,ivx,1) = - ng * q(ivx) - nc + evroe(1,ivy,1) = - ng * q(ivy) + evroe(1,ivz,1) = - ng * q(ivz) + evroe(1,ien,1) = ng + + evroe(1,idn,2) = - q(ivy) + + evroe(1,idn,3) = - q(ivz) + + evroe(1,idn,4) = 1.0d+00 - ng * vv + evroe(1,ivx,4) = nd * q(ivx) + evroe(1,ivy,4) = nd * q(ivy) + evroe(1,ivz,4) = nd * q(ivz) + evroe(1,ien,4) = - nd + + evroe(1,idn,5) = nh - nw + evroe(1,ivx,5) = - ng * q(ivx) + nc + evroe(1,ivy,5) = - ng * q(ivy) + evroe(1,ivz,5) = - ng * q(ivz) + evroe(1,ien,5) = ng + +! update the varying elements of the matrix of right eigenvectors +! + evroe(2,1,ivx) = q(ivx) - cc + evroe(2,1,ivy) = q(ivy) + evroe(2,1,ivz) = q(ivz) + evroe(2,1,ien) = q(ien) - vc + + evroe(2,2,ien) = q(ivy) + + evroe(2,3,ien) = q(ivz) + + evroe(2,4,ivx) = q(ivx) + evroe(2,4,ivy) = q(ivy) + evroe(2,4,ivz) = q(ivz) + evroe(2,4,ien) = vh + + evroe(2,5,ivx) = q(ivx) + cc + evroe(2,5,ivy) = q(ivy) + evroe(2,5,ivz) = q(ivz) + evroe(2,5,ien) = q(ien) + vc + +! copy matrices of eigenvectors +! + l(1:nv,1:nv) = evroe(1,1:nv,1:nv) + r(1:nv,1:nv) = evroe(2,1:nv,1:nv) + +!------------------------------------------------------------------------------- +! + end subroutine esystem_roe_hd_adi +! !******************************************************************************* ! ! ISOTHERMAL MAGNETOHYDRODYNAMIC EQUATIONS diff --git a/src/schemes.F90 b/src/schemes.F90 index 0e82819..5a81fa2 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -202,6 +202,16 @@ module schemes ! riemann => riemann_hd_adi_hllc + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_hd_adi_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -1528,6 +1538,155 @@ module schemes ! !=============================================================================== ! +! subroutine RIEMANN_HD_ADI_ROE: +! ----------------------------- +! +! Subroutine solves one dimensional Riemann problem using +! the Roe's method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 +! +!=============================================================================== +! + subroutine riemann_hd_adi_roe(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ien + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ipr) = ((ul(ien,i) + ql(ipr,i)) / sdl & + + (ur(ien,i) + qr(ipr,i)) / sdr) / sds + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + else + do p = nv, 1, - 1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_hd_adi_roe +! +!=============================================================================== +! !***** ISOTHERMAL MAGNETOHYDRODYNAMICS ***** ! !=============================================================================== From 42c738c78d8b88bd4e74525886a456c954143545 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Mar 2014 13:55:09 -0300 Subject: [PATCH 5/7] EQUATIONS, SCHEMES: Implement Roe's solver for isothermal MHD. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 289 +++++++++++++++++++++++++++++++++++++++++++++- src/schemes.F90 | 193 ++++++++++++++++++++++++++++++- 2 files changed, 474 insertions(+), 8 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 7927385..e553f31 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -348,10 +348,11 @@ module equations ! set pointers to the subroutines ! - prim2cons => prim2cons_mhd_iso - cons2prim => cons2prim_mhd_iso - fluxspeed => fluxspeed_mhd_iso - maxspeed => maxspeed_mhd_iso + prim2cons => prim2cons_mhd_iso + cons2prim => cons2prim_mhd_iso + fluxspeed => fluxspeed_mhd_iso + maxspeed => maxspeed_mhd_iso + eigensystem_roe => esystem_roe_mhd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -941,6 +942,8 @@ module equations ! ! Arguments: ! +! x - ratio of the perpendicular magnetic field component difference +! y - ratio of the density ! q - the intermediate Roe state vector; ! c - the vector of eigenvalues; ! r - the matrix of right eigenvectors; @@ -958,7 +961,7 @@ module equations ! !=============================================================================== ! - subroutine esystem_roe_hd_iso(q, c, r, l) + subroutine esystem_roe_hd_iso(x, y, q, c, r, l) ! local variables are not implicit by default ! @@ -966,6 +969,7 @@ module equations ! subroutine arguments ! + real(kind=8) , intent(in) :: x, y real(kind=8), dimension(nv) , intent(in) :: q real(kind=8), dimension(nv) , intent(inout) :: c real(kind=8), dimension(nv,nv), intent(inout) :: l, r @@ -1346,6 +1350,8 @@ module equations ! ! Arguments: ! +! x - ratio of the perpendicular magnetic field component difference +! y - ratio of the density ! q - the intermediate Roe state vector; ! c - the vector of eigenvalues; ! r - the matrix of right eigenvectors; @@ -1363,7 +1369,7 @@ module equations ! !=============================================================================== ! - subroutine esystem_roe_hd_adi(q, c, r, l) + subroutine esystem_roe_hd_adi(x, y, q, c, r, l) ! local variables are not implicit by default ! @@ -1371,6 +1377,7 @@ module equations ! subroutine arguments ! + real(kind=8) , intent(in) :: x, y real(kind=8), dimension(nv) , intent(in) :: q real(kind=8), dimension(nv) , intent(inout) :: c real(kind=8), dimension(nv,nv), intent(inout) :: l, r @@ -1800,6 +1807,276 @@ module equations ! end function maxspeed_mhd_iso ! +!=============================================================================== +! +! subroutine ESYSTEM_ROE_MHD_ISO: +! ------------------------------ +! +! Subroutine computes eigenvalues and eigenvectors for a given set of +! equations and input variables. +! +! Arguments: +! +! x - ratio of the perpendicular magnetic field component difference +! y - ratio of the density +! q - the intermediate Roe state vector; +! c - the vector of eigenvalues; +! r - the matrix of right eigenvectors; +! l - the matrix of left eigenvectors; +! +! References: +! +! [1] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! [2] Balsara, D. S. +! "Linearized Formulation of the Riemann Problem for Adiabatic and +! Isothermal Magnetohydrodynamics", +! The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131 +! +!=============================================================================== +! + subroutine esystem_roe_mhd_iso(x, y, q, c, r, l) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8) , intent(in) :: x, y + real(kind=8), dimension(nv) , intent(in) :: q + real(kind=8), dimension(nv) , intent(inout) :: c + real(kind=8), dimension(nv,nv), intent(inout) :: l, r + +! saved variables +! + logical , save :: first = .true. + +! local variables +! + real(kind=8) :: di, btsq, bt_starsq, casq, twid_csq + real(kind=8) :: ct2, tsum, tdif, cf2_cs2, cfsq, cf, cssq, cs, ca + real(kind=8) :: bt, bt_star, bet2, bet3, bet2_star, bet3_star, bet_starsq + real(kind=8) :: alpha_f, alpha_s + real(kind=8) :: sqrtd, s, twid_c, qf, qs, af_prime, as_prime + real(kind=8) :: norm, cff, css, af, as, afpb, aspb, q2_star, q3_star, vqstr +! +!------------------------------------------------------------------------------- +! +! prepare the internal arrays at the first run +! + if (first) then + +! reset all elements +! + evroe(:, : ,:) = 0.0d+00 + +! unset the first execution flag +! + first = .false. + + end if ! first execution + +! prepare coefficients for eigenvalues +! + di = 1.0d+00 / q(idn) + casq = q(ibx) * q(ibx) * di + ca = sqrt(casq) + btsq = q(iby) * q(iby) + q(ibz) * q(ibz) + bt_starsq = btsq * y + twid_csq = csnd2 + x + ct2 = bt_starsq * di + tsum = casq + ct2 + twid_csq + tdif = casq + ct2 - twid_csq + cf2_cs2 = sqrt(tdif * tdif + 4.0d+00 * twid_csq * ct2) + cfsq = 0.5d+00 * (tsum + cf2_cs2) + cf = sqrt(cfsq) + cssq = twid_csq * casq / cfsq + cs = sqrt(cssq) + +! prepare eigenvalues +! + c(1) = q(ivx) - cf + c(2) = q(ivx) - ca + c(3) = q(ivx) - cs + c(4) = q(ivx) + c(5) = q(ivx) + cs + c(6) = q(ivx) + ca + c(7) = q(ivx) + cf + c(8) = c(7) + +! calculate the eigenvectors only if the waves propagate in both direction +! + if (c(1) >= 0.0d+00) return + if (c(7) <= 0.0d+00) return + +! prepare remaining coefficients for eigenvectors +! + bt = sqrt(btsq) + bt_star = sqrt(bt_starsq) + if (bt == 0.0d+00) then + bet2 = 1.0d+00 + bet3 = 0.0d+00 + else + bet2 = q(iby) / bt + bet3 = q(ibz) / bt + end if + bet2_star = bet2 / sqrt(y) + bet3_star = bet3 / sqrt(y) + bet_starsq = bet2_star * bet2_star + bet3_star * bet3_star + + if ((cfsq - cssq) == 0.0d+00) then + alpha_f = 1.0d+00 + alpha_s = 0.0d+00 + else if ((twid_csq - cssq) <= 0.0d+00) then + alpha_f = 0.0d+00 + alpha_s = 1.0d+00 + else if ((cfsq - twid_csq) <= 0.0d+00) then + alpha_f = 1.0d+00 + alpha_s = 0.0d+00 + else + alpha_f = sqrt((twid_csq - cssq) / (cfsq - cssq)) + alpha_s = sqrt((cfsq - twid_csq) / (cfsq - cssq)) + end if + + sqrtd = sqrt(q(idn)) + s = sign(1.0d+00, q(ibx)) + twid_c = sqrt(twid_csq) + qf = cf * alpha_f * s + qs = cs * alpha_s * s + af_prime = twid_c * alpha_f / sqrtd + as_prime = twid_c * alpha_s / sqrtd + +! update the varying elements of the matrix of right eigenvectors +! +! left-going fast wave +! + evroe(2,1,idn) = alpha_f + evroe(2,1,ivx) = alpha_f * c(1) + evroe(2,1,ivy) = alpha_f * q(ivy) + qs * bet2_star + evroe(2,1,ivz) = alpha_f * q(ivz) + qs * bet3_star + evroe(2,1,iby) = as_prime * bet2_star + evroe(2,1,ibz) = as_prime * bet3_star + +! left-going Alfvèn wave +! + evroe(2,2,ivy) = - bet3 + evroe(2,2,ivz) = bet2 + evroe(2,2,iby) = - bet3 * s / sqrtd + evroe(2,2,ibz) = bet2 * s / sqrtd + +! left-going slow wave +! + evroe(2,3,idn) = alpha_s + evroe(2,3,ivx) = alpha_s * c(3) + evroe(2,3,ivy) = alpha_s * q(ivy) - qf * bet2_star + evroe(2,3,ivz) = alpha_s * q(ivz) - qf * bet3_star + evroe(2,3,iby) = - af_prime * bet2_star + evroe(2,3,ibz) = - af_prime * bet3_star + +! right-going slow wave +! + evroe(2,5,idn) = alpha_s + evroe(2,5,ivx) = alpha_s * c(5) + evroe(2,5,ivy) = alpha_s * q(ivy) + qf * bet2_star + evroe(2,5,ivz) = alpha_s * q(ivz) + qf * bet3_star + evroe(2,5,iby) = evroe(2,3,iby) + evroe(2,5,ibz) = evroe(2,3,ibz) + +! right-going Alfvèn wave +! + evroe(2,6,ivy) = bet3 + evroe(2,6,ivz) = - bet2 + evroe(2,6,iby) = evroe(2,2,iby) + evroe(2,6,ibz) = evroe(2,2,ibz) + +! right-going fast wave +! + evroe(2,7,idn) = alpha_f + evroe(2,7,ivx) = alpha_f * c(7) + evroe(2,7,ivy) = alpha_f * q(ivy) - qs * bet2_star + evroe(2,7,ivz) = alpha_f * q(ivz) - qs * bet3_star + evroe(2,7,iby) = evroe(2,1,iby) + evroe(2,7,ibz) = evroe(2,1,ibz) + +! update the varying elements of the matrix of left eigenvectors +! + norm = 0.5d+00 / twid_csq + cff = norm * alpha_f * cf + css = norm * alpha_s * cs + qf = qf * norm + qs = qs * norm + af = norm * af_prime * q(idn) + as = norm * as_prime * q(idn) + afpb = norm * af_prime * bt_star + aspb = norm * as_prime * bt_star + + q2_star = bet2_star / bet_starsq + q3_star = bet3_star / bet_starsq + vqstr = q(ivy) * q2_star + q(ivz) * q3_star + +! left-going fast wave +! + evroe(1,idn,1) = cff * c(7) - qs * vqstr - aspb + evroe(1,ivx,1) = - cff + evroe(1,ivy,1) = qs * q2_star + evroe(1,ivz,1) = qs * q3_star + evroe(1,iby,1) = as * q2_star + evroe(1,ibz,1) = as * q3_star + +! left-going Alfvèn wave +! + evroe(1,idn,2) = 0.5d+00 * (q(ivy) * bet3 - q(ivz) * bet2) + evroe(1,ivy,2) = - 0.5d+00 * bet3 + evroe(1,ivz,2) = 0.5d+00 * bet2 + evroe(1,iby,2) = - 0.5d+00 * sqrtd * bet3 * s + evroe(1,ibz,2) = 0.5d+00 * sqrtd * bet2 * s + +! left-going slow wave +! + evroe(1,idn,3) = css * c(5) + qf * vqstr + afpb + evroe(1,ivx,3) = - css + evroe(1,ivy,3) = - qf * q2_star + evroe(1,ivz,3) = - qf * q3_star + evroe(1,iby,3) = - af * q2_star + evroe(1,ibz,3) = - af * q3_star + +! right-going slow wave +! + evroe(1,idn,5) = - css * c(3) - qf * vqstr + afpb + evroe(1,ivx,5) = css + evroe(1,ivy,5) = - evroe(1,ivy,3) + evroe(1,ivz,5) = - evroe(1,ivz,3) + evroe(1,iby,5) = evroe(1,iby,3) + evroe(1,ibz,5) = evroe(1,ibz,3) + +! right-going Alfvèn wave +! + evroe(1,idn,6) = - evroe(1,idn,2) + evroe(1,ivy,6) = - evroe(1,ivy,2) + evroe(1,ivz,6) = - evroe(1,ivz,2) + evroe(1,iby,6) = evroe(1,iby,2) + evroe(1,ibz,6) = evroe(1,ibz,2) + +! right-going fast wave +! + evroe(1,idn,7) = - cff * c(1) + qs * vqstr - aspb + evroe(1,ivx,7) = cff + evroe(1,ivy,7) = - evroe(1,ivy,1) + evroe(1,ivz,7) = - evroe(1,ivz,1) + evroe(1,iby,7) = evroe(1,iby,1) + evroe(1,ibz,7) = evroe(1,ibz,1) + +! copy matrices of eigenvectors +! + l(1:nv,1:nv) = evroe(1,1:nv,1:nv) + r(1:nv,1:nv) = evroe(2,1:nv,1:nv) + +!------------------------------------------------------------------------------- +! + end subroutine esystem_roe_mhd_iso +! !******************************************************************************* ! ! ADIABATIC MAGNETOHYDRODYNAMIC EQUATIONS diff --git a/src/schemes.F90 b/src/schemes.F90 index 5a81fa2..606eef6 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -267,6 +267,16 @@ module schemes ! riemann => riemann_mhd_iso_hlldm + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_mhd_iso_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -936,7 +946,7 @@ module schemes ! obtain eigenvalues and eigenvectors ! - call eigensystem_roe(qi(:), ci(:), ri(:,:), li(:,:)) + call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:)) ! return upwind fluxes ! @@ -1632,7 +1642,7 @@ module schemes ! obtain eigenvalues and eigenvectors ! - call eigensystem_roe(qi(:), ci(:), ri(:,:), li(:,:)) + call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:)) ! return upwind fluxes ! @@ -2829,6 +2839,185 @@ module schemes ! !=============================================================================== ! +! subroutine RIEMANN_MHD_ISO_ROE: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Roe's method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 +! +!=============================================================================== +! + subroutine riemann_mhd_iso_roe(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp + use equations , only : imx, imy, imz + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx, yy + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ibx) = ql(ibx,i) + qi(iby) = (sdr * ql(iby,i) + sdl * qr(iby,i)) / sds + qi(ibz) = (sdr * ql(ibz,i) + sdl * qr(ibz,i)) / sds + qi(ibp) = ql(ibp,i) + +! prepare coefficients +! + xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 & + + (ql(ibz,i) - qr(ibz,i))**2) / sds**2 + yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn) + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! prepare fluxes which do not change across the states +! + f(ibx,i) = fl(ibx,i) + f(ibp,i) = fl(ibp,i) + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) + f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) + f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) + f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) + f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) + f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + else + do p = nv, 1, -1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_iso_roe +! +!=============================================================================== +! !***** ADIABATIC MAGNETOHYDRODYNAMICS ***** ! !=============================================================================== From 88edaec837621984b2e2764d74fbe0c174d6676d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Mar 2014 14:19:31 -0300 Subject: [PATCH 6/7] EQUATIONS, SCHEMES: Implement Roe's solver for adiabatic MHD. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 332 +++++++++++++++++++++++++++++++++++++++++++++- src/schemes.F90 | 200 ++++++++++++++++++++++++++++ 2 files changed, 528 insertions(+), 4 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index e553f31..bde35e8 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -371,10 +371,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_mhd_adi - cons2prim => cons2prim_mhd_adi - fluxspeed => fluxspeed_mhd_adi - maxspeed => maxspeed_mhd_adi + prim2cons => prim2cons_mhd_adi + cons2prim => cons2prim_mhd_adi + fluxspeed => fluxspeed_mhd_adi + maxspeed => maxspeed_mhd_adi + eigensystem_roe => esystem_roe_mhd_adi case default @@ -2404,6 +2405,329 @@ module equations !------------------------------------------------------------------------------- ! end function maxspeed_mhd_adi +! +!=============================================================================== +! +! subroutine ESYSTEM_ROE_MHD_ADI: +! ------------------------------ +! +! Subroutine computes eigenvalues and eigenvectors for a given set of +! equations and input variables. +! +! Arguments: +! +! x - ratio of the perpendicular magnetic field component difference +! y - ratio of the density +! q - the intermediate Roe state vector; +! c - the vector of eigenvalues; +! r - the matrix of right eigenvectors; +! l - the matrix of left eigenvectors; +! +! References: +! +! [1] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! [2] Balsara, D. S. +! "Linearized Formulation of the Riemann Problem for Adiabatic and +! Isothermal Magnetohydrodynamics", +! The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131 +! +!=============================================================================== +! + subroutine esystem_roe_mhd_adi(x, y, q, c, r, l) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8) , intent(in) :: x, y + real(kind=8), dimension(nv) , intent(in) :: q + real(kind=8), dimension(nv) , intent(inout) :: c + real(kind=8), dimension(nv,nv), intent(inout) :: l, r + +! saved variables +! + logical , save :: first = .true. + real(kind=8), save :: gammam2 + +! local variables +! + real(kind=8) :: di, vsq, btsq, bt_starsq, casq, hp, twid_asq + real(kind=8) :: ct2, tsum, tdif, cf2_cs2, cfsq, cf, cssq, cs, ca, bt + real(kind=8) :: bt_star, bet2, bet3, bet2_star, bet3_star, bet_starsq, vbet + real(kind=8) :: alpha_f, alpha_s, af_prime, as_prime + real(kind=8) :: sqrtd, isqrtd, s, twid_a, qf, qs, afpbb, aspbb + real(kind=8) :: qa, qb, qc, qd + real(kind=8) :: norm, cff, css, af, as, afpb, aspb + real(kind=8) :: q2_star, q3_star, vqstr + +! local parameters +! + real(kind=8), parameter :: eps = epsilon(di) +! +!------------------------------------------------------------------------------- +! +! prepare the internal arrays at the first run +! + if (first) then + +! prepare coefficients +! + gammam2 = gamma - 2.0d+00 + +! reset all elements +! + evroe(:, : ,:) = 0.0d+00 + +! initiate the matrix of right eigenvectors +! + evroe(2,4,idn) = 1.0d+00 + +! unset the first execution flag +! + first = .false. + + end if ! first execution + +! prepare coefficients for eigenvalues +! + di = 1.0d+00 / q(idn) + casq = q(ibx) * q(ibx) * di + ca = sqrt(casq) + vsq = q(ivx) * q(ivx) + q(ivy) * q(ivy) + q(ivz) * q(ivz) + btsq = q(iby) * q(iby) + q(ibz) * q(ibz) + bt_starsq = (gammam1 - gammam2 * y) * btsq + hp = q(ien) - (casq + btsq * di) + twid_asq = max(eps, (gammam1 * (hp - 0.5d+00 * vsq) - gammam2 * x)) + ct2 = bt_starsq * di + tsum = casq + ct2 + twid_asq + tdif = casq + ct2 - twid_asq + cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_asq * ct2)) + cfsq = 0.5d+00 * (tsum + cf2_cs2) + cf = sqrt(cfsq) + cssq = twid_asq * casq / cfsq + cs = sqrt(cssq) + +! prepare eigenvalues +! + c(1) = q(ivx) - cf + c(2) = q(ivx) - ca + c(3) = q(ivx) - cs + c(4) = q(ivx) + c(5) = q(ivx) + c(6) = q(ivx) + cs + c(7) = q(ivx) + ca + c(8) = q(ivx) + cf + c(9) = c(8) + +! calculate the eigenvectors only if the waves propagate in both direction +! + if (c(1) >= 0.0d+00) return + if (c(8) <= 0.0d+00) return + +! prepare remaining coefficients for eigenvectors +! + bt = sqrt(btsq) + bt_star = sqrt(bt_starsq) + if (bt == 0.0d+00) then + bet2 = 1.0d+00 + bet3 = 0.0d+00 + else + bet2 = q(iby) / bt + bet3 = q(ibz) / bt + end if + bet2_star = bet2 / sqrt(gammam1 - gammam2 * y) + bet3_star = bet3 / sqrt(gammam1 - gammam2 * y) + bet_starsq = bet2_star * bet2_star + bet3_star * bet3_star + vbet = q(ivy) * bet2_star + q(ivz) * bet3_star + + if ( (cfsq - cssq) == 0.0d+00 ) then + alpha_f = 1.0d+00 + alpha_s = 0.0d+00 + else if ( (twid_asq - cssq) <= 0.0d+00 ) then + alpha_f = 0.0d+00 + alpha_s = 1.0d+00 + else if ( (cfsq - twid_asq) <= 0.0d+00 ) then + alpha_f = 1.0d+00 + alpha_s = 0.0d+00 + else + alpha_f = sqrt((twid_asq - cssq) / (cfsq - cssq)) + alpha_s = sqrt((cfsq - twid_asq) / (cfsq - cssq)) + end if + + sqrtd = sqrt(q(idn)) + isqrtd = 1.0d+00 / sqrtd + s = sign(1.0d+00, q(ibx)) + twid_a = sqrt(twid_asq) + qf = cf * alpha_f * s + qs = cs * alpha_s * s + af_prime = twid_a * alpha_f * isqrtd + as_prime = twid_a * alpha_s * isqrtd + afpbb = af_prime * bt_star * bet_starsq + aspbb = as_prime * bt_star * bet_starsq + +! update the varying elements of the matrix of right eigenvectors +! + evroe(2,1,idn) = alpha_f + evroe(2,3,idn) = alpha_s + evroe(2,6,idn) = alpha_s + evroe(2,8,idn) = alpha_f + + evroe(2,1,ivx) = alpha_f * c(1) + evroe(2,3,ivx) = alpha_s * c(3) + evroe(2,4,ivx) = q(ivx) + evroe(2,6,ivx) = alpha_s * c(6) + evroe(2,8,ivx) = alpha_f * c(8) + + qa = alpha_f * q(ivy) + qb = alpha_s * q(ivy) + qc = qs * bet2_star + qd = qf * bet2_star + + evroe(2,1,ivy) = qa + qc + evroe(2,2,ivy) = - bet3 + evroe(2,3,ivy) = qb - qd + evroe(2,4,ivy) = q(ivy) + evroe(2,6,ivy) = qb + qd + evroe(2,7,ivy) = bet3 + evroe(2,8,ivy) = qa - qc + + qa = alpha_f * q(ivz) + qb = alpha_s * q(ivz) + qc = qs * bet3_star + qd = qf * bet3_star + + evroe(2,1,ivz) = qa + qc + evroe(2,2,ivz) = bet2 + evroe(2,3,ivz) = qb - qd + evroe(2,4,ivz) = q(ivz) + evroe(2,6,ivz) = qb + qd + evroe(2,7,ivz) = - bet2 + evroe(2,8,ivz) = qa - qc + + evroe(2,1,ipr) = alpha_f * (hp - q(ivx) * cf) + qs * vbet + aspbb + evroe(2,2,ipr) = -(q(ivy) * bet3 - q(ivz) * bet2) + evroe(2,3,ipr) = alpha_s * (hp - q(ivx) * cs) - qf * vbet - afpbb + evroe(2,4,ipr) = 0.5d+00 * vsq + gammam2 * x / gammam1 + evroe(2,6,ipr) = alpha_s * (hp + q(ivx) * cs) + qf * vbet - afpbb + evroe(2,7,ipr) = - evroe(2,2,ipr) + evroe(2,8,ipr) = alpha_f * (hp + q(ivx) * cf) - qs * vbet + aspbb + + evroe(2,1,iby) = as_prime * bet2_star + evroe(2,2,iby) = - bet3 * s * isqrtd + evroe(2,3,iby) = - af_prime * bet2_star + evroe(2,6,iby) = evroe(2,3,iby) + evroe(2,7,iby) = evroe(2,2,iby) + evroe(2,8,iby) = evroe(2,1,iby) + + evroe(2,1,ibz) = as_prime * bet3_star + evroe(2,2,ibz) = bet2 * s * isqrtd + evroe(2,3,ibz) = - af_prime * bet3_star + evroe(2,6,ibz) = evroe(2,3,ibz) + evroe(2,7,ibz) = evroe(2,2,ibz) + evroe(2,8,ibz) = evroe(2,1,ibz) + +! update the varying elements of the matrix of left eigenvectors +! + norm = 0.5d+00 / twid_asq + cff = norm * alpha_f * cf + css = norm * alpha_s * cs + qf = qf * norm + qs = qs * norm + af = norm * af_prime * q(idn) + as = norm * as_prime * q(idn) + afpb = norm * af_prime * bt_star + aspb = norm * as_prime * bt_star + + norm = norm * gammam1 + alpha_f = alpha_f * norm + alpha_s = alpha_s * norm + q2_star = bet2_star / bet_starsq + q3_star = bet3_star / bet_starsq + vqstr = (q(ivy) * q2_star + q(ivz) * q3_star) + norm = 2.0d+00 * norm + +! left-going fast wave +! + evroe(1,idn,1) = alpha_f * (vsq - hp) + cff * (cf + q(ivx)) & + - qs * vqstr - aspb + evroe(1,ivx,1) = - alpha_f * q(ivx) - cff + evroe(1,ivy,1) = - alpha_f * q(ivy) + qs * q2_star + evroe(1,ivz,1) = - alpha_f * q(ivz) + qs * q3_star + evroe(1,ipr,1) = alpha_f + evroe(1,iby,1) = as * q2_star - alpha_f * q(iby) + evroe(1,ibz,1) = as * q3_star - alpha_f * q(ibz) + +! left-going Alfvèn wave +! + evroe(1,idn,2) = 0.5d+00 * (q(ivy) * bet3 - q(ivz) * bet2) + evroe(1,ivy,2) = - 0.5d+00 * bet3 + evroe(1,ivz,2) = 0.5d+00 * bet2 + evroe(1,iby,2) = - 0.5d+00 * sqrtd * bet3 * s + evroe(1,ibz,2) = 0.5d+00 * sqrtd * bet2 * s + +! left-going slow wave +! + evroe(1,idn,3) = alpha_s * (vsq - hp) + css * (cs + q(ivx)) & + + qf * vqstr + afpb + evroe(1,ivx,3) = - alpha_s * q(ivx) - css + evroe(1,ivy,3) = - alpha_s * q(ivy) - qf * q2_star + evroe(1,ivz,3) = - alpha_s * q(ivz) - qf * q3_star + evroe(1,ipr,3) = alpha_s + evroe(1,iby,3) = - af * q2_star - alpha_s * q(iby) + evroe(1,ibz,3) = - af * q3_star - alpha_s * q(ibz) + +! entropy wave +! + evroe(1,idn,4) = 1.0d+00 - norm * (0.5d+00 * vsq - gammam2 * x / gammam1) + evroe(1,ivx,4) = norm * q(ivx) + evroe(1,ivy,4) = norm * q(ivy) + evroe(1,ivz,4) = norm * q(ivz) + evroe(1,ipr,4) = - norm + evroe(1,iby,4) = norm * q(iby) + evroe(1,ibz,4) = norm * q(ibz) + +! right-going slow wave +! + evroe(1,idn,6) = alpha_s * (vsq - hp) + css * (cs - q(ivx)) & + - qf * vqstr + afpb + evroe(1,ivx,6) = - alpha_s * q(ivx) + css + evroe(1,ivy,6) = - alpha_s * q(ivy) + qf * q2_star + evroe(1,ivz,6) = - alpha_s * q(ivz) + qf * q3_star + evroe(1,ipr,6) = alpha_s + evroe(1,iby,6) = evroe(1,iby,3) + evroe(1,ibz,6) = evroe(1,ibz,3) + +! right-going Alfvèn wave +! + evroe(1,idn,7) = - evroe(1,idn,2) + evroe(1,ivy,7) = - evroe(1,ivy,2) + evroe(1,ivz,7) = - evroe(1,ivz,2) + evroe(1,iby,7) = evroe(1,iby,2) + evroe(1,ibz,7) = evroe(1,ibz,2) + +! right-going fast wave +! + evroe(1,idn,8) = alpha_f * (vsq - hp) + cff * (cf - q(ivx)) & + + qs * vqstr - aspb + evroe(1,ivx,8) = - alpha_f * q(ivx) + cff + evroe(1,ivy,8) = - alpha_f * q(ivy) - qs * q2_star + evroe(1,ivz,8) = - alpha_f * q(ivz) - qs * q3_star + evroe(1,ipr,8) = alpha_f + evroe(1,iby,8) = evroe(1,iby,1) + evroe(1,ibz,8) = evroe(1,ibz,1) + +! copy matrices of eigenvectors +! + l(1:nv,1:nv) = evroe(1,1:nv,1:nv) + r(1:nv,1:nv) = evroe(2,1:nv,1:nv) + +!------------------------------------------------------------------------------- +! + end subroutine esystem_roe_mhd_adi !=============================================================================== ! diff --git a/src/schemes.F90 b/src/schemes.F90 index 606eef6..d1d53a5 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -322,6 +322,16 @@ module schemes ! riemann => riemann_mhd_adi_hlld + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_mhd_adi_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -4283,6 +4293,196 @@ module schemes !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_adi_hlld +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ADI_ROE: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Roe's method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Stone, J. M. & Gardiner, T. A., +! "ATHENA: A New Code for Astrophysical MHD", +! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 +! +!=============================================================================== +! + subroutine riemann_mhd_adi_roe(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: pml, pmr + real(kind=8) :: xx, yy + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare magnetic pressures +! + pml = 0.5d+00 * sum(ql(ibx:ibz,i)**2) + pmr = 0.5d+00 * sum(qr(ibx:ibz,i)**2) + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl & + + (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds + qi(ibx) = ql(ibx,i) + qi(iby) = (sdr * ql(iby,i) + sdl * qr(iby,i)) / sds + qi(ibz) = (sdr * ql(ibz,i) + sdl * qr(ibz,i)) / sds + qi(ibp) = ql(ibp,i) + +! prepare coefficients +! + xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 & + + (ql(ibz,i) - qr(ibz,i))**2) / sds**2 + yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn) + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! prepare fluxes which do not change across the states +! + f(ibx,i) = fl(ibx,i) + f(ibp,i) = fl(ibp,i) + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) + f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) + f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) + f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) + f(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i)) + f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) + f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(ien,i) = f(ien,i) + xx * ri(p,ien) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + else + do p = nv, 1, -1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(ien,i) = f(ien,i) + xx * ri(p,ien) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_adi_roe !=============================================================================== ! From be40a6113c58d95224e7b0f8b2b30d7414de005f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 19 Mar 2014 19:04:42 -0300 Subject: [PATCH 7/7] IO: Store ghost cells in the snapshots by default. This is just in case we would like to use interpolation to reconstruct the full domain resolution. Signed-off-by: Grzegorz Kowal --- src/io.F90 | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 7c95c6c..cecf5a3 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -61,17 +61,17 @@ module io ! (in hours); ! hsnap - the problem time interval for regular snapshot storing; ! - character(len=255), save :: respath = "./" - character , save :: ftype = "p" - integer , save :: nrest = -1 - integer(kind=4) , save :: irest = 1 - integer(kind=4) , save :: isnap = -1 - real , save :: hrest = 6.0d+00 - real , save :: hsnap = 1.0d+00 + character(len=255), save :: respath = "./" + character , save :: ftype = "p" + integer , save :: nrest = -1 + integer(kind=4) , save :: irest = 1 + integer(kind=4) , save :: isnap = -1 + real , save :: hrest = 6.0d+00 + real , save :: hsnap = 1.0d+00 ! flags to determine the way of data writing ! - logical , save :: with_ghosts = .false. + logical , save :: with_ghosts = .true. ! local variables to store the number of processors and maximum level read from ! the restart file @@ -141,7 +141,7 @@ module io ! local variables ! - character(len=255) :: ghosts = "off" + character(len=255) :: ghosts = "on" integer :: dd, hh, mm, ss ! !------------------------------------------------------------------------------- @@ -171,15 +171,15 @@ module io ! get the flag determining if the ghost cells are stored ! - call get_parameter_string ("include_ghosts", ghosts ) + call get_parameter_string ("include_ghosts" , ghosts ) ! check ghost cell storing flag ! select case(trim(ghosts)) - case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - with_ghosts = .true. - case default + case ("off", "OFF", "n", "N", "false", "FALSE", "no", "NO") with_ghosts = .false. + case default + with_ghosts = .true. end select ! return the run number @@ -193,6 +193,11 @@ module io "snapshot type", "primitive variables" if (ftype == 'c') write (*,"(4x,a13,10x,'=',1x,a)") & "snapshot type", "conservative variables" + if (with_ghosts) then + write (*,"(4x,a21,2x,'=',1x,a)") "with ghosts cells ", "on" + else + write (*,"(4x,a21,2x,'=',1x,a)") "with ghosts cells ", "off" + end if write (*,"(4x,a21,2x,'=',1x,e9.2)") "snapshot interval ", hsnap if (hrest > 0.0d+00) then dd = int(hrest / 2.4d+01)