!!******************************************************************************
!!
!! 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-2021 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
implicit none
! 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_solver_iface(ql, qr, ul, ur, fl, fr, cl, cr, f)
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr
real(kind=8), dimension(:,:), intent(in) :: cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
end subroutine
end interface
! pointer to the Riemann solver
!
procedure(riemann_solver_iface), pointer, save :: riemann_solver => 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)
use equations , only : eqsys, eos
use parameters, only : get_parameter
implicit none
logical, intent(in) :: verbose
integer, intent(out) :: status
character(len=64) :: solver = "HLL"
character(len=64) :: statev = "primitive"
character(len=64) :: flux = "off"
!-------------------------------------------------------------------------------
!
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")
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann_solver => riemann_hll
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann_solver => 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")
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann_solver => riemann_hll
case("hllc", "HLLC")
! set the solver name
!
name_sol = "HLLC"
! set pointers to subroutines
!
riemann_solver => riemann_hd_hllc
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann_solver => 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")
! select the Riemann solver
!
select case(trim(solver))
case("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann_solver => riemann_hll
case ("hlld", "HLLD")
! set the solver name
!
name_sol = "HLLD"
! set the solver pointer
!
riemann_solver => riemann_mhd_iso_hlld
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann_solver => 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'" // &
", 'roe'."
end if
status = 1
end select
case("adi", "ADI", "adiabatic", "ADIABATIC")
! select the Riemann solver
!
select case(trim(solver))
case ("hll", "HLL")
! set the solver name
!
name_sol = "HLL"
! set pointers to subroutines
!
riemann_solver => riemann_hll
case ("hllc", "HLLC")
! set the solver name
!
name_sol = "HLLC"
! set the solver pointer
!
riemann_solver => riemann_mhd_hllc
case ("hlld", "HLLD")
! set the solver name
!
name_sol = "HLLD"
! set the solver pointer
!
riemann_solver => riemann_mhd_adi_hlld
case("roe", "ROE")
! set the solver name
!
name_sol = "ROE"
! set pointers to subroutines
!
riemann_solver => 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")
! 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_solver => riemann_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_solver => riemann_srhd_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")
! 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_solver => riemann_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_solver => riemann_srmhd_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
!-------------------------------------------------------------------------------
!
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)
implicit none
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
status = 0
nullify(riemann_solver)
!-------------------------------------------------------------------------------
!
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)
use helpers, only : print_section, print_parameter
implicit none
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
!
!===============================================================================
!
! subroutine UPDATE_FLUX:
! ----------------------
!
! Subroutine solves the Riemann problem along each direction and calculates
! the numerical fluxes, which are used later to calculate the conserved
! variable increments.
!
! Arguments:
!
! dx - the spatial step;
! q - the array of primitive variables;
! s - the state arrays for primitive variables;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux(dx, q, s, f)
use coordinates, only : nn => bcells, nbl, neu
use equations , only : relativistic
use equations , only : nf, ivars, ivx, ivz
implicit none
real(kind=8), dimension(:) , intent(in) :: dx
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
real(kind=8), dimension(:,:,:,:,:,:), intent(inout) :: s
real(kind=8), dimension(:,:,:,:,:) , intent(out) :: f
integer :: n, i, j, k = 1, l, p
real(kind=8) :: vm
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!-------------------------------------------------------------------------------
!
! in the relativistic case, apply the reconstruction on variables
! using the four-velocities if requested
!
if (relativistic .and. states_4vec) then
f(:,:,:,:,1) = q(:,:,:,:)
#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))
f(ivx:ivz,i,j,k,1) = f(ivx:ivz,i,j,k,1) / vm
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
call reconstruct_interfaces(dx(:), f(:,:,:,:,1), s(:,:,:,:,:,:))
! convert the states' 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(s(ivx:ivz,i,j,k,l,p)**2))
s(ivx:ivz,i,j,k,l,p) = s(ivx: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
call reconstruct_interfaces(dx(:), q(:,:,:,:), s(:,:,:,:,:,:))
end if
f(:,:,:,:,:) = 0.0d+00
! 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
!
do n = 1, nf
qi(n,:,1:2) = s(ivars(1,n),:,j,k,1:2,1)
end do
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call numerical_flux(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
do n = 1, nf
f(ivars(1,n),:,j,k,1) = fi(n,:)
end do
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
!
do n = 1, nf
qi(n,:,1:2) = s(ivars(2,n),i,:,k,1:2,2)
end do
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call numerical_flux(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
do n = 1, nf
f(ivars(2,n),i,:,k,2) = fi(n,:)
end do
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
!
do n = 1, nf
qi(n,:,1:2) = s(ivars(3,n),i,j,:,1:2,3)
end do
! call one dimensional Riemann solver in order to obtain numerical fluxes
!
call numerical_flux(qi(:,:,:), fi(:,:))
! update the array of fluxes
!
do n = 1, nf
f(ivars(3,n),i,j,:,3) = fi(n,:)
end do
end do ! i = nbl, neu
end do ! j = nbl, neu
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine update_flux
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine RECONSTRUCT_INTERFACES:
! ---------------------------------
!
! Subroutine reconstructs the Riemann states along all directions.
!
! 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
use equations , only : positive
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
!
integer :: p
!
!-------------------------------------------------------------------------------
!
! interpolate interfaces for all flux variables
!
do p = 1, nf
call interfaces(positive(p), dx(1:NDIMS), q (p,:,:,:), &
qi(p,:,:,:,1:2,1:NDIMS))
end do
!-------------------------------------------------------------------------------
!
end subroutine reconstruct_interfaces
!
!===============================================================================
!
! subroutine NUMERICAL_FLUX:
! -------------------------
!
! Subroutine prepares Riemann states and calls the selected Riemann solver
! in order to get the numerical flux. If requested, the resulting flux
! is corrected by higher order correction terms.
!
! Arguments:
!
! q - the array of primitive variables at the Riemann states;
! f - the output array of fluxes;
!
!===============================================================================
!
subroutine numerical_flux(q, f)
use equations, only : magnetized, cmax, ibx, ibp
use equations, only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:,:), intent(inout) :: q
real(kind=8), dimension(:,:) , intent(out) :: f
integer :: i
real(kind=8) :: bx, bp
real(kind=8), dimension(size(q,1),size(q,2)) :: ql, qr, ul, ur, fl, fr
real(kind=8), dimension( 2,size(q,2)) :: cl, cr
!-------------------------------------------------------------------------------
!
! copy the state vectors
!
ql(:,:) = q(:,:,1)
qr(:,:) = q(:,:,2)
! obtain the state values for Bx and Psi for the GLM-MHD equations
!
if (magnetized) then
do i = 1, size(q, 2)
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
end if
! 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(:,:))
! get the HLL fluxes
!
call riemann_solver(ql(:,:), qr(:,:), ul(:,:), ur(:,:), &
fl(:,:), fr(:,:), cl(:,:), cr(:,:), f(:,:))
! higher order flux corrections
!
call higher_order_flux_correction(f)
!-------------------------------------------------------------------------------
!
end subroutine numerical_flux
!
!===============================================================================
!
! subroutine RIEMANN_HLL:
! ----------------------
!
! Subroutine solves one dimensional general Riemann problem using
! the Harten-Lax-van Leer (HLL) method.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_hll(ql, qr, ul, ur, fl, fr, cl, cr, f)
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
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(size(ql,1)) :: wl, wr
!
!-------------------------------------------------------------------------------
!
! get the length of vector
!
n = size(ql, 2)
! 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(:) / srml - sr * wl(:) / srml
end if ! sl < 0 < sr
end do ! i = 1, n
!-------------------------------------------------------------------------------
!
end subroutine riemann_hll
!
!===============================================================================
!
! subroutine RIEMANN_HD_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:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_hd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use equations, only : idn, ivy, ivz, imx, imy, imz, ien
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i
real(kind=8) :: dn, pr
real(kind=8) :: sl, sr, sm, slmm, srmm
! local arrays to store the states
!
real(kind=8), dimension(size(ql,1)) :: wl, wr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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) / dn - 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
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_MHD_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:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_mhd_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use equations, only : idn, ivy, ivz, imx, imy, imz, ien
use equations, only : ibx, iby, ibz, ibp
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i
real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
! local arrays to store the states
!
real(kind=8), dimension(size(ql,1)) :: wl, wr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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
! Bx and square of Bx, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! 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
! the total pressure, constant across the contact discontinuity
!
pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2
! separate the cases when Bₓ = 0 and Bₓ ≠ 0
!
if (b2 > 0.0d+00) then
! 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(:) / srml - 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
! 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
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sm = 0
end if ! Bₓ = 0
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_SRHD_HLLC:
! ----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
! by Mignone & Bodo.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use algebra , only : quadratic
use equations, only : idn, imx, imy, imz, ien
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i, nr
real(kind=8) :: sl, sr, srml, sm
real(kind=8) :: pr, dv
! local arrays to store the states
!
real(kind=8), dimension(size(ql,1)) :: wl, wr, uh, us, fh
real(kind=8), dimension(3) :: a
real(kind=8), dimension(2) :: x
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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
!-------------------------------------------------------------------------------
!
end subroutine riemann_srhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_SRMHD_HLLC:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
! by Mignone & Bodo.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_hllc(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use algebra , only : quadratic
use equations, only : idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i, nr
real(kind=8) :: sl, sr, srml, sm
real(kind=8) :: 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(size(ql,1)) :: wl, wr, uh, us, fh
real(kind=8), dimension(3) :: a, vs
real(kind=8), dimension(2) :: x
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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 of Lorentz factor
!
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
!-------------------------------------------------------------------------------
!
end subroutine riemann_srmhd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_ISO_MHD_HLLD:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the isothermal HLLD
! method described by Mignone.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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_hlld(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use equations, only : idn, imx, imy, imz, ibx, iby, ibz, ibp
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm, sdif
real(kind=8) :: bx, b2, dn, dvl, dvr, ca
! local arrays to store the states
!
real(kind=8), dimension(size(ql,1)) :: wl, wr, wcl, wcr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! 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
! speed differences
!
slmm = sl - sm
srmm = sr - sm
! calculate density of the intermediate states
!
dn = dn / srml
! left and right Alfvén speeds
!
ca = abs(bx) / sqrt(dn)
sml = sm - ca
smr = sm + ca
! 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 (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracies
! 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) = (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) = 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) = (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*
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! 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) = (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) = 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) = (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) = 0.5d+00 * (sml * wcr(:) / ca - smr * wcl(:) / ca)
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
end if ! sl* < 0 < sr*
else ! some degeneracies
! separate degeneracies
!
if (sml > sl .and. 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) = (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
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the left-going Alfvén wave
!
wcl(:) = (sml - sl) * ui(:) + wl(:)
! the flux corresponding to the middle intermediate state
!
sdif = sr - sml
f(:,i) = sml * wr(:) / sdif - sr * wcl(:) / sdif
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
end if ! sl* < 0
else if (smr < sr .and. 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) = (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
! separate cases for non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! vector of the right-going Alfvén wave
!
wcr(:) = (smr - sr) * ui(:) + wr(:)
! the flux corresponding to the middle intermediate state
!
sdif = smr - sl
f(:,i) = sl * wcr(:) / sdif - smr * wl(:) / sdif
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
end if ! sr* > 0
else ! sl* < sl & sr* > sr
! both outer states are degenerate so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sl* < sl & sr* > sr
end if ! degeneracies
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_hlld
!
!===============================================================================
!
! subroutine RIEMANN_ADI_MHD_HLLD:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the adiabatic HLLD
! method described by Miyoshi & Kusano.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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(ql, qr, ul, ur, fl, fr, cl, cr, f)
! include external procedures
!
use equations, only : idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
! local variables
!
integer :: i
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
real(kind=8) :: dn, bx, 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(size(ql,1)) :: wl, wr, wcl, wcr, ui
!
!-------------------------------------------------------------------------------
!
! iterate over all points
!
do i = 1, size(ql, 2)
! 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
! square of Bₓ, i.e. Bₓ²
!
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
! 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
! the total pressure, constant across the contact discontinuity and Alfvén waves
!
pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2
! left and right Alfvén speeds
!
cal = abs(bx) / sqrt(dnl)
car = abs(bx) / sqrt(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 (sml > sl .and. smr < sr .and. 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
! 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*
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! 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 = abs(bx) * (sqrt(dnr) + sqrt(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 a simple HLL formula;
!
f(:,i) = sml * wcr(:) / dv - smr * wcl(:) / dv
end if ! sm = 0
else ! Bx = 0
if (sm > 0.0d+00) then
! primitive variables for the outer left intermediate state
!
vy = slmm * wl(imy) / dvl
vz = slmm * wl(imz) / dvl
by = wl(idn) * wl(iby) / dvl
bz = wl(idn) * wl(ibz) / dvl
vb = 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) = 0.0d+00
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt) / slmm
! the outer left intermediate flux
!
f(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then
! primitive variables for the outer right intermediate state
!
vy = ( srmm * wr(imy)) / dvr
vz = ( srmm * wr(imz)) / dvr
by = (wr(idn) * wr(iby)) / dvr
bz = (wr(idn) * wr(ibz)) / dvr
vb = 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) = 0.0d+00
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt) / srmm
! the outer right intermediate flux
!
f(:,i) = sr * ui(:) - wr(:)
else ! Sm = 0
! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
end if ! Bx = 0
end if ! sl* < 0 < sr*
else ! some degeneracies
! separate degeneracies
!
if (sml > sl .and. 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
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! 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(:) / srmm - sr * wcr(:) / srmm
end if ! sm < 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
end if ! sl* < 0
else if (smr < sr .and. 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
! separate cases with non-zero and zero Bx
!
if (b2 > 0.0d+00) then
! 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(:) / slmm - sl * wcl(:) / slmm
end if ! sm > 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
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(:) / srml - sr * wl(:) / srml
end if ! sl* < sl & sr* > sr
end if ! deneneracies
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_hlld
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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(ql, qr, ul, ur, fl, fr, cl, cr, f)
use equations, only : idn, ivx, ivz
use equations, only : eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
integer :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8), dimension(size(ql,1)) :: qi, ci, al
real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
nf = size(ql, 1)
! iterate over all points
!
do i = 1, size(ql, 2)
! 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx: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
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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(ql, qr, ul, ur, fl, fr, cl, cr, f)
use equations, only : idn, ivx, ivz, ipr, ien
use equations, only : eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
integer :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8), dimension(size(ql,1)) :: qi, ci, al
real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
nf = size(ql, 1)
! iterate over all points
!
do i = 1, size(ql, 2)
! 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx: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
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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(ql, qr, ul, ur, fl, fr, cl, cr, f)
use equations, only : idn, ivx, ivz, imx, imy, imz, ibx, iby, ibz, ibp
use equations, only : eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
integer :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
real(kind=8), dimension(size(ql,1)) :: qi, ci, al
real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
nf = size(ql, 1)
! iterate over all points
!
do i = 1, size(ql, 2)
! 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
qi(ibx) = ql(ibx,i)
qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby: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
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the Riemann states;
! ul, ur - the conservative variables of the Riemann states;
! fl, fr - the physical fluxes of the Riemann states;
! f - the Riemann average flux;
!
! 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(ql, qr, ul, ur, fl, fr, cl, cr, f)
use equations, only : idn, ivx, ivz, imx, imy, imz, ipr, ien
use equations, only : ibx, iby, ibz, ibp
use equations, only : eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr, ul, ur, fl, fr, cl, cr
real(kind=8), dimension(:,:), intent(out) :: f
integer :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
real(kind=8) :: pml, pmr
real(kind=8), dimension(size(ql,1)) :: qi, ci, al
real(kind=8), dimension(size(ql,1),size(ql,1)) :: li, ri
!-------------------------------------------------------------------------------
!
! get the number of fluxes
!
nf = size(ql, 1)
! iterate over all points
!
do i = 1, size(ql, 2)
! 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx: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:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby: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
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_roe
!
!===============================================================================
!
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
! ---------------------------------------
!
! Subroutine adds higher order corrections to the numerical fluxes.
!
! Arguments:
!
! f - the vector 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
!===============================================================================
!
end module schemes