4096 lines
111 KiB
Fortran
4096 lines
111 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
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_adi_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index
|
||
use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x
|
||
real(kind=8), dimension(5,5), save :: rm
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
|
||
|
||
integer :: i, p
|
||
real(kind=8) :: bl, br, btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra
|
||
real(kind=8) :: csn, ent
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm, lv
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adi_m1 / adiabatic_index
|
||
|
||
rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||
rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
! various parameters
|
||
!
|
||
bl = ql(idn,i) / ql(ipr,i)
|
||
br = qr(idn,i) / qr(ipr,i)
|
||
|
||
! the intermediate state averages
|
||
!
|
||
btl = lmean(bl, br)
|
||
bta = amean(bl, br)
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
pra = dna / bta
|
||
prl = dnl / btl
|
||
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
ent = 1.0d+00 / adi_m1x / btl + vva
|
||
|
||
! the difference in entropy variables between states
|
||
!
|
||
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
|
||
+ (log(br) - log(bl)) / adi_m1 &
|
||
- 5.0d-01 * (br * sum(qr(ivx:ivz,i)**2) &
|
||
- bl * sum(ql(ivx:ivz,i)**2))
|
||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||
v(ipr) = bl - br
|
||
|
||
! the eigenvalues
|
||
!
|
||
csn = sqrt(adiabatic_index * pra / dnl)
|
||
|
||
! the diagonal matrix of eigenvalues
|
||
!
|
||
lm(1) = vxa + csn
|
||
lm(2) = vxa
|
||
lm(3) = vxa
|
||
lm(4) = vxa
|
||
lm(5) = vxa - csn
|
||
|
||
! the scaling diagonal matrix
|
||
!
|
||
tm(1) = 5.0d-01 * dnl / adiabatic_index
|
||
tm(2) = dnl * adi_m1x
|
||
tm(3) = pra
|
||
tm(4) = pra
|
||
tm(5) = tm(1)
|
||
|
||
! the average of the right eigenvector matrix
|
||
!
|
||
rm(2,1) = vxa + csn
|
||
rm(3,1) = vya
|
||
rm(4,1) = vza
|
||
rm(2,2) = vxa
|
||
rm(3,2) = vya
|
||
rm(4,2) = vza
|
||
rm(2,5) = vxa - csn
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,:) = [ ent + vxa * csn, vva, vya, vza, ent - vxa * csn ]
|
||
|
||
! KEPEC flux
|
||
!
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa + pra
|
||
fi(imy,i) = fi(idn,i) * vya
|
||
fi(imz,i) = fi(idn,i) * vza
|
||
fi(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
|
||
|
||
! KEPES correction
|
||
!
|
||
lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
|
||
if (vxa >= 0.0d+00) then
|
||
do p = 1, nf
|
||
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
|
||
end do
|
||
else
|
||
do p = nf, 1, -1
|
||
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
|
||
end do
|
||
end if
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_adi_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! RELATIVISTIC ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_SRHD_HLLC:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
|
||
! by Mignone & Bodo.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A. & Bodo, G.
|
||
! "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society,
|
||
! 2005, Volume 364, Pages 126-136
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_srhd_hllc(ql, qr, fi)
|
||
|
||
use algebra , only : quadratic
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ien
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, nr
|
||
real(kind=8) :: sl, sr, srml, sm
|
||
real(kind=8) :: pr, dv
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
|
||
real(kind=8), dimension(3) :: a
|
||
real(kind=8), dimension(2) :: x
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
uh(:) = ( wr(:) - wl(:)) / srml
|
||
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
wl(ien) = wl(ien) + wl(idn)
|
||
wr(ien) = wr(ien) + wr(idn)
|
||
|
||
! prepare the quadratic coefficients (eq. 18 in [1])
|
||
!
|
||
a(1) = uh(imx)
|
||
a(2) = - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn)
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! calculate total pressure (eq. 17 in [1])
|
||
!
|
||
pr = fh(imx) - (fh(ien) + fh(idn)) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pr <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 16 in [1])
|
||
!
|
||
dv = sl - sm
|
||
us(idn) = wl(idn) / dv
|
||
us(imy) = wl(imy) / dv
|
||
us(imz) = wl(imz) / dv
|
||
us(ien) = (wl(ien) + pr * sm) / dv
|
||
us(imx) = (us(ien) + pr) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 16 in [1])
|
||
!
|
||
dv = sr - sm
|
||
us(idn) = wr(idn) / dv
|
||
us(imy) = wr(imy) / dv
|
||
us(imz) = wr(imz) / dv
|
||
us(ien) = (wr(ien) + pr * sm) / dv
|
||
us(imx) = (us(ien) + pr) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity and all fluxes
|
||
! except the parallel momentum one are zero
|
||
!
|
||
fi(idn,i) = 0.0d+00
|
||
fi(imx,i) = pr
|
||
fi(imy,i) = 0.0d+00
|
||
fi(imz,i) = 0.0d+00
|
||
fi(ien,i) = 0.0d+00
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_srhd_hllc
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ISOTHERMAL MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_ISO_MHD_HLLD:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the modified
|
||
! isothermal HLLD method by Mignone. The difference is in allowing a density
|
||
! jumps across the Alfvén waves, due to the inner intermediate states being
|
||
! an average of the states between two Alfvén waves, which also include
|
||
! two slow magnetosonic waves. Moreover, due to the different left and right
|
||
! intermediate state densities the left and right Alfvén speeds might be
|
||
! different too.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A.,
|
||
! "A simple and accurate Riemann solver for isothermal MHD",
|
||
! Journal of Computational Physics, 2007, 225, pp. 1427-1441,
|
||
! https://doi.org/10.1016/j.jcp.2007.01.033
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_hlld(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||
real(kind=8) :: bx, b2, dvl, dvr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, uil, uir
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
! constant advection speed across the intermediate states
|
||
!
|
||
sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn))
|
||
|
||
! speed differences
|
||
!
|
||
srml = sr - sl
|
||
slmm = sl - sm
|
||
srmm = sr - sm
|
||
|
||
! set variables which do not change across the Riemann fan
|
||
!
|
||
uil(idn) = wl(idn) / slmm
|
||
uir(idn) = wr(idn) / srmm
|
||
uil(imx) = uil(idn) * sm
|
||
uir(imx) = uir(idn) * sm
|
||
uil(ibx) = bx
|
||
uir(ibx) = bx
|
||
uil(ibp) = ql(ibp,i)
|
||
uir(ibp) = qr(ibp,i)
|
||
|
||
! Alfvén speeds
|
||
!
|
||
sml = sm - abs(bx) / sqrt(uil(idn))
|
||
smr = sm + abs(bx) / sqrt(uir(idn))
|
||
|
||
! calculate denominators in order to check degeneracy
|
||
!
|
||
dvl = slmm * wl(idn) - b2
|
||
dvr = srmm * wr(idn) - b2
|
||
|
||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||
!
|
||
if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then
|
||
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
fi(:,i) = sl * uil(:) - wl(:)
|
||
|
||
else if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
fi(:,i) = sr * uir(:) - wr(:)
|
||
|
||
else ! sl* <= 0 <= sr*
|
||
|
||
if (smr > sml) then
|
||
|
||
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||
|
||
dvl = smr - sml
|
||
fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dvl
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sl* <= 0 <= sr*
|
||
|
||
else if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr
|
||
|
||
uil(imy) = uil(idn) * (slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
uil(imz) = uil(idn) * (slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
fi(:,i) = sl * uil(:) - wl(:)
|
||
|
||
else if (sr > sml) then
|
||
|
||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||
|
||
dvl = sr - sml
|
||
fi(:,i) = (sml * wr(:) - sr * wcl(:)) / dvl
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sl* < 0
|
||
|
||
else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl
|
||
|
||
uir(imy) = uir(idn) * (srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
uir(imz) = uir(idn) * (srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
fi(:,i) = sr * uir(:) - wr(:)
|
||
|
||
else if (smr > sl) then
|
||
|
||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||
|
||
dvr = smr - sl
|
||
fi(:,i) = (sl * wcr(:) - smr * wl(:)) / dvr
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sr* > 0
|
||
|
||
else ! sl* < sl & sr* > sr
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_hlld
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ISO_ROE:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Stone, J. M. & Gardiner, T. A.,
|
||
! "ATHENA: A New Code for Astrophysical MHD",
|
||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed, eigensystem_roe
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, p
|
||
real(kind=8) :: sdl, sdr, sds
|
||
real(kind=8) :: xx, yy
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: qi, ci, al
|
||
real(kind=8), dimension(nf,nf) :: li, ri
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
! calculate variables for the Roe intermediate state
|
||
!
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
! prepare the Roe intermediate state vector (eq. 11.60 in [2])
|
||
!
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
! prepare coefficients
|
||
!
|
||
xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 &
|
||
+ (ql(ibz,i) - qr(ibz,i))**2) / sds**2
|
||
yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
! obtain eigenvalues and eigenvectors
|
||
!
|
||
call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))
|
||
|
||
! return upwind fluxes
|
||
!
|
||
if (ci(1) >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (ci(nf) <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else
|
||
|
||
! prepare fluxes which do not change across the states
|
||
!
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
! calculate wave amplitudes α = L.ΔU
|
||
!
|
||
al(:) = 0.0d+00
|
||
do p = 1, nf
|
||
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
|
||
end do
|
||
|
||
! calculate the flux average
|
||
!
|
||
fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
|
||
fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
|
||
fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
|
||
fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
|
||
fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
|
||
fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))
|
||
|
||
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
|
||
!
|
||
if (qi(ivx) >= 0.0d+00) then
|
||
do p = 1, nf
|
||
xx = - 0.5d+00 * al(p) * abs(ci(p))
|
||
|
||
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
|
||
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
|
||
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
|
||
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
|
||
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
|
||
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
|
||
end do
|
||
else
|
||
do p = nf, 1, -1
|
||
xx = - 0.5d+00 * al(p) * abs(ci(p))
|
||
|
||
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
|
||
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
|
||
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
|
||
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
|
||
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
|
||
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
|
||
end do
|
||
end if
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ISO_KEPES:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine solves one dimensional isothermal magnetohydrodynamic
|
||
! Riemann problem using the entropy stable KEPES method. The method is
|
||
! a modification of the method described in [1] for the isothermal equation
|
||
! of state.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
|
||
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
|
||
! field divergence diminishing ideal magnetohydrodynamics equations",
|
||
! Journal of Computational Physics, 2018, 364, pp. 420-467,
|
||
! https://doi.org/10.1016/j.jcp.2018.03.002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, csnd, csnd2, ch => cglm
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
ibx, iby, ibz, ibp
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8), dimension(8,8), save :: rm
|
||
!$omp threadprivate(first, rm)
|
||
|
||
integer :: i, p
|
||
real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
|
||
real(kind=8) :: dnl, sqd, bp
|
||
real(kind=8) :: bb2, bp2, b1, b2, b3, bb
|
||
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
|
||
real(kind=8) :: alf2, als2, alf, als
|
||
real(kind=8) :: f1, f2, f3, sgnb1, v2l, v2r, b2l, b2r
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm, lv
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
rm(:,:) = 0.0d+00
|
||
rm(5,4) = 1.0d+00
|
||
rm(8,4) = 1.0d+00
|
||
rm(5,5) = 1.0d+00
|
||
rm(8,5) = - 1.0d+00
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
! various parameters
|
||
!
|
||
v2l = sum(ql(ivx:ivz,i)**2)
|
||
v2r = sum(qr(ivx:ivz,i)**2)
|
||
b2l = sum(ql(ibx:ibz,i)**2)
|
||
b2r = sum(qr(ibx:ibz,i)**2)
|
||
|
||
! the difference in entropy variables between the states
|
||
!
|
||
v(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) &
|
||
- 5.0d-01 * (v2r - v2l)
|
||
v(ivx) = qr(ivx,i) - ql(ivx,i)
|
||
v(ivy) = qr(ivy,i) - ql(ivy,i)
|
||
v(ivz) = qr(ivz,i) - ql(ivz,i)
|
||
v(ibx) = qr(ibx,i) - ql(ibx,i)
|
||
v(iby) = qr(iby,i) - ql(iby,i)
|
||
v(ibz) = qr(ibz,i) - ql(ibz,i)
|
||
v(ibp) = qr(ibp,i) - ql(ibp,i)
|
||
|
||
! the intermediate state averages
|
||
!
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
bxa = amean(ql(ibx,i), qr(ibx,i))
|
||
bya = amean(ql(iby,i), qr(iby,i))
|
||
bza = amean(ql(ibz,i), qr(ibz,i))
|
||
bpa = amean(ql(ibp,i), qr(ibp,i))
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
pta = csnd2 * dna + 5.0d-01 * amean(b2l, b2r)
|
||
|
||
! the eigenvalues and related parameters
|
||
!
|
||
sqd = sqrt(dnl)
|
||
b1 = bxa / sqd
|
||
b2 = bya / sqd
|
||
b3 = bza / sqd
|
||
bp2 = b2**2 + b3**2
|
||
bp = sqrt(bp2)
|
||
bb2 = b1**2 + bp2
|
||
bb = sqrt(bb2)
|
||
sgnb1 = sign(1.0d+00, b1)
|
||
|
||
if (bp > 0.0d+00) then
|
||
x2 = b2 / bp
|
||
x3 = b3 / bp
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
ca2 = b1**2
|
||
ca = sqrt(ca2)
|
||
f1 = csnd2 + bb2
|
||
f3 = max(0.0d+00, f1**2 - 4.0d+00 * csnd2 * ca2)
|
||
f2 = sqrt(f3)
|
||
cf2 = 5.0d-01 * (f1 + f2)
|
||
cf = sqrt(cf2)
|
||
cs2 = 5.0d-01 * (f1 - f2)
|
||
cs = sqrt(cs2)
|
||
|
||
f2 = cf2 - cs2
|
||
if (f2 > 0.0d+00) then
|
||
alf2 = max(0.0d+00, csnd2 - cs2) / f2
|
||
als2 = max(0.0d+00, cf2 - csnd2) / f2
|
||
else
|
||
alf2 = 0.0d+00
|
||
als2 = 0.0d+00
|
||
end if
|
||
alf = sqrt(alf2)
|
||
als = sqrt(als2)
|
||
|
||
! === the average of the right eigenvector matrix ===
|
||
!
|
||
! right (1) and left (9) fast magnetosonic eigenvectors
|
||
!
|
||
f1 = dnl * alf
|
||
f2 = dnl * als * cs * sgnb1
|
||
f3 = als * csnd * sqd
|
||
rm(1,1) = f1
|
||
rm(2,1) = f1 * (vxa + cf)
|
||
rm(3,1) = f1 * vya - f2 * x2
|
||
rm(4,1) = f1 * vza - f2 * x3
|
||
rm(6,1) = f3 * x2
|
||
rm(7,1) = f3 * x3
|
||
|
||
rm(1,8) = f1
|
||
rm(2,8) = f1 * (vxa - cf)
|
||
rm(3,8) = f1 * vya + f2 * x2
|
||
rm(4,8) = f1 * vza + f2 * x3
|
||
rm(6,8) = rm(6,1)
|
||
rm(7,8) = rm(7,1)
|
||
|
||
! right (2) and left (8) Alfvénic eigenvectors
|
||
!
|
||
f1 = dnl * sqrt(dna)
|
||
rm(3,2) = f1 * x3
|
||
rm(4,2) = - f1 * x2
|
||
rm(6,2) = - dnl * x3
|
||
rm(7,2) = dnl * x2
|
||
|
||
rm(3,7) = - rm(3,2)
|
||
rm(4,7) = - rm(4,2)
|
||
rm(6,7) = rm(6,2)
|
||
rm(7,7) = rm(7,2)
|
||
|
||
! right (3) and left (7) slow magnetosonic eigenvectors
|
||
!
|
||
f1 = dnl * als
|
||
f2 = dnl * alf * cf * sgnb1
|
||
f3 = - alf * csnd * sqd
|
||
rm(1,3) = f1
|
||
rm(2,3) = f1 * (vxa + cs)
|
||
rm(3,3) = f1 * vya + f2 * x2
|
||
rm(4,3) = f1 * vza + f2 * x3
|
||
rm(6,3) = f3 * x2
|
||
rm(7,3) = f3 * x3
|
||
|
||
rm(1,6) = f1
|
||
rm(2,6) = f1 * (vxa - cs)
|
||
rm(3,6) = f1 * vya - f2 * x2
|
||
rm(4,6) = f1 * vza - f2 * x3
|
||
rm(6,6) = rm(6,3)
|
||
rm(7,6) = rm(7,3)
|
||
|
||
! the eigenvalue diagonal matrix
|
||
!
|
||
lm(1) = vxa + cf
|
||
lm(2) = vxa + ca
|
||
lm(3) = vxa + cs
|
||
lm(4) = vxa + ch
|
||
lm(5) = vxa - ch
|
||
lm(6) = vxa - cs
|
||
lm(7) = vxa - ca
|
||
lm(8) = vxa - cf
|
||
|
||
! the scaling diagonal matrix
|
||
!
|
||
tm(1) = 2.5d-01 / (dnl * csnd2)
|
||
tm(2) = 2.5d-01 / dnl**2
|
||
tm(3) = tm(1)
|
||
tm(4) = 2.5d-01
|
||
tm(5) = tm(4)
|
||
tm(6) = tm(1)
|
||
tm(7) = tm(2)
|
||
tm(8) = tm(1)
|
||
|
||
! KEPEC flux
|
||
!
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
|
||
fi(imy,i) = fi(idn,i) * vya - bxa * bya
|
||
fi(imz,i) = fi(idn,i) * vza - bxa * bza
|
||
fi(ibx,i) = ch * bpa
|
||
fi(iby,i) = vxa * bya - bxa * vya
|
||
fi(ibz,i) = vxa * bza - bxa * vza
|
||
fi(ibp,i) = ch * bxa
|
||
|
||
! KEPES correction
|
||
!
|
||
lv = abs(lm) * tm * matmul(v, rm)
|
||
if (vxa >= 0.0d+00) then
|
||
do p = 1, nf
|
||
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
|
||
end do
|
||
else
|
||
do p = nf, 1, -1
|
||
fi(:,i) = fi(:,i) - lv(p) * rm(:,p)
|
||
end do
|
||
end if
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_HLLC:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the HLLC method,
|
||
! by Toro. In the HLLC method the tangential components of the velocity are
|
||
! discontinuous actoss the contact dicontinuity.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Toro, E. F., Spruce, M., & Speares, W.
|
||
! "Restoration of the contact surface in the HLL-Riemann solver",
|
||
! Shock Waves, 1994, Volume 4, Issue 1, pp. 25-34
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_hllc(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivy, ivz, imx, imy, imz, ien, &
|
||
ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2
|
||
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, ui
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
dn = wr(idn) - wl(idn)
|
||
sm = (wr(imx) - wl(imx)) / dn
|
||
|
||
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
|
||
|
||
if (b2 > 0.0d+00) then
|
||
|
||
vy = (wr(imy) - wl(imy)) / dn
|
||
vz = (wr(imz) - wl(imz)) / dn
|
||
by = (wr(iby) - wl(iby)) / srml
|
||
bz = (wr(ibz) - wl(ibz)) / srml
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
slmm = sl - sm
|
||
|
||
ui(idn) = wl(idn) / slmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * vy
|
||
ui(imz) = ui(idn) * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ql(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
srmm = sr - sm
|
||
|
||
ui(idn) = wr(idn) / srmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * vy
|
||
ui(imz) = ui(idn) * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = qr(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
|
||
! to the HLL one
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sm = 0
|
||
|
||
else ! Bₓ = 0
|
||
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
slmm = sl - sm
|
||
|
||
ui(idn) = wl(idn) / slmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * ql(ivy,i)
|
||
ui(imz) = ui(idn) * ql(ivz,i)
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = wl(iby) / slmm
|
||
ui(ibz) = wl(ibz) / slmm
|
||
ui(ibp) = ql(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt) / slmm
|
||
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
srmm = sr - sm
|
||
|
||
ui(idn) = wr(idn) / srmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * qr(ivy,i)
|
||
ui(imz) = ui(idn) * qr(ivz,i)
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = wr(iby) / srmm
|
||
ui(ibz) = wr(ibz) / srmm
|
||
ui(ibp) = qr(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt) / srmm
|
||
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
|
||
! to the HLL one
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sm = 0
|
||
|
||
end if ! Bₓ = 0
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_hllc
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_ADI_MHD_HLLD:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the adiabatic HLLD
|
||
! method described by Miyoshi & Kusano.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Miyoshi, T. & Kusano, K.,
|
||
! "A multi-state HLL approximate Riemann solver for ideal
|
||
! magnetohydrodynamics",
|
||
! Journal of Computational Physics, 2005, 208, pp. 315-344
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_hlld(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
|
||
real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb
|
||
real(kind=8) :: dnl, dnr, cal, car, sml, smr
|
||
real(kind=8) :: dv, dvl, dvr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, ui
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then ! sl ≥ 0
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then ! sr ≤ 0
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
dn = wr(idn) - wl(idn)
|
||
sm = (wr(imx) - wl(imx)) / dn
|
||
|
||
slmm = sl - sm
|
||
srmm = sr - sm
|
||
|
||
dnl = wl(idn) / slmm
|
||
dnr = wr(idn) / srmm
|
||
|
||
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
|
||
|
||
cal = abs(bx) / sqrt(dnl)
|
||
car = abs(bx) / sqrt(dnr)
|
||
sml = sm - cal
|
||
smr = sm + car
|
||
|
||
! calculate division factors (also used to determine degeneracies)
|
||
!
|
||
dvl = slmm * wl(idn) - b2
|
||
dvr = srmm * wr(idn) - b2
|
||
|
||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||
!
|
||
if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracy
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sl* < 0 < sr*
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||
|
||
! prepare constant primitive variables of the intermediate states
|
||
!
|
||
dv = abs(bx) * (sqrt(dnr) + sqrt(dnl))
|
||
vy = (wcr(imy) - wcl(imy)) / dv
|
||
vz = (wcr(imz) - wcl(imz)) / dv
|
||
dv = car + cal
|
||
by = (wcr(iby) - wcl(iby)) / dv
|
||
bz = (wcr(ibz) - wcl(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
! conservative variables for the inmost left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
|
||
|
||
! the inmost left intermediate flux
|
||
!
|
||
fi(:,i) = sml * ui(:) - wcl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
! conservative variables for the inmost right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
|
||
|
||
! the inmost right intermediate flux
|
||
!
|
||
fi(:,i) = smr * ui(:) - wcr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
|
||
! the flux can be averaged from the Alfvén waves using a simple HLL formula;
|
||
!
|
||
fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dv
|
||
|
||
end if ! sm = 0
|
||
|
||
else ! Bx = 0
|
||
|
||
if (sm > 0.0d+00) then
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = slmm * wl(imy) / dvl
|
||
vz = slmm * wl(imz) / dvl
|
||
by = wl(idn) * wl(iby) / dvl
|
||
bz = wl(idn) * wl(ibz) / dvl
|
||
vb = vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt) / slmm
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy)) / dvr
|
||
vz = ( srmm * wr(imz)) / dvr
|
||
by = (wr(idn) * wr(iby)) / dvr
|
||
bz = (wr(idn) * wr(ibz)) / dvr
|
||
vb = vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt) / srmm
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! Sm = 0
|
||
|
||
! since Bx = 0 and Sm = 0, then revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! Bx = 0
|
||
|
||
end if ! sl* < 0 < sr*
|
||
|
||
else ! some degeneracies
|
||
|
||
! separate degeneracies
|
||
!
|
||
if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else ! sl* <= 0
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||
|
||
! primitive variables for the inner left intermediate state
|
||
!
|
||
dv = srmm * dnr + cal * dnl
|
||
vy = (wr(imy) - wcl(imy)) / dv
|
||
vz = (wr(imz) - wcl(imz)) / dv
|
||
dv = sr - sml
|
||
by = (wr(iby) - wcl(iby)) / dv
|
||
bz = (wr(ibz) - wcl(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the inner left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm >= 0.0d+00) then ! sm >= 0
|
||
|
||
! the inner left intermediate flux
|
||
!
|
||
fi(:,i) = sml * ui(:) - wcl(:)
|
||
|
||
else ! sm < 0
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcr(:) = (sm - sml) * ui(:) + wcl(:)
|
||
|
||
! calculate the average flux over the right inner intermediate state
|
||
!
|
||
fi(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm
|
||
|
||
end if ! sm < 0
|
||
|
||
else ! bx = 0
|
||
|
||
! no Alfvén wave, so revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sl* < 0
|
||
|
||
else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sr* > 0
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||
|
||
! primitive variables for the inner right intermediate state
|
||
!
|
||
dv = slmm * dnl - car * dnr
|
||
vy = (wl(imy) - wcr(imy)) / dv
|
||
vz = (wl(imz) - wcr(imz)) / dv
|
||
dv = sl - smr
|
||
by = (wl(iby) - wcr(iby)) / dv
|
||
bz = (wl(ibz) - wcr(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the inner left intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm <= 0.0d+00) then ! sm <= 0
|
||
|
||
! the inner right intermediate flux
|
||
!
|
||
fi(:,i) = smr * ui(:) - wcr(:)
|
||
|
||
else ! sm > 0
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcl(:) = (sm - smr) * ui(:) + wcr(:)
|
||
|
||
! calculate the average flux over the left inner intermediate state
|
||
!
|
||
fi(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm
|
||
|
||
end if ! sm > 0
|
||
|
||
else ! bx = 0
|
||
|
||
! no Alfvén wave, so revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sr* > 0
|
||
|
||
else ! sl* < sl & sr* > sr
|
||
|
||
! so far we revert to HLL flux in the case of degeneracies
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sl* < sl & sr* > sr
|
||
|
||
end if ! deneneracies
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_hlld
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ADI_ROE:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Stone, J. M. & Gardiner, T. A.,
|
||
! "ATHENA: A New Code for Astrophysical MHD",
|
||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ipr, ien, &
|
||
ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed, eigensystem_roe
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, p
|
||
real(kind=8) :: sdl, sdr, sds
|
||
real(kind=8) :: xx, yy
|
||
real(kind=8) :: pml, pmr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: qi, ci, al
|
||
real(kind=8), dimension(nf,nf) :: li, ri
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
! calculate variables for the Roe intermediate state
|
||
!
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
! prepare magnetic pressures
|
||
!
|
||
pml = 0.5d+00 * sum(ql(ibx:ibz,i)**2)
|
||
pmr = 0.5d+00 * sum(qr(ibx:ibz,i)**2)
|
||
|
||
! prepare the Roe intermediate state vector (eq. 11.60 in [2])
|
||
!
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds
|
||
qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl &
|
||
+ (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby:ibz,i) / sds
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
! prepare coefficients
|
||
!
|
||
xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 &
|
||
+ (ql(ibz,i) - qr(ibz,i))**2) / sds**2
|
||
yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
! obtain eigenvalues and eigenvectors
|
||
!
|
||
call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:))
|
||
|
||
! return upwind fluxes
|
||
!
|
||
if (ci(1) >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (ci(nf) <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else
|
||
|
||
! prepare fluxes which do not change across the states
|
||
!
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
! calculate wave amplitudes α = L.ΔU
|
||
!
|
||
al(:) = 0.0d+00
|
||
do p = 1, nf
|
||
al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i))
|
||
end do
|
||
|
||
! calculate the flux average
|
||
!
|
||
fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i))
|
||
fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i))
|
||
fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i))
|
||
fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i))
|
||
fi(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i))
|
||
fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i))
|
||
fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i))
|
||
|
||
! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K)
|
||
!
|
||
if (qi(ivx) >= 0.0d+00) then
|
||
do p = 1, nf
|
||
xx = - 0.5d+00 * al(p) * abs(ci(p))
|
||
|
||
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
|
||
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
|
||
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
|
||
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
|
||
fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
|
||
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
|
||
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
|
||
end do
|
||
else
|
||
do p = nf, 1, -1
|
||
xx = - 0.5d+00 * al(p) * abs(ci(p))
|
||
|
||
fi(idn,i) = fi(idn,i) + xx * ri(p,idn)
|
||
fi(imx,i) = fi(imx,i) + xx * ri(p,imx)
|
||
fi(imy,i) = fi(imy,i) + xx * ri(p,imy)
|
||
fi(imz,i) = fi(imz,i) + xx * ri(p,imz)
|
||
fi(ien,i) = fi(ien,i) + xx * ri(p,ien)
|
||
fi(iby,i) = fi(iby,i) + xx * ri(p,iby)
|
||
fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz)
|
||
end do
|
||
end if
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ADI_KEPES:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine solves one dimensional adiabatic magnetohydrodynamic
|
||
! Riemann problem using the entropy stable KEPES method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
|
||
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
|
||
! field divergence diminishing ideal magnetohydrodynamics equations",
|
||
! Journal of Computational Physics, 2018, 364, pp. 420-467,
|
||
! https://doi.org/10.1016/j.jcp.2018.03.002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index, ch => cglm
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
ibx, iby, ibz, ibp, ipr, ien
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x
|
||
real(kind=8), dimension(9,9), save :: rm
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
|
||
|
||
integer :: i
|
||
real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa
|
||
real(kind=8) :: dnl, prl, pta, vva, br, bl, bp
|
||
real(kind=8) :: btl, bta, eka, ema, ub2, uba
|
||
real(kind=8) :: aa2, al2, ab2, bb2, bp2, aa, al, ab, b1, b2, b3, bb
|
||
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
|
||
real(kind=8) :: alf2, als2, alf, als, wps, wms, wpf, wmf
|
||
real(kind=8) :: f1, f2, f3, f4, sgnb1, v2l, v2r, b2l, b2r
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adi_m1 / adiabatic_index
|
||
|
||
rm(:,:) = 0.0d+00
|
||
rm(6,4) = 1.0d+00
|
||
rm(9,4) = 1.0d+00
|
||
rm(1,5) = 1.0d+00
|
||
rm(6,6) = 1.0d+00
|
||
rm(9,6) = - 1.0d+00
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
! the difference in entropy variables between the states
|
||
!
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
|
||
b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
|
||
|
||
bl = ql(idn,i) / ql(ipr,i)
|
||
br = qr(idn,i) / qr(ipr,i)
|
||
|
||
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
|
||
+ (log(br) - log(bl)) / adi_m1 &
|
||
- 5.0d-01 * (br * v2r - bl * v2l)
|
||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||
v(ipr) = bl - br
|
||
v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
|
||
v(iby) = br * qr(iby,i) - bl * ql(iby,i)
|
||
v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
|
||
v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)
|
||
|
||
! the intermediate state averages
|
||
!
|
||
btl = lmean(bl, br)
|
||
bta = amean(bl, br)
|
||
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
bxa = amean(ql(ibx,i), qr(ibx,i))
|
||
bya = amean(ql(iby,i), qr(iby,i))
|
||
bza = amean(ql(ibz,i), qr(ibz,i))
|
||
bpa = amean(ql(ibp,i), qr(ibp,i))
|
||
pra = dna / bta
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
prl = lmean(ql(ipr,i), qr(ipr,i))
|
||
eka = 5.0d-01 * amean(v2l, v2r)
|
||
ema = 5.0d-01 * amean(b2l, b2r)
|
||
pta = pra + ema
|
||
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r)
|
||
uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), &
|
||
sum(qr(ivx:ivz,i) * qr(ibx:ibz,i)))
|
||
|
||
! the eigenvalues and related parameters
|
||
!
|
||
aa2 = adiabatic_index * pra / dnl
|
||
aa = sqrt(aa2)
|
||
al2 = adiabatic_index * prl / dnl
|
||
al = sqrt(al2)
|
||
ab2 = adiabatic_index / bta
|
||
ab = sqrt(ab2)
|
||
|
||
b1 = bxa / sqrt(dnl)
|
||
b2 = bya / sqrt(dnl)
|
||
b3 = bza / sqrt(dnl)
|
||
bp2 = b2**2 + b3**2
|
||
bp = sqrt(bp2)
|
||
bb2 = b1**2 + bp2
|
||
bb = sqrt(bb2)
|
||
sgnb1 = sign(1.0d+00, b1)
|
||
|
||
if (bp > 0.0d+00) then
|
||
x2 = b2 / bp
|
||
x3 = b3 / bp
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
ca2 = b1**2
|
||
ca = sqrt(ca2)
|
||
f1 = aa2 + bb2
|
||
f3 = max(0.0d+00, f1**2 - 4.0d+00 * aa2 * ca2)
|
||
f2 = sqrt(f3)
|
||
cf2 = 5.0d-01 * (f1 + f2)
|
||
cf = sqrt(cf2)
|
||
cs2 = 5.0d-01 * (f1 - f2)
|
||
cs = sqrt(cs2)
|
||
|
||
f2 = cf2 - cs2
|
||
if (f2 > 0.0d+00) then
|
||
alf2 = max(0.0d+00, aa2 - cs2) / f2
|
||
als2 = max(0.0d+00, cf2 - aa2) / f2
|
||
else
|
||
alf2 = 0.0d+00
|
||
als2 = 0.0d+00
|
||
end if
|
||
alf = sqrt(alf2)
|
||
als = sqrt(als2)
|
||
|
||
f3 = vva + al2 / adi_m1
|
||
f4 = sgnb1 * (vya * x2 + vza * x3)
|
||
f1 = dnl * (als * f3 - alf * ab * bp)
|
||
f2 = dnl * (als * cs * vxa + alf * cf * f4)
|
||
wps = f1 + f2
|
||
wms = f1 - f2
|
||
f1 = dnl * (alf * f3 + als * ab * bp)
|
||
f2 = dnl * (alf * cf * vxa - als * cs * f4)
|
||
wpf = f1 + f2
|
||
wmf = f1 - f2
|
||
|
||
! the scaling diagonal matrix
|
||
!
|
||
f1 = 5.0d-01 / bta
|
||
tm(1) = 5.0d-01 / (adiabatic_index * dnl)
|
||
tm(2) = f1 / dnl**2
|
||
tm(3) = tm(1)
|
||
tm(4) = f1
|
||
tm(5) = dnl * adi_m1x
|
||
tm(6) = tm(4)
|
||
tm(7) = tm(1)
|
||
tm(8) = tm(2)
|
||
tm(9) = tm(1)
|
||
|
||
! === the average of the right eigenvector matrix ===
|
||
!
|
||
! right (1) and left (9) fast magnetosonic eigenvectors
|
||
!
|
||
f1 = dnl * alf
|
||
f2 = dnl * als * cs * sgnb1
|
||
f3 = als * ab * sqrt(dnl)
|
||
rm(1,1) = f1
|
||
rm(2,1) = f1 * (vxa + cf)
|
||
rm(3,1) = f1 * vya - f2 * x2
|
||
rm(4,1) = f1 * vza - f2 * x3
|
||
rm(5,1) = wpf
|
||
rm(7,1) = f3 * x2
|
||
rm(8,1) = f3 * x3
|
||
|
||
rm(1,9) = f1
|
||
rm(2,9) = f1 * (vxa - cf)
|
||
rm(3,9) = f1 * vya + f2 * x2
|
||
rm(4,9) = f1 * vza + f2 * x3
|
||
rm(5,9) = wmf
|
||
rm(7,9) = rm(7,1)
|
||
rm(8,9) = rm(8,1)
|
||
|
||
! right (2) and left (8) Alfvénic eigenvectors
|
||
!
|
||
f1 = dnl * sqrt(dna)
|
||
rm(3,2) = f1 * x3
|
||
rm(4,2) = - f1 * x2
|
||
rm(5,2) = - f1 * (x2 * vza - x3 * vya)
|
||
rm(7,2) = - dnl * x3
|
||
rm(8,2) = dnl * x2
|
||
|
||
rm(3,8) = - rm(3,2)
|
||
rm(4,8) = - rm(4,2)
|
||
rm(5,8) = - rm(5,2)
|
||
rm(7,8) = rm(7,2)
|
||
rm(8,8) = rm(8,2)
|
||
|
||
! right (3) and left (7) slow magnetosonic eigenvectors
|
||
!
|
||
f1 = dnl * als
|
||
f2 = dnl * alf * cf * sgnb1
|
||
f3 = - alf * ab * sqrt(dnl)
|
||
rm(1,3) = f1
|
||
rm(2,3) = f1 * (vxa + cs)
|
||
rm(3,3) = f1 * vya + f2 * x2
|
||
rm(4,3) = f1 * vza + f2 * x3
|
||
rm(5,3) = wps
|
||
rm(7,3) = f3 * x2
|
||
rm(8,3) = f3 * x3
|
||
|
||
rm(1,7) = f1
|
||
rm(2,7) = f1 * (vxa - cs)
|
||
rm(3,7) = f1 * vya - f2 * x2
|
||
rm(4,7) = f1 * vza - f2 * x3
|
||
rm(5,7) = wms
|
||
rm(7,7) = rm(7,3)
|
||
rm(8,7) = rm(8,3)
|
||
|
||
! right (4) and left (6) GLM eigenvectors
|
||
!
|
||
rm(5,4) = bxa + bpa
|
||
rm(5,6) = bxa - bpa
|
||
|
||
! entropy (5) eigenvectors
|
||
!
|
||
rm(2,5) = vxa
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,5) = vva
|
||
|
||
! the eigenvalue diagonal matrix
|
||
!
|
||
lm(1) = vxa + cf
|
||
lm(2) = vxa + ca
|
||
lm(3) = vxa + cs
|
||
lm(4) = vxa + ch
|
||
lm(5) = vxa
|
||
lm(6) = vxa - ch
|
||
lm(7) = vxa - cs
|
||
lm(8) = vxa - ca
|
||
lm(9) = vxa - cf
|
||
|
||
! KEPEC flux
|
||
!
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
|
||
fi(imy,i) = fi(idn,i) * vya - bxa * bya
|
||
fi(imz,i) = fi(idn,i) * vza - bxa * bza
|
||
fi(ibx,i) = ch * bpa
|
||
fi(iby,i) = vxa * bya - bxa * vya
|
||
fi(ibz,i) = vxa * bza - bxa * vza
|
||
fi(ibp,i) = ch * bxa
|
||
fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
|
||
+ (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) &
|
||
+ (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) &
|
||
+ fi(ibp,i) * bpa - ub2 + bxa * uba &
|
||
- ch * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
|
||
|
||
! KEPES correction
|
||
!
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! RELATIVISTIC ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_SRMHD_HLLC:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
|
||
! by Mignone & Bodo.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A. & Bodo, G.
|
||
! "An HLLC Riemann solver for relativistic flows - II.
|
||
! Magnetohydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society,
|
||
! 2006, Volume 368, Pages 1040-1054
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_srmhd_hllc(ql, qr, fi)
|
||
|
||
use algebra , only : quadratic
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, nr
|
||
real(kind=8) :: sl, sr, srml, sm
|
||
real(kind=8) :: bx2
|
||
real(kind=8) :: pt, dv, fc, uu, ff, uf
|
||
real(kind=8) :: vv, vb, gi
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
|
||
real(kind=8), dimension(3) :: a, vs
|
||
real(kind=8), dimension(2) :: x
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
uh(:) = ( wr(:) - wl(:)) / srml
|
||
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
wl(ien) = wl(ien) + wl(idn)
|
||
wr(ien) = wr(ien) + wr(idn)
|
||
|
||
! calculate Bₓ²
|
||
!
|
||
bx2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
! calculate the contact dicontinuity speed solving eq. (41)
|
||
!
|
||
if (bx2 >= 1.0d-16) then
|
||
|
||
! prepare the quadratic coefficients, (eq. 42 in [1])
|
||
!
|
||
uu = sum(uh(iby:ibz) * uh(iby:ibz))
|
||
ff = sum(fh(iby:ibz) * fh(iby:ibz))
|
||
uf = sum(uh(iby:ibz) * fh(iby:ibz))
|
||
|
||
a(1) = uh(imx) - uf
|
||
a(2) = uu + ff - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn) - uf
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! substitute magnetic field components from the HLL state (eq. 37 in [1])
|
||
!
|
||
us(ibx) = ql(ibx,i)
|
||
us(iby) = uh(iby)
|
||
us(ibz) = uh(ibz)
|
||
us(ibp) = ql(ibp,i)
|
||
|
||
! calculate velocity components without Bₓ normalization (eq. 38 in [1])
|
||
!
|
||
vs(1) = sm
|
||
vs(2) = (us(iby) * sm - fh(iby)) / us(ibx)
|
||
vs(3) = (us(ibz) * sm - fh(ibz)) / us(ibx)
|
||
|
||
! calculate v² and v.B
|
||
!
|
||
vv = sum(vs(1:3) * vs( 1: 3))
|
||
vb = sum(vs(1:3) * us(ibx:ibz))
|
||
|
||
! calculate inverse of Lorentz factor
|
||
!
|
||
gi = 1.0d+00 - vv
|
||
|
||
! calculate total pressure (eq. 40 in [1])
|
||
!
|
||
pt = fh(imx) + gi * bx2 - (fh(ien) + fh(idn) - us(ibx) * vb) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pt <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sl - sm
|
||
fc = (sl - ql(ivx,i)) / dv
|
||
us(idn) = fc * ul(idn,i)
|
||
us(imy) = (sl * ul(imy,i) - fl(imy,i) &
|
||
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
|
||
us(imz) = (sl * ul(imz,i) - fl(imz,i) &
|
||
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
|
||
us(ien) = (sl * (ul(ien,i) + ul(idn,i)) &
|
||
- (fl(ien,i) + fl(idn,i)) &
|
||
+ pt * sm - us(ibx) * vb) / dv - us(idn)
|
||
|
||
! calculate normal component of momentum
|
||
!
|
||
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sr - sm
|
||
fc = (sr - qr(ivx,i)) / dv
|
||
us(idn) = fc * ur(idn,i)
|
||
us(imy) = (sr * ur(imy,i) - fr(imy,i) &
|
||
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
|
||
us(imz) = (sr * ur(imz,i) - fr(imz,i) &
|
||
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
|
||
us(ien) = (sr * (ur(ien,i) + ur(idn,i)) &
|
||
- (fr(ien,i) + fr(idn,i)) &
|
||
+ pt * sm - us(ibx) * vb) / dv - us(idn)
|
||
|
||
! calculate normal component of momentum
|
||
!
|
||
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity so use HLL flux
|
||
!
|
||
fi(:,i) = fh(:)
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
else ! Bx² > 0
|
||
|
||
! prepare the quadratic coefficients (eq. 47 in [1])
|
||
!
|
||
a(1) = uh(imx)
|
||
a(2) = - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn)
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! calculate total pressure (eq. 48 in [1])
|
||
!
|
||
pt = fh(imx) - (fh(ien) + fh(idn)) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pt <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sl - sm
|
||
us(idn) = wl(idn) / dv
|
||
us(imy) = wl(imy) / dv
|
||
us(imz) = wl(imz) / dv
|
||
us(ien) = (wl(ien) + pt * sm) / dv
|
||
us(imx) = (us(ien) + pt) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
us(ibx) = 0.0d+00
|
||
us(iby) = wl(iby) / dv
|
||
us(ibz) = wl(ibz) / dv
|
||
us(ibp) = ql(ibp,i)
|
||
|
||
! calculate the flux (eq. 27 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sr - sm
|
||
us(idn) = wr(idn) / dv
|
||
us(imy) = wr(imy) / dv
|
||
us(imz) = wr(imz) / dv
|
||
us(ien) = (wr(ien) + pt * sm) / dv
|
||
us(imx) = (us(ien) + pt) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
us(ibx) = 0.0d+00
|
||
us(iby) = wr(iby) / dv
|
||
us(ibz) = wr(ibz) / dv
|
||
us(ibp) = qr(ibp,i)
|
||
|
||
! calculate the flux (eq. 27 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity so use HLL flux
|
||
!
|
||
fi(:,i) = fh(:)
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
end if ! Bx² == 0 or > 0
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_srmhd_hllc
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! SUPPORTING SUBROUTINES/FUNCTIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
|
||
! ---------------------------------------
|
||
!
|
||
! Subroutine adds higher order corrections to the numerical fluxes.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! f - the vector of numerical fluxes;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine higher_order_flux_correction(f)
|
||
|
||
! include external procedures
|
||
!
|
||
use interpolations, only : order
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:,:), intent(inout) :: f
|
||
|
||
! local variables
|
||
!
|
||
integer :: n
|
||
|
||
! local arrays to store the states
|
||
!
|
||
real(kind=8), dimension(size(f,1),size(f,2)) :: f2, f4
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! depending on the scheme order calculate and add higher order corrections,
|
||
! if required
|
||
!
|
||
if (.not. high_order_flux .or. order <= 2) return
|
||
|
||
! get the length of vector
|
||
!
|
||
n = size(f, 2)
|
||
|
||
! 2nd order correction
|
||
!
|
||
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
|
||
f2(:,1 ) = 0.0d+00
|
||
f2(:, n ) = 0.0d+00
|
||
|
||
! 4th order correction
|
||
!
|
||
if (order > 4) then
|
||
|
||
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
|
||
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
|
||
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
|
||
f4(:,1 ) = 0.0d+00
|
||
f4(:,2 ) = 0.0d+00
|
||
f4(:, n-1) = 0.0d+00
|
||
f4(:, n ) = 0.0d+00
|
||
|
||
! correct the flux for 4th order
|
||
!
|
||
f(:,:) = f(:,:) + f4(:,:)
|
||
|
||
end if
|
||
|
||
! correct the flux for 2nd order
|
||
!
|
||
f(:,:) = f(:,:) + f2(:,:)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine higher_order_flux_correction
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LMEAN:
|
||
! --------------
|
||
!
|
||
! Function calculates the logarithmic mean using the optimized algorithm
|
||
! by Ranocha et al.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! l, r - the values to calculate the mean of;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ranocha et al.,
|
||
! "Efficient implementation of modern entropy stable and kinetic energy
|
||
! preserving discontinuous Galerkin methods for conservation laws",
|
||
! http://arxiv.org/abs/2112.10517v1
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function lmean(l, r)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: l, r
|
||
|
||
real(kind=8) :: u, d, s
|
||
|
||
real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, &
|
||
c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
d = r - l
|
||
s = r + l
|
||
u = (d / s)**2
|
||
if (u < 1.0d-04) then
|
||
lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
|
||
else
|
||
if (d >= 0.0d+00) then
|
||
s = log(r / l)
|
||
else
|
||
s = log(l / r)
|
||
end if
|
||
lmean = abs(d) / s
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function lmean
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function AMEAN:
|
||
! --------------
|
||
!
|
||
! Function calculate the arithmetic mean.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! l, r - the values to calculate the mean of;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function amean(l, r)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: l, r
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
amean = 5.0d-01 * (l + r)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function amean
|
||
|
||
!===============================================================================
|
||
!
|
||
end module schemes
|