gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/schemes.F90
Grzegorz Kowal 1acb7cca04 SCHEMES: Add a separate subroutine for higher order flux corrections. 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2020-03-03 10:44:04 -03:00

6439 lines
169 KiB
Fortran
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2008-2020 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: SCHEMES
!!
!! The module provides and interface to numerical schemes, subroutines to
!! calculate variable increment and one dimensional Riemann solvers.
!!
!! If you implement a new set of equations, you have to add at least one
!! corresponding update_flux subroutine, and one Riemann solver.
!!
!!******************************************************************************
!
module schemes
#ifdef PROFILE
! import external subroutines
!
use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */
! module variables are not implicit by default
!
implicit none
#ifdef PROFILE
! timer indices
!
integer , save :: imi, imf, ims, imr
#endif /* PROFILE */
! Rieman solver and state vectors names
!
character(len=255) , save :: name_sol = ""
character(len=255) , save :: name_sts = "primitive"
! 4-vector reconstruction flag
!
logical , save :: states_4vec = .false.
! high order fluxes
!
logical , save :: high_order_flux = .false.
! interfaces for procedure pointers
!
abstract interface
subroutine update_flux_iface(dx, q, f)
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
end subroutine
subroutine riemann_iface(q, f)
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
end subroutine
end interface
! pointer to the flux update procedure
!
procedure(update_flux_iface), pointer, save :: update_flux => null()
! pointer to the Riemann solver
!
procedure(riemann_iface) , pointer, save :: riemann => null()
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_schemes, finalize_schemes, print_schemes
public :: update_flux
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!!
!!*** PUBLIC SUBROUTINES *****************************************************
!!
!===============================================================================
!
! subroutine INITIALIZE_SCHEMES:
! -----------------------------
!
! Subroutine initiate the module by setting module parameters and subroutine
! pointers.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine initialize_schemes(verbose, status)
! include external procedures and variables
!
use equations , only : eqsys, eos
use parameters, only : get_parameter
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
integer, intent(out) :: status
! local variables
!
character(len=64) :: solver = "HLL"
character(len=64) :: statev = "primitive"
character(len=64) :: flux = "off"
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('schemes:: initialization' , imi)
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
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! get the Riemann solver
!
call get_parameter("riemann_solver" , solver)
call get_parameter("state_variables", statev)
! depending on the system of equations initialize the module variables
!
select case(trim(eqsys))
!--- HYDRODYNAMICS ---
!
case("hd", "HD", "hydro", "HYDRO", "hydrodynamic", "HYDRODYNAMIC")
! depending on the equation of state complete the initialization
!
select case(trim(eos))
case("iso", "ISO", "isothermal", "ISOTHERMAL")
! set pointers to subroutines
!
update_flux => update_flux_hd_iso
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_hd_iso_hll
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann => riemann_hd_iso_roe
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for isothermal HD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'roe'."
end if
status = 1
end select
case("adi", "ADI", "adiabatic", "ADIABATIC")
! set pointers to subroutines
!
update_flux => update_flux_hd_adi
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_hd_adi_hll
case("hllc", "HLLC")
! set the solver name
!
name_sol = "HLLC"
! set pointers to subroutines
!
riemann => riemann_hd_adi_hllc
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann => riemann_hd_adi_roe
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for adiabatic HD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc', 'roe'."
end if
status = 1
end select
end select
!--- MAGNETOHYDRODYNAMICS ---
!
case("mhd", "MHD", "magnetohydrodynamic", "MAGNETOHYDRODYNAMIC")
! depending on the equation of state complete the initialization
!
select case(trim(eos))
case("iso", "ISO", "isothermal", "ISOTHERMAL")
! set pointers to the subroutines
!
update_flux => update_flux_mhd_iso
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_mhd_iso_hll
case ("hlld", "HLLD")
! set the solver name
!
name_sol = "HLLD"
! set the solver pointer
!
riemann => riemann_mhd_iso_hlld
case ("hlldm", "hlld-m", "HLLDM", "HLLD-M")
! set the solver name
!
name_sol = "HLLD-M"
! set the solver pointer
!
riemann => riemann_mhd_iso_hlldm
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann => riemann_mhd_iso_roe
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for isothermal MHD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hlld'" // &
", 'hlld-m', 'roe'."
end if
status = 1
end select
case("adi", "ADI", "adiabatic", "ADIABATIC")
! set pointers to subroutines
!
update_flux => update_flux_mhd_adi
! select the Riemann solver
!
select case(trim(solver))
case ("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_mhd_adi_hll
case ("hllc", "HLLC")
! set the solver name
!
name_sol = "HLLC"
! set the solver pointer
!
riemann => riemann_mhd_adi_hllc
case ("hlld", "HLLD")
! set the solver name
!
name_sol = "HLLD"
! set the solver pointer
!
riemann => riemann_mhd_adi_hlld
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann => riemann_mhd_adi_roe
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for adiabatic MHD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'" // &
", 'hlld', 'roe'."
end if
status = 1
end select
end select
!--- SPECIAL RELATIVITY HYDRODYNAMICS ---
!
case("srhd", "SRHD")
! depending on the equation of state complete the initialization
!
select case(trim(eos))
case("adi", "ADI", "adiabatic", "ADIABATIC")
! set pointers to subroutines
!
update_flux => update_flux_srhd_adi
! select the state reconstruction method
!
select case(trim(statev))
case("4vec", "4-vector", "4VEC", "4-VECTOR")
! set the state reconstruction name
!
name_sts = "4-vector"
! set 4-vector reconstruction flag
!
states_4vec = .true.
! in the case of state variables, revert to primitive
!
case default
! set the state reconstruction name
!
name_sts = "primitive"
end select
! select the Riemann solver
!
select case(trim(solver))
case ("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_srhd_adi_hll
case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
! set the solver name
!
name_sol = "HLLC (Mignone & Bodo 2005)"
! set pointers to subroutines
!
riemann => riemann_srhd_adi_hllc
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for adiabatic SR-HD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'."
end if
status = 1
end select
end select
!--- SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS ---
!
case("srmhd", "SRMHD")
! depending on the equation of state complete the initialization
!
select case(trim(eos))
case("adi", "ADI", "adiabatic", "ADIABATIC")
! set pointers to subroutines
!
update_flux => update_flux_srmhd_adi
! select the state reconstruction method
!
select case(trim(statev))
case("4vec", "4-vector", "4VEC", "4-VECTOR")
! set the state reconstruction name
!
name_sts = "4-vector"
! set 4-vector reconstruction flag
!
states_4vec = .true.
! in the case of state variables, revert to primitive
!
case default
! set the state reconstruction name
!
name_sts = "primitive"
end select
! select the Riemann solver
!
select case(trim(solver))
case ("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann => riemann_srmhd_adi_hll
case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
! set the solver name
!
name_sol = "HLLC (Mignone & Bodo 2006)"
! set pointers to subroutines
!
riemann => riemann_srmhd_adi_hllc
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected Riemann solver implemented " // &
"for adiabatic SR-MHD: '" // trim(solver) // "'."
write(*,"(1x,a)") "Available Riemann solvers: 'hll', 'hllc'."
end if
status = 1
end select
end select
end select
! flag for higher order flux correction
!
call get_parameter("high_order_flux", flux)
select case(trim(flux))
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
high_order_flux = .true.
case default
high_order_flux = .false.
end select
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine initialize_schemes
!
!===============================================================================
!
! subroutine FINALIZE_SCHEMES:
! ---------------------------
!
! Subroutine releases memory used by the module.
!
! Arguments:
!
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine finalize_schemes(status)
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! nullify procedure pointers
!
nullify(update_flux)
nullify(riemann)
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine finalize_schemes
!
!===============================================================================
!
! subroutine PRINT_SCHEMES:
! ------------------------
!
! Subroutine prints module parameters and settings.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
!
!===============================================================================
!
subroutine print_schemes(verbose)
! import external procedures and variables
!
use helpers, only : print_section, print_parameter
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
!
!-------------------------------------------------------------------------------
!
if (verbose) then
call print_section(verbose, "Schemes")
call print_parameter(verbose, "Riemann solver" , name_sol)
call print_parameter(verbose, "state variables", name_sts)
if (high_order_flux) then
call print_parameter(verbose, "high order flux correction", "on" )
else
call print_parameter(verbose, "high order flux correction", "off")
end if
end if
!-------------------------------------------------------------------------------
!
end subroutine print_schemes
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine RECONSTRUCT_INTERFACES:
! ---------------------------------
!
! Subroutine reconstructs the Riemann states using 3D reconstruction methods.
!
! Arguments:
!
! dx - the spatial step;
! q - the array of primitive variables;
! qi - the array of reconstructed states (2 in each direction);
!
!===============================================================================
!
subroutine reconstruct_interfaces(dx, q, qi)
! include external procedures
!
use equations , only : nf, idn, ipr
use interpolations , only : interfaces
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:,:), intent(out) :: qi
! local variables
!
logical :: positive
integer :: p
!
!-------------------------------------------------------------------------------
!
! iterate over all flux variables
!
do p = 1, nf
! determine if the variable is positive
!
positive = (p == idn .or. p == ipr)
! interpolate interfaces
!
call interfaces(positive, dx(1:NDIMS), q (p,:,:,:), &
qi(p,:,:,:,1:2,1:NDIMS))
end do ! p = 1, nf
!-------------------------------------------------------------------------------
!
end subroutine reconstruct_interfaces
!
!===============================================================================
!
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
! ---------------------------------------
!
! Subroutine adds higher order correctione to the numerical fluxes.
!
! Arguments:
!
! f - the vecotr of numerical fluxes;
!
!===============================================================================
!
subroutine higher_order_flux_correction(f)
! include external procedures
!
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(inout) :: f
! local variables
!
integer :: n
! local arrays to store the states
!
real(kind=8), dimension(size(f,1),size(f,2)) :: f2, f4
!
!-------------------------------------------------------------------------------
!
! depending on the scheme order calculate and add higher order corrections,
! if required
!
if (.not. high_order_flux .or. order <= 2) return
! get the length of vector
!
n = size(f, 2)
! 2nd order correction
!
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
! 4th order correction
!
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
! correct the flux for 4th order
!
f(:,:) = f(:,:) + f4(:,:)
end if
! correct the flux for 2nd order
!
f(:,:) = f(:,:) + f2(:,:)
!-------------------------------------------------------------------------------
!
end subroutine higher_order_flux_correction
!
!===============================================================================
!
!***** ISOTHERMAL HYDRODYNAMICS *****
!
!===============================================================================
!
! subroutine UPDATE_FLUX_HD_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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_hd_iso(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
! local variables are not implicit by default
!
implicit none
! input arguments
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy states to directional lines with proper vector component ordering
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#ifdef PROFILE
! stop accounting time for flux update
!
call stop_timer(imf)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_flux_hd_iso
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_HLL:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate the fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_iso_hll
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:))
! 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(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:))
! return upwind fluxes
!
if (ci(1) >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
else
! calculate wave amplitudes α = L.ΔU
!
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
!
f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nf
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
else
do p = nf, 1, - 1
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
end if
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_iso_roe
!
!===============================================================================
!
!***** 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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_hd_adi(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
! local variables are not implicit by default
!
implicit none
! input arguments
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy states to directional lines with proper vector component ordering
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
f(1,ien,:,j,k) = fi(ien,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
f(2,ien,i,:,k) = fi(ien,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
f(3,ien,i,j,:) = fi(ien,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#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:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all position
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate the inverse of speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for the Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_hll
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_HLLC:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using the HLLC method,
! by Toro. In the HLLC method the tangential components of the velocity are
! discontinuous actoss the contact dicontinuity.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Toro, E. F., Spruce, M., & Speares, W.
! "Restoration of the contact surface in the HLL-Riemann solver",
! Shock Waves, 1994, Volume 4, Issue 1, pp. 25-34
!
!===============================================================================
!
subroutine riemann_hd_adi_hllc(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, sm
real(kind=8) :: slmm, srmm
real(kind=8) :: dn, pr
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, ui
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the speed of contact discontinuity
!
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! calculate the pressure of the intermediate state
!
pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! the left speed difference
!
slmm = sl - sm
! conservative variables for the left intermediate state
!
ui(idn) = wl(idn) / slmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * ql(ivy,i)
ui(imz) = ui(idn) * ql(ivz,i)
ui(ien) = (wl(ien) + sm * pr) / slmm
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! the right speed difference
!
srmm = sr - sm
! conservative variables for the right intermediate state
!
ui(idn) = wr(idn) / srmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * qr(ivy,i)
ui(imz) = ui(idn) * qr(ivz,i)
ui(ien) = (wr(ien) + sm * pr) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! intermediate flux is constant across the contact discontinuity and all except
! the parallel momentum flux are zero
!
f(idn,i) = 0.0d+00
f(imx,i) = - wl(imx)
f(imy,i) = 0.0d+00
f(imz,i) = 0.0d+00
f(ien,i) = 0.0d+00
end if ! sm = 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_hllc
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, ien
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:))
! 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(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:))
! return upwind fluxes
!
if (ci(1) >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (ci(nf) <= 0.0d+00) then
f(:,i) = fr(:,i)
else
! calculate wave amplitudes α = L.ΔU
!
al(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
end do
! calculate the flux average
!
f(:,i) = 0.5d+00 * (fl(:,i) + fr(:,i))
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
!
if (qi(ivx) >= 0.0d+00) then
do p = 1, nf
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
else
do p = nf, 1, - 1
f(:,i) = f(:,i) - 0.5d+00 * al(p) * abs(ci(p)) * ri(p,:)
end do
end if
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_roe
!
!===============================================================================
!
!***** 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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_mhd_iso(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
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
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
f(1,ibx,:,j,k) = fi(ibx,:)
f(1,iby,:,j,k) = fi(iby,:)
f(1,ibz,:,j,k) = fi(ibz,:)
f(1,ibp,:,j,k) = fi(ibp,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
f(2,ibx,i,:,k) = fi(ibz,:)
f(2,iby,i,:,k) = fi(ibx,:)
f(2,ibz,i,:,k) = fi(iby,:)
f(2,ibp,i,:,k) = fi(ibp,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
f(3,ibx,i,j,:) = fi(iby,:)
f(3,iby,i,j,:) = fi(ibz,:)
f(3,ibz,i,j,:) = fi(ibx,:)
f(3,ibp,i,j,:) = fi(ibp,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#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.
!
! Arguments:
!
! q - 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(q, f)
! include external variables and procedures
!
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
real(kind=8) :: bx, bp
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate the fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_hll
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_HLLD:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
! method with correct state averaging and degeneracies treatement.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Kowal, G.,
! "An isothermal MHD Riemann solver with correct state averaging and
! degeneracies treatement",
! Journal of Computational Physics, in preparation
!
!===============================================================================
!
subroutine riemann_mhd_iso_hlld(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the advection speed in the intermediate states
!
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! speed differences
!
slmm = sl - sm
srmm = sr - sm
! left and right state densities
!
dnl = wl(idn) / slmm
dnr = wr(idn) / srmm
! apply full HLLD solver only if the Alfvén wave is strong enough
!
if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
! left and right Alfvén speeds
!
sml = sm - sqrt(b2 / dnl)
smr = sm + sqrt(b2 / dnr)
! calculate denominators in order to check degeneracy
!
dvl = slmm * wl(idn) - b2
dvr = srmm * wr(idn) - b2
! check degeneracy Sl* -> Sl or Sr* -> Sr
!
if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (smr <= 0.0d+00) then ! sr* ≤ 0
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sl* < 0 < sr*
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
end if ! sl* < 0 < sr*
else ! one degeneracy
! separate degeneracies
!
if (dvl > 0.0d+00) then ! sr* > sr
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else ! sl* < 0
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
end if ! sl* < 0
else if (dvr > 0.0d+00) then ! sl* < sl
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! choose the correct state depending on the speed signs
!
if (smr <= 0.0d+00) then ! sr* ≤ 0
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sr* > 0
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
end if ! sr* > 0
else ! sl* < sl & sr* > sr
! both outer states are degenerate so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl* < sl & sr* > sr
end if ! one degeneracy
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! the Alfvén wave is too weak, so ignore it in order to not introduce any
! numerical instabilities; since Bₓ is not zero, the perpendicular components
! of velocity and magnetic field are continuous; we have only one intermediate
! state, therefore simply apply the HLL solver
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
else ! Bₓ² = 0
! take the right flux depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! conservative variables for the left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = wl(imy) / slmm
ui(imz) = wl(imz) / slmm
ui(ibx) = 0.0d+00
ui(iby) = wl(iby) / slmm
ui(ibz) = wl(ibz) / slmm
ui(ibp) = ul(ibp,i)
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! conservative variables for the right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = wr(imy) / srmm
ui(imz) = wr(imz) / srmm
ui(ibx) = 0.0d+00
ui(iby) = wr(iby) / srmm
ui(ibz) = wr(ibz) / srmm
ui(ibp) = ur(ibp,i)
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! the intermediate flux; since the advection speed is zero, perpendicular
! components do not change, so revert it to HLL
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
end if ! Bₓ² = 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_hlld
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_HLLDM:
! --------------------------------
!
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
! method described by Mignone.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Mignone, A.,
! "A simple and accurate Riemann solver for isothermal MHD",
! Journal of Computational Physics, 2007, 225, pp. 1427-1441
!
!===============================================================================
!
subroutine riemann_mhd_iso_hlldm(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr, ca, ca2
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the advection speed in the intermediate states
!
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! speed differences
!
slmm = sl - sm
srmm = sr - sm
! calculate density of the intermediate states
!
dn = dn / srml
! get the square of the Alfvén speed
!
ca2 = b2 / dn
! apply full HLLD solver only if the Alfvén wave is strong enough
!
if (ca2 > 1.0d-08 * max(sl**2, sr**2)) then ! Bₓ² > ε
! left and right Alfvén speeds
!
ca = sqrt(ca2)
sml = sm - ca
smr = sm + ca
! calculate denominators in order to check degeneracy
!
dvl = slmm * slmm * dn - b2
dvr = srmm * srmm * dn - b2
! check degeneracy Sl* -> Sl or Sr* -> Sr
!
if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! conservative variables for the outer left intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (smr <= 0.0d+00) then ! sr* ≤ 0
! conservative variables for the outer right intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sl* < 0 < sr*
! conservative variables for the outer left intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! conservative variables for the outer right intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
end if ! sl* < 0 < sr*
else ! one degeneracy
! separate degeneracies
!
if (dvl > 0.0d+00) then ! sr* > sr
! conservative variables for the outer left intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl
ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl
ui(ibx) = bx
ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl
ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl
ui(ibp) = ul(ibp,i)
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else ! sl* < 0
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
end if ! sl* < 0
else if (dvr > 0.0d+00) then ! sl* < sl
! conservative variables for the outer right intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr
ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr
ui(ibx) = bx
ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr
ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr
ui(ibp) = ur(ibp,i)
! choose the correct state depending on the speed signs
!
if (smr <= 0.0d+00) then ! sr* ≤ 0
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sr* > 0
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
end if ! sr* > 0
else ! sl* < sl & sr* > sr
! both outer states are degenerate so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl* < sl & sr* > sr
end if ! one degeneracy
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! the Alfvén wave is very weak, so ignore it in order to not introduce any
! numerical instabilities; since Bₓ is not zero, the perpendicular components
! of velocity and magnetic field are continuous; we have only one intermediate
! state, therefore simply apply the HLL solver
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
else ! Bₓ² = 0
! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and
! still we can have a discontinuity in the perpendicular components
!
! take the right flux depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! conservative variables for the left intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = wl(imy) / slmm
ui(imz) = wl(imz) / slmm
ui(ibx) = 0.0d+00
ui(iby) = wl(iby) / slmm
ui(ibz) = wl(ibz) / slmm
ui(ibp) = ul(ibp,i)
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! conservative variables for the right intermediate state
!
ui(idn) = dn
ui(imx) = dn * sm
ui(imy) = wr(imy) / srmm
ui(imz) = wr(imz) / srmm
ui(ibx) = 0.0d+00
ui(iby) = wr(iby) / srmm
ui(ibz) = wr(ibz) / srmm
ui(ibp) = ur(ibp,i)
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! both states are equal so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
end if ! Bₓ² = 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_hlldm
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp
use equations , only : imx, imy, imz
use equations , only : cmax
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: bx, bp
real(kind=8) :: xx, yy
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:))
! 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(nf) <= 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(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (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, nf
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 = nf, 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
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_roe
!
!===============================================================================
!
!***** 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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_mhd_adi(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
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
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
f(1,ibx,:,j,k) = fi(ibx,:)
f(1,iby,:,j,k) = fi(iby,:)
f(1,ibz,:,j,k) = fi(ibz,:)
f(1,ibp,:,j,k) = fi(ibp,:)
f(1,ien,:,j,k) = fi(ien,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
f(2,ibx,i,:,k) = fi(ibz,:)
f(2,iby,i,:,k) = fi(ibx,:)
f(2,ibz,i,:,k) = fi(iby,:)
f(2,ibp,i,:,k) = fi(ibp,:)
f(2,ien,i,:,k) = fi(ien,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
f(3,ibx,i,j,:) = fi(iby,:)
f(3,iby,i,j,:) = fi(ibz,:)
f(3,ibz,i,j,:) = fi(ibx,:)
f(3,ibp,i,j,:) = fi(ibp,:)
f(3,ien,i,j,:) = fi(ien,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#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:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
real(kind=8) :: bx, bp
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate the fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_hll
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_HLLC:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the HLLC
! method by Gurski or Li. The HLLC and HLLC-C differ by definitions of
! the tangential components of the velocity and magnetic field.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Gurski, K. F.,
! "An HLLC-Type Approximate Riemann Solver for Ideal
! Magnetohydrodynamics",
! SIAM Journal on Scientific Computing, 2004, Volume 25, Issue 6,
! pp. 2165–2187
! [2] Li, S.,
! "An HLLC Riemann solver for magneto-hydrodynamics",
! Journal of Computational Physics, 2005, Volume 203, Issue 1,
! pp. 344-357
!
!===============================================================================
!
subroutine riemann_mhd_adi_hllc(q, f)
! include external procedures
!
use equations , only : nf
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 : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
real(kind=8) :: dn, bx, bp, b2, pt, vy, vz, by, bz, vb
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, ui
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLLC flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the speed of contact discontinuity
!
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! square of Bx, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! separate the cases when Bₓ = 0 and Bₓ ≠ 0
!
if (b2 > 0.0d+00) then
! the total pressure, constant across the contact discontinuity
!
pt = 0.5d+00 * ((wl(idn) * sm - wl(imx)) &
+ (wr(idn) * sm - wr(imx))) + b2
! constant intermediate state tangential components of velocity and
! magnetic field
!
vy = (wr(imy) - wl(imy)) / dn
vz = (wr(imz) - wl(imz)) / dn
by = (wr(iby) - wl(iby)) / srml
bz = (wr(ibz) - wl(ibz)) / srml
! scalar product of velocity and magnetic field for the intermediate states
!
vb = sm * bx + vy * by + vz * bz
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! the left speed difference
!
slmm = sl - sm
! conservative variables for the left intermediate state
!
ui(idn) = wl(idn) / slmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * vy
ui(imz) = ui(idn) * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ql(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! the right speed difference
!
srmm = sr - sm
! conservative variables for the right intermediate state
!
ui(idn) = wr(idn) / srmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * vy
ui(imz) = ui(idn) * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = qr(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
else ! Bₓ = 0
! the total pressure, constant across the contact discontinuity
!
pt = 0.5d+00 * ((wl(idn) * sm - wl(imx)) &
+ (wr(idn) * sm - wr(imx)))
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! the left speed difference
!
slmm = sl - sm
! conservative variables for the left intermediate state
!
ui(idn) = wl(idn) / slmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * ql(ivy,i)
ui(imz) = ui(idn) * ql(ivz,i)
ui(ibx) = 0.0d+00
ui(iby) = wl(iby) / slmm
ui(ibz) = wl(ibz) / slmm
ui(ibp) = ql(ibp,i)
ui(ien) = (wl(ien) + sm * pt) / slmm
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! the right speed difference
!
srmm = sr - sm
! conservative variables for the right intermediate state
!
ui(idn) = wr(idn) / srmm
ui(imx) = ui(idn) * sm
ui(imy) = ui(idn) * qr(ivy,i)
ui(imz) = ui(idn) * qr(ivz,i)
ui(ibx) = 0.0d+00
ui(iby) = wr(iby) / srmm
ui(ibz) = wr(ibz) / srmm
ui(ibp) = qr(ibp,i)
ui(ien) = (wr(ien) + sm * pt) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! intermediate flux is constant across the contact discontinuity and all except
! the parallel momentum flux are zero
!
f(idn,i) = 0.0d+00
f(imx,i) = - 0.5d+00 * (wl(imx) + wr(imx))
f(imy,i) = 0.0d+00
f(imz,i) = 0.0d+00
f(ibx,i) = fl(ibx,i)
f(iby,i) = 0.0d+00
f(ibz,i) = 0.0d+00
f(ibp,i) = fl(ibp,i)
f(ien,i) = 0.0d+00
end if ! sm = 0
end if ! Bₓ = 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_hllc
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_HLLD:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the adiabatic HLLD
! method described by Miyoshi & Kusano.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Miyoshi, T. & Kusano, K.,
! "A multi-state HLL approximate Riemann solver for ideal
! magnetohydrodynamics",
! Journal of Computational Physics, 2005, 208, pp. 315-344
!
!===============================================================================
!
subroutine riemann_mhd_adi_hlld(q, f)
! include external procedures
!
use equations , only : nf
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 : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
real(kind=8) :: dn, bx, bp, b2, pt, vy, vz, by, bz, vb
real(kind=8) :: dnl, dnr, cal, car, sml, smr
real(kind=8) :: dv, dvl, dvr
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr, wcl, wcr, ui
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all points
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLLD flux
!
if (sl >= 0.0d+00) then ! sl ≥ 0
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then ! sr ≤ 0
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! the speed of contact discontinuity
!
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
! speed differences
!
slmm = sl - sm
srmm = sr - sm
! left and right state densities
!
dnl = wl(idn) / slmm
dnr = wr(idn) / srmm
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! the total pressure, constant across the contact discontinuity and Alfvén waves
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
! if Alfvén wave is strong enough, apply the full HLLD solver, otherwise
! revert to the HLLC one
!
if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε
! left and right Alfvén speeds
!
cal = sqrt(b2 / dnl)
car = sqrt(b2 / dnr)
sml = sm - cal
smr = sm + car
! calculate division factors (also used to determine degeneracies)
!
dvl = slmm * wl(idn) - b2
dvr = srmm * wr(idn) - b2
! check degeneracy Sl* -> Sl or Sr* -> Sr
!
if (dvl > 0.0d+00 .and. dvr > 0.0d+00) then ! no degeneracy
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! primitive variables for the outer left intermediate state
!
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (smr <= 0.0d+00) then ! sr* ≤ 0
! primitive variables for the outer right intermediate state
!
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sl* < 0 < sr*
! primitive variables for the outer left intermediate state
!
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! primitive variables for the outer right intermediate state
!
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! prepare constant primitive variables of the intermediate states
!
dv = car * dnr + cal * dnl
vy = (wcr(imy) - wcl(imy)) / dv
vz = (wcr(imz) - wcl(imz)) / dv
dv = car + cal
by = (wcr(iby) - wcl(iby)) / dv
bz = (wcr(ibz) - wcl(ibz)) / dv
vb = sm * bx + vy * by + vz * bz
! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! conservative variables for the inmost left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
! the inmost left intermediate flux
!
f(:,i) = sml * ui(:) - wcl(:)
else if (sm < 0.0d+00) then ! sm < 0
! conservative variables for the inmost right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
! the inmost right intermediate flux
!
f(:,i) = smr * ui(:) - wcr(:)
else ! sm = 0
! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
! the flux can be averaged from the Alfvén waves using the simple HLL formula
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
end if ! sm = 0
end if ! sl* < 0 < sr*
else ! one degeneracy
! separate degeneracies
!
if (dvl > 0.0d+00) then ! sr* > sr
! primitive variables for the outer left intermediate state
!
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! choose the correct state depending on the speed signs
!
if (sml >= 0.0d+00) then ! sl* ≥ 0
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else ! sl* < 0
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! primitive variables for the inner left intermediate state
!
dv = srmm * dnr + cal * dnl
vy = (wr(imy) - wcl(imy)) / dv
vz = (wr(imz) - wcl(imz)) / dv
dv = sr - sml
by = (wr(iby) - wcl(iby)) / dv
bz = (wr(ibz) - wcl(ibz)) / dv
vb = sm * bx + vy * by + vz * bz
! conservative variables for the inner left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
if (sm >= 0.0d+00) then ! sm ≥ 0
! the inner left intermediate flux
!
f(:,i) = sml * ui(:) - wcl(:)
else ! sm < 0
! vector of the left-going Alfvén wave
!
wcr(:) = (sm - sml) * ui(:) + wcl(:)
! calculate the average flux over the right inner intermediate state
!
f(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm
end if ! sm < 0
end if ! sl* < 0
else if (dvr > 0.0d+00) then ! sl* < sl
! primitive variables for the outer right intermediate state
!
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
vb = sm * bx + vy * by + vz * bz
! conservative variables for the outer right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! choose the correct state depending on the speed signs
!
if (smr <= 0.0d+00) then ! sr* ≤ 0
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sr* > 0
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! primitive variables for the inner right intermediate state
!
dv = slmm * dnl - car * dnr
vy = (wl(imy) - wcr(imy)) / dv
vz = (wl(imz) - wcr(imz)) / dv
dv = sl - smr
by = (wl(iby) - wcr(iby)) / dv
bz = (wl(ibz) - wcr(ibz)) / dv
vb = sm * bx + vy * by + vz * bz
! conservative variables for the inner left intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
! choose the correct state depending on the sign of contact discontinuity
! advection speed
!
if (sm <= 0.0d+00) then ! sm ≤ 0
! the inner right intermediate flux
!
f(:,i) = smr * ui(:) - wcr(:)
else ! sm > 0
! vector of the right-going Alfvén wave
!
wcl(:) = (sm - smr) * ui(:) + wcr(:)
! calculate the average flux over the left inner intermediate state
!
f(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm
end if ! sm > 0
end if ! sr* > 0
else ! sl* < sl & sr* > sr
! so far we revert to HLL flux in the case of degeneracies
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl* < sl & sr* > sr
end if ! one degeneracy
else if (b2 > 0.0d+00) then ! Bₓ² > 0
! constant intermediate state tangential components of velocity and
! magnetic field
!
vy = (wr(imy) - wl(imy)) / dn
vz = (wr(imz) - wl(imz)) / dn
by = (wr(iby) - wl(iby)) / srml
bz = (wr(ibz) - wl(ibz)) / srml
! scalar product of velocity and magnetic field for the intermediate states
!
vb = sm * bx + vy * by + vz * bz
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! conservative variables for the left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = dnl * vy
ui(imz) = dnl * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! conservative variables for the right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = dnr * vy
ui(imz) = dnr * vz
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
else ! Bₓ² = 0
! separate intermediate states depending on the sign of the advection speed
!
if (sm > 0.0d+00) then ! sm > 0
! conservative variables for the left intermediate state
!
ui(idn) = dnl
ui(imx) = dnl * sm
ui(imy) = wl(imy) / slmm
ui(imz) = wl(imz) / slmm
ui(ibx) = 0.0d+00
ui(iby) = wl(iby) / slmm
ui(ibz) = wl(ibz) / slmm
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt) / slmm
! the left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
! conservative variables for the right intermediate state
!
ui(idn) = dnr
ui(imx) = dnr * sm
ui(imy) = wr(imy) / srmm
ui(imz) = wr(imz) / srmm
ui(ibx) = 0.0d+00
ui(iby) = wr(iby) / srmm
ui(ibz) = wr(ibz) / srmm
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt) / srmm
! the right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
end if ! Bₓ² = 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_hlld
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! q - 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(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : imx, imy, imz, ien
use equations , only : cmax
use equations , only : prim2cons, fluxspeed, eigensystem_roe
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, p, i
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: bx, bp
real(kind=8) :: pml, pmr
real(kind=8) :: xx, yy
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate corresponding conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at the states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:))
! 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(nf) <= 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(:) = 0.0d+00
do p = 1, nf
al(:) = al(:) + li(p,:) * (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, nf
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 = nf, 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
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_roe
!
!===============================================================================
!
!***** ADIABATIC SPECIAL RELATIVITY HYDRODYNAMICS *****
!
!===============================================================================
!
! subroutine UPDATE_FLUX_SRHD_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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_srhd_adi(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
! local variables are not implicit by default
!
implicit none
! input arguments
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1, l, p
real(kind=8) :: vm
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn) :: qq
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1) :: qq
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! apply reconstruction on variables using 4-vector or 3-vector
!
if (states_4vec) then
! convert velocities to four-velocities for physical reconstruction
!
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
qq(idn,i,j,k) = q(idn,i,j,k) / vm
qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
qq(ipr,i,j,k) = q(ipr,i,j,k)
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! convert state four-velocities back to velocities
!
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
do l = 1, 2
do p = 1, NDIMS
vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
end do ! p = 1, ndims
end do ! l = 1, 2
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
else
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
end if
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
f(1,ien,:,j,k) = fi(ien,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
f(2,ien,i,:,k) = fi(ien,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
f(3,ien,i,j,:) = fi(ien,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#ifdef PROFILE
! stop accounting time for flux update
!
call stop_timer(imf)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_flux_srhd_adi
!
!===============================================================================
!
! subroutine RIEMANN_SRHD_ADI_HLL:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! q - 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_srhd_adi_hll(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : ivx
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all position
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate the inverse of speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for the Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_srhd_adi_hll
!
!===============================================================================
!
! subroutine RIEMANN_SRHD_ADI_HLLC:
! --------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
! by Mignone & Bodo.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Mignone, A. & Bodo, G.
! "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
! Monthly Notices of the Royal Astronomical Society,
! 2005, Volume 364, Pages 126-136
!
!===============================================================================
!
subroutine riemann_srhd_adi_hllc(q, f)
! include external procedures
!
use algebra , only : quadratic
use equations , only : nf
use equations , only : ivx, idn, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i, nr
real(kind=8) :: sl, sr, srml, sm
real(kind=8) :: pr, dv, fc
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: uh, us, fh, wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
real(kind=8), dimension(3) :: a
real(kind=8), dimension(2) :: x
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all position
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate the inverse of speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
! correct the energy waves
!
wl(ien) = wl(ien) + wl(idn)
wr(ien) = wr(ien) + wr(idn)
! prepare the quadratic coefficients (eq. 18 in [1])
!
a(1) = uh(imx)
a(2) = - (fh(imx) + uh(ien) + uh(idn))
a(3) = fh(ien) + fh(idn)
! solve the quadratic equation
!
nr = quadratic(a(1:3), x(1:2))
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
!
sm = x(1)
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(:,i) = fh(:)
else
! calculate total pressure (eq. 17 in [1])
!
pr = fh(imx) - (fh(ien) + fh(idn)) * sm
! if the pressure is negative, use the HLL flux
!
if (pr <= 0.0d+00) then
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
if (sm > 0.0d+00) then
! calculate the conserved variable vector (eqs. 16 in [1])
!
dv = sl - sm
us(idn) = wl(idn) / dv
us(imy) = wl(imy) / dv
us(imz) = wl(imz) / dv
us(ien) = (wl(ien) + pr * sm) / dv
us(imx) = (us(ien) + pr) * sm
us(ien) = us(ien) - us(idn)
! calculate the flux (eq. 14 in [1])
!
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
! calculate the conserved variable vector (eqs. 16 in [1])
!
dv = sr - sm
us(idn) = wr(idn) / dv
us(imy) = wr(imy) / dv
us(imz) = wr(imz) / dv
us(ien) = (wr(ien) + pr * sm) / dv
us(imx) = (us(ien) + pr) * sm
us(ien) = us(ien) - us(idn)
! calculate the flux (eq. 14 in [1])
!
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity and all fluxes
! except the parallel momentum one are zero
!
f(idn,i) = 0.0d+00
f(imx,i) = pr
f(imy,i) = 0.0d+00
f(imz,i) = 0.0d+00
f(ien,i) = 0.0d+00
end if ! sm == 0
end if ! p* < 0
end if ! sl < sm < sr
end if ! nr < 1
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for the Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_srhd_adi_hllc
!
!===============================================================================
!
!***** ADIABATIC SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS *****
!
!===============================================================================
!
! subroutine UPDATE_FLUX_SRMHD_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:
!
! dx - the spatial step;
! q - the array of primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux_srmhd_adi(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : nf
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
!
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
! local variables
!
integer :: i, j, k = 1, l, p
real(kind=8) :: vm
! local temporary arrays
!
#if NDIMS == 3
real(kind=8), dimension(nf,nn,nn,nn) :: qq
real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,nn, 1) :: qq
real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
#endif /* NDIMS == 3 */
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for flux update
!
call start_timer(imf)
#endif /* PROFILE */
! initialize fluxes
!
f(:,:,:,:,:) = 0.0d+00
! apply reconstruction on variables using 4-vector or 3-vector
!
if (states_4vec) then
! convert velocities to four-velocities for physical reconstruction
!
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
qq(idn,i,j,k) = q(idn,i,j,k) / vm
qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
qq(ipr,i,j,k) = q(ipr,i,j,k)
qq(ibx,i,j,k) = q(ibx,i,j,k)
qq(iby,i,j,k) = q(iby,i,j,k)
qq(ibz,i,j,k) = q(ibz,i,j,k)
qq(ibp,i,j,k) = q(ibp,i,j,k)
qq(ipr,i,j,k) = q(ipr,i,j,k)
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
! convert state four-velocities back to velocities
!
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
do l = 1, 2
do p = 1, NDIMS
vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
end do ! p = 1, ndims
end do ! l = 1, 2
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
else
! reconstruct interfaces
!
call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
end if
! calculate the flux along the X-direction
!
#if NDIMS == 3
do k = nbl, neu
#endif /* NDIMS == 3 */
do j = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(1,idn,:,j,k) = fi(idn,:)
f(1,imx,:,j,k) = fi(imx,:)
f(1,imy,:,j,k) = fi(imy,:)
f(1,imz,:,j,k) = fi(imz,:)
f(1,ibx,:,j,k) = fi(ibx,:)
f(1,iby,:,j,k) = fi(iby,:)
f(1,ibz,:,j,k) = fi(ibz,:)
f(1,ibp,:,j,k) = fi(ibp,:)
f(1,ien,:,j,k) = fi(ien,:)
end do ! j = nbl, neu
! calculate the flux along the Y direction
!
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(2,idn,i,:,k) = fi(idn,:)
f(2,imx,i,:,k) = fi(imz,:)
f(2,imy,i,:,k) = fi(imx,:)
f(2,imz,i,:,k) = fi(imy,:)
f(2,ibx,i,:,k) = fi(ibz,:)
f(2,iby,i,:,k) = fi(ibx,:)
f(2,ibz,i,:,k) = fi(iby,:)
f(2,ibp,i,:,k) = fi(ibp,:)
f(2,ien,i,:,k) = fi(ien,:)
end do ! i = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
#if NDIMS == 3
! calculate the flux along the Z direction
!
do j = nbl, neu
do i = nbl, neu
! copy directional variable vectors to pass to the one dimensional solver
!
qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call riemann(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
f(3,idn,i,j,:) = fi(idn,:)
f(3,imx,i,j,:) = fi(imy,:)
f(3,imy,i,j,:) = fi(imz,:)
f(3,imz,i,j,:) = fi(imx,:)
f(3,ibx,i,j,:) = fi(iby,:)
f(3,iby,i,j,:) = fi(ibz,:)
f(3,ibz,i,j,:) = fi(ibx,:)
f(3,ibp,i,j,:) = fi(ibp,:)
f(3,ien,i,j,:) = fi(ien,:)
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
#ifdef PROFILE
! stop accounting time for flux update
!
call stop_timer(imf)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_flux_srmhd_adi
!
!===============================================================================
!
! subroutine RIEMANN_SRMHD_ADI_HLL:
! --------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! q - 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_srmhd_adi_hll(q, f)
! include external procedures
!
use equations , only : nf
use equations , only : ivx, ibx, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i
real(kind=8) :: sl, sr, srml
real(kind=8) :: bx, bp
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all position
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate the inverse of speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for the Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_srmhd_adi_hll
!
!===============================================================================
!
! subroutine RIEMANN_SRMHD_ADI_HLLC:
! ---------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
! by Mignone & Bodo.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
! References:
!
! [1] Mignone, A. & Bodo, G.
! "An HLLC Riemann solver for relativistic flows - II.
! Magnetohydrodynamics",
! Monthly Notices of the Royal Astronomical Society,
! 2006, Volume 368, Pages 1040-1054
!
!===============================================================================
!
subroutine riemann_srmhd_adi_hllc(q, f)
! include external procedures
!
use algebra , only : quadratic
use equations , only : nf
use equations , only : ivx, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
use equations , only : cmax
use equations , only : prim2cons, fluxspeed
use interpolations, only : order
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
! local variables
!
integer :: n, i, nr
real(kind=8) :: sl, sr, srml, sm
real(kind=8) :: bx, bp, bx2
real(kind=8) :: pt, dv, fc, uu, ff, uf
real(kind=8) :: vv, vb, gi
! local arrays to store the states
!
real(kind=8), dimension(nf,size(q,2)) :: ql, qr, ul, ur, fl, fr, f2, f4
real(kind=8), dimension(nf) :: uh, us, fh, wl, wr
real(kind=8), dimension(2,size(q,2)) :: cl, cr
real(kind=8), dimension(3) :: a, vs
real(kind=8), dimension(2) :: x
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! get the length of vector
!
n = size(q, 2)
! copy state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! 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
! calculate the conserved variables of the left and right states
!
call prim2cons(ql(:,:), ul(:,:))
call prim2cons(qr(:,:), ur(:,:))
! calculate the physical fluxes and speeds at both states
!
call fluxspeed(ql(:,:), ul(:,:), fl(:,:), cl(:,:))
call fluxspeed(qr(:,:), ur(:,:), fr(:,:), cr(:,:))
! iterate over all position
!
do i = 1, n
! estimate the minimum and maximum speeds
!
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
! calculate the HLL flux
!
if (sl >= 0.0d+00) then
f(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
f(:,i) = fr(:,i)
else ! sl < 0 < sr
! calculate the inverse of speed difference
!
srml = sr - sl
! calculate vectors of the left and right-going waves
!
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! calculate fluxes for the intermediate state
!
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
! correct the energy waves
!
wl(ien) = wl(ien) + wl(idn)
wr(ien) = wr(ien) + wr(idn)
! calculate Bₓ²
!
bx2 = ql(ibx,i) * qr(ibx,i)
! calculate the contact dicontinuity speed solving eq. (41)
!
if (bx2 >= 1.0d-16) then
! prepare the quadratic coefficients, (eq. 42 in [1])
!
uu = sum(uh(iby:ibz) * uh(iby:ibz))
ff = sum(fh(iby:ibz) * fh(iby:ibz))
uf = sum(uh(iby:ibz) * fh(iby:ibz))
a(1) = uh(imx) - uf
a(2) = uu + ff - (fh(imx) + uh(ien) + uh(idn))
a(3) = fh(ien) + fh(idn) - uf
! solve the quadratic equation
!
nr = quadratic(a(1:3), x(1:2))
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
!
sm = x(1)
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(:,i) = fh(:)
else
! substitute magnetic field components from the HLL state (eq. 37 in [1])
!
us(ibx) = ql(ibx,i)
us(iby) = uh(iby)
us(ibz) = uh(ibz)
us(ibp) = ql(ibp,i)
! calculate velocity components without Bₓ normalization (eq. 38 in [1])
!
vs(1) = sm
vs(2) = (us(iby) * sm - fh(iby)) / us(ibx)
vs(3) = (us(ibz) * sm - fh(ibz)) / us(ibx)
! calculate v² and v.B
!
vv = sum(vs(1:3) * vs( 1: 3))
vb = sum(vs(1:3) * us(ibx:ibz))
! calculate inverse gamma
!
gi = 1.0d+00 - vv
! calculate total pressure (eq. 40 in [1])
!
pt = fh(imx) + gi * bx2 - (fh(ien) + fh(idn) - us(ibx) * vb) * sm
! if the pressure is negative, use the HLL flux
!
if (pt <= 0.0d+00) then
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
if (sm > 0.0d+00) then
! calculate the conserved variable vector (eqs. 43-46 in [1])
!
dv = sl - sm
fc = (sl - ql(ivx,i)) / dv
us(idn) = fc * ul(idn,i)
us(imy) = (sl * ul(imy,i) - fl(imy,i) &
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
us(imz) = (sl * ul(imz,i) - fl(imz,i) &
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
us(ien) = (sl * (ul(ien,i) + ul(idn,i)) &
- (fl(ien,i) + fl(idn,i)) &
+ pt * sm - us(ibx) * vb) / dv - us(idn)
! calculate normal component of momentum
!
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
! calculate the flux (eq. 14 in [1])
!
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
! calculate the conserved variable vector (eqs. 43-46 in [1])
!
dv = sr - sm
fc = (sr - qr(ivx,i)) / dv
us(idn) = fc * ur(idn,i)
us(imy) = (sr * ur(imy,i) - fr(imy,i) &
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
us(imz) = (sr * ur(imz,i) - fr(imz,i) &
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
us(ien) = (sr * (ur(ien,i) + ur(idn,i)) &
- (fr(ien,i) + fr(idn,i)) &
+ pt * sm - us(ibx) * vb) / dv - us(idn)
! calculate normal component of momentum
!
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
! calculate the flux (eq. 14 in [1])
!
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
f(:,i) = fh(:)
end if ! sm == 0
end if ! p* < 0
end if ! sl < sm < sr
end if ! nr < 1
else ! Bx² > 0
! prepare the quadratic coefficients (eq. 47 in [1])
!
a(1) = uh(imx)
a(2) = - (fh(imx) + uh(ien) + uh(idn))
a(3) = fh(ien) + fh(idn)
! solve the quadratic equation
!
nr = quadratic(a(1:3), x(1:2))
! if Δ < 0, just use the HLL flux
!
if (nr < 1) then
f(:,i) = fh(:)
else
! get the contact dicontinuity speed
!
sm = x(1)
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
!
if ((sm <= sl) .or. (sm >= sr)) then
f(:,i) = fh(:)
else
! calculate total pressure (eq. 48 in [1])
!
pt = fh(imx) - (fh(ien) + fh(idn)) * sm
! if the pressure is negative, use the HLL flux
!
if (pt <= 0.0d+00) then
f(:,i) = fh(:)
else
! depending in the sign of the contact dicontinuity speed, calculate the proper
! state and corresponding flux
!
if (sm > 0.0d+00) then
! calculate the conserved variable vector (eqs. 43-46 in [1])
!
dv = sl - sm
us(idn) = wl(idn) / dv
us(imy) = wl(imy) / dv
us(imz) = wl(imz) / dv
us(ien) = (wl(ien) + pt * sm) / dv
us(imx) = (us(ien) + pt) * sm
us(ien) = us(ien) - us(idn)
us(ibx) = 0.0d+00
us(iby) = wl(iby) / dv
us(ibz) = wl(ibz) / dv
us(ibp) = ql(ibp,i)
! calculate the flux (eq. 27 in [1])
!
f(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
else if (sm < 0.0d+00) then
! calculate the conserved variable vector (eqs. 43-46 in [1])
!
dv = sr - sm
us(idn) = wr(idn) / dv
us(imy) = wr(imy) / dv
us(imz) = wr(imz) / dv
us(ien) = (wr(ien) + pt * sm) / dv
us(imx) = (us(ien) + pt) * sm
us(ien) = us(ien) - us(idn)
us(ibx) = 0.0d+00
us(iby) = wr(iby) / dv
us(ibz) = wr(ibz) / dv
us(ibp) = qr(ibp,i)
! calculate the flux (eq. 27 in [1])
!
f(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
f(:,i) = fh(:)
end if ! sm == 0
end if ! p* < 0
end if ! sl < sm < sr
end if ! nr < 1
end if ! Bx² == 0 or > 0
end if ! sl < 0 < sr
end do ! i = 1, n
! calculate and add high-order corrections, if required
!
if (high_order_flux .and. order > 2) then
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
f2(:,1 ) = 0.0d+00
f2(:, n ) = 0.0d+00
if (order > 4) then
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
f4(:,1 ) = 0.0d+00
f4(:,2 ) = 0.0d+00
f4(:, n-1) = 0.0d+00
f4(:, n ) = 0.0d+00
f(:,:) = f(:,:) + f4(:,:)
end if
f(:,:) = f(:,:) + f2(:,:)
end if
#ifdef PROFILE
! stop accounting time for the Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine riemann_srmhd_adi_hllc
!===============================================================================
!
end module schemes
Reference in New Issue View Git Blame Copy Permalink