Grzegorz Kowal
fb25489487
The tearing test problem was numerically unstable in the adiabatic case. Setting bty to 1 in the case of br beeing zero solved the issue. Make it consistent for the isothermal ROE solver as well. Additionaly, one typo was fixed in the calculation of the left eigenvectors. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
4145 lines
115 KiB
Fortran
4145 lines
115 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-2023 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 helpers , only : print_message, uppercase
|
||
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"
|
||
|
||
character(len=256) :: msg
|
||
|
||
character(len=*), parameter :: fmt = "('The selected Riemann solver " // &
|
||
"is not implemented for the ',a,': ',a,'.'" // &
|
||
",a,4x,'Available Riemann solvers are: ', a, '.')"
|
||
character(len=*), parameter :: loc = 'SCHEMES::initialize_schemes()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
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")
|
||
|
||
select case(trim(eos))
|
||
|
||
case("iso", "ISO", "isothermal", "ISOTHERMAL")
|
||
|
||
select case(trim(solver))
|
||
|
||
case("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case("roe", "ROE")
|
||
name_sol = "ROE"
|
||
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(msg,fmt) 'isothermal hydrodynamics', &
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'ROE', 'KEPES'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
select case(trim(solver))
|
||
|
||
case("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case("hllc", "HLLC")
|
||
name_sol = "HLLC"
|
||
riemann_solver => riemann_hd_hllc
|
||
|
||
case("roe", "ROE")
|
||
name_sol = "ROE"
|
||
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(msg,fmt) 'adiabatic hydrodynamics', &
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'HLLC', 'ROE', 'KEPES'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
end select
|
||
|
||
!--- MAGNETOHYDRODYNAMICS ---
|
||
!
|
||
case("mhd", "MHD", "magnetohydrodynamic", "MAGNETOHYDRODYNAMIC")
|
||
|
||
select case(trim(eos))
|
||
|
||
case("iso", "ISO", "isothermal", "ISOTHERMAL")
|
||
|
||
select case(trim(solver))
|
||
|
||
case("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case ("hlld", "HLLD")
|
||
name_sol = "HLLD"
|
||
riemann_solver => riemann_mhd_iso_hlld
|
||
|
||
case("roe", "ROE")
|
||
name_sol = "ROE"
|
||
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(msg,fmt) 'isothermal magnetohydrodynamics', &
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'HLLD', 'ROE', 'KEPES'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
select case(trim(solver))
|
||
|
||
case ("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case ("hllc", "HLLC")
|
||
name_sol = "HLLC"
|
||
riemann_solver => riemann_mhd_hllc
|
||
|
||
case ("hlld", "HLLD")
|
||
name_sol = "HLLD"
|
||
riemann_solver => riemann_mhd_adi_hlld
|
||
|
||
case("roe", "ROE")
|
||
name_sol = "ROE"
|
||
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(msg,fmt) 'adiabatic magnetohydrodynamics', &
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'HLLC', 'HLLD', 'ROE', 'KEPES'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
end select
|
||
|
||
!--- SPECIAL RELATIVITY HYDRODYNAMICS ---
|
||
!
|
||
case("srhd", "SRHD")
|
||
|
||
select case(trim(eos))
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
select case(trim(statev))
|
||
|
||
case("4vec", "4-vector", "4VEC", "4-VECTOR")
|
||
name_sts = "4-vector"
|
||
states_4vec = .true.
|
||
|
||
case default
|
||
name_sts = "primitive"
|
||
|
||
end select
|
||
|
||
select case(trim(solver))
|
||
|
||
case ("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
|
||
name_sol = "HLLC (Mignone & Bodo 2005)"
|
||
riemann_solver => riemann_srhd_hllc
|
||
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(msg,fmt) 'adiabatic special relativity hydrodynamics', &
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'HLLC'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
end select
|
||
|
||
!--- SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS ---
|
||
!
|
||
case("srmhd", "SRMHD")
|
||
|
||
select case(trim(eos))
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
select case(trim(statev))
|
||
|
||
case("4vec", "4-vector", "4VEC", "4-VECTOR")
|
||
name_sts = "4-vector"
|
||
states_4vec = .true.
|
||
|
||
case default
|
||
name_sts = "primitive"
|
||
|
||
end select
|
||
|
||
select case(trim(solver))
|
||
|
||
case ("hll", "HLL")
|
||
name_sol = "HLL"
|
||
riemann_solver => riemann_hll
|
||
|
||
case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M")
|
||
name_sol = "HLLC (Mignone & Bodo 2006)"
|
||
riemann_solver => riemann_srmhd_hllc
|
||
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(msg,fmt) 'adiabatic special relativity magnetohydrodynamics',&
|
||
"'" // uppercase(trim(solver)) // "'", &
|
||
new_line('a'), &
|
||
"'HLL', 'HLLC'"
|
||
call print_message(loc, msg)
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
end select
|
||
|
||
end select
|
||
|
||
! flag for higher order flux correction
|
||
!
|
||
call get_parameter("high_order_flux", flux)
|
||
|
||
select case(trim(flux))
|
||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||
high_order_flux = .true.
|
||
case default
|
||
high_order_flux = .false.
|
||
end select
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine initialize_schemes
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FINALIZE_SCHEMES:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine releases memory used by the module.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine finalize_schemes(status)
|
||
|
||
implicit none
|
||
|
||
integer, intent(out) :: status
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
nullify(riemann_solver)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine finalize_schemes
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRINT_SCHEMES:
|
||
! ------------------------
|
||
!
|
||
! Subroutine prints module parameters and settings.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine print_schemes(verbose)
|
||
|
||
use helpers, only : print_section, print_parameter
|
||
|
||
implicit none
|
||
|
||
logical, intent(in) :: verbose
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (verbose) then
|
||
call print_section(verbose, "Schemes")
|
||
call print_parameter(verbose, "Riemann solver" , name_sol)
|
||
call print_parameter(verbose, "state variables", name_sts)
|
||
if (high_order_flux) then
|
||
call print_parameter(verbose, "high order flux correction", "on" )
|
||
else
|
||
call print_parameter(verbose, "high order flux correction", "off")
|
||
end if
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine print_schemes
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine UPDATE_FLUX:
|
||
! ----------------------
|
||
!
|
||
! Subroutine solves the Riemann problem along each direction and calculates
|
||
! the numerical fluxes, which are used later to calculate the conserved
|
||
! variable increments.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! dx - the spatial step;
|
||
! q - the array of primitive variables;
|
||
! s - the state arrays for primitive variables;
|
||
! f - the array of numerical fluxes;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine update_flux(dx, q, s, f)
|
||
|
||
use coordinates, only : nn => bcells, nbl, neu
|
||
use equations , only : relativistic
|
||
use equations , only : nf, ivars, ivx, ivz
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:) , intent(in) :: dx
|
||
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:,:), intent(inout) :: s
|
||
real(kind=8), dimension(:,:,:,:,:) , intent(out) :: f
|
||
|
||
integer :: n, i, j, k, l, p
|
||
real(kind=8) :: vm
|
||
|
||
real(kind=8), dimension(nf,nn,2) :: qi
|
||
real(kind=8), dimension(nf,nn) :: fi
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
|
||
! in the relativistic case, apply the reconstruction on variables
|
||
! using the four-velocities if requested
|
||
!
|
||
if (relativistic .and. states_4vec) then
|
||
|
||
f(:,:,:,:,1) = q(:,:,:,:)
|
||
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
|
||
|
||
f(ivx:ivz,i,j,k,1) = f(ivx:ivz,i,j,k,1) / vm
|
||
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
call reconstruct_interfaces(dx(:), f(:,:,:,:,1), s(:,:,:,:,:,:))
|
||
|
||
! convert the states' four-velocities back to velocities
|
||
!
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
do l = 1, 2
|
||
do p = 1, NDIMS
|
||
|
||
vm = sqrt(1.0d+00 + sum(s(ivx:ivz,i,j,k,l,p)**2))
|
||
|
||
s(ivx:ivz,i,j,k,l,p) = s(ivx:ivz,i,j,k,l,p) / vm
|
||
|
||
end do ! p = 1, ndims
|
||
end do ! l = 1, 2
|
||
end do ! i = 1, nn
|
||
end do ! j = 1, nn
|
||
#if NDIMS == 3
|
||
end do ! k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
|
||
else
|
||
|
||
call reconstruct_interfaces(dx(:), q(:,:,:,:), s(:,:,:,:,:,:))
|
||
|
||
end if
|
||
|
||
f(:,:,:,:,:) = 0.0d+00
|
||
|
||
! calculate the flux along the X-direction
|
||
!
|
||
#if NDIMS == 3
|
||
do k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
do j = nbl, neu
|
||
|
||
! copy states to directional lines with proper vector component ordering
|
||
!
|
||
do n = 1, nf
|
||
qi(n,:,1:2) = s(ivars(1,n),:,j,k,1:2,1)
|
||
end do
|
||
|
||
! call one dimensional Riemann solver in order to obtain numerical fluxes
|
||
!
|
||
call numerical_flux(qi(:,:,:), fi(:,:))
|
||
|
||
! update the array of fluxes
|
||
!
|
||
do n = 1, nf
|
||
f(ivars(1,n),:,j,k,1) = fi(n,:)
|
||
end do
|
||
|
||
end do ! j = nbl, neu
|
||
|
||
! calculate the flux along the Y direction
|
||
!
|
||
do i = nbl, neu
|
||
|
||
! copy directional variable vectors to pass to the one dimensional solver
|
||
!
|
||
do n = 1, nf
|
||
qi(n,:,1:2) = s(ivars(2,n),i,:,k,1:2,2)
|
||
end do
|
||
|
||
! call one dimensional Riemann solver in order to obtain numerical fluxes
|
||
!
|
||
call numerical_flux(qi(:,:,:), fi(:,:))
|
||
|
||
! update the array of fluxes
|
||
!
|
||
do n = 1, nf
|
||
f(ivars(2,n),i,:,k,2) = fi(n,:)
|
||
end do
|
||
|
||
end do ! i = nbl, neu
|
||
#if NDIMS == 3
|
||
end do ! k = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
#if NDIMS == 3
|
||
! calculate the flux along the Z direction
|
||
!
|
||
do j = nbl, neu
|
||
do i = nbl, neu
|
||
|
||
! copy directional variable vectors to pass to the one dimensional solver
|
||
!
|
||
do n = 1, nf
|
||
qi(n,:,1:2) = s(ivars(3,n),i,j,:,1:2,3)
|
||
end do
|
||
|
||
! call one dimensional Riemann solver in order to obtain numerical fluxes
|
||
!
|
||
call numerical_flux(qi(:,:,:), fi(:,:))
|
||
|
||
! update the array of fluxes
|
||
!
|
||
do n = 1, nf
|
||
f(ivars(3,n),i,j,:,3) = fi(n,:)
|
||
end do
|
||
|
||
end do ! i = nbl, neu
|
||
end do ! j = nbl, neu
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine update_flux
|
||
!
|
||
!===============================================================================
|
||
!!
|
||
!!*** PRIVATE SUBROUTINES ****************************************************
|
||
!!
|
||
!===============================================================================
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RECONSTRUCT_INTERFACES:
|
||
! ---------------------------------
|
||
!
|
||
! Subroutine reconstructs the Riemann states along all directions.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! dx - the spatial step;
|
||
! q - the array of primitive variables;
|
||
! qi - the array of reconstructed states (2 in each direction);
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine reconstruct_interfaces(dx, q, qi)
|
||
|
||
! include external procedures
|
||
!
|
||
use equations , only : nf
|
||
use equations , only : positive
|
||
use interpolations, only : interfaces
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:) , intent(in) :: dx
|
||
real(kind=8), dimension(:,:,:,:) , intent(in) :: q
|
||
real(kind=8), dimension(:,:,:,:,:,:), intent(out) :: qi
|
||
|
||
! local variables
|
||
!
|
||
integer :: p
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! interpolate interfaces for all flux variables
|
||
!
|
||
do p = 1, nf
|
||
call interfaces(positive(p), dx(1:NDIMS), q (p,:,:,:), &
|
||
qi(p,:,:,:,1:2,1:NDIMS))
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine reconstruct_interfaces
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NUMERICAL_FLUX:
|
||
! -------------------------
|
||
!
|
||
! Subroutine prepares Riemann states and calls the selected Riemann solver
|
||
! in order to get the numerical flux. If requested, the resulting flux
|
||
! is corrected by higher order correction terms.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the array of primitive variables at the Riemann states;
|
||
! f - the output array of fluxes;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine numerical_flux(q, f)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, ibx, ibp
|
||
use equations , only : magnetized, cmax
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:), intent(inout) :: q
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
|
||
integer :: i
|
||
real(kind=8) :: bx, bp
|
||
|
||
real(kind=8), dimension(nf,nn) :: ql, qr
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (kepes) then
|
||
call riemann_solver(q(:,:,1), q(:,:,2), f(:,:))
|
||
else
|
||
|
||
! copy the state vectors
|
||
!
|
||
ql(:,:) = q(:,:,1)
|
||
qr(:,:) = q(:,:,2)
|
||
|
||
! obtain the state values for Bx and Psi for the GLM-MHD equations
|
||
!
|
||
if (magnetized) then
|
||
do i = 1, size(q, 2)
|
||
|
||
bx = 5.0d-01 * ((qr(ibx,i) + ql(ibx,i)) &
|
||
- (qr(ibp,i) - ql(ibp,i)) / cmax)
|
||
bp = 5.0d-01 * ((qr(ibp,i) + ql(ibp,i)) &
|
||
- (qr(ibx,i) - ql(ibx,i)) * cmax)
|
||
|
||
ql(ibx,i) = bx
|
||
qr(ibx,i) = bx
|
||
ql(ibp,i) = bp
|
||
qr(ibp,i) = bp
|
||
|
||
end do
|
||
end if
|
||
|
||
! get the numerical fluxes
|
||
!
|
||
call riemann_solver(ql(:,:), qr(:,:), f(:,:))
|
||
|
||
! higher order flux corrections
|
||
!
|
||
call higher_order_flux_correction(f)
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine numerical_flux
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! GENERIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HLL:
|
||
! ----------------------
|
||
!
|
||
! Subroutine solves one dimensional general Riemann problem using
|
||
! the Harten-Lax-van Leer (HLL) method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Harten, A., Lax, P. D., Van Leer, B.
|
||
! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic
|
||
! Conservation Laws",
|
||
! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61,
|
||
! https://doi.org/10.1137/1025002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hll(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sl, sr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
wl = sl * ul(:,i) - fl(:,i)
|
||
wr = sr * ur(:,i) - fr(:,i)
|
||
|
||
fi(:,i) = (sl * wr - sr * wl) / (sr - sl)
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hll
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ISOTHERMAL HYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HD_ISO_ROE:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Roe, P. L.
|
||
! "Approximate Riemann Solvers, Parameter Vectors, and Difference
|
||
! Schemes",
|
||
! Journal of Computational Physics, 1981, 43, pp. 357-372
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_iso_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, csnd
|
||
use equations , only : idn, ivx, ivy, ivz
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sdl, sdr, sds, cs, ch
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qs, lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: chi
|
||
real(kind=8), dimension(4,4), save :: rvec, lvec
|
||
!$omp threadprivate(first, chi, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
chi = 5.0d-01 / csnd
|
||
|
||
rvec(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||
rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
|
||
lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
lvec(:,2) = [ chi , 0.0d+00, 0.0d+00, - chi ]
|
||
lvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
lvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl)
|
||
call fluxspeed(qr, ur, fr)
|
||
|
||
do i = 1, nn
|
||
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
qs(idn) = sdl * sdr
|
||
qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
|
||
|
||
cs = sign(csnd, qs(ivx))
|
||
ch = sign(chi , qs(ivx))
|
||
|
||
lm(1) = qs(ivx) + cs
|
||
lm(2) = qs(ivx)
|
||
lm(3) = qs(ivx)
|
||
lm(4) = qs(ivx) - cs
|
||
|
||
rvec(1,2) = lm(1)
|
||
rvec(4,2) = lm(4)
|
||
rvec(1,3) = qs(ivy)
|
||
rvec(4,3) = qs(ivy)
|
||
rvec(1,4) = qs(ivz)
|
||
rvec(4,4) = qs(ivz)
|
||
|
||
lvec(1,1) = - ch * lm(4)
|
||
lvec(2,1) = - qs(ivy)
|
||
lvec(3,1) = - qs(ivz)
|
||
lvec(4,1) = ch * lm(1)
|
||
lvec(1,2) = ch
|
||
lvec(4,2) = - ch
|
||
|
||
al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
|
||
df = matmul(al, rvec)
|
||
fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_iso_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HD_ISO_KEPES:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional isothermal hydrodynamic Riemann problem
|
||
! using the entropy stable KEPES method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J.,
|
||
! "Entropy stable numerical approximations for the isothermal and
|
||
! polytropic Euler equations",
|
||
! BIT Numerical Mathematics (2020) 60:791–824,
|
||
! https://doi.org/10.1007/s10543-019-00789-w
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_iso_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, csnd, csnd2
|
||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: v2l, v2r, cs
|
||
real(kind=8) :: dnl, dna, vxa, vya, vza
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8), dimension(4,4), save :: rm
|
||
!$omp threadprivate(first, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
rm(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rm(:,2) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
rm(:,3) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||
rm(:,4) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
|
||
v(idn) = 5.0d-01 * (v2l - v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
v(ivx) = qr(ivx,i) - ql(ivx,i)
|
||
v(ivy) = qr(ivy,i) - ql(ivy,i)
|
||
v(ivz) = qr(ivz,i) - ql(ivz,i)
|
||
|
||
cs = sign(csnd, vxa)
|
||
|
||
lm(:) = [ vxa + cs, vxa, vxa, vxa - cs ]
|
||
|
||
tm(1) = 5.0d-01 * dnl / csnd2
|
||
tm(2) = dnl
|
||
tm(3) = dnl
|
||
tm(4) = tm(1)
|
||
|
||
rm(2:4,1) = [ lm(1), vya, vza ]
|
||
rm(2:4,4) = [ lm(4), vya, vza ]
|
||
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa + csnd2 * dna
|
||
fi(imy,i) = fi(idn,i) * vya
|
||
fi(imz,i) = fi(idn,i) * vza
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_iso_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HD_HLLC:
|
||
! --------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the HLLC
|
||
! method by Gurski or Li. The HLLC and HLLC-C differ by definitions of
|
||
! the tangential components of the velocity and magnetic field.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Toro, E. F.,
|
||
! "The HLLC Riemann solver",
|
||
! Shock Waves, 2019, 29, 8, 1065-1082,
|
||
! https://doi.org/10.1007/s00193-019-00912-4
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_hllc(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf
|
||
use equations , only : idn, ivy, ivz, imx, imy, imz, ien
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: dn, pr
|
||
real(kind=8) :: sl, sr, sm, slmm, srmm
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, ui
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
wl = sl * ul(:,i) - fl(:,i)
|
||
wr = sr * ur(:,i) - fr(:,i)
|
||
|
||
dn = wr(idn) - wl(idn)
|
||
sm = (wr(imx) - wl(imx)) / dn
|
||
|
||
pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
|
||
|
||
if (sm > 0.0d+00) then
|
||
|
||
slmm = sl - sm
|
||
|
||
ui(idn) = wl(idn) / slmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * ql(ivy,i)
|
||
ui(imz) = ui(idn) * ql(ivz,i)
|
||
ui(ien) = (wl(ien) + sm * pr) / slmm
|
||
|
||
fi(:,i) = sl * ui - wl
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
srmm = sr - sm
|
||
|
||
ui(idn) = wr(idn) / srmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * qr(ivy,i)
|
||
ui(imz) = ui(idn) * qr(ivz,i)
|
||
ui(ien) = (wr(ien) + sm * pr) / srmm
|
||
|
||
fi(:,i) = sr * ui - wr
|
||
|
||
else
|
||
|
||
fi(:,i) = (sl * wr - sr * wl) / (sr - sl)
|
||
|
||
end if
|
||
|
||
end if
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_hllc
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HD_ADI_ROE:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Roe, P. L.
|
||
! "Approximate Riemann Solvers, Parameter Vectors, and Difference
|
||
! Schemes",
|
||
! Journal of Computational Physics, 1981, 43, pp. 357-372
|
||
! https://doi.org/10.1016/0021-9991(81)90128-5
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_adi_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index
|
||
use equations , only : idn, ivx, ivy, ivz, ipr, ien
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sdl, sdr, sds, vv, vh, c2, cs, vc, fa, fb, fc, f1, f2, f3
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qs, lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x
|
||
real(kind=8), dimension(5,5), save :: rvec, lvec
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adiabatic_index / adi_m1
|
||
|
||
rvec(:,1) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||
rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
rvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
|
||
lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
lvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
lvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
lvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
lvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl)
|
||
call fluxspeed(qr, ur, fr)
|
||
|
||
do i = 1, nn
|
||
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
qs(idn) = sdl * sdr
|
||
qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
|
||
qs(ien) = ((ul(ien,i) + ql(ipr,i)) / sdl &
|
||
+ (ur(ien,i) + qr(ipr,i)) / sdr) / sds ! enthalpy
|
||
|
||
vv = sum(qs(ivx:ivz) * qs(ivx:ivz))
|
||
vh = 5.0d-01 * vv
|
||
c2 = adi_m1 * (qs(ien) - vh)
|
||
cs = sign(sqrt(c2), qs(ivx))
|
||
fa = adi_m1 / c2
|
||
fb = 5.0d-01 * fa
|
||
fc = 5.0d-01 / cs
|
||
vc = cs * qs(ivx)
|
||
f1 = fb * vh
|
||
f2 = fc * qs(ivx)
|
||
f3 = - fb * qs(ivx)
|
||
|
||
lm(1) = qs(ivx) + cs
|
||
lm(2) = qs(ivx)
|
||
lm(3) = qs(ivx)
|
||
lm(4) = qs(ivx)
|
||
lm(5) = qs(ivx) - cs
|
||
|
||
rvec(1,2) = lm(1)
|
||
rvec(1,3) = qs(ivy)
|
||
rvec(1,4) = qs(ivz)
|
||
rvec(1,5) = qs(ien) + vc
|
||
rvec(2,2) = qs(ivx)
|
||
rvec(2,3) = qs(ivy)
|
||
rvec(2,4) = qs(ivz)
|
||
rvec(2,5) = vh
|
||
rvec(3,5) = qs(ivy)
|
||
rvec(4,5) = qs(ivz)
|
||
rvec(5,2) = lm(5)
|
||
rvec(5,3) = qs(ivy)
|
||
rvec(5,4) = qs(ivz)
|
||
rvec(5,5) = qs(ien) - vc
|
||
|
||
lvec(1,1) = f1 - f2
|
||
lvec(1,2) = f3 + fc
|
||
lvec(1,3) = - fb * qs(ivy)
|
||
lvec(1,4) = - fb * qs(ivz)
|
||
lvec(1,5) = fb
|
||
lvec(2,1) = - fa * vh + 1.0d+00
|
||
lvec(2,2) = fa * qs(ivx)
|
||
lvec(2,3) = fa * qs(ivy)
|
||
lvec(2,4) = fa * qs(ivz)
|
||
lvec(2,5) = - fa
|
||
lvec(3,1) = - qs(ivy)
|
||
lvec(4,1) = - qs(ivz)
|
||
lvec(5,1) = f1 + f2
|
||
lvec(5,2) = f3 - fc
|
||
lvec(5,3) = lvec(1,3)
|
||
lvec(5,4) = lvec(1,4)
|
||
lvec(5,5) = fb
|
||
|
||
al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i))
|
||
df = matmul(al, rvec)
|
||
fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df)
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_adi_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_HD_ADI_KEPES:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional adiabatic hydrodynamic Riemann problem
|
||
! using the entropy stable KEPES method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Winters, A. R., Derigs, D., Gassner, G. J., Walch, S.
|
||
! "A uniquely defined entropy stable matrix dissipation operator
|
||
! for high Mach number ideal MHD and compressible Euler
|
||
! simulations",
|
||
! Journal of Computational Physics, 2017, 332, pp. 274-289
|
||
! https://doi.org/10.1016/j.jcp.2016.12.006
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_hd_adi_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index
|
||
use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: bl, br, v2l, v2r, cs, vc
|
||
real(kind=8) :: btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra, ent
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x
|
||
real(kind=8), dimension(5,5), save :: rm
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adi_m1 / adiabatic_index
|
||
|
||
rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||
rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||
rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||
rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
bl = ql(idn,i) / ql(ipr,i)
|
||
br = qr(idn,i) / qr(ipr,i)
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
|
||
btl = lmean(bl, br)
|
||
bta = amean(bl, br)
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
pra = dna / bta
|
||
prl = dnl / btl
|
||
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
ent = 1.0d+00 / adi_m1x / btl + vva
|
||
|
||
v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
if (br > bl) then
|
||
v(idn) = v(idn) + log(br / bl) / adi_m1
|
||
else if (bl > br) then
|
||
v(idn) = v(idn) - log(bl / br) / adi_m1
|
||
end if
|
||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||
v(ipr) = bl - br
|
||
|
||
cs = sign(sqrt(adiabatic_index * pra / dnl), vxa)
|
||
vc = vxa * cs
|
||
|
||
lm(1) = vxa + cs
|
||
lm(2) = vxa
|
||
lm(3) = vxa
|
||
lm(4) = vxa
|
||
lm(5) = vxa - cs
|
||
|
||
tm(1) = 5.0d-01 * dnl / adiabatic_index
|
||
tm(2) = dnl * adi_m1x
|
||
tm(3) = pra
|
||
tm(4) = pra
|
||
tm(5) = tm(1)
|
||
|
||
rm(2,1) = lm(1)
|
||
rm(3,1) = vya
|
||
rm(4,1) = vza
|
||
rm(2,2) = vxa
|
||
rm(3,2) = vya
|
||
rm(4,2) = vza
|
||
rm(2,5) = lm(5)
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,:) = [ ent + vc, vva, vya, vza, ent - vc ]
|
||
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa + pra
|
||
fi(imy,i) = fi(idn,i) * vya
|
||
fi(imz,i) = fi(idn,i) * vza
|
||
fi(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_hd_adi_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! RELATIVISTIC ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_SRHD_HLLC:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
|
||
! by Mignone & Bodo.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A. & Bodo, G.
|
||
! "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society,
|
||
! 2005, Volume 364, Pages 126-136
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_srhd_hllc(ql, qr, fi)
|
||
|
||
use algebra , only : quadratic
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ien
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, nr
|
||
real(kind=8) :: sl, sr, srml, sm
|
||
real(kind=8) :: pr, dv
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
|
||
real(kind=8), dimension(3) :: a
|
||
real(kind=8), dimension(2) :: x
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
uh(:) = ( wr(:) - wl(:)) / srml
|
||
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
wl(ien) = wl(ien) + wl(idn)
|
||
wr(ien) = wr(ien) + wr(idn)
|
||
|
||
! prepare the quadratic coefficients (eq. 18 in [1])
|
||
!
|
||
a(1) = uh(imx)
|
||
a(2) = - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn)
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! calculate total pressure (eq. 17 in [1])
|
||
!
|
||
pr = fh(imx) - (fh(ien) + fh(idn)) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pr <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 16 in [1])
|
||
!
|
||
dv = sl - sm
|
||
us(idn) = wl(idn) / dv
|
||
us(imy) = wl(imy) / dv
|
||
us(imz) = wl(imz) / dv
|
||
us(ien) = (wl(ien) + pr * sm) / dv
|
||
us(imx) = (us(ien) + pr) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 16 in [1])
|
||
!
|
||
dv = sr - sm
|
||
us(idn) = wr(idn) / dv
|
||
us(imy) = wr(imy) / dv
|
||
us(imz) = wr(imz) / dv
|
||
us(ien) = (wr(ien) + pr * sm) / dv
|
||
us(imx) = (us(ien) + pr) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity and all fluxes
|
||
! except the parallel momentum one are zero
|
||
!
|
||
fi(idn,i) = 0.0d+00
|
||
fi(imx,i) = pr
|
||
fi(imy,i) = 0.0d+00
|
||
fi(imz,i) = 0.0d+00
|
||
fi(ien,i) = 0.0d+00
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_srhd_hllc
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ISOTHERMAL MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_ISO_MHD_HLLD:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the modified
|
||
! isothermal HLLD method by Mignone. The difference is in allowing a density
|
||
! jumps across the Alfvén waves, due to the inner intermediate states being
|
||
! an average of the states between two Alfvén waves, which also include
|
||
! two slow magnetosonic waves. Moreover, due to the different left and right
|
||
! intermediate state densities the left and right Alfvén speeds might be
|
||
! different too.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A.,
|
||
! "A simple and accurate Riemann solver for isothermal MHD",
|
||
! Journal of Computational Physics, 2007, 225, pp. 1427-1441,
|
||
! https://doi.org/10.1016/j.jcp.2007.01.033
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_hlld(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
|
||
real(kind=8) :: bx, b2, dvl, dvr, dnl, dnr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, uil, uir
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn))
|
||
|
||
srml = sr - sl
|
||
slmm = sl - sm
|
||
srmm = sr - sm
|
||
|
||
uil(idn) = wl(idn) / slmm
|
||
uir(idn) = wr(idn) / srmm
|
||
uil(imx) = uil(idn) * sm
|
||
uir(imx) = uir(idn) * sm
|
||
uil(ibx) = bx
|
||
uir(ibx) = bx
|
||
uil(ibp) = ql(ibp,i)
|
||
uir(ibp) = qr(ibp,i)
|
||
|
||
sml = sm - abs(bx) / sqrt(uil(idn))
|
||
smr = sm + abs(bx) / sqrt(uir(idn))
|
||
|
||
dvl = slmm * wl(idn) - b2
|
||
dvr = srmm * wr(idn) - b2
|
||
|
||
dnr = uir(idn) / dvr
|
||
dnl = uil(idn) / dvl
|
||
|
||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||
!
|
||
if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then
|
||
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
fi(:,i) = sl * uil(:) - wl(:)
|
||
|
||
else if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
fi(:,i) = sr * uir(:) - wr(:)
|
||
|
||
else ! sl* <= 0 <= sr*
|
||
|
||
if (smr > sml) then
|
||
|
||
uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||
|
||
dvl = smr - sml
|
||
fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dvl
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sl* <= 0 <= sr*
|
||
|
||
else if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr
|
||
|
||
uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby))
|
||
uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz))
|
||
uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
fi(:,i) = sl * uil(:) - wl(:)
|
||
|
||
else if (sr > sml) then
|
||
|
||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||
|
||
dvl = sr - sml
|
||
fi(:,i) = (sml * wr(:) - sr * wcl(:)) / dvl
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sl* < 0
|
||
|
||
else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl
|
||
|
||
uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby))
|
||
uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz))
|
||
uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
|
||
if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
fi(:,i) = sr * uir(:) - wr(:)
|
||
|
||
else if (smr > sl) then
|
||
|
||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||
|
||
dvr = smr - sl
|
||
fi(:,i) = (sl * wcr(:) - smr * wl(:)) / dvr
|
||
|
||
else ! Bx = 0 -> Sₘ = 0
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sr* > 0
|
||
|
||
else ! sl* < sl & sr* > sr
|
||
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_hlld
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ISO_ROE:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Stone, J. M. & Gardiner, T. A.,
|
||
! "ATHENA: A New Code for Astrophysical MHD",
|
||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, csnd2
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2
|
||
real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als
|
||
real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff
|
||
real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb
|
||
real(kind=8) :: f1, f2, f3, f4
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qi
|
||
real(kind=8), dimension(6 ) :: lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
integer , dimension(6) , save :: ivs
|
||
real(kind=8), dimension(6,6), save :: rvec, lvec
|
||
!$omp threadprivate(first, ivs, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
|
||
rvec(:,:) = 0.0d+00
|
||
lvec(:,:) = 0.0d+00
|
||
|
||
ivs = [ idn, imx, imy, imz, iby, ibz ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl)
|
||
call fluxspeed(qr, ur, fr)
|
||
|
||
do i = 1, nn
|
||
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
|
||
qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
f1 = ql(iby,i) - qr(iby,i)
|
||
f2 = ql(ibz,i) - qr(ibz,i)
|
||
xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
|
||
yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
|
||
br = sqrt(br2)
|
||
if (br2 > 0.0d+00) then
|
||
bty = qi(iby) / br
|
||
btz = qi(ibz) / br
|
||
else
|
||
bty = 1.0d+00
|
||
btz = 0.0d+00
|
||
end if
|
||
|
||
cc2 = csnd2 + xx
|
||
ca2 = qi(ibx) * qi(ibx) / qi(idn)
|
||
f1 = br2 / qi(idn)
|
||
f2 = ca2 + f1
|
||
f3 = f2 - cc2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * cc2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + cc2 + f4)
|
||
cs2 = cc2 * ca2 / cf2
|
||
|
||
if (cs2 < cc2 .and. cc2 < cf2) then
|
||
f1 = (cc2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (cc2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cc = sqrt(cc2)
|
||
ca = sqrt(ca2)
|
||
cf = sign(sqrt(cf2), qi(ivx))
|
||
cs = sign(sqrt(cs2), qi(ivx))
|
||
|
||
lm(1) = qi(ivx) + cf
|
||
lm(2) = qi(ivx) + ca
|
||
lm(3) = qi(ivx) + cs
|
||
lm(4) = qi(ivx) - cs
|
||
lm(5) = qi(ivx) - ca
|
||
lm(6) = qi(ivx) - cf
|
||
|
||
sqrtd = sqrt(qi(idn))
|
||
f1 = sqrt(yy)
|
||
qf = cf * sign(alf, qi(ibx)) / f1
|
||
qs = cs * sign(als, qi(ibx)) / f1
|
||
f1 = cc / (f1 * sqrtd)
|
||
af_prime = alf * f1
|
||
as_prime = als * f1
|
||
|
||
rvec(1,1) = alf
|
||
rvec(1,2) = alf * lm(1)
|
||
rvec(1,3) = alf * qi(ivy) - qs * bty
|
||
rvec(1,4) = alf * qi(ivz) - qs * btz
|
||
rvec(1,5) = as_prime * bty
|
||
rvec(1,6) = as_prime * btz
|
||
|
||
rvec(6,1) = alf
|
||
rvec(6,2) = alf * lm(6)
|
||
rvec(6,3) = alf * qi(ivy) + qs * bty
|
||
rvec(6,4) = alf * qi(ivz) + qs * btz
|
||
rvec(6,5) = rvec(1,5)
|
||
rvec(6,6) = rvec(1,6)
|
||
|
||
rvec(3,1) = als
|
||
rvec(3,2) = als * lm(3)
|
||
rvec(3,3) = als * qi(ivy) + qf * bty
|
||
rvec(3,4) = als * qi(ivz) + qf * btz
|
||
rvec(3,5) = - af_prime * bty
|
||
rvec(3,6) = - af_prime * btz
|
||
|
||
rvec(4,1) = als
|
||
rvec(4,2) = als * lm(4)
|
||
rvec(4,3) = als * qi(ivy) - qf * bty
|
||
rvec(4,4) = als * qi(ivz) - qf * btz
|
||
rvec(4,5) = rvec(3,5)
|
||
rvec(4,6) = rvec(3,6)
|
||
|
||
f1 = sign(1.0d+00 / sqrtd, qi(ibx))
|
||
rvec(2,3) = btz
|
||
rvec(2,4) = - bty
|
||
rvec(2,5) = - btz * f1
|
||
rvec(2,6) = bty * f1
|
||
|
||
rvec(5,3) = - rvec(2,3)
|
||
rvec(5,4) = - rvec(2,4)
|
||
rvec(5,5) = rvec(2,5)
|
||
rvec(5,6) = rvec(2,6)
|
||
|
||
norm = 2.0d+00 * cc2
|
||
cff = alf * cf / norm
|
||
css = als * cs / norm
|
||
norm = norm / yy
|
||
qf = qf / norm
|
||
qs = qs / norm
|
||
f1 = qi(idn) / norm
|
||
af = af_prime * f1
|
||
as = as_prime * f1
|
||
f1 = br / norm
|
||
afpb = af_prime * f1
|
||
aspb = as_prime * f1
|
||
vqstr = (qi(ivy) * bty + qi(ivz) * btz)
|
||
|
||
f1 = qs * vqstr
|
||
lvec(1,1) = - cff * lm(6) + f1 - aspb
|
||
lvec(1,2) = cff
|
||
lvec(1,3) = - qs * bty
|
||
lvec(1,4) = - qs * btz
|
||
lvec(1,5) = as * bty
|
||
lvec(1,6) = as * btz
|
||
|
||
lvec(6,1) = cff * lm(1) - f1 - aspb
|
||
lvec(6,2) = - lvec(1,2)
|
||
lvec(6,3) = - lvec(1,3)
|
||
lvec(6,4) = - lvec(1,4)
|
||
lvec(6,5) = lvec(1,5)
|
||
lvec(6,6) = lvec(1,6)
|
||
|
||
f1 = qf * vqstr
|
||
lvec(3,1) = - css * lm(4) - f1 + afpb
|
||
lvec(3,2) = css
|
||
lvec(3,3) = qf * bty
|
||
lvec(3,4) = qf * btz
|
||
lvec(3,5) = - af * bty
|
||
lvec(3,6) = - af * btz
|
||
|
||
lvec(4,1) = css * lm(3) + f1 + afpb
|
||
lvec(4,2) = - lvec(3,2)
|
||
lvec(4,3) = - lvec(3,3)
|
||
lvec(4,4) = - lvec(3,4)
|
||
lvec(4,5) = lvec(3,5)
|
||
lvec(4,6) = lvec(3,6)
|
||
|
||
f1 = sign(5.0d-01 * sqrtd, qi(ibx))
|
||
lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty)
|
||
lvec(2,3) = 5.0d-01 * btz
|
||
lvec(2,4) = - 5.0d-01 * bty
|
||
lvec(2,5) = - f1 * btz
|
||
lvec(2,6) = f1 * bty
|
||
|
||
lvec(5,1) = - lvec(2,1)
|
||
lvec(5,3) = - lvec(2,3)
|
||
lvec(5,4) = - lvec(2,4)
|
||
lvec(5,5) = lvec(2,5)
|
||
lvec(5,6) = lvec(2,6)
|
||
|
||
al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
|
||
df = matmul(al, rvec)
|
||
fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ISO_KEPES:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine solves one dimensional isothermal magnetohydrodynamic
|
||
! Riemann problem using the entropy stable KEPES method. The method is
|
||
! a modification of the method described in [1] for the isothermal equation
|
||
! of state.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
|
||
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
|
||
! field divergence diminishing ideal magnetohydrodynamics equations",
|
||
! Journal of Computational Physics, 2018, 364, pp. 420-467,
|
||
! https://doi.org/10.1016/j.jcp.2018.03.002
|
||
! [2] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J.,
|
||
! "Entropy stable numerical approximations for the isothermal and
|
||
! polytropic Euler equations",
|
||
! BIT Numerical Mathematics (2020) 60:791–824,
|
||
! https://doi.org/10.1007/s10543-019-00789-w
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_iso_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, csnd, csnd2, ch => cglm
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
ibx, iby, ibz, ibp
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
|
||
real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq
|
||
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
|
||
real(kind=8) :: alf, als, f1, f2, f3, f4
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8), dimension(8,8), save :: rm
|
||
!$omp threadprivate(first, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
rm(:,:) = 0.0d+00
|
||
rm(5,4) = 1.0d+00
|
||
rm(8,4) = 1.0d+00
|
||
rm(5,5) = 1.0d+00
|
||
rm(8,5) = - 1.0d+00
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
|
||
b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
|
||
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
bxa = amean(ql(ibx,i), qr(ibx,i))
|
||
bya = amean(ql(iby,i), qr(iby,i))
|
||
bza = amean(ql(ibz,i), qr(ibz,i))
|
||
bpa = amean(ql(ibp,i), qr(ibp,i))
|
||
pta = 5.0d-01 * amean(b2l, b2r) + csnd2 * dna
|
||
|
||
v(idn) = 5.0d-01 * (v2l - v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
v(ivx) = qr(ivx,i) - ql(ivx,i)
|
||
v(ivy) = qr(ivy,i) - ql(ivy,i)
|
||
v(ivz) = qr(ivz,i) - ql(ivz,i)
|
||
v(ibx) = qr(ibx,i) - ql(ibx,i)
|
||
v(iby) = qr(iby,i) - ql(iby,i)
|
||
v(ibz) = qr(ibz,i) - ql(ibz,i)
|
||
v(ibp) = qr(ibp,i) - ql(ibp,i)
|
||
|
||
bp2 = bya * bya + bza * bza
|
||
if (bp2 > 0.0d+00) then
|
||
f1 = sqrt(bp2)
|
||
x2 = bya / f1
|
||
x3 = bza / f1
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
ca2 = bxa * bxa / dnl
|
||
f1 = bp2 / dnl
|
||
f2 = ca2 + f1
|
||
f3 = f2 - csnd2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * csnd2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + csnd2 + f4)
|
||
cs2 = csnd2 * ca2 / cf2
|
||
|
||
if (cs2 < csnd2 .and. csnd2 < cf2) then
|
||
f1 = (csnd2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (csnd2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), vxa)
|
||
ca = sqrt(ca2)
|
||
cs = sign(sqrt(cs2), vxa)
|
||
|
||
lm(1) = vxa + cf
|
||
lm(2) = vxa + ca
|
||
lm(3) = vxa + cs
|
||
lm(4) = vxa + ch
|
||
lm(5) = vxa - ch
|
||
lm(6) = vxa - cs
|
||
lm(7) = vxa - ca
|
||
lm(8) = vxa - cf
|
||
|
||
tm(1) = 5.0d-01 / (dnl * csnd2)
|
||
tm(2) = 5.0d-01 / dnl**2
|
||
tm(3) = tm(1)
|
||
tm(4) = 5.0d-01
|
||
tm(5) = tm(4)
|
||
tm(6) = tm(1)
|
||
tm(7) = tm(2)
|
||
tm(8) = tm(1)
|
||
|
||
das = dnl * als
|
||
daf = dnl * alf
|
||
csq = csnd * sqrt(dnl)
|
||
|
||
f1 = daf
|
||
f2 = sign(das, bxa) * cs
|
||
f3 = als * csq
|
||
rm(1,1) = f1
|
||
rm(2,1) = f1 * (vxa + cf)
|
||
rm(3,1) = f1 * vya - f2 * x2
|
||
rm(4,1) = f1 * vza - f2 * x3
|
||
rm(6,1) = f3 * x2
|
||
rm(7,1) = f3 * x3
|
||
|
||
rm(1,8) = f1
|
||
rm(2,8) = f1 * (vxa - cf)
|
||
rm(3,8) = f1 * vya + f2 * x2
|
||
rm(4,8) = f1 * vza + f2 * x3
|
||
rm(6,8) = rm(6,1)
|
||
rm(7,8) = rm(7,1)
|
||
|
||
f1 = dnl * sqrt(dna)
|
||
rm(3,2) = f1 * x3
|
||
rm(4,2) = - f1 * x2
|
||
rm(6,2) = - dnl * x3
|
||
rm(7,2) = dnl * x2
|
||
|
||
rm(3,7) = - rm(3,2)
|
||
rm(4,7) = - rm(4,2)
|
||
rm(6,7) = rm(6,2)
|
||
rm(7,7) = rm(7,2)
|
||
|
||
f1 = das
|
||
f2 = sign(daf, bxa) * cf
|
||
f3 = - alf * csq
|
||
rm(1,3) = f1
|
||
rm(2,3) = f1 * (vxa + cs)
|
||
rm(3,3) = f1 * vya + f2 * x2
|
||
rm(4,3) = f1 * vza + f2 * x3
|
||
rm(6,3) = f3 * x2
|
||
rm(7,3) = f3 * x3
|
||
|
||
rm(1,6) = f1
|
||
rm(2,6) = f1 * (vxa - cs)
|
||
rm(3,6) = f1 * vya - f2 * x2
|
||
rm(4,6) = f1 * vza - f2 * x3
|
||
rm(6,6) = rm(6,3)
|
||
rm(7,6) = rm(7,3)
|
||
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
|
||
fi(imy,i) = fi(idn,i) * vya - bxa * bya
|
||
fi(imz,i) = fi(idn,i) * vza - bxa * bza
|
||
fi(ibx,i) = ch * bpa
|
||
fi(iby,i) = vxa * bya - bxa * vya
|
||
fi(ibz,i) = vxa * bza - bxa * vza
|
||
fi(ibp,i) = ch * bxa
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_iso_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_HLLC:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the HLLC method,
|
||
! by Toro. In the HLLC method the tangential components of the velocity are
|
||
! discontinuous actoss the contact dicontinuity.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Toro, E. F., Spruce, M., & Speares, W.
|
||
! "Restoration of the contact surface in the HLL-Riemann solver",
|
||
! Shock Waves, 1994, Volume 4, Issue 1, pp. 25-34
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_hllc(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivy, ivz, imx, imy, imz, ien, &
|
||
ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2
|
||
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, ui
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
dn = wr(idn) - wl(idn)
|
||
sm = (wr(imx) - wl(imx)) / dn
|
||
|
||
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
|
||
|
||
if (b2 > 0.0d+00) then
|
||
|
||
vy = (wr(imy) - wl(imy)) / dn
|
||
vz = (wr(imz) - wl(imz)) / dn
|
||
by = (wr(iby) - wl(iby)) / srml
|
||
bz = (wr(ibz) - wl(ibz)) / srml
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
slmm = sl - sm
|
||
|
||
ui(idn) = wl(idn) / slmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * vy
|
||
ui(imz) = ui(idn) * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ql(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
srmm = sr - sm
|
||
|
||
ui(idn) = wr(idn) / srmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * vy
|
||
ui(imz) = ui(idn) * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = qr(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
|
||
! to the HLL one
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sm = 0
|
||
|
||
else ! Bₓ = 0
|
||
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
slmm = sl - sm
|
||
|
||
ui(idn) = wl(idn) / slmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * ql(ivy,i)
|
||
ui(imz) = ui(idn) * ql(ivz,i)
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = wl(iby) / slmm
|
||
ui(ibz) = wl(ibz) / slmm
|
||
ui(ibp) = ql(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt) / slmm
|
||
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
srmm = sr - sm
|
||
|
||
ui(idn) = wr(idn) / srmm
|
||
ui(imx) = ui(idn) * sm
|
||
ui(imy) = ui(idn) * qr(ivy,i)
|
||
ui(imz) = ui(idn) * qr(ivz,i)
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = wr(iby) / srmm
|
||
ui(ibz) = wr(ibz) / srmm
|
||
ui(ibp) = qr(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt) / srmm
|
||
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
|
||
! to the HLL one
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sm = 0
|
||
|
||
end if ! Bₓ = 0
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_hllc
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_ADI_MHD_HLLD:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using the adiabatic HLLD
|
||
! method described by Miyoshi & Kusano.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Miyoshi, T. & Kusano, K.,
|
||
! "A multi-state HLL approximate Riemann solver for ideal
|
||
! magnetohydrodynamics",
|
||
! Journal of Computational Physics, 2005, 208, pp. 315-344
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_hlld(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sl, sr, sm, srml, slmm, srmm
|
||
real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb
|
||
real(kind=8) :: dnl, dnr, cal, car, sml, smr
|
||
real(kind=8) :: dv, dvl, dvr
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, ui
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then ! sl ≥ 0
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then ! sr ≤ 0
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
bx = ql(ibx,i)
|
||
b2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
dn = wr(idn) - wl(idn)
|
||
sm = (wr(imx) - wl(imx)) / dn
|
||
|
||
slmm = sl - sm
|
||
srmm = sr - sm
|
||
|
||
dnl = wl(idn) / slmm
|
||
dnr = wr(idn) / srmm
|
||
|
||
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
|
||
|
||
cal = abs(bx) / sqrt(dnl)
|
||
car = abs(bx) / sqrt(dnr)
|
||
sml = sm - cal
|
||
smr = sm + car
|
||
|
||
! calculate division factors (also used to determine degeneracies)
|
||
!
|
||
dvl = slmm * wl(idn) - b2
|
||
dvr = srmm * wr(idn) - b2
|
||
|
||
! check degeneracy Sl* -> Sl or Sr* -> Sr
|
||
!
|
||
if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracy
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sl* < 0 < sr*
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||
|
||
! prepare constant primitive variables of the intermediate states
|
||
!
|
||
dv = abs(bx) * (sqrt(dnr) + sqrt(dnl))
|
||
vy = (wcr(imy) - wcl(imy)) / dv
|
||
vz = (wcr(imz) - wcl(imz)) / dv
|
||
dv = car + cal
|
||
by = (wcr(iby) - wcl(iby)) / dv
|
||
bz = (wcr(ibz) - wcl(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm > 0.0d+00) then ! sm > 0
|
||
|
||
! conservative variables for the inmost left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
|
||
|
||
! the inmost left intermediate flux
|
||
!
|
||
fi(:,i) = sml * ui(:) - wcl(:)
|
||
|
||
else if (sm < 0.0d+00) then ! sm < 0
|
||
|
||
! conservative variables for the inmost right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
|
||
|
||
! the inmost right intermediate flux
|
||
!
|
||
fi(:,i) = smr * ui(:) - wcr(:)
|
||
|
||
else ! sm = 0
|
||
|
||
! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
|
||
! the flux can be averaged from the Alfvén waves using a simple HLL formula;
|
||
!
|
||
fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dv
|
||
|
||
end if ! sm = 0
|
||
|
||
else ! Bx = 0
|
||
|
||
if (sm > 0.0d+00) then
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = slmm * wl(imy) / dvl
|
||
vz = slmm * wl(imz) / dvl
|
||
by = wl(idn) * wl(iby) / dvl
|
||
bz = wl(idn) * wl(ibz) / dvl
|
||
vb = vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt) / slmm
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy)) / dvr
|
||
vz = ( srmm * wr(imz)) / dvr
|
||
by = (wr(idn) * wr(iby)) / dvr
|
||
bz = (wr(idn) * wr(ibz)) / dvr
|
||
vb = vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = 0.0d+00
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt) / srmm
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! Sm = 0
|
||
|
||
! since Bx = 0 and Sm = 0, then revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! Bx = 0
|
||
|
||
end if ! sl* < 0 < sr*
|
||
|
||
else ! some degeneracies
|
||
|
||
! separate degeneracies
|
||
!
|
||
if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr
|
||
|
||
! primitive variables for the outer left intermediate state
|
||
!
|
||
vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl
|
||
vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl
|
||
by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl
|
||
bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (sml > 0.0d+00) then ! sl* > 0
|
||
|
||
! the outer left intermediate flux
|
||
!
|
||
fi(:,i) = sl * ui(:) - wl(:)
|
||
|
||
else ! sl* <= 0
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcl(:) = (sml - sl) * ui(:) + wl(:)
|
||
|
||
! primitive variables for the inner left intermediate state
|
||
!
|
||
dv = srmm * dnr + cal * dnl
|
||
vy = (wr(imy) - wcl(imy)) / dv
|
||
vz = (wr(imz) - wcl(imz)) / dv
|
||
dv = sr - sml
|
||
by = (wr(iby) - wcl(iby)) / dv
|
||
bz = (wr(ibz) - wcl(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the inner left intermediate state
|
||
!
|
||
ui(idn) = dnl
|
||
ui(imx) = dnl * sm
|
||
ui(imy) = dnl * vy
|
||
ui(imz) = dnl * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ul(ibp,i)
|
||
ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm >= 0.0d+00) then ! sm >= 0
|
||
|
||
! the inner left intermediate flux
|
||
!
|
||
fi(:,i) = sml * ui(:) - wcl(:)
|
||
|
||
else ! sm < 0
|
||
|
||
! vector of the left-going Alfvén wave
|
||
!
|
||
wcr(:) = (sm - sml) * ui(:) + wcl(:)
|
||
|
||
! calculate the average flux over the right inner intermediate state
|
||
!
|
||
fi(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm
|
||
|
||
end if ! sm < 0
|
||
|
||
else ! bx = 0
|
||
|
||
! no Alfvén wave, so revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sl* < 0
|
||
|
||
else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl
|
||
|
||
! primitive variables for the outer right intermediate state
|
||
!
|
||
vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr
|
||
vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr
|
||
by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr
|
||
bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the outer right intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm
|
||
|
||
! choose the correct state depending on the speed signs
|
||
!
|
||
if (smr < 0.0d+00) then ! sr* < 0
|
||
|
||
! the outer right intermediate flux
|
||
!
|
||
fi(:,i) = sr * ui(:) - wr(:)
|
||
|
||
else ! sr* > 0
|
||
|
||
! separate cases with non-zero and zero Bx
|
||
!
|
||
if (b2 > 0.0d+00) then
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcr(:) = (smr - sr) * ui(:) + wr(:)
|
||
|
||
! primitive variables for the inner right intermediate state
|
||
!
|
||
dv = slmm * dnl - car * dnr
|
||
vy = (wl(imy) - wcr(imy)) / dv
|
||
vz = (wl(imz) - wcr(imz)) / dv
|
||
dv = sl - smr
|
||
by = (wl(iby) - wcr(iby)) / dv
|
||
bz = (wl(ibz) - wcr(ibz)) / dv
|
||
vb = sm * bx + vy * by + vz * bz
|
||
|
||
! conservative variables for the inner left intermediate state
|
||
!
|
||
ui(idn) = dnr
|
||
ui(imx) = dnr * sm
|
||
ui(imy) = dnr * vy
|
||
ui(imz) = dnr * vz
|
||
ui(ibx) = bx
|
||
ui(iby) = by
|
||
ui(ibz) = bz
|
||
ui(ibp) = ur(ibp,i)
|
||
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
|
||
|
||
! choose the correct state depending on the sign of contact discontinuity
|
||
! advection speed
|
||
!
|
||
if (sm <= 0.0d+00) then ! sm <= 0
|
||
|
||
! the inner right intermediate flux
|
||
!
|
||
fi(:,i) = smr * ui(:) - wcr(:)
|
||
|
||
else ! sm > 0
|
||
|
||
! vector of the right-going Alfvén wave
|
||
!
|
||
wcl(:) = (sm - smr) * ui(:) + wcr(:)
|
||
|
||
! calculate the average flux over the left inner intermediate state
|
||
!
|
||
fi(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm
|
||
|
||
end if ! sm > 0
|
||
|
||
else ! bx = 0
|
||
|
||
! no Alfvén wave, so revert to the HLL flux
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if
|
||
|
||
end if ! sr* > 0
|
||
|
||
else ! sl* < sl & sr* > sr
|
||
|
||
! so far we revert to HLL flux in the case of degeneracies
|
||
!
|
||
fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
end if ! sl* < sl & sr* > sr
|
||
|
||
end if ! deneneracies
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_hlld
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ADI_ROE:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Roe's method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Stone, J. M. & Gardiner, T. A.,
|
||
! "ATHENA: A New Code for Astrophysical MHD",
|
||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||
! [2] Toro, E. F.,
|
||
! "Riemann Solvers and Numerical Methods for Fluid dynamics",
|
||
! Springer-Verlag Berlin Heidelberg, 2009
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_roe(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index
|
||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien, &
|
||
ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: sdl, sdr, sds
|
||
real(kind=8) :: xx, yy, pml, pmr
|
||
|
||
real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs
|
||
real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty
|
||
real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm
|
||
real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb
|
||
real(kind=8) :: f1, f2, f3, f4, f5
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension(nf ) :: qi
|
||
real(kind=8), dimension(7 ) :: lm, al, df
|
||
|
||
logical , save :: first = .true.
|
||
integer , dimension(7) , save :: ivs
|
||
real(kind=8) , save :: adi_m1, adi_m2, adi_m2d1
|
||
real(kind=8), dimension(7,7), save :: rvec, lvec
|
||
!$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m2 = adiabatic_index - 2.0d+00
|
||
adi_m2d1 = adi_m2 / adi_m1
|
||
|
||
rvec(:,:) = 0.0d+00
|
||
lvec(:,:) = 0.0d+00
|
||
rvec(4,1) = 1.0d+00
|
||
|
||
ivs = [ idn, imx, imy, imz, ipr, iby, ibz ]
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl)
|
||
call fluxspeed(qr, ur, fr)
|
||
|
||
do i = 1, nn
|
||
|
||
pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
|
||
pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
|
||
sdl = sqrt(ql(idn,i))
|
||
sdr = sqrt(qr(idn,i))
|
||
sds = sdl + sdr
|
||
|
||
qi(idn) = sdl * sdr
|
||
qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds
|
||
qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl &
|
||
+ (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds
|
||
qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds
|
||
qi(ibx) = ql(ibx,i)
|
||
qi(ibp) = ql(ibp,i)
|
||
|
||
f1 = qr(iby,i) - ql(iby,i)
|
||
f2 = qr(ibz,i) - ql(ibz,i)
|
||
xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds)
|
||
yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn)
|
||
|
||
br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz)
|
||
if (br2 > 0.0d+00) then
|
||
br = sqrt(br2)
|
||
bty = qi(iby) / br
|
||
btz = qi(ibz) / br
|
||
else
|
||
br = 0.0d+00
|
||
bty = 1.0d+00
|
||
btz = 0.0d+00
|
||
end if
|
||
|
||
v2 = sum(qi(ivx:ivz) * qi(ivx:ivz))
|
||
v2h = 5.0d-01 * v2
|
||
ayy = adi_m1 - adi_m2 * yy
|
||
sqty = sqrt(ayy)
|
||
ca2 = qi(ibx) * qi(ibx)
|
||
hp = qi(ien) - (ca2 + br2) / qi(idn)
|
||
ca2 = ca2 / qi(idn)
|
||
f1 = adi_m1 * (hp - v2h)
|
||
f2 = adi_m2 * xx
|
||
if (f1 > f2) then
|
||
cc2 = f1 - f2
|
||
cc = sqrt(cc2)
|
||
f1 = ayy * br2 / qi(idn)
|
||
f2 = ca2 + f1
|
||
f3 = f2 + cc2
|
||
f4 = f2 - cc2
|
||
f5 = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1)
|
||
cf2 = 5.0d-01 * (f3 + f5)
|
||
cs2 = cc2 * ca2 / cf2
|
||
else
|
||
cf2 = ca2 + ayy * br2 / qi(idn)
|
||
cc2 = 0.0d+00
|
||
cc = 0.0d+00
|
||
cs2 = 0.0d+00
|
||
end if
|
||
|
||
if (cs2 < cc2 .and. cc2 < cf2) then
|
||
f1 = (cc2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (cc2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), qi(ivx))
|
||
ca = sign(sqrt(ca2), qi(ivx))
|
||
cs = sign(sqrt(cs2), qi(ivx))
|
||
|
||
lm(1) = qi(ivx) + cf
|
||
lm(2) = qi(ivx) + ca
|
||
lm(3) = qi(ivx) + cs
|
||
lm(4) = qi(ivx)
|
||
lm(5) = qi(ivx) - cs
|
||
lm(6) = qi(ivx) - ca
|
||
lm(7) = qi(ivx) - cf
|
||
|
||
vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty
|
||
|
||
sqrtd = sqrt(qi(idn))
|
||
qf = cf * sign(alf, qi(ibx))
|
||
qs = cs * sign(als, qi(ibx))
|
||
f1 = cc / sqrtd
|
||
af_prime = f1 * alf
|
||
as_prime = f1 * als
|
||
f1 = br / sqty
|
||
afpbb = af_prime * f1
|
||
aspbb = as_prime * f1
|
||
|
||
f1 = qs / sqty
|
||
f2 = qs * vbet
|
||
f3 = as_prime / sqty
|
||
rvec(1,1) = alf
|
||
rvec(1,2) = alf * lm(1)
|
||
rvec(1,3) = alf * qi(ivy) - f1 * bty
|
||
rvec(1,4) = alf * qi(ivz) - f1 * btz
|
||
rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb
|
||
rvec(1,6) = f3 * bty
|
||
rvec(1,7) = f3 * btz
|
||
|
||
rvec(7,1) = alf
|
||
rvec(7,2) = alf * lm(7)
|
||
rvec(7,3) = alf * qi(ivy) + f1 * bty
|
||
rvec(7,4) = alf * qi(ivz) + f1 * btz
|
||
rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb
|
||
rvec(7,6) = rvec(1,6)
|
||
rvec(7,7) = rvec(1,7)
|
||
|
||
f1 = qf / sqty
|
||
f2 = qf * vbet
|
||
f3 = - af_prime / sqty
|
||
rvec(3,1) = als
|
||
rvec(3,2) = als * lm(3)
|
||
rvec(3,3) = als * qi(ivy) + f1 * bty
|
||
rvec(3,4) = als * qi(ivz) + f1 * btz
|
||
rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb
|
||
rvec(3,6) = f3 * bty
|
||
rvec(3,7) = f3 * btz
|
||
|
||
rvec(5,1) = als
|
||
rvec(5,2) = als * lm(5)
|
||
rvec(5,3) = als * qi(ivy) - f1 * bty
|
||
rvec(5,4) = als * qi(ivz) - f1 * btz
|
||
rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb
|
||
rvec(5,6) = rvec(3,6)
|
||
rvec(5,7) = rvec(3,7)
|
||
|
||
f1 = sign(1.0d+00, qi(ivx))
|
||
f2 = sign(1.0d+00 / sqrtd, qi(ibx))
|
||
rvec(2,3) = f1 * btz
|
||
rvec(2,4) = - f1 * bty
|
||
rvec(2,5) = f1 * (qi(ivy) * btz - qi(ivz) * bty)
|
||
rvec(2,6) = - f2 * btz
|
||
rvec(2,7) = f2 * bty
|
||
|
||
rvec(6,3) = - rvec(2,3)
|
||
rvec(6,4) = - rvec(2,4)
|
||
rvec(6,5) = - rvec(2,5)
|
||
rvec(6,6) = rvec(2,6)
|
||
rvec(6,7) = rvec(2,7)
|
||
|
||
rvec(4,2) = qi(ivx)
|
||
rvec(4,3) = qi(ivy)
|
||
rvec(4,4) = qi(ivz)
|
||
rvec(4,5) = v2h + adi_m2d1 * xx
|
||
|
||
norm = 2.0d+00 * cc2
|
||
f1 = sqty * br / norm
|
||
afpb = af_prime * f1
|
||
aspb = as_prime * f1
|
||
sqty = sqty / norm
|
||
vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty
|
||
|
||
f1 = adi_m1 * alf / norm
|
||
f2 = qs * sqty
|
||
f3 = as_prime * qi(idn) * sqty
|
||
f4 = alf * cf / norm
|
||
f5 = qs * vqstr
|
||
lvec(1,1) = f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb
|
||
lvec(1,2) = - f1 * qi(ivx) + f4
|
||
lvec(1,3) = - f1 * qi(ivy) - f2 * bty
|
||
lvec(1,4) = - f1 * qi(ivz) - f2 * btz
|
||
lvec(1,5) = f1
|
||
lvec(1,6) = f3 * bty - f1 * qi(iby)
|
||
lvec(1,7) = f3 * btz - f1 * qi(ibz)
|
||
|
||
lvec(7,1) = f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb
|
||
lvec(7,2) = - f1 * qi(ivx) - f4
|
||
lvec(7,3) = - f1 * qi(ivy) + f2 * bty
|
||
lvec(7,4) = - f1 * qi(ivz) + f2 * btz
|
||
lvec(7,5) = f1
|
||
lvec(7,6) = lvec(1,6)
|
||
lvec(7,7) = lvec(1,7)
|
||
|
||
f1 = adi_m1 * als / norm
|
||
f2 = qf * sqty
|
||
f3 = af_prime * qi(idn) * sqty
|
||
f4 = als * cs / norm
|
||
f5 = qf * vqstr
|
||
lvec(3,1) = f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb
|
||
lvec(3,2) = - f1 * qi(ivx) + f4
|
||
lvec(3,3) = - f1 * qi(ivy) + f2 * bty
|
||
lvec(3,4) = - f1 * qi(ivz) + f2 * btz
|
||
lvec(3,5) = f1
|
||
lvec(3,6) = - f3 * bty - f1 * qi(iby)
|
||
lvec(3,7) = - f3 * btz - f1 * qi(ibz)
|
||
|
||
lvec(5,1) = f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb
|
||
lvec(5,2) = - f1 * qi(ivx) - f4
|
||
lvec(5,3) = - f1 * qi(ivy) - f2 * bty
|
||
lvec(5,4) = - f1 * qi(ivz) - f2 * btz
|
||
lvec(5,5) = f1
|
||
lvec(5,6) = lvec(3,6)
|
||
lvec(5,7) = lvec(3,7)
|
||
|
||
f1 = sign(5.0d-01, qi(ivx)) * bty
|
||
f2 = sign(5.0d-01, qi(ivx)) * btz
|
||
f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx))
|
||
lvec(2,1) = qi(ivz) * f1 - qi(ivy) * f2
|
||
lvec(2,3) = f2
|
||
lvec(2,4) = - f1
|
||
lvec(2,6) = - f2 * f3
|
||
lvec(2,7) = f1 * f3
|
||
|
||
lvec(6,1) = - lvec(2,1)
|
||
lvec(6,3) = - lvec(2,3)
|
||
lvec(6,4) = - lvec(2,4)
|
||
lvec(6,6) = lvec(2,6)
|
||
lvec(6,7) = lvec(2,7)
|
||
|
||
f1 = 1.0d+00 / cc2
|
||
lvec(4,1) = 1.0d+00 - f1 * (v2h - adi_m2d1 * xx)
|
||
lvec(4,2) = f1 * qi(ivx)
|
||
lvec(4,3) = f1 * qi(ivy)
|
||
lvec(4,4) = f1 * qi(ivz)
|
||
lvec(4,5) = - f1
|
||
lvec(4,6) = f1 * qi(iby)
|
||
lvec(4,7) = f1 * qi(ibz)
|
||
|
||
al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i))
|
||
df = matmul(al, rvec)
|
||
fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:))
|
||
fi(ibx,i) = fl(ibx,i)
|
||
fi(ibp,i) = fl(ibp,i)
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_roe
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_MHD_ADI_KEPES:
|
||
! --------------------------------
|
||
!
|
||
! Subroutine solves one dimensional adiabatic magnetohydrodynamic
|
||
! Riemann problem using the entropy stable KEPES method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M.,
|
||
! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic
|
||
! field divergence diminishing ideal magnetohydrodynamics equations",
|
||
! Journal of Computational Physics, 2018, 364, pp. 420-467,
|
||
! https://doi.org/10.1016/j.jcp.2018.03.002
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_mhd_adi_kepes(ql, qr, fi)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, adiabatic_index, cglm
|
||
use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, &
|
||
ibx, iby, ibz, ibp, ipr, ien
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i
|
||
real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa, dnl, prl, pta
|
||
real(kind=8) :: v2l, v2r, b2l, b2r, bl, br, bta, btl, bp2, eka, ema, vva
|
||
real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, aa2, x2, x3, ub2, uba
|
||
real(kind=8) :: alf, als, das, daf, csq, wps, wms, wpf, wmf
|
||
real(kind=8) :: f1, f2, f3, f4, f5
|
||
|
||
real(kind=8), dimension(nf) :: v, lm, tm
|
||
|
||
logical , save :: first = .true.
|
||
real(kind=8) , save :: adi_m1, adi_m1x, adi_m1xi
|
||
real(kind=8), dimension(9,9), save :: rm
|
||
!$omp threadprivate(first, adi_m1, adi_m1x, adi_m1xi, rm)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (first) then
|
||
adi_m1 = adiabatic_index - 1.0d+00
|
||
adi_m1x = adi_m1 / adiabatic_index
|
||
adi_m1xi = adiabatic_index / adi_m1
|
||
|
||
rm(:,:) = 0.0d+00
|
||
rm(1,5) = 1.0d+00
|
||
rm(6,4) = 1.0d+00
|
||
rm(6,6) = 1.0d+00
|
||
|
||
first = .false.
|
||
end if
|
||
|
||
do i = 1, nn
|
||
|
||
v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i))
|
||
v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i))
|
||
b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i))
|
||
bl = ql(idn,i) / ql(ipr,i)
|
||
br = qr(idn,i) / qr(ipr,i)
|
||
|
||
btl = lmean(bl, br)
|
||
bta = amean(bl, br)
|
||
dnl = lmean(ql(idn,i), qr(idn,i))
|
||
dna = amean(ql(idn,i), qr(idn,i))
|
||
vxa = amean(ql(ivx,i), qr(ivx,i))
|
||
vya = amean(ql(ivy,i), qr(ivy,i))
|
||
vza = amean(ql(ivz,i), qr(ivz,i))
|
||
bxa = amean(ql(ibx,i), qr(ibx,i))
|
||
bya = amean(ql(iby,i), qr(iby,i))
|
||
bza = amean(ql(ibz,i), qr(ibz,i))
|
||
bpa = amean(ql(ibp,i), qr(ibp,i))
|
||
prl = lmean(ql(ipr,i), qr(ipr,i))
|
||
pra = dna / bta
|
||
eka = 5.0d-01 * amean(v2l, v2r)
|
||
ema = 5.0d-01 * amean(b2l, b2r)
|
||
pta = pra + ema
|
||
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
|
||
ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r)
|
||
uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), &
|
||
sum(qr(ivx:ivz,i) * qr(ibx:ibz,i)))
|
||
|
||
v(idn) = 5.0d-01 * (bl * v2l - br * v2r)
|
||
if (qr(idn,i) > ql(idn,i)) then
|
||
v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i))
|
||
else if (ql(idn,i) > qr(idn,i)) then
|
||
v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i))
|
||
end if
|
||
if (br > bl) then
|
||
v(idn) = v(idn) + log(br / bl) / adi_m1
|
||
else if (bl > br) then
|
||
v(idn) = v(idn) - log(bl / br) / adi_m1
|
||
end if
|
||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||
v(ipr) = bl - br
|
||
v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
|
||
v(iby) = br * qr(iby,i) - bl * ql(iby,i)
|
||
v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
|
||
v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)
|
||
|
||
bp2 = bya * bya + bza * bza
|
||
if (bp2 > 0.0d+00) then
|
||
f1 = sqrt(bp2)
|
||
x2 = bya / f1
|
||
x3 = bza / f1
|
||
else
|
||
x2 = 0.0d+00
|
||
x3 = 0.0d+00
|
||
end if
|
||
|
||
aa2 = adiabatic_index * pra / dnl
|
||
ca2 = bxa * bxa / dnl
|
||
f1 = bp2 / dnl
|
||
f2 = ca2 + f1
|
||
f3 = f2 - aa2
|
||
f4 = sqrt(f3 * f3 + 4.0d+00 * aa2 * f1)
|
||
cf2 = 5.0d-01 * (f2 + aa2 + f4)
|
||
cs2 = aa2 * ca2 / cf2
|
||
|
||
if (cs2 < aa2 .and. aa2 < cf2) then
|
||
f1 = (aa2 - cs2) / (cf2 - cs2)
|
||
alf = sqrt(f1)
|
||
als = sqrt(1.0d+00 - f1)
|
||
else if (aa2 >= cf2) then
|
||
alf = 1.0d+00
|
||
als = 0.0d+00
|
||
else
|
||
alf = 0.0d+00
|
||
als = 1.0d+00
|
||
end if
|
||
|
||
cf = sign(sqrt(cf2), vxa)
|
||
cs = sign(sqrt(cs2), vxa)
|
||
ca = sign(sqrt(ca2), vxa)
|
||
ch = sign( cglm, vxa)
|
||
|
||
lm(1) = vxa + cf
|
||
lm(2) = vxa + ca
|
||
lm(3) = vxa + cs
|
||
lm(4) = vxa + ch
|
||
lm(5) = vxa
|
||
lm(6) = vxa - ch
|
||
lm(7) = vxa - cs
|
||
lm(8) = vxa - ca
|
||
lm(9) = vxa - cf
|
||
|
||
f1 = 5.0d-01 / bta
|
||
f2 = adiabatic_index * dnl
|
||
tm(1) = 5.0d-01 / f2
|
||
tm(2) = f1 / dnl**2
|
||
tm(3) = tm(1)
|
||
tm(4) = f1
|
||
tm(5) = dnl * adi_m1x
|
||
tm(6) = tm(4)
|
||
tm(7) = tm(1)
|
||
tm(8) = tm(2)
|
||
tm(9) = tm(1)
|
||
|
||
das = dnl * als
|
||
daf = dnl * alf
|
||
csq = sqrt(f2 / bta)
|
||
|
||
f3 = vva + adi_m1xi * prl / dnl
|
||
f4 = vya * x2 + vza * x3
|
||
f5 = sqrt(adiabatic_index * bp2 / (dnl * bta))
|
||
f1 = dnl * (als * f3 - alf * f5)
|
||
f2 = dnl * (als * cs * vxa + sign(alf, bxa) * cf * f4)
|
||
wps = f1 + f2
|
||
wms = f1 - f2
|
||
f1 = dnl * (alf * f3 + als * f5)
|
||
f2 = dnl * (alf * cf * vxa - sign(als, bxa) * cs * f4)
|
||
wpf = f1 + f2
|
||
wmf = f1 - f2
|
||
|
||
f1 = daf
|
||
f2 = sign(das, bxa) * cs
|
||
f3 = als * csq
|
||
rm(1,1) = f1
|
||
rm(2,1) = f1 * lm(1)
|
||
rm(3,1) = f1 * vya - f2 * x2
|
||
rm(4,1) = f1 * vza - f2 * x3
|
||
rm(5,1) = wpf
|
||
rm(7,1) = f3 * x2
|
||
rm(8,1) = f3 * x3
|
||
|
||
rm(1,9) = f1
|
||
rm(2,9) = f1 * lm(9)
|
||
rm(3,9) = f1 * vya + f2 * x2
|
||
rm(4,9) = f1 * vza + f2 * x3
|
||
rm(5,9) = wmf
|
||
rm(7,9) = rm(7,1)
|
||
rm(8,9) = rm(8,1)
|
||
|
||
f1 = sign(dnl * sqrt(dna), vxa)
|
||
rm(3,2) = f1 * x3
|
||
rm(4,2) = - f1 * x2
|
||
rm(5,2) = - f1 * (x2 * vza - x3 * vya)
|
||
rm(7,2) = - dnl * x3
|
||
rm(8,2) = dnl * x2
|
||
|
||
rm(3,8) = - rm(3,2)
|
||
rm(4,8) = - rm(4,2)
|
||
rm(5,8) = - rm(5,2)
|
||
rm(7,8) = rm(7,2)
|
||
rm(8,8) = rm(8,2)
|
||
|
||
f1 = das
|
||
f2 = sign(daf, bxa) * cf
|
||
f3 = - alf * csq
|
||
rm(1,3) = f1
|
||
rm(2,3) = f1 * lm(3)
|
||
rm(3,3) = f1 * vya + f2 * x2
|
||
rm(4,3) = f1 * vza + f2 * x3
|
||
rm(5,3) = wps
|
||
rm(7,3) = f3 * x2
|
||
rm(8,3) = f3 * x3
|
||
|
||
rm(1,7) = f1
|
||
rm(2,7) = f1 * lm(7)
|
||
rm(3,7) = f1 * vya - f2 * x2
|
||
rm(4,7) = f1 * vza - f2 * x3
|
||
rm(5,7) = wms
|
||
rm(7,7) = rm(7,3)
|
||
rm(8,7) = rm(8,3)
|
||
|
||
f1 = sign(1.0d+00, vxa)
|
||
rm(5,4) = bxa + bpa * f1
|
||
rm(5,6) = bxa - bpa * f1
|
||
rm(9,4) = f1
|
||
rm(9,6) = - f1
|
||
|
||
rm(2,5) = vxa
|
||
rm(3,5) = vya
|
||
rm(4,5) = vza
|
||
rm(5,5) = vva
|
||
|
||
fi(idn,i) = dnl * vxa
|
||
fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta
|
||
fi(imy,i) = fi(idn,i) * vya - bxa * bya
|
||
fi(imz,i) = fi(idn,i) * vza - bxa * bza
|
||
fi(ibx,i) = cglm * bpa
|
||
fi(iby,i) = vxa * bya - bxa * vya
|
||
fi(ibz,i) = vxa * bza - bxa * vza
|
||
fi(ibp,i) = cglm * bxa
|
||
fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
|
||
+ (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) &
|
||
+ (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) &
|
||
+ fi(ibp,i) * bpa - ub2 + bxa * uba &
|
||
- cglm * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
|
||
|
||
fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_mhd_adi_kepes
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! RELATIVISTIC ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine RIEMANN_SRMHD_HLLC:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine solves one dimensional Riemann problem using
|
||
! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC)
|
||
! by Mignone & Bodo.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! ql, qr - the primitive variables of the left and right Riemann states;
|
||
! fi - the numerical flux at the cell interface;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A. & Bodo, G.
|
||
! "An HLLC Riemann solver for relativistic flows - II.
|
||
! Magnetohydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society,
|
||
! 2006, Volume 368, Pages 1040-1054
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine riemann_srmhd_hllc(ql, qr, fi)
|
||
|
||
use algebra , only : quadratic
|
||
use coordinates, only : nn => bcells
|
||
use equations , only : nf, idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp
|
||
use equations , only : prim2cons, fluxspeed
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:), intent(in) :: ql, qr
|
||
real(kind=8), dimension(:,:), intent(out) :: fi
|
||
|
||
integer :: i, nr
|
||
real(kind=8) :: sl, sr, srml, sm
|
||
real(kind=8) :: bx2
|
||
real(kind=8) :: pt, dv, fc, uu, ff, uf
|
||
real(kind=8) :: vv, vb, gi
|
||
|
||
real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr
|
||
real(kind=8), dimension( 2,nn) :: cl, cr
|
||
real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh
|
||
real(kind=8), dimension(3) :: a, vs
|
||
real(kind=8), dimension(2) :: x
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
call prim2cons(ql, ul)
|
||
call prim2cons(qr, ur)
|
||
|
||
call fluxspeed(ql, ul, fl, cl)
|
||
call fluxspeed(qr, ur, fr, cr)
|
||
|
||
do i = 1, nn
|
||
|
||
sl = min(cl(1,i), cr(1,i))
|
||
sr = max(cl(2,i), cr(2,i))
|
||
|
||
if (sl >= 0.0d+00) then
|
||
|
||
fi(:,i) = fl(:,i)
|
||
|
||
else if (sr <= 0.0d+00) then
|
||
|
||
fi(:,i) = fr(:,i)
|
||
|
||
else ! sl < 0 < sr
|
||
|
||
srml = sr - sl
|
||
|
||
wl(:) = sl * ul(:,i) - fl(:,i)
|
||
wr(:) = sr * ur(:,i) - fr(:,i)
|
||
|
||
uh(:) = ( wr(:) - wl(:)) / srml
|
||
fh(:) = (sl * wr(:) - sr * wl(:)) / srml
|
||
|
||
wl(ien) = wl(ien) + wl(idn)
|
||
wr(ien) = wr(ien) + wr(idn)
|
||
|
||
! calculate Bₓ²
|
||
!
|
||
bx2 = ql(ibx,i) * qr(ibx,i)
|
||
|
||
! calculate the contact dicontinuity speed solving eq. (41)
|
||
!
|
||
if (bx2 >= 1.0d-16) then
|
||
|
||
! prepare the quadratic coefficients, (eq. 42 in [1])
|
||
!
|
||
uu = sum(uh(iby:ibz) * uh(iby:ibz))
|
||
ff = sum(fh(iby:ibz) * fh(iby:ibz))
|
||
uf = sum(uh(iby:ibz) * fh(iby:ibz))
|
||
|
||
a(1) = uh(imx) - uf
|
||
a(2) = uu + ff - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn) - uf
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! substitute magnetic field components from the HLL state (eq. 37 in [1])
|
||
!
|
||
us(ibx) = ql(ibx,i)
|
||
us(iby) = uh(iby)
|
||
us(ibz) = uh(ibz)
|
||
us(ibp) = ql(ibp,i)
|
||
|
||
! calculate velocity components without Bₓ normalization (eq. 38 in [1])
|
||
!
|
||
vs(1) = sm
|
||
vs(2) = (us(iby) * sm - fh(iby)) / us(ibx)
|
||
vs(3) = (us(ibz) * sm - fh(ibz)) / us(ibx)
|
||
|
||
! calculate v² and v.B
|
||
!
|
||
vv = sum(vs(1:3) * vs( 1: 3))
|
||
vb = sum(vs(1:3) * us(ibx:ibz))
|
||
|
||
! calculate inverse of Lorentz factor
|
||
!
|
||
gi = 1.0d+00 - vv
|
||
|
||
! calculate total pressure (eq. 40 in [1])
|
||
!
|
||
pt = fh(imx) + gi * bx2 - (fh(ien) + fh(idn) - us(ibx) * vb) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pt <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sl - sm
|
||
fc = (sl - ql(ivx,i)) / dv
|
||
us(idn) = fc * ul(idn,i)
|
||
us(imy) = (sl * ul(imy,i) - fl(imy,i) &
|
||
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
|
||
us(imz) = (sl * ul(imz,i) - fl(imz,i) &
|
||
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
|
||
us(ien) = (sl * (ul(ien,i) + ul(idn,i)) &
|
||
- (fl(ien,i) + fl(idn,i)) &
|
||
+ pt * sm - us(ibx) * vb) / dv - us(idn)
|
||
|
||
! calculate normal component of momentum
|
||
!
|
||
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sr - sm
|
||
fc = (sr - qr(ivx,i)) / dv
|
||
us(idn) = fc * ur(idn,i)
|
||
us(imy) = (sr * ur(imy,i) - fr(imy,i) &
|
||
- us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv
|
||
us(imz) = (sr * ur(imz,i) - fr(imz,i) &
|
||
- us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv
|
||
us(ien) = (sr * (ur(ien,i) + ur(idn,i)) &
|
||
- (fr(ien,i) + fr(idn,i)) &
|
||
+ pt * sm - us(ibx) * vb) / dv - us(idn)
|
||
|
||
! calculate normal component of momentum
|
||
!
|
||
us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb
|
||
|
||
! calculate the flux (eq. 14 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity so use HLL flux
|
||
!
|
||
fi(:,i) = fh(:)
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
else ! Bx² > 0
|
||
|
||
! prepare the quadratic coefficients (eq. 47 in [1])
|
||
!
|
||
a(1) = uh(imx)
|
||
a(2) = - (fh(imx) + uh(ien) + uh(idn))
|
||
a(3) = fh(ien) + fh(idn)
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! if Δ < 0, just use the HLL flux
|
||
!
|
||
if (nr < 1) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! get the contact dicontinuity speed
|
||
!
|
||
sm = x(1)
|
||
|
||
! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux
|
||
!
|
||
if ((sm <= sl) .or. (sm >= sr)) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! calculate total pressure (eq. 48 in [1])
|
||
!
|
||
pt = fh(imx) - (fh(ien) + fh(idn)) * sm
|
||
|
||
! if the pressure is negative, use the HLL flux
|
||
!
|
||
if (pt <= 0.0d+00) then
|
||
fi(:,i) = fh(:)
|
||
else
|
||
|
||
! depending in the sign of the contact dicontinuity speed, calculate the proper
|
||
! state and corresponding flux
|
||
!
|
||
if (sm > 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sl - sm
|
||
us(idn) = wl(idn) / dv
|
||
us(imy) = wl(imy) / dv
|
||
us(imz) = wl(imz) / dv
|
||
us(ien) = (wl(ien) + pt * sm) / dv
|
||
us(imx) = (us(ien) + pt) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
us(ibx) = 0.0d+00
|
||
us(iby) = wl(iby) / dv
|
||
us(ibz) = wl(ibz) / dv
|
||
us(ibp) = ql(ibp,i)
|
||
|
||
! calculate the flux (eq. 27 in [1])
|
||
!
|
||
fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i))
|
||
|
||
else if (sm < 0.0d+00) then
|
||
|
||
! calculate the conserved variable vector (eqs. 43-46 in [1])
|
||
!
|
||
dv = sr - sm
|
||
us(idn) = wr(idn) / dv
|
||
us(imy) = wr(imy) / dv
|
||
us(imz) = wr(imz) / dv
|
||
us(ien) = (wr(ien) + pt * sm) / dv
|
||
us(imx) = (us(ien) + pt) * sm
|
||
us(ien) = us(ien) - us(idn)
|
||
us(ibx) = 0.0d+00
|
||
us(iby) = wr(iby) / dv
|
||
us(ibz) = wr(ibz) / dv
|
||
us(ibp) = qr(ibp,i)
|
||
|
||
! calculate the flux (eq. 27 in [1])
|
||
!
|
||
fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i))
|
||
|
||
else
|
||
|
||
! intermediate flux is constant across the contact discontinuity so use HLL flux
|
||
!
|
||
fi(:,i) = fh(:)
|
||
|
||
end if ! sm == 0
|
||
|
||
end if ! p* < 0
|
||
|
||
end if ! sl < sm < sr
|
||
|
||
end if ! nr < 1
|
||
|
||
end if ! Bx² == 0 or > 0
|
||
|
||
end if ! sl < 0 < sr
|
||
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine riemann_srmhd_hllc
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! SUPPORTING SUBROUTINES/FUNCTIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine HIGHER_ORDER_FLUX_CORRECTION:
|
||
! ---------------------------------------
|
||
!
|
||
! Subroutine adds higher order corrections to the numerical fluxes.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! f - the vector of numerical fluxes;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine higher_order_flux_correction(f)
|
||
|
||
! include external procedures
|
||
!
|
||
use interpolations, only : order
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
real(kind=8), dimension(:,:), intent(inout) :: f
|
||
|
||
! local variables
|
||
!
|
||
integer :: n
|
||
|
||
! local arrays to store the states
|
||
!
|
||
real(kind=8), dimension(size(f,1),size(f,2)) :: f2, f4
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! depending on the scheme order calculate and add higher order corrections,
|
||
! if required
|
||
!
|
||
if (.not. high_order_flux .or. order <= 2) return
|
||
|
||
! get the length of vector
|
||
!
|
||
n = size(f, 2)
|
||
|
||
! 2nd order correction
|
||
!
|
||
f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01
|
||
f2(:,1 ) = 0.0d+00
|
||
f2(:, n ) = 0.0d+00
|
||
|
||
! 4th order correction
|
||
!
|
||
if (order > 4) then
|
||
|
||
f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) &
|
||
- 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) &
|
||
+ 6.0d+00 * f(:,3:n-2)) / 6.4d+02
|
||
f4(:,1 ) = 0.0d+00
|
||
f4(:,2 ) = 0.0d+00
|
||
f4(:, n-1) = 0.0d+00
|
||
f4(:, n ) = 0.0d+00
|
||
|
||
! correct the flux for 4th order
|
||
!
|
||
f(:,:) = f(:,:) + f4(:,:)
|
||
|
||
end if
|
||
|
||
! correct the flux for 2nd order
|
||
!
|
||
f(:,:) = f(:,:) + f2(:,:)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine higher_order_flux_correction
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function LMEAN:
|
||
! --------------
|
||
!
|
||
! Function calculates the logarithmic mean using the optimized algorithm
|
||
! by Ranocha et al.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! l, r - the values to calculate the mean of;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Ranocha et al.,
|
||
! "Efficient implementation of modern entropy stable and kinetic energy
|
||
! preserving discontinuous Galerkin methods for conservation laws",
|
||
! http://arxiv.org/abs/2112.10517v1
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function lmean(l, r)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: l, r
|
||
|
||
real(kind=8) :: u, d, s
|
||
|
||
real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, &
|
||
c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
d = r - l
|
||
s = r + l
|
||
u = (d / s)**2
|
||
if (u < 1.0d-04) then
|
||
lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
|
||
else
|
||
if (d >= 0.0d+00) then
|
||
s = log(r / l)
|
||
else
|
||
s = log(l / r)
|
||
end if
|
||
lmean = abs(d) / s
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function lmean
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! function AMEAN:
|
||
! --------------
|
||
!
|
||
! Function calculate the arithmetic mean.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! l, r - the values to calculate the mean of;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
real(kind=8) function amean(l, r)
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: l, r
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
amean = 5.0d-01 * (l + r)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end function amean
|
||
|
||
!===============================================================================
|
||
!
|
||
end module schemes
|