4101 lines
111 KiB
Fortran
Raw Normal View History

!!******************************************************************************
!!
!! 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.
2012-08-01 17:25:49 -03:00
!!
!!******************************************************************************
!
module schemes
implicit none
! Rieman solver and state vectors names
!
character(len=255) , save :: name_sol = ""
character(len=255) , save :: name_sts = "primitive"
! KEPES solver
!
logical , save :: kepes = .false.
! 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, fi)
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
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("kepes", "KEPES")
name_sol = "KEPES"
kepes = .true.
riemann_solver => riemann_hd_iso_kepes
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("kepes", "KEPES")
name_sol = "KEPES"
kepes = .true.
riemann_solver => riemann_hd_adi_kepes
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("kepes", "KEPES")
name_sol = "KEPES"
kepes = .true.
riemann_solver => riemann_mhd_iso_kepes
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("kepes", "KEPES")
name_sol = "KEPES"
kepes = .true.
riemann_solver => riemann_mhd_adi_kepes
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, l, p
real(kind=8) :: vm
real(kind=8), dimension(nf,nn,2) :: qi
real(kind=8), dimension(nf,nn) :: fi
!-------------------------------------------------------------------------------
!
#if NDIMS == 2
k = 1
#endif /* NDIMS == 2 */
! 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 coordinates, only : nn => bcells
use equations , only : nf, ibx, ibp
use equations , only : magnetized, cmax
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(nf,nn) :: ql, qr
!-------------------------------------------------------------------------------
!
if (kepes) then
call riemann_solver(q(:,:,1), q(:,:,2), f(:,:))
else
! 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 = 5.0d-01 * ((qr(ibx,i) + ql(ibx,i)) &
- (qr(ibp,i) - ql(ibp,i)) / cmax)
bp = 5.0d-01 * ((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
! get the numerical fluxes
!
call riemann_solver(ql(:,:), qr(:,:), f(:,:))
! higher order flux corrections
!
call higher_order_flux_correction(f)
end if
!-------------------------------------------------------------------------------
!
end subroutine numerical_flux
!
!*******************************************************************************
!
! GENERIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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,
! https://doi.org/10.1137/1025002
!
!===============================================================================
!
subroutine riemann_hll(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: sl, sr
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
wl = sl * ul(:,i) - fl(:,i)
wr = sr * ur(:,i) - fr(:,i)
fi(:,i) = (sl * wr - sr * wl) / (sr - sl)
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hll
!
!*******************************************************************************
!
! ISOTHERMAL HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use coordinates, only : nn => bcells
use equations , only : nf, csnd
use equations , only : idn, ivx, ivy, ivz
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: sdl, sdr, sds, cs, ch
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension(nf ) :: qs, lm, al, df
logical , save :: first = .true.
real(kind=8) , save :: chi
real(kind=8), dimension(4,4), save :: rvec, lvec
!$omp threadprivate(first, chi, rvec, lvec)
!-------------------------------------------------------------------------------
!
if (first) then
chi = 5.0d-01 / csnd
rvec(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
rvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
rvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,2) = [ chi , 0.0d+00, 0.0d+00, - chi ]
lvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
first = .false.
end if
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl)
call fluxspeed(qr, ur, fr)
do i = 1, nn
sdl = sqrt(ql(idn,i))
sdr = sqrt(qr(idn,i))
sds = sdl + sdr
qs(idn) = sdl * sdr
qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
cs = sign(csnd, qs(ivx))
ch = sign(chi , qs(ivx))
lm(1) = qs(ivx) + cs
lm(2) = qs(ivx)
lm(3) = qs(ivx)
lm(4) = qs(ivx) - cs
rvec(1,2) = lm(1)
rvec(4,2) = lm(4)
rvec(1,3) = qs(ivy)
rvec(4,3) = qs(ivy)
rvec(1,4) = qs(ivz)
rvec(4,4) = qs(ivz)
lvec(1,1) = - ch * lm(4)
lvec(2,1) = - qs(ivy)
lvec(3,1) = - qs(ivz)
lvec(4,1) = ch * lm(1)
lvec(1,2) = ch
lvec(4,2) = - ch
al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
df = matmul(al, rvec)
fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_iso_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ISO_KEPES:
! -------------------------------
!
! Subroutine solves one dimensional isothermal hydrodynamic Riemann problem
! using the entropy stable KEPES method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J.,
! "Entropy stable numerical approximations for the isothermal and
! polytropic Euler equations",
! BIT Numerical Mathematics (2020) 60:791824,
! https://doi.org/10.1007/s10543-019-00789-w
!
!===============================================================================
!
subroutine riemann_hd_iso_kepes(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, csnd, csnd2
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: v2l, v2r
real(kind=8) :: dnl, dna, vxa, vya, vza
real(kind=8), dimension(nf) :: v, lm, tm, dm
logical , save :: first = .true.
real(kind=8), dimension(4,4), save :: rm
!$omp threadprivate(first, rm)
!-------------------------------------------------------------------------------
!
if (first) then
rm(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
rm(:,2) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
rm(:,3) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
rm(:,4) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
first = .false.
end if
do i = 1, nn
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
dnl = lmean(ql(idn,i), qr(idn,i))
dna = amean(ql(idn,i), qr(idn,i))
vxa = amean(ql(ivx,i), qr(ivx,i))
vya = amean(ql(ivy,i), qr(ivy,i))
vza = amean(ql(ivz,i), qr(ivz,i))
v(idn) = 5.0d-01 * (v2l - v2r)
if (qr(idn,i) > ql(idn,i)) then
v(idn) = v(idn) + csnd2 * (log(qr(idn,i) / ql(idn,i)))
else if (ql(idn,i) > qr(idn,i)) then
v(idn) = v(idn) - csnd2 * (log(ql(idn,i) / qr(idn,i)))
end if
v(ivx) = qr(ivx,i) - ql(ivx,i)
v(ivy) = qr(ivy,i) - ql(ivy,i)
v(ivz) = qr(ivz,i) - ql(ivz,i)
lm(:) = [ vxa + csnd, vxa, vxa, vxa - csnd ]
tm(1) = 5.0d-01 * dnl / csnd2
tm(2) = dnl
tm(3) = dnl
tm(4) = tm(1)
rm(2:4,1) = [ lm(1), vya, vza ]
rm(2:4,4) = [ lm(4), vya, vza ]
fi(idn,i) = dnl * vxa
fi(imx,i) = fi(idn,i) * vxa + csnd2 * dna
fi(imy,i) = fi(idn,i) * vya
fi(imz,i) = fi(idn,i) * vza
dm = (abs(lm) * tm) * matmul(v, rm)
if (vxa >= 0.0d+00) then
dm = matmul(rm , dm )
else
dm = matmul(rm(:,nf:1:-1), dm(nf:1:-1))
end if
fi(:,i) = fi(:,i) - 5.0d-01 * dm(:)
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_iso_kepes
!
!*******************************************************************************
!
! ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Toro, E. F.,
! "The HLLC Riemann solver",
! Shock Waves, 2019, 29, 8, 1065-1082,
! https://doi.org/10.1007/s00193-019-00912-4
!
!===============================================================================
!
subroutine riemann_hd_hllc(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf
use equations , only : idn, ivy, ivz, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: dn, pr
real(kind=8) :: sl, sr, sm, slmm, srmm
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, ui
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
wl = sl * ul(:,i) - fl(:,i)
wr = sr * ur(:,i) - fr(:,i)
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
if (sm > 0.0d+00) then
slmm = sl - sm
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
fi(:,i) = sl * ui - wl
else if (sm < 0.0d+00) then
srmm = sr - sm
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
fi(:,i) = sr * ui - wr
else
fi(:,i) = (sl * wr - sr * wl) / (sr - sl)
end if
end if
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_hllc
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_ROE:
! -----------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Roe, P. L.
! "Approximate Riemann Solvers, Parameter Vectors, and Difference
! Schemes",
! Journal of Computational Physics, 1981, 43, pp. 357-372
! https://doi.org/10.1016/0021-9991(81)90128-5
! [2] Toro, E. F.,
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
! Springer-Verlag Berlin Heidelberg, 2009
!
!===============================================================================
!
subroutine riemann_hd_adi_roe(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, adiabatic_index
use equations , only : idn, ivx, ivy, ivz, ipr, ien
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: sdl, sdr, sds, vv, vh, c2, cs, vc, fa, fb, fc, f1, f2, f3
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension(nf ) :: qs, lm, al, df
logical , save :: first = .true.
real(kind=8) , save :: adi_m1, adi_m1x
real(kind=8), dimension(5,5), save :: rvec, lvec
!$omp threadprivate(first, adi_m1, adi_m1x, rvec, lvec)
!-------------------------------------------------------------------------------
!
if (first) then
adi_m1 = adiabatic_index - 1.0d+00
adi_m1x = adiabatic_index / adi_m1
rvec(:,1) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
rvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
rvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
rvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
lvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
lvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
first = .false.
end if
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl)
call fluxspeed(qr, ur, fr)
do i = 1, nn
sdl = sqrt(ql(idn,i))
sdr = sqrt(qr(idn,i))
sds = sdl + sdr
qs(idn) = sdl * sdr
qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
qs(ien) = ((ul(ien,i) + ql(ipr,i)) / sdl &
+ (ur(ien,i) + qr(ipr,i)) / sdr) / sds ! enthalpy
vv = sum(qs(ivx:ivz) * qs(ivx:ivz))
vh = 5.0d-01 * vv
c2 = adi_m1 * (qs(ien) - vh)
cs = sign(sqrt(c2), qs(ivx))
fa = adi_m1 / c2
fb = 5.0d-01 * fa
fc = 5.0d-01 / cs
vc = cs * qs(ivx)
f1 = fb * vh
f2 = fc * qs(ivx)
f3 = - fb * qs(ivx)
lm(1) = qs(ivx) + cs
lm(2) = qs(ivx)
lm(3) = qs(ivx)
lm(4) = qs(ivx)
lm(5) = qs(ivx) - cs
rvec(1,2) = lm(1)
rvec(1,3) = qs(ivy)
rvec(1,4) = qs(ivz)
rvec(1,5) = qs(ien) + vc
rvec(2,2) = qs(ivx)
rvec(2,3) = qs(ivy)
rvec(2,4) = qs(ivz)
rvec(2,5) = vh
rvec(3,5) = qs(ivy)
rvec(4,5) = qs(ivz)
rvec(5,2) = lm(5)
rvec(5,3) = qs(ivy)
rvec(5,4) = qs(ivz)
rvec(5,5) = qs(ien) - vc
lvec(1,1) = f1 - f2
lvec(1,2) = f3 + fc
lvec(1,3) = - fb * qs(ivy)
lvec(1,4) = - fb * qs(ivz)
lvec(1,5) = fb
lvec(2,1) = - fa * vh + 1.0d+00
lvec(2,2) = fa * qs(ivx)
lvec(2,3) = fa * qs(ivy)
lvec(2,4) = fa * qs(ivz)
lvec(2,5) = - fa
lvec(3,1) = - qs(ivy)
lvec(4,1) = - qs(ivz)
lvec(5,1) = f1 + f2
lvec(5,2) = f3 - fc
lvec(5,3) = lvec(1,3)
lvec(5,4) = lvec(1,4)
lvec(5,5) = fb
al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
df = matmul(al, rvec)
fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_HD_ADI_KEPES:
! -------------------------------
!
! Subroutine solves one dimensional adiabatic hydrodynamic Riemann problem
! using the entropy stable KEPES method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Winters, A. R., Derigs, D., Gassner, G. J., Walch, S.
! "A uniquely defined entropy stable matrix dissipation operator
! for high Mach number ideal MHD and compressible Euler
! simulations",
! Journal of Computational Physics, 2017, 332, pp. 274-289
!
!===============================================================================
!
subroutine riemann_hd_adi_kepes(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, adiabatic_index
use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
logical , save :: first = .true.
real(kind=8) , save :: adi_m1, adi_m1x
real(kind=8), dimension(5,5), save :: rm
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
integer :: i, p
real(kind=8) :: bl, br, btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra
real(kind=8) :: csn, ent
real(kind=8), dimension(nf) :: v, lm, tm, lv
!-------------------------------------------------------------------------------
!
if (first) then
adi_m1 = adiabatic_index - 1.0d+00
adi_m1x = adi_m1 / adiabatic_index
rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
first = .false.
end if
do i = 1, nn
! various parameters
!
bl = ql(idn,i) / ql(ipr,i)
br = qr(idn,i) / qr(ipr,i)
! the intermediate state averages
!
btl = lmean(bl, br)
bta = amean(bl, br)
dnl = lmean(ql(idn,i), qr(idn,i))
dna = amean(ql(idn,i), qr(idn,i))
vxa = amean(ql(ivx,i), qr(ivx,i))
vya = amean(ql(ivy,i), qr(ivy,i))
vza = amean(ql(ivz,i), qr(ivz,i))
pra = dna / bta
prl = dnl / btl
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
ent = 1.0d+00 / adi_m1x / btl + vva
! the difference in entropy variables between states
!
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
+ (log(br) - log(bl)) / adi_m1 &
- 5.0d-01 * (br * sum(qr(ivx:ivz,i)**2) &
- bl * sum(ql(ivx:ivz,i)**2))
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
v(ipr) = bl - br
! the eigenvalues
!
csn = sqrt(adiabatic_index * pra / dnl)
! the diagonal matrix of eigenvalues
!
lm(1) = vxa + csn
lm(2) = vxa
lm(3) = vxa
lm(4) = vxa
lm(5) = vxa - csn
! the scaling diagonal matrix
!
tm(1) = 5.0d-01 * dnl / adiabatic_index
tm(2) = dnl * adi_m1x
tm(3) = pra
tm(4) = pra
tm(5) = tm(1)
! the average of the right eigenvector matrix
!
rm(2,1) = vxa + csn
rm(3,1) = vya
rm(4,1) = vza
rm(2,2) = vxa
rm(3,2) = vya
rm(4,2) = vza
rm(2,5) = vxa - csn
rm(3,5) = vya
rm(4,5) = vza
rm(5,:) = [ ent + vxa * csn, vva, vya, vza, ent - vxa * csn ]
! KEPEC flux
!
fi(idn,i) = dnl * vxa
fi(imx,i) = fi(idn,i) * vxa + pra
fi(imy,i) = fi(idn,i) * vya
fi(imz,i) = fi(idn,i) * vza
fi(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
! KEPES correction
!
lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
if (vxa >= 0.0d+00) then
do p = 1, nf
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
end do
else
do p = nf, 1, -1
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
end do
end if
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_hd_adi_kepes
!
!*******************************************************************************
!
! RELATIVISTIC ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use algebra , only : quadratic
use coordinates, only : nn => bcells
use equations , only : nf, idn, imx, imy, imz, ien
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i, nr
real(kind=8) :: sl, sr, srml, sm
real(kind=8) :: pr, dv
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
real(kind=8), dimension(3) :: a
real(kind=8), dimension(2) :: x
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
srml = sr - sl
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
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
fi(:,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
fi(:,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
fi(:,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])
!
fi(:,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])
!
fi(:,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
!
fi(idn,i) = 0.0d+00
fi(imx,i) = pr
fi(imy,i) = 0.0d+00
fi(imz,i) = 0.0d+00
fi(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
!
!*******************************************************************************
!
! ISOTHERMAL MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! subroutine RIEMANN_ISO_MHD_HLLD:
! -------------------------------
!
! Subroutine solves one dimensional Riemann problem using the modified
! isothermal HLLD method by Mignone. The difference is in allowing a density
! jumps across the Alfvén waves, due to the inner intermediate states being
! an average of the states between two Alfvén waves, which also include
! two slow magnetosonic waves. Moreover, due to the different left and right
! intermediate state densities the left and right Alfvén speeds might be
! different too.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Mignone, A.,
! "A simple and accurate Riemann solver for isothermal MHD",
! Journal of Computational Physics, 2007, 225, pp. 1427-1441,
! https://doi.org/10.1016/j.jcp.2007.01.033
!
!===============================================================================
!
subroutine riemann_mhd_iso_hlld(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, idn, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
real(kind=8) :: bx, b2, dvl, dvr
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, uil, uir
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
! constant advection speed across the intermediate states
!
sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn))
! speed differences
!
srml = sr - sl
slmm = sl - sm
srmm = sr - sm
! set variables which do not change across the Riemann fan
!
uil(idn) = wl(idn) / slmm
uir(idn) = wr(idn) / srmm
uil(imx) = uil(idn) * sm
uir(imx) = uir(idn) * sm
uil(ibx) = bx
uir(ibx) = bx
uil(ibp) = ql(ibp,i)
uir(ibp) = qr(ibp,i)
! Alfvén speeds
!
sml = sm - abs(bx) / sqrt(uil(idn))
smr = sm + abs(bx) / sqrt(uir(idn))
! 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
if (sml > 0.0d+00) then ! sl* > 0
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
fi(:,i) = sl * uil(:) - wl(:)
else if (smr < 0.0d+00) then ! sr* < 0
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
fi(:,i) = sr * uir(:) - wr(:)
else ! sl* <= 0 <= sr*
if (smr > sml) then
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
wcl(:) = (sml - sl) * uil(:) + wl(:)
wcr(:) = (smr - sr) * uir(:) + wr(:)
dvl = smr - sml
fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dvl
else ! Bx = 0 -> Sₘ = 0
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sl* <= 0 <= sr*
else if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
if (sml > 0.0d+00) then ! sl* > 0
fi(:,i) = sl * uil(:) - wl(:)
else if (sr > sml) then
wcl(:) = (sml - sl) * uil(:) + wl(:)
dvl = sr - sml
fi(:,i) = (sml * wr(:) - sr * wcl(:)) / dvl
else ! Bx = 0 -> Sₘ = 0
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sl* < 0
else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
if (smr < 0.0d+00) then ! sr* < 0
fi(:,i) = sr * uir(:) - wr(:)
else if (smr > sl) then
wcr(:) = (smr - sr) * uir(:) + wr(:)
dvr = smr - sl
fi(:,i) = (sl * wcr(:) - smr * wl(:)) / dvr
else ! Bx = 0 -> Sₘ = 0
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sr* > 0
else ! sl* < sl & sr* > sr
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_hlld
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ISO_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use coordinates, only : nn => bcells
use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed, eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
! 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
fi(:,i) = fl(:,i)
else if (ci(nf) <= 0.0d+00) then
fi(:,i) = fr(:,i)
else
! prepare fluxes which do not change across the states
!
fi(ibx,i) = fl(ibx,i)
fi(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
!
fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
fi(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))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
end do
else
do p = nf, 1, -1
xx = - 0.5d+00 * al(p) * abs(ci(p))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(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_ISO_KEPES:
! --------------------------------
!
! Subroutine solves one dimensional isothermal magnetohydrodynamic
! Riemann problem using the entropy stable KEPES method. The method is
! a modification of the method described in [1] for the isothermal equation
! of state.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
! field divergence diminishing ideal magnetohydrodynamics equations",
! Journal of Computational Physics, 2018, 364, pp. 420-467,
! https://doi.org/10.1016/j.jcp.2018.03.002
!
!===============================================================================
!
subroutine riemann_mhd_iso_kepes(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, csnd, csnd2, ch => cglm
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
ibx, iby, ibz, ibp
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
logical , save :: first = .true.
real(kind=8), dimension(8,8), save :: rm
!$omp threadprivate(first, rm)
integer :: i, p
real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
real(kind=8) :: dnl, sqd, bp
real(kind=8) :: bb2, bp2, b1, b2, b3, bb
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
real(kind=8) :: alf2, als2, alf, als
real(kind=8) :: f1, f2, f3, sgnb1, v2l, v2r, b2l, b2r
real(kind=8), dimension(nf) :: v, lm, tm, lv
!-------------------------------------------------------------------------------
!
if (first) then
rm(:,:) = 0.0d+00
rm(5,4) = 1.0d+00
rm(8,4) = 1.0d+00
rm(5,5) = 1.0d+00
rm(8,5) = - 1.0d+00
first = .false.
end if
do i = 1, nn
! various parameters
!
v2l = sum(ql(ivx:ivz,i)**2)
v2r = sum(qr(ivx:ivz,i)**2)
b2l = sum(ql(ibx:ibz,i)**2)
b2r = sum(qr(ibx:ibz,i)**2)
! the difference in entropy variables between the states
!
v(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) &
- 5.0d-01 * (v2r - v2l)
v(ivx) = qr(ivx,i) - ql(ivx,i)
v(ivy) = qr(ivy,i) - ql(ivy,i)
v(ivz) = qr(ivz,i) - ql(ivz,i)
v(ibx) = qr(ibx,i) - ql(ibx,i)
v(iby) = qr(iby,i) - ql(iby,i)
v(ibz) = qr(ibz,i) - ql(ibz,i)
v(ibp) = qr(ibp,i) - ql(ibp,i)
! the intermediate state averages
!
dna = amean(ql(idn,i), qr(idn,i))
vxa = amean(ql(ivx,i), qr(ivx,i))
vya = amean(ql(ivy,i), qr(ivy,i))
vza = amean(ql(ivz,i), qr(ivz,i))
bxa = amean(ql(ibx,i), qr(ibx,i))
bya = amean(ql(iby,i), qr(iby,i))
bza = amean(ql(ibz,i), qr(ibz,i))
bpa = amean(ql(ibp,i), qr(ibp,i))
dnl = lmean(ql(idn,i), qr(idn,i))
pta = csnd2 * dna + 5.0d-01 * amean(b2l, b2r)
! the eigenvalues and related parameters
!
sqd = sqrt(dnl)
b1 = bxa / sqd
b2 = bya / sqd
b3 = bza / sqd
bp2 = b2**2 + b3**2
bp = sqrt(bp2)
bb2 = b1**2 + bp2
bb = sqrt(bb2)
sgnb1 = sign(1.0d+00, b1)
if (bp > 0.0d+00) then
x2 = b2 / bp
x3 = b3 / bp
else
x2 = 0.0d+00
x3 = 0.0d+00
end if
ca2 = b1**2
ca = sqrt(ca2)
f1 = csnd2 + bb2
f3 = max(0.0d+00, f1**2 - 4.0d+00 * csnd2 * ca2)
f2 = sqrt(f3)
cf2 = 5.0d-01 * (f1 + f2)
cf = sqrt(cf2)
cs2 = 5.0d-01 * (f1 - f2)
cs = sqrt(cs2)
f2 = cf2 - cs2
if (f2 > 0.0d+00) then
alf2 = max(0.0d+00, csnd2 - cs2) / f2
als2 = max(0.0d+00, cf2 - csnd2) / f2
else
alf2 = 0.0d+00
als2 = 0.0d+00
end if
alf = sqrt(alf2)
als = sqrt(als2)
! === the average of the right eigenvector matrix ===
!
! right (1) and left (9) fast magnetosonic eigenvectors
!
f1 = dnl * alf
f2 = dnl * als * cs * sgnb1
f3 = als * csnd * sqd
rm(1,1) = f1
rm(2,1) = f1 * (vxa + cf)
rm(3,1) = f1 * vya - f2 * x2
rm(4,1) = f1 * vza - f2 * x3
rm(6,1) = f3 * x2
rm(7,1) = f3 * x3
rm(1,8) = f1
rm(2,8) = f1 * (vxa - cf)
rm(3,8) = f1 * vya + f2 * x2
rm(4,8) = f1 * vza + f2 * x3
rm(6,8) = rm(6,1)
rm(7,8) = rm(7,1)
! right (2) and left (8) Alfvénic eigenvectors
!
f1 = dnl * sqrt(dna)
rm(3,2) = f1 * x3
rm(4,2) = - f1 * x2
rm(6,2) = - dnl * x3
rm(7,2) = dnl * x2
rm(3,7) = - rm(3,2)
rm(4,7) = - rm(4,2)
rm(6,7) = rm(6,2)
rm(7,7) = rm(7,2)
! right (3) and left (7) slow magnetosonic eigenvectors
!
f1 = dnl * als
f2 = dnl * alf * cf * sgnb1
f3 = - alf * csnd * sqd
rm(1,3) = f1
rm(2,3) = f1 * (vxa + cs)
rm(3,3) = f1 * vya + f2 * x2
rm(4,3) = f1 * vza + f2 * x3
rm(6,3) = f3 * x2
rm(7,3) = f3 * x3
rm(1,6) = f1
rm(2,6) = f1 * (vxa - cs)
rm(3,6) = f1 * vya - f2 * x2
rm(4,6) = f1 * vza - f2 * x3
rm(6,6) = rm(6,3)
rm(7,6) = rm(7,3)
! the eigenvalue diagonal matrix
!
lm(1) = vxa + cf
lm(2) = vxa + ca
lm(3) = vxa + cs
lm(4) = vxa + ch
lm(5) = vxa - ch
lm(6) = vxa - cs
lm(7) = vxa - ca
lm(8) = vxa - cf
! the scaling diagonal matrix
!
tm(1) = 2.5d-01 / (dnl * csnd2)
tm(2) = 2.5d-01 / dnl**2
tm(3) = tm(1)
tm(4) = 2.5d-01
tm(5) = tm(4)
tm(6) = tm(1)
tm(7) = tm(2)
tm(8) = tm(1)
! KEPEC flux
!
fi(idn,i) = dnl * vxa
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
fi(imy,i) = fi(idn,i) * vya - bxa * bya
fi(imz,i) = fi(idn,i) * vza - bxa * bza
fi(ibx,i) = ch * bpa
fi(iby,i) = vxa * bya - bxa * vya
fi(ibz,i) = vxa * bza - bxa * vza
fi(ibp,i) = ch * bxa
! KEPES correction
!
lv = abs(lm) * tm * matmul(v, rm)
if (vxa >= 0.0d+00) then
do p = 1, nf
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
end do
else
do p = nf, 1, -1
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
end do
end if
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_iso_kepes
!
!*******************************************************************************
!
! ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use coordinates, only : nn => bcells
use equations , only : nf, idn, ivy, ivz, imx, imy, imz, ien, &
ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i
real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, ui
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
srml = sr - sl
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
if (b2 > 0.0d+00) then
vy = (wr(imy) - wl(imy)) / dn
vz = (wr(imz) - wl(imz)) / dn
by = (wr(iby) - wl(iby)) / srml
bz = (wr(ibz) - wl(ibz)) / srml
vb = sm * bx + vy * by + vz * bz
if (sm > 0.0d+00) then ! sm > 0
slmm = sl - sm
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
fi(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
srmm = sr - sm
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
fi(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
end if ! sm = 0
else ! Bₓ = 0
if (sm > 0.0d+00) then ! sm > 0
slmm = sl - sm
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
fi(:,i) = sl * ui(:) - wl(:)
else if (sm < 0.0d+00) then ! sm < 0
srmm = sr - sm
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
fi(:,i) = sr * ui(:) - wr(:)
else ! sm = 0
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
fi(:,i) = (sl * wr(:) - 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_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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use coordinates, only : nn => bcells
use equations , only : nf, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
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
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, ui
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then ! sl ≥ 0
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then ! sr ≤ 0
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
bx = ql(ibx,i)
b2 = ql(ibx,i) * qr(ibx,i)
srml = sr - sl
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
dn = wr(idn) - wl(idn)
sm = (wr(imx) - wl(imx)) / dn
slmm = sl - sm
srmm = sr - sm
dnl = wl(idn) / slmm
dnr = wr(idn) / srmm
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
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
!
fi(:,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
!
fi(:,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
!
fi(:,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
!
fi(:,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;
!
fi(:,i) = (sml * wcr(:) - 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
!
fi(:,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
!
fi(:,i) = sr * ui(:) - wr(:)
else ! Sm = 0
! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
fi(:,i) = (sl * wr(:) - 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
!
fi(:,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
!
fi(:,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
!
fi(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm
end if ! sm < 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
fi(:,i) = (sl * wr(:) - 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
!
fi(:,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
!
fi(:,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
!
fi(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm
end if ! sm > 0
else ! bx = 0
! no Alfvén wave, so revert to the HLL flux
!
fi(:,i) = (sl * wr(:) - 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
!
fi(:,i) = (sl * wr(:) - 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_MHD_ADI_ROE:
! ------------------------------
!
! Subroutine solves one dimensional Riemann problem using
! the Roe's method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use coordinates, only : nn => bcells
use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ipr, ien, &
ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed, eigensystem_roe
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
integer :: i, p
real(kind=8) :: sdl, sdr, sds
real(kind=8) :: xx, yy
real(kind=8) :: pml, pmr
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: qi, ci, al
real(kind=8), dimension(nf,nf) :: li, ri
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
! 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
fi(:,i) = fl(:,i)
else if (ci(nf) <= 0.0d+00) then
fi(:,i) = fr(:,i)
else
! prepare fluxes which do not change across the states
!
fi(ibx,i) = fl(ibx,i)
fi(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
!
fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
fi(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i))
fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
fi(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))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
end do
else
do p = nf, 1, -1
xx = - 0.5d+00 * al(p) * abs(ci(p))
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
end do
end if
end if ! sl < 0 < sr
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_roe
!
!===============================================================================
!
! subroutine RIEMANN_MHD_ADI_KEPES:
! --------------------------------
!
! Subroutine solves one dimensional adiabatic magnetohydrodynamic
! Riemann problem using the entropy stable KEPES method.
!
! Arguments:
!
! ql, qr - the primitive variables of the left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! References:
!
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
! field divergence diminishing ideal magnetohydrodynamics equations",
! Journal of Computational Physics, 2018, 364, pp. 420-467,
! https://doi.org/10.1016/j.jcp.2018.03.002
!
!===============================================================================
!
subroutine riemann_mhd_adi_kepes(ql, qr, fi)
use coordinates, only : nn => bcells
use equations , only : nf, adiabatic_index, ch => cglm
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
ibx, iby, ibz, ibp, ipr, ien
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
logical , save :: first = .true.
real(kind=8) , save :: adi_m1, adi_m1x
real(kind=8), dimension(9,9), save :: rm
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
integer :: i
real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa
real(kind=8) :: dnl, prl, pta, vva, br, bl, bp
real(kind=8) :: btl, bta, eka, ema, ub2, uba
real(kind=8) :: aa2, al2, ab2, bb2, bp2, aa, al, ab, b1, b2, b3, bb
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
real(kind=8) :: alf2, als2, alf, als, wps, wms, wpf, wmf
real(kind=8) :: f1, f2, f3, f4, sgnb1, v2l, v2r, b2l, b2r
real(kind=8), dimension(nf) :: v, lm, tm
!-------------------------------------------------------------------------------
!
if (first) then
adi_m1 = adiabatic_index - 1.0d+00
adi_m1x = adi_m1 / adiabatic_index
rm(:,:) = 0.0d+00
rm(6,4) = 1.0d+00
rm(9,4) = 1.0d+00
rm(1,5) = 1.0d+00
rm(6,6) = 1.0d+00
rm(9,6) = - 1.0d+00
first = .false.
end if
do i = 1, nn
! the difference in entropy variables between the states
!
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
bl = ql(idn,i) / ql(ipr,i)
br = qr(idn,i) / qr(ipr,i)
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
+ (log(br) - log(bl)) / adi_m1 &
- 5.0d-01 * (br * v2r - bl * v2l)
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
v(ipr) = bl - br
v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
v(iby) = br * qr(iby,i) - bl * ql(iby,i)
v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)
! the intermediate state averages
!
btl = lmean(bl, br)
bta = amean(bl, br)
dna = amean(ql(idn,i), qr(idn,i))
vxa = amean(ql(ivx,i), qr(ivx,i))
vya = amean(ql(ivy,i), qr(ivy,i))
vza = amean(ql(ivz,i), qr(ivz,i))
bxa = amean(ql(ibx,i), qr(ibx,i))
bya = amean(ql(iby,i), qr(iby,i))
bza = amean(ql(ibz,i), qr(ibz,i))
bpa = amean(ql(ibp,i), qr(ibp,i))
pra = dna / bta
dnl = lmean(ql(idn,i), qr(idn,i))
prl = lmean(ql(ipr,i), qr(ipr,i))
eka = 5.0d-01 * amean(v2l, v2r)
ema = 5.0d-01 * amean(b2l, b2r)
pta = pra + ema
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r)
uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), &
sum(qr(ivx:ivz,i) * qr(ibx:ibz,i)))
! the eigenvalues and related parameters
!
aa2 = adiabatic_index * pra / dnl
aa = sqrt(aa2)
al2 = adiabatic_index * prl / dnl
al = sqrt(al2)
ab2 = adiabatic_index / bta
ab = sqrt(ab2)
b1 = bxa / sqrt(dnl)
b2 = bya / sqrt(dnl)
b3 = bza / sqrt(dnl)
bp2 = b2**2 + b3**2
bp = sqrt(bp2)
bb2 = b1**2 + bp2
bb = sqrt(bb2)
sgnb1 = sign(1.0d+00, b1)
if (bp > 0.0d+00) then
x2 = b2 / bp
x3 = b3 / bp
else
x2 = 0.0d+00
x3 = 0.0d+00
end if
ca2 = b1**2
ca = sqrt(ca2)
f1 = aa2 + bb2
f3 = max(0.0d+00, f1**2 - 4.0d+00 * aa2 * ca2)
f2 = sqrt(f3)
cf2 = 5.0d-01 * (f1 + f2)
cf = sqrt(cf2)
cs2 = 5.0d-01 * (f1 - f2)
cs = sqrt(cs2)
f2 = cf2 - cs2
if (f2 > 0.0d+00) then
alf2 = max(0.0d+00, aa2 - cs2) / f2
als2 = max(0.0d+00, cf2 - aa2) / f2
else
alf2 = 0.0d+00
als2 = 0.0d+00
end if
alf = sqrt(alf2)
als = sqrt(als2)
f3 = vva + al2 / adi_m1
f4 = sgnb1 * (vya * x2 + vza * x3)
f1 = dnl * (als * f3 - alf * ab * bp)
f2 = dnl * (als * cs * vxa + alf * cf * f4)
wps = f1 + f2
wms = f1 - f2
f1 = dnl * (alf * f3 + als * ab * bp)
f2 = dnl * (alf * cf * vxa - als * cs * f4)
wpf = f1 + f2
wmf = f1 - f2
! the scaling diagonal matrix
!
f1 = 5.0d-01 / bta
tm(1) = 5.0d-01 / (adiabatic_index * dnl)
tm(2) = f1 / dnl**2
tm(3) = tm(1)
tm(4) = f1
tm(5) = dnl * adi_m1x
tm(6) = tm(4)
tm(7) = tm(1)
tm(8) = tm(2)
tm(9) = tm(1)
! === the average of the right eigenvector matrix ===
!
! right (1) and left (9) fast magnetosonic eigenvectors
!
f1 = dnl * alf
f2 = dnl * als * cs * sgnb1
f3 = als * ab * sqrt(dnl)
rm(1,1) = f1
rm(2,1) = f1 * (vxa + cf)
rm(3,1) = f1 * vya - f2 * x2
rm(4,1) = f1 * vza - f2 * x3
rm(5,1) = wpf
rm(7,1) = f3 * x2
rm(8,1) = f3 * x3
rm(1,9) = f1
rm(2,9) = f1 * (vxa - cf)
rm(3,9) = f1 * vya + f2 * x2
rm(4,9) = f1 * vza + f2 * x3
rm(5,9) = wmf
rm(7,9) = rm(7,1)
rm(8,9) = rm(8,1)
! right (2) and left (8) Alfvénic eigenvectors
!
f1 = dnl * sqrt(dna)
rm(3,2) = f1 * x3
rm(4,2) = - f1 * x2
rm(5,2) = - f1 * (x2 * vza - x3 * vya)
rm(7,2) = - dnl * x3
rm(8,2) = dnl * x2
rm(3,8) = - rm(3,2)
rm(4,8) = - rm(4,2)
rm(5,8) = - rm(5,2)
rm(7,8) = rm(7,2)
rm(8,8) = rm(8,2)
! right (3) and left (7) slow magnetosonic eigenvectors
!
f1 = dnl * als
f2 = dnl * alf * cf * sgnb1
f3 = - alf * ab * sqrt(dnl)
rm(1,3) = f1
rm(2,3) = f1 * (vxa + cs)
rm(3,3) = f1 * vya + f2 * x2
rm(4,3) = f1 * vza + f2 * x3
rm(5,3) = wps
rm(7,3) = f3 * x2
rm(8,3) = f3 * x3
rm(1,7) = f1
rm(2,7) = f1 * (vxa - cs)
rm(3,7) = f1 * vya - f2 * x2
rm(4,7) = f1 * vza - f2 * x3
rm(5,7) = wms
rm(7,7) = rm(7,3)
rm(8,7) = rm(8,3)
! right (4) and left (6) GLM eigenvectors
!
rm(5,4) = bxa + bpa
rm(5,6) = bxa - bpa
! entropy (5) eigenvectors
!
rm(2,5) = vxa
rm(3,5) = vya
rm(4,5) = vza
rm(5,5) = vva
! the eigenvalue diagonal matrix
!
lm(1) = vxa + cf
lm(2) = vxa + ca
lm(3) = vxa + cs
lm(4) = vxa + ch
lm(5) = vxa
lm(6) = vxa - ch
lm(7) = vxa - cs
lm(8) = vxa - ca
lm(9) = vxa - cf
! KEPEC flux
!
fi(idn,i) = dnl * vxa
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
fi(imy,i) = fi(idn,i) * vya - bxa * bya
fi(imz,i) = fi(idn,i) * vza - bxa * bza
fi(ibx,i) = ch * bpa
fi(iby,i) = vxa * bya - bxa * vya
fi(ibz,i) = vxa * bza - bxa * vza
fi(ibp,i) = ch * bxa
fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
+ (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) &
+ (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) &
+ fi(ibp,i) * bpa - ub2 + bxa * uba &
- ch * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
! KEPES correction
!
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
end do
!-------------------------------------------------------------------------------
!
end subroutine riemann_mhd_adi_kepes
!
!*******************************************************************************
!
! RELATIVISTIC ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
!
!*******************************************************************************
!
!===============================================================================
!
! 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 left and right Riemann states;
! fi - the numerical flux at the cell interface;
!
! 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, fi)
use algebra , only : quadratic
use coordinates, only : nn => bcells
use equations , only : nf, idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp
use equations , only : prim2cons, fluxspeed
implicit none
real(kind=8), dimension(:,:), intent(in) :: ql, qr
real(kind=8), dimension(:,:), intent(out) :: fi
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
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
real(kind=8), dimension( 2,nn) :: cl, cr
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
real(kind=8), dimension(3) :: a, vs
real(kind=8), dimension(2) :: x
!-------------------------------------------------------------------------------
!
call prim2cons(ql, ul)
call prim2cons(qr, ur)
call fluxspeed(ql, ul, fl, cl)
call fluxspeed(qr, ur, fr, cr)
do i = 1, nn
sl = min(cl(1,i), cr(1,i))
sr = max(cl(2,i), cr(2,i))
if (sl >= 0.0d+00) then
fi(:,i) = fl(:,i)
else if (sr <= 0.0d+00) then
fi(:,i) = fr(:,i)
else ! sl < 0 < sr
srml = sr - sl
wl(:) = sl * ul(:,i) - fl(:,i)
wr(:) = sr * ur(:,i) - fr(:,i)
uh(:) = ( wr(:) - wl(:)) / srml
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
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
fi(:,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
fi(:,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
fi(:,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])
!
fi(:,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])
!
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
fi(:,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
fi(:,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
fi(:,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
fi(:,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])
!
fi(:,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])
!
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
else
! intermediate flux is constant across the contact discontinuity so use HLL flux
!
fi(:,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
!
!*******************************************************************************
!
! SUPPORTING SUBROUTINES/FUNCTIONS
!
!*******************************************************************************
!
!===============================================================================
!
! 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
!
!===============================================================================
!
! function LMEAN:
! --------------
!
! Function calculates the logarithmic mean using the optimized algorithm
! by Ranocha et al.
!
! Arguments:
!
! l, r - the values to calculate the mean of;
!
! References:
!
! [1] Ranocha et al.,
! "Efficient implementation of modern entropy stable and kinetic energy
! preserving discontinuous Galerkin methods for conservation laws",
! http://arxiv.org/abs/2112.10517v1
!
!===============================================================================
!
real(kind=8) function lmean(l, r)
implicit none
real(kind=8), intent(in) :: l, r
real(kind=8) :: u, d, s
real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, &
c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00
!-------------------------------------------------------------------------------
!
d = r - l
s = r + l
u = (d / s)**2
if (u < 1.0d-04) then
lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
else
if (d >= 0.0d+00) then
s = log(r / l)
else
s = log(l / r)
end if
lmean = abs(d) / s
end if
!-------------------------------------------------------------------------------
!
end function lmean
!
!===============================================================================
!
! function AMEAN:
! --------------
!
! Function calculate the arithmetic mean.
!
! Arguments:
!
! l, r - the values to calculate the mean of;
!
!===============================================================================
!
real(kind=8) function amean(l, r)
implicit none
real(kind=8), intent(in) :: l, r
!-------------------------------------------------------------------------------
!
amean = 5.0d-01 * (l + r)
!-------------------------------------------------------------------------------
!
end function amean
2008-12-08 21:04:20 -06:00
!===============================================================================
!
end module schemes