From 0cf24819920992f0039f884d91e95fef22fde82a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 6 Jan 2022 22:44:21 -0300 Subject: [PATCH] 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 --- sources/equations.F90 | 417 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 386 insertions(+), 31 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index eed5ec1..b74c3d9 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -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: ! ----------------------------------- !