gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/schemes.F90
Grzegorz Kowal 7fdac2f584 SCHEMES: Make HD ROE solvers numerically symmetric. 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2020-11-17 17:52:22 -03:00

3540 lines
91 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_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)
! 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")
! 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
#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(riemann_solver)
#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 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;
! f - the array of numerical fluxes;
!
!===============================================================================
!
subroutine update_flux(dx, q, f)
! include external variables
!
use coordinates, only : nn => bcells, nbl, neu
use equations , only : relativistic
use equations , only : nf, ivars, ivx, ivz
! 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 :: n, 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 (relativistic .and. 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( : ,i,j,k) = q( : ,i,j,k)
qq(ivx:ivz,i,j,k) = q(ivx:ivz,i,j,k) / vm
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(ivx:ivz,i,j,k,l,p) = qs(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
! 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 states to directional lines with proper vector component ordering
!
do n = 1, nf
qi(n,:,1:2) = qs(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(1,ivars(1,n),:,j,k) = 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) = qs(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(2,ivars(2,n),i,:,k) = 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) = qs(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(3,ivars(3,n),i,j,:) = fi(n,:)
end do
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
!
!===============================================================================
!
! 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)
! include external procedures
!
use equations, only : magnetized, cmax, ibx, ibp
use equations, only : prim2cons, fluxspeed
! 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 :: i
real(kind=8) :: bx, bp
! local arrays to store the states
!
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
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for Riemann solver
!
call start_timer(imr)
#endif /* PROFILE */
! copy 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)
#ifdef PROFILE
! stop accounting time for Riemann solver
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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) 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) 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) 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) 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) 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) 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)
! include external procedures
!
use equations, only : idn, ivx, ivy, ivz
use equations, only : eigensystem_roe
! 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 :: nf, i, p
real(kind=8) :: sdl, sdr, sds
! local arrays to store the states
!
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)
! include external procedures
!
use equations, only : idn, ivx, ivy, ivz, ipr, ien
use equations, only : eigensystem_roe
! 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 :: nf, i, p
real(kind=8) :: sdl, sdr, sds
! local arrays to store the states
!
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)
! include external procedures
!
use equations, only : idn, ivx, ivy, ivz, imx, imy, imz, ibx, iby, ibz, ibp
use equations, only : eigensystem_roe
! 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 :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
! local arrays to store the states
!
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)
! include external procedures
!
use equations, only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
use equations, only : ibx, iby, ibz, ibp
use equations, only : eigensystem_roe
! 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 :: nf, i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
real(kind=8) :: pml, pmr
! local arrays to store the states
!
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
Reference in New Issue View Git Blame Copy Permalink