EQUATIONS: Add subroutines get_maximum_speeds().

These subroutines determine the maximum physical and characteristic
speeds, and the maximum eigenvalue in the system.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-01-06 22:44:21 -03:00
parent 226998e80e
commit 0cf2481992

View File

@ -75,6 +75,10 @@ module equations
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) :: maxspeed
end function
subroutine get_maximum_speeds_iface(qq, um, cm, lm)
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out):: um, cm, lm
end subroutine
subroutine esystem_roe_iface(x, y, q, c, r, l)
real(kind=8) , intent(in) :: x, y
real(kind=8), dimension(:) , intent(in) :: q
@ -101,6 +105,8 @@ module equations
!
procedure(maxspeed_iface) , pointer, save :: maxspeed => null()
procedure(get_maximum_speeds_iface), pointer, save :: get_maximum_speeds => null()
! pointer to the Roe eigensystem procedure
!
procedure(esystem_roe_iface), pointer, save :: eigensystem_roe => null()
@ -164,9 +170,10 @@ module equations
!
real(kind=8), save :: csnd = 1.0d+00, csnd2 = 1.0d+00
! maximum speed in the system
! maximum speeds in the system
!
real(kind=8), save :: cmax = 0.0d+00, cmax2 = 0.0d+00
real(kind=8), save :: umax = 0.0d+00, amax = 0.0d+00, cglm = 0.0d+00
! the lower limits for density and pressure to be treated as physical
!
@ -222,13 +229,13 @@ module equations
public :: prim2cons_mhd_adi, fluxspeed_mhd_adi
public :: prim2cons_srhd_adi, fluxspeed_srhd_adi
public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi
public :: maxspeed, reset_maxspeed, get_maxspeed
public :: maxspeed, reset_maxspeed, get_maximum_speeds
public :: eigensystem_roe
public :: update_primitive_variables
public :: fix_unphysical_cells, correct_unphysical_states
public :: adiabatic_index, relativistic, magnetized
public :: csnd, csnd2
public :: cmax, cmax2
public :: cmax, cmax2, umax, amax, cglm
public :: nv, nf, ns
public :: inx, iny, inz
public :: idn, ivx, ivy, ivz, imx, imy, imz
@ -361,11 +368,12 @@ module equations
! set pointers to subroutines
!
prim2cons => prim2cons_hd_iso
cons2prim => cons2prim_hd_iso
fluxspeed => fluxspeed_hd_iso
maxspeed => maxspeed_hd_iso
eigensystem_roe => esystem_roe_hd_iso
prim2cons => prim2cons_hd_iso
cons2prim => cons2prim_hd_iso
fluxspeed => fluxspeed_hd_iso
maxspeed => maxspeed_hd_iso
eigensystem_roe => esystem_roe_hd_iso
get_maximum_speeds => get_maximum_speeds_hd_iso
case("adi", "ADI", "adiabatic", "ADIABATIC")
@ -401,11 +409,12 @@ module equations
! set pointers to subroutines
!
prim2cons => prim2cons_hd_adi
cons2prim => cons2prim_hd_adi
fluxspeed => fluxspeed_hd_adi
maxspeed => maxspeed_hd_adi
eigensystem_roe => esystem_roe_hd_adi
prim2cons => prim2cons_hd_adi
cons2prim => cons2prim_hd_adi
fluxspeed => fluxspeed_hd_adi
maxspeed => maxspeed_hd_adi
eigensystem_roe => esystem_roe_hd_adi
get_maximum_speeds => get_maximum_speeds_hd_adi
! warn about the unimplemented equation of state
!
@ -527,11 +536,12 @@ module equations
! set pointers to the subroutines
!
prim2cons => prim2cons_mhd_iso
cons2prim => cons2prim_mhd_iso
fluxspeed => fluxspeed_mhd_iso
maxspeed => maxspeed_mhd_iso
eigensystem_roe => esystem_roe_mhd_iso
prim2cons => prim2cons_mhd_iso
cons2prim => cons2prim_mhd_iso
fluxspeed => fluxspeed_mhd_iso
maxspeed => maxspeed_mhd_iso
eigensystem_roe => esystem_roe_mhd_iso
get_maximum_speeds => get_maximum_speeds_mhd_iso
case("adi", "ADI", "adiabatic", "ADIABATIC")
@ -574,11 +584,12 @@ module equations
! set pointers to subroutines
!
prim2cons => prim2cons_mhd_adi
cons2prim => cons2prim_mhd_adi
fluxspeed => fluxspeed_mhd_adi
maxspeed => maxspeed_mhd_adi
eigensystem_roe => esystem_roe_mhd_adi
prim2cons => prim2cons_mhd_adi
cons2prim => cons2prim_mhd_adi
fluxspeed => fluxspeed_mhd_adi
maxspeed => maxspeed_mhd_adi
eigensystem_roe => esystem_roe_mhd_adi
get_maximum_speeds => get_maximum_speeds_mhd_adi
case default
@ -703,10 +714,11 @@ module equations
! set pointers to subroutines
!
prim2cons => prim2cons_srhd_adi
cons2prim => cons2prim_srhd_adi
fluxspeed => fluxspeed_srhd_adi
maxspeed => maxspeed_srhd_adi
prim2cons => prim2cons_srhd_adi
cons2prim => cons2prim_srhd_adi
fluxspeed => fluxspeed_srhd_adi
maxspeed => maxspeed_srhd_adi
get_maximum_speeds => get_maximum_speeds_srhd_adi
! warn about the unimplemented equation of state
!
@ -884,10 +896,11 @@ module equations
! set pointers to subroutines
!
prim2cons => prim2cons_srmhd_adi
cons2prim => cons2prim_srmhd_adi
fluxspeed => fluxspeed_srmhd_adi
maxspeed => maxspeed_srmhd_adi
prim2cons => prim2cons_srmhd_adi
cons2prim => cons2prim_srmhd_adi
fluxspeed => fluxspeed_srmhd_adi
maxspeed => maxspeed_srmhd_adi
get_maximum_speeds => get_maximum_speeds_srmhd_adi
! warn about the unimplemented equation of state
!
@ -1881,6 +1894,61 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_HD_ISO:
! ------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_hd_iso(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
integer :: i, j, k = 1
real(kind=8) :: vl, vu
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = csnd
lm = 0.0d+00
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
vl = minval(qq(ivx:ivz,i,j,k))
vu = maxval(qq(ivx:ivz,i,j,k))
vm = max(vm, abs(vl), abs(vu))
lm = max(lm, abs(vl - csnd), abs(vu + csnd))
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_hd_iso
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_HD_ISO:
! -----------------------------
!
@ -2255,6 +2323,64 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_HD_ADI:
! ------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_hd_adi(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
integer :: i, j, k = 1
real(kind=8) :: vl, vu, cc
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = 0.0d+00
lm = 0.0d+00
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
vl = minval(qq(ivx:ivz,i,j,k))
vu = maxval(qq(ivx:ivz,i,j,k))
vm = max(vm, abs(vl), abs(vu))
cc = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k))
cm = max(cm, cc)
lm = max(lm, abs(vl - cc), abs(vu + cc))
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_hd_adi
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_HD_ADI:
! -----------------------------
!
@ -2682,6 +2808,69 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_MHD_ISO:
! -------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_mhd_iso(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
integer :: i, j, k = 1
real(kind=8) :: vl, vu, cc, xx
real(kind=8), dimension(3) :: bb
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = 0.0d+00
lm = 0.0d+00
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
vl = minval(qq(ivx:ivz,i,j,k))
vu = maxval(qq(ivx:ivz,i,j,k))
vm = max(vm, abs(vl), abs(vu))
bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
xx = csnd2 + bb(1) + bb(2) + bb(3)
cc = max(0.0d+00, xx**2 - 4.0d+00 * csnd2 * bb(1))
cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
cm = max(cm, cc)
lm = max(lm, abs(vl - cc), abs(vu + cc))
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_mhd_iso
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_MHD_ISO:
! ------------------------------
!
@ -3236,6 +3425,70 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_MHD_ADI:
! -------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_mhd_adi(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
integer :: i, j, k = 1
real(kind=8) :: vl, vu, aa, cc, xx
real(kind=8), dimension(3) :: bb
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = 0.0d+00
lm = 0.0d+00
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
vl = minval(qq(ivx:ivz,i,j,k))
vu = maxval(qq(ivx:ivz,i,j,k))
vm = max(vm, abs(vl), abs(vu))
aa = adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)
bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
xx = aa + bb(1) + bb(2) + bb(3)
cc = max(0.0d+00, xx**2 - 4.0d+00 * aa * bb(1))
cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
cm = max(cm, cc)
lm = max(lm, abs(vl - cc), abs(vu + cc))
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_mhd_adi
!
!===============================================================================
!
! subroutine ESYSTEM_ROE_MHD_ADI:
! ------------------------------
!
@ -3860,6 +4113,72 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_SRHD_ADI:
! --------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_srhd_adi(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
integer :: i, j, k = 1
real(kind=8) :: vl, vu, vv, ww, aa, cc, ss, fc
real(kind=8), dimension(3) :: bb
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = 0.0d+00
lm = 0.0d+00
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
vl = minval(qq(ivx:ivz,i,j,k))
vu = maxval(qq(ivx:ivz,i,j,k))
vm = max(vm, abs(vl), abs(vu))
vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi
aa = adiabatic_index * qq(ipr,i,j,k) / ww
ss = aa * (1.0d+00 - vv) / (1.0d+00 - aa)
fc = 1.0d+00 + ss
cc = sqrt(ss * (fc - vv))
cm = max(cm, cc)
lm = max(lm, abs(vl - cc), abs(vu + cc))
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_srhd_adi
!
!===============================================================================
!
! subroutine NR_FUNCTION_SRHD_ADI_1D:
! ----------------------------------
!
@ -5064,6 +5383,42 @@ module equations
!
!===============================================================================
!
! subroutine GET_MAXIMUM_SPEEDS_SRMHD_ADI:
! ---------------------------------------
!
! Subroutine determines the maximum characteristic speed and eigenvalue
! in the input array of primitive variables.
!
! Arguments:
!
! qq - the input array of primitive variables;
! vm - the maximum physical speed;
! cm - the maximum characteristic speed;
! lm - the maximum eigenvalue;
!
!===============================================================================
!
subroutine get_maximum_speeds_srmhd_adi(qq, vm, cm, lm)
use coordinates, only : nb, ne
implicit none
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
real(kind=8) , intent(out) :: vm, cm, lm
!-------------------------------------------------------------------------------
!
vm = 0.0d+00
cm = 1.0d+00
lm = 1.0d+00
!-------------------------------------------------------------------------------
!
end subroutine get_maximum_speeds_srmhd_adi
!
!===============================================================================
!
! subroutine NR_VELOCITY_SRMHD_ADI_1D:
! -----------------------------------
!