4283 lines
117 KiB
Fortran
4283 lines
117 KiB
Fortran
!!******************************************************************************
|
||
!!
|
||
!! This file is part of the AMUN source code, a program to perform
|
||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||
!! adaptive mesh.
|
||
!!
|
||
!! Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
|
||
!!
|
||
!! This program is free software: you can redistribute it and/or modify
|
||
!! it under the terms of the GNU General Public License as published by
|
||
!! the Free Software Foundation, either version 3 of the License, or
|
||
!! (at your option) any later version.
|
||
!!
|
||
!! This program is distributed in the hope that it will be useful,
|
||
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
!! GNU General Public License for more details.
|
||
!!
|
||
!! You should have received a copy of the GNU General Public License
|
||
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
!!
|
||
!!******************************************************************************
|
||
!!
|
||
!! module: SCHEMES
|
||
!!
|
||
!! The module provides and interface to numerical schemes, subroutines to
|
||
!! calculate variable increment and one dimensional Riemann solvers.
|
||
!!
|
||
!! If you implement a new set of equations, you have to add at least one
|
||
!! corresponding update_flux subroutine, and one Riemann solver.
|
||
!!
|
||
!!******************************************************************************
|
||
!
|
||
module schemes
|
||
|
||
implicit none
|
||
|
||
! Rieman solver and state vectors names
|
||
!
|
||
character(len=255) , save :: name_sol = ""
|
||
character(len=255) , save :: name_sts = "primitive"
|
||
|
||
! 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:791–824,
|
||
! 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, cs
|
||
real(kind=8) :: dnl, dna, vxa, vya, vza
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
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)
|
||
|
||
cs = sign(csnd, vxa)
|
||
|
||
lm(:) = [ vxa + cs, vxa, vxa, vxa - cs ]
|
||
|
||
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
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
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
|
||
! https://doi.org/10.1016/j.jcp.2016.12.006
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
|
||
integer :: i
|
||
real(kind=8) :: bl, br, v2l, v2r, cs, vc
|
||
real(kind=8) :: btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra, ent
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
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)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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
|
||
|
||
bl = ql(idn,i) / ql(ipr,i)
|
||
br = qr(idn,i) / qr(ipr,i)
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
|
||
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
|
||
|
||
v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
if (br > bl) then
|
||
v(idn) = v(idn) + log(br / bl) / adi_m1
|
||
else if (bl > br) then
|
||
v(idn) = v(idn) - log(bl / br) / adi_m1
|
||
end if
|
||
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
|
||
|
||
cs = sign(sqrt(adiabatic_index * pra / dnl), vxa)
|
||
vc = vxa * cs
|
||
|
||
lm(1) = vxa + cs
|
||
lm(2) = vxa
|
||
lm(3) = vxa
|
||
lm(4) = vxa
|
||
lm(5) = vxa - cs
|
||
|
||
tm(1) = 5.0d-01 * dnl / adiabatic_index
|
||
tm(2) = dnl * adi_m1x
|
||
tm(3) = pra
|
||
tm(4) = pra
|
||
tm(5) = tm(1)
|
||
|
||
rm(2,1) = lm(1)
|
||
rm(3,1) = vya
|
||
rm(4,1) = vza
|
||
rm(2,2) = vxa
|
||
rm(3,2) = vya
|
||
rm(4,2) = vza
|
||
rm(2,5) = lm(5)
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,:) = [ ent + vc, vva, vya, vza, ent - vc ]
|
||
|
||
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
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
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, dnl, dnr
|
||
|
||
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
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn))
|
||
|
||
srml = sr - sl
|
||
slmm = sl - sm
|
||
srmm = sr - sm
|
||
|
||
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)
|
||
|
||
sml = sm - abs(bx) / sqrt(uil(idn))
|
||
smr = sm + abs(bx) / sqrt(uir(idn))
|
||
|
||
dvl = slmm * wl(idn) - b2
|
||
dvr = srmm * wr(idn) - b2
|
||
|
||
dnr = uir(idn) / dvr
|
||
dnl = uil(idn) / dvl
|
||
|
||
! 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) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
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) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
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) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
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) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
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) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
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
|
||
|
||
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, csnd2
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
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) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2
|
||
real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als
|
||
real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff
|
||
real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb
|
||
real(kind=8) :: f1, f2, f3, f4
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qi
|
||
real(kind=8), dimension(6 ) :: lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
integer , dimension(6) , save :: ivs
|
||
real(kind=8), dimension(6,6), save :: rvec, lvec
|
||
!$omp threadprivate(first, ivs, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
|
||
rvec(:,:) = 0.0d+00
|
||
lvec(:,:) = 0.0d+00
|
||
|
||
ivs = [ idn, imx, imy, imz, iby, ibz ]
|
||
|
||
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
|
||
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
|
||
qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
f1 = ql(iby,i) - qr(iby,i)
|
||
f2 = ql(ibz,i) - qr(ibz,i)
|
||
xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
|
||
yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
|
||
br = sqrt(br2)
|
||
if (br2 > 0.0d+00) then
|
||
bty = qi(iby) / br
|
||
btz = qi(ibz) / br
|
||
else
|
||
bty = 0.0d+00
|
||
btz = 0.0d+00
|
||
end if
|
||
|
||
cc2 = csnd2 + xx
|
||
ca2 = qi(ibx) * qi(ibx) / qi(idn)
|
||
f1 = br2 / qi(idn)
|
||
f2 = ca2 + f1
|
||
f3 = f2 - cc2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * cc2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + cc2 + f4)
|
||
cs2 = cc2 * ca2 / cf2
|
||
|
||
if (cs2 < cc2 .and. cc2 < cf2) then
|
||
f1 = (cc2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (cc2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cc = sqrt(cc2)
|
||
ca = sqrt(ca2)
|
||
cf = sign(sqrt(cf2), qi(ivx))
|
||
cs = sign(sqrt(cs2), qi(ivx))
|
||
|
||
lm(1) = qi(ivx) + cf
|
||
lm(2) = qi(ivx) + ca
|
||
lm(3) = qi(ivx) + cs
|
||
lm(4) = qi(ivx) - cs
|
||
lm(5) = qi(ivx) - ca
|
||
lm(6) = qi(ivx) - cf
|
||
|
||
sqrtd = sqrt(qi(idn))
|
||
f1 = sqrt(yy)
|
||
qf = cf * sign(alf, qi(ibx)) / f1
|
||
qs = cs * sign(als, qi(ibx)) / f1
|
||
f1 = cc / (f1 * sqrtd)
|
||
af_prime = alf * f1
|
||
as_prime = als * f1
|
||
|
||
rvec(1,1) = alf
|
||
rvec(1,2) = alf * lm(1)
|
||
rvec(1,3) = alf * qi(ivy) - qs * bty
|
||
rvec(1,4) = alf * qi(ivz) - qs * btz
|
||
rvec(1,5) = as_prime * bty
|
||
rvec(1,6) = as_prime * btz
|
||
|
||
rvec(6,1) = alf
|
||
rvec(6,2) = alf * lm(6)
|
||
rvec(6,3) = alf * qi(ivy) + qs * bty
|
||
rvec(6,4) = alf * qi(ivz) + qs * btz
|
||
rvec(6,5) = rvec(1,5)
|
||
rvec(6,6) = rvec(1,6)
|
||
|
||
rvec(3,1) = als
|
||
rvec(3,2) = als * lm(3)
|
||
rvec(3,3) = als * qi(ivy) + qf * bty
|
||
rvec(3,4) = als * qi(ivz) + qf * btz
|
||
rvec(3,5) = - af_prime * bty
|
||
rvec(3,6) = - af_prime * btz
|
||
|
||
rvec(4,1) = als
|
||
rvec(4,2) = als * lm(4)
|
||
rvec(4,3) = als * qi(ivy) - qf * bty
|
||
rvec(4,4) = als * qi(ivz) - qf * btz
|
||
rvec(4,5) = rvec(3,5)
|
||
rvec(4,6) = rvec(3,6)
|
||
|
||
f1 = sign(1.0d+00 / sqrtd, qi(ibx))
|
||
rvec(2,3) = btz
|
||
rvec(2,4) = - bty
|
||
rvec(2,5) = - btz * f1
|
||
rvec(2,6) = bty * f1
|
||
|
||
rvec(5,3) = - rvec(2,3)
|
||
rvec(5,4) = - rvec(2,4)
|
||
rvec(5,5) = rvec(2,5)
|
||
rvec(5,6) = rvec(2,6)
|
||
|
||
norm = 2.0d+00 * cc2
|
||
cff = alf * cf / norm
|
||
css = als * cs / norm
|
||
norm = norm / yy
|
||
qf = qf / norm
|
||
qs = qs / norm
|
||
f1 = qi(idn) / norm
|
||
af = af_prime * f1
|
||
as = as_prime * f1
|
||
f1 = br / norm
|
||
afpb = af_prime * f1
|
||
aspb = as_prime * f1
|
||
vqstr = (qi(ivy) * bty + qi(ivz) * btz)
|
||
|
||
f1 = qs * vqstr
|
||
lvec(1,1) = - cff * lm(6) + f1 - aspb
|
||
lvec(1,2) = cff
|
||
lvec(1,3) = - qs * bty
|
||
lvec(1,4) = - qs * btz
|
||
lvec(1,5) = as * bty
|
||
lvec(1,6) = as * btz
|
||
|
||
lvec(6,1) = cff * lm(1) - f1 - aspb
|
||
lvec(6,2) = - lvec(1,2)
|
||
lvec(6,3) = - lvec(1,3)
|
||
lvec(6,4) = - lvec(1,4)
|
||
lvec(6,5) = lvec(1,5)
|
||
lvec(6,6) = lvec(1,6)
|
||
|
||
f1 = qf * vqstr
|
||
lvec(3,1) = - css * lm(4) - f1 + afpb
|
||
lvec(3,2) = css
|
||
lvec(3,3) = qf * bty
|
||
lvec(3,4) = qf * btz
|
||
lvec(3,5) = - af * bty
|
||
lvec(3,6) = - af * btz
|
||
|
||
lvec(4,1) = css * lm(3) + f1 + afpb
|
||
lvec(4,2) = - lvec(3,2)
|
||
lvec(4,3) = - lvec(3,3)
|
||
lvec(4,4) = - lvec(3,4)
|
||
lvec(4,5) = lvec(3,5)
|
||
lvec(4,6) = lvec(3,6)
|
||
|
||
f1 = sign(5.0d-01 * sqrtd, qi(ibx))
|
||
lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty)
|
||
lvec(2,3) = 5.0d-01 * btz
|
||
lvec(2,4) = - 5.0d-01 * bty
|
||
lvec(2,5) = - f1 * btz
|
||
lvec(2,6) = f1 * bty
|
||
|
||
lvec(5,1) = - lvec(2,1)
|
||
lvec(5,3) = - lvec(2,3)
|
||
lvec(5,4) = - lvec(2,4)
|
||
lvec(5,5) = lvec(2,5)
|
||
lvec(5,6) = lvec(2,6)
|
||
|
||
al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
|
||
df = matmul(al, rvec)
|
||
fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
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
|
||
! [2] 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:791–824,
|
||
! https://doi.org/10.1007/s10543-019-00789-w
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
|
||
integer :: i
|
||
real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
|
||
real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq
|
||
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
|
||
real(kind=8) :: alf, als, f1, f2, f3, f4
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8), dimension(8,8), save :: rm
|
||
!$omp threadprivate(first, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
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
|
||
|
||
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))
|
||
|
||
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))
|
||
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))
|
||
pta = 5.0d-01 * amean(b2l, b2r) + csnd2 * dna
|
||
|
||
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)
|
||
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)
|
||
|
||
bp2 = bya * bya + bza * bza
|
||
if (bp2 > 0.0d+00) then
|
||
f1 = sqrt(bp2)
|
||
x2 = bya / f1
|
||
x3 = bza / f1
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
ca2 = bxa * bxa / dnl
|
||
f1 = bp2 / dnl
|
||
f2 = ca2 + f1
|
||
f3 = f2 - csnd2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * csnd2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + csnd2 + f4)
|
||
cs2 = csnd2 * ca2 / cf2
|
||
|
||
if (cs2 < csnd2 .and. csnd2 < cf2) then
|
||
f1 = (csnd2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (csnd2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), vxa)
|
||
ca = sqrt(ca2)
|
||
cs = sign(sqrt(cs2), vxa)
|
||
|
||
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
|
||
|
||
tm(1) = 5.0d-01 / (dnl * csnd2)
|
||
tm(2) = 5.0d-01 / dnl**2
|
||
tm(3) = tm(1)
|
||
tm(4) = 5.0d-01
|
||
tm(5) = tm(4)
|
||
tm(6) = tm(1)
|
||
tm(7) = tm(2)
|
||
tm(8) = tm(1)
|
||
|
||
das = dnl * als
|
||
daf = dnl * alf
|
||
csq = csnd * sqrt(dnl)
|
||
|
||
f1 = daf
|
||
f2 = sign(das, bxa) * cs
|
||
f3 = als * csq
|
||
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)
|
||
|
||
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)
|
||
|
||
f1 = das
|
||
f2 = sign(daf, bxa) * cf
|
||
f3 = - alf * csq
|
||
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)
|
||
|
||
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(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
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, adiabatic_index
|
||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, 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) :: sdl, sdr, sds
|
||
real(kind=8) :: xx, yy, pml, pmr
|
||
|
||
real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs
|
||
real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty
|
||
real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm
|
||
real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb
|
||
real(kind=8) :: f1, f2, f3, f4, f5
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qi
|
||
real(kind=8), dimension(7 ) :: lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
integer , dimension(7) , save :: ivs
|
||
real(kind=8) , save :: adi_m1, adi_m2, adi_m2d1
|
||
real(kind=8), dimension(7,7), save :: rvec, lvec
|
||
!$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m2 = adiabatic_index - 2.0d+00
|
||
adi_m2d1 = adi_m2 / adi_m1
|
||
|
||
rvec(:,:) = 0.0d+00
|
||
lvec(:,:) = 0.0d+00
|
||
rvec(4,1) = 1.0d+00
|
||
|
||
ivs = [ idn, imx, imy, imz, ipr, iby, ibz ]
|
||
|
||
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
|
||
|
||
pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
|
||
pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + 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(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
f1 = qr(iby,i) - ql(iby,i)
|
||
f2 = qr(ibz,i) - ql(ibz,i)
|
||
xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
|
||
yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
|
||
if (br2 > 0.0d+00) then
|
||
br = sqrt(br2)
|
||
bty = qi(iby) / br
|
||
btz = qi(ibz) / br
|
||
else
|
||
br = 0.0d+00
|
||
bty = 0.0d+00
|
||
btz = 0.0d+00
|
||
end if
|
||
|
||
v2 = sum(qi(ivx:ivz) * qi(ivx:ivz))
|
||
v2h = 5.0d-01 * v2
|
||
ayy = adi_m1 - adi_m2 * yy
|
||
sqty = sqrt(ayy)
|
||
ca2 = qi(ibx) * qi(ibx)
|
||
hp = qi(ien) - (ca2 + br2) / qi(idn)
|
||
ca2 = ca2 / qi(idn)
|
||
f1 = adi_m1 * (hp - v2h)
|
||
f2 = adi_m2 * xx
|
||
if (f1 > f2) then
|
||
cc2 = f1 - f2
|
||
cc = sqrt(cc2)
|
||
f1 = ayy * br2 / qi(idn)
|
||
f2 = ca2 + f1
|
||
f3 = f2 + cc2
|
||
f4 = f2 - cc2
|
||
f5 = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1)
|
||
cf2 = 5.0d-01 * (f3 + f5)
|
||
cs2 = cc2 * ca2 / cf2
|
||
else
|
||
cf2 = ca2 + ayy * br2 / qi(idn)
|
||
cc2 = 0.0d+00
|
||
cc = 0.0d+00
|
||
cs2 = 0.0d+00
|
||
end if
|
||
|
||
if (cs2 < cc2 .and. cc2 < cf2) then
|
||
f1 = (cc2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (cc2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), qi(ivx))
|
||
ca = sign(sqrt(ca2), qi(ivx))
|
||
cs = sign(sqrt(cs2), qi(ivx))
|
||
|
||
lm(1) = qi(ivx) + cf
|
||
lm(2) = qi(ivx) + ca
|
||
lm(3) = qi(ivx) + cs
|
||
lm(4) = qi(ivx)
|
||
lm(5) = qi(ivx) - cs
|
||
lm(6) = qi(ivx) - ca
|
||
lm(7) = qi(ivx) - cf
|
||
|
||
vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty
|
||
|
||
sqrtd = sqrt(qi(idn))
|
||
qf = cf * sign(alf, qi(ibx))
|
||
qs = cs * sign(als, qi(ibx))
|
||
f1 = cc / sqrtd
|
||
af_prime = f1 * alf
|
||
as_prime = f1 * als
|
||
f1 = br / sqty
|
||
afpbb = af_prime * f1
|
||
aspbb = as_prime * f1
|
||
|
||
f1 = qs / sqty
|
||
f2 = qs * vbet
|
||
f3 = as_prime / sqty
|
||
rvec(1,1) = alf
|
||
rvec(1,2) = alf * lm(1)
|
||
rvec(1,3) = alf * qi(ivy) - f1 * bty
|
||
rvec(1,4) = alf * qi(ivz) - f1 * btz
|
||
rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb
|
||
rvec(1,6) = f3 * bty
|
||
rvec(1,7) = f3 * btz
|
||
|
||
rvec(7,1) = alf
|
||
rvec(7,2) = alf * lm(7)
|
||
rvec(7,3) = alf * qi(ivy) + f1 * bty
|
||
rvec(7,4) = alf * qi(ivz) + f1 * btz
|
||
rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb
|
||
rvec(7,6) = rvec(1,6)
|
||
rvec(7,7) = rvec(1,7)
|
||
|
||
f1 = qf / sqty
|
||
f2 = qf * vbet
|
||
f3 = - af_prime / sqty
|
||
rvec(3,1) = als
|
||
rvec(3,2) = als * lm(3)
|
||
rvec(3,3) = als * qi(ivy) + f1 * bty
|
||
rvec(3,4) = als * qi(ivz) + f1 * btz
|
||
rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb
|
||
rvec(3,6) = f3 * bty
|
||
rvec(3,7) = f3 * btz
|
||
|
||
rvec(5,1) = als
|
||
rvec(5,2) = als * lm(5)
|
||
rvec(5,3) = als * qi(ivy) - f1 * bty
|
||
rvec(5,4) = als * qi(ivz) - f1 * btz
|
||
rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb
|
||
rvec(5,6) = rvec(3,6)
|
||
rvec(5,7) = rvec(3,7)
|
||
|
||
f1 = sign(1.0d+00, qi(ivx))
|
||
f2 = sign(1.0d+00 / sqrtd, qi(ibx))
|
||
rvec(2,3) = f1 * btz
|
||
rvec(2,4) = - f1 * bty
|
||
rvec(2,5) = f1 * (qi(ivy) * btz - qi(ivz) * bty)
|
||
rvec(2,6) = - f2 * btz
|
||
rvec(2,7) = f2 * bty
|
||
|
||
rvec(6,3) = - rvec(2,3)
|
||
rvec(6,4) = - rvec(2,4)
|
||
rvec(6,5) = - rvec(2,5)
|
||
rvec(6,6) = rvec(2,6)
|
||
rvec(6,7) = rvec(2,7)
|
||
|
||
rvec(4,2) = qi(ivx)
|
||
rvec(4,3) = qi(ivy)
|
||
rvec(4,4) = qi(ivz)
|
||
rvec(4,5) = v2h + adi_m2d1 * xx
|
||
|
||
norm = 2.0d+00 * cc2
|
||
f1 = sqty * br / norm
|
||
afpb = af_prime * f1
|
||
aspb = as_prime * f1
|
||
sqty = sqty / norm
|
||
vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty
|
||
|
||
f1 = adi_m1 * alf / norm
|
||
f2 = qs * sqty
|
||
f3 = as_prime * qi(idn) * sqty
|
||
f4 = alf * cf / norm
|
||
f5 = qs * vqstr
|
||
lvec(1,1) = f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb
|
||
lvec(1,5) = f1
|
||
lvec(1,2) = - f1 * qi(ivx) + f4
|
||
lvec(1,3) = - f1 * qi(ivy) - f2 * bty
|
||
lvec(1,4) = - f1 * qi(ivz) - f2 * btz
|
||
lvec(1,6) = f3 * bty - f1 * qi(iby)
|
||
lvec(1,7) = f3 * btz - f1 * qi(ibz)
|
||
|
||
lvec(7,1) = f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb
|
||
lvec(7,5) = f1
|
||
lvec(7,2) = - f1 * qi(ivx) - f4
|
||
lvec(7,3) = - f1 * qi(ivy) + f2 * bty
|
||
lvec(7,4) = - f1 * qi(ivz) + f2 * btz
|
||
lvec(7,6) = lvec(1,6)
|
||
lvec(7,7) = lvec(1,7)
|
||
|
||
f1 = adi_m1 * als / norm
|
||
f2 = qf * sqty
|
||
f3 = af_prime * qi(idn) * sqty
|
||
f4 = als * cs / norm
|
||
f5 = qf * vqstr
|
||
lvec(3,1) = f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb
|
||
lvec(3,5) = f1
|
||
lvec(3,2) = - f1 * qi(ivx) + f4
|
||
lvec(3,3) = - f1 * qi(ivy) + f2 * bty
|
||
lvec(3,4) = - f1 * qi(ivz) + f2 * btz
|
||
lvec(3,6) = - f3 * bty - f1 * qi(iby)
|
||
lvec(3,6) = - f3 * btz - f1 * qi(ibz)
|
||
|
||
lvec(5,1) = f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb
|
||
lvec(5,5) = f1
|
||
lvec(5,2) = - f1 * qi(ivx) - f4
|
||
lvec(5,3) = - f1 * qi(ivy) - f2 * bty
|
||
lvec(5,4) = - f1 * qi(ivz) - f2 * btz
|
||
lvec(5,6) = lvec(3,6)
|
||
lvec(5,7) = lvec(3,7)
|
||
|
||
f1 = sign(5.0d-01, qi(ivx)) * bty
|
||
f2 = sign(5.0d-01, qi(ivx)) * btz
|
||
f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx))
|
||
lvec(2,1) = qi(ivz) * f1 - qi(ivy) * f2
|
||
lvec(2,3) = f2
|
||
lvec(2,4) = - f1
|
||
lvec(2,6) = - f2 * f3
|
||
lvec(2,7) = f1 * f3
|
||
|
||
lvec(6,1) = - lvec(2,1)
|
||
lvec(6,3) = - lvec(2,3)
|
||
lvec(6,4) = - lvec(2,4)
|
||
lvec(6,6) = lvec(2,6)
|
||
lvec(6,7) = lvec(2,7)
|
||
|
||
f1 = 1.0d+00 / cc2
|
||
lvec(4,1) = 1.0d+00 - f1 * (v2h - adi_m2d1 * xx)
|
||
lvec(4,5) = - f1
|
||
lvec(4,2) = f1 * qi(ivx)
|
||
lvec(4,3) = f1 * qi(ivy)
|
||
lvec(4,4) = f1 * qi(ivz)
|
||
lvec(4,6) = f1 * qi(iby)
|
||
lvec(4,7) = f1 * qi(ibz)
|
||
|
||
al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
|
||
df = matmul(al, rvec)
|
||
fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
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, 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
|
||
|
||
integer :: i
|
||
real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa, dnl, prl, pta
|
||
real(kind=8) :: v2l, v2r, b2l, b2r, bl, br, bta, btl, bp2, eka, ema, vva
|
||
real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, aa2, x2, x3, ub2, uba
|
||
real(kind=8) :: alf, als, das, daf, csq, wps, wms, wpf, wmf
|
||
real(kind=8) :: f1, f2, f3, f4, f5
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x, adi_m1xi
|
||
real(kind=8), dimension(9,9), save :: rm
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, adi_m1xi, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adi_m1 / adiabatic_index
|
||
adi_m1xi = adiabatic_index / adi_m1
|
||
|
||
rm(:,:) = 0.0d+00
|
||
rm(1,5) = 1.0d+00
|
||
rm(6,4) = 1.0d+00
|
||
rm(6,6) = 1.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))
|
||
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)
|
||
|
||
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))
|
||
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))
|
||
prl = lmean(ql(ipr,i), qr(ipr,i))
|
||
pra = dna / bta
|
||
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)))
|
||
|
||
v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
if (br > bl) then
|
||
v(idn) = v(idn) + log(br / bl) / adi_m1
|
||
else if (bl > br) then
|
||
v(idn) = v(idn) - log(bl / br) / adi_m1
|
||
end if
|
||
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)
|
||
|
||
bp2 = bya * bya + bza * bza
|
||
if (bp2 > 0.0d+00) then
|
||
f1 = sqrt(bp2)
|
||
x2 = bya / f1
|
||
x3 = bza / f1
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
aa2 = adiabatic_index * pra / dnl
|
||
ca2 = bxa * bxa / dnl
|
||
f1 = bp2 / dnl
|
||
f2 = ca2 + f1
|
||
f3 = f2 - aa2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * aa2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + aa2 + f4)
|
||
cs2 = aa2 * ca2 / cf2
|
||
|
||
if (cs2 < aa2 .and. aa2 < cf2) then
|
||
f1 = (aa2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (aa2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), vxa)
|
||
cs = sign(sqrt(cs2), vxa)
|
||
ca = sign(sqrt(ca2), vxa)
|
||
ch = sign( cglm, vxa)
|
||
|
||
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
|
||
|
||
f1 = 5.0d-01 / bta
|
||
f2 = adiabatic_index * dnl
|
||
tm(1) = 5.0d-01 / f2
|
||
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)
|
||
|
||
das = dnl * als
|
||
daf = dnl * alf
|
||
csq = sqrt(f2 / bta)
|
||
|
||
f3 = vva + adi_m1xi * prl / dnl
|
||
f4 = vya * x2 + vza * x3
|
||
f5 = sqrt(adiabatic_index * bp2 / (dnl * bta))
|
||
f1 = dnl * (als * f3 - alf * f5)
|
||
f2 = dnl * (als * cs * vxa + sign(alf, bxa) * cf * f4)
|
||
wps = f1 + f2
|
||
wms = f1 - f2
|
||
f1 = dnl * (alf * f3 + als * f5)
|
||
f2 = dnl * (alf * cf * vxa - sign(als, bxa) * cs * f4)
|
||
wpf = f1 + f2
|
||
wmf = f1 - f2
|
||
|
||
f1 = daf
|
||
f2 = sign(das, bxa) * cs
|
||
f3 = als * csq
|
||
rm(1,1) = f1
|
||
rm(2,1) = f1 * lm(1)
|
||
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 * lm(9)
|
||
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)
|
||
|
||
f1 = sign(dnl * sqrt(dna), vxa)
|
||
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)
|
||
|
||
f1 = das
|
||
f2 = sign(daf, bxa) * cf
|
||
f3 = - alf * csq
|
||
rm(1,3) = f1
|
||
rm(2,3) = f1 * lm(3)
|
||
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 * lm(7)
|
||
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)
|
||
|
||
f1 = sign(1.0d+00, vxa)
|
||
rm(5,4) = bxa + bpa * f1
|
||
rm(5,6) = bxa - bpa * f1
|
||
rm(9,4) = f1
|
||
rm(9,6) = - f1
|
||
|
||
rm(2,5) = vxa
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,5) = vva
|
||
|
||
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) = cglm * bpa
|
||
fi(iby,i) = vxa * bya - bxa * vya
|
||
fi(ibz,i) = vxa * bza - bxa * vza
|
||
fi(ibp,i) = cglm * 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 &
|
||
- cglm * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
|
||
|
||
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
|
||
|
||
!===============================================================================
|
||
!
|
||
end module schemes
|