Merge branch 'master' into reconnection
This commit is contained in:
commit
650b5fe8f2
@ -215,11 +215,12 @@ class Amun:
|
||||
if all(v in self.variables for v in ['dens', 'pres']):
|
||||
self.variables['temperature'] = 'temp'
|
||||
if self.attributes['eqsys'] in [ 'hd', 'mhd' ] \
|
||||
and all(v in self.variables for v in ['dens','velo']):
|
||||
and all(v in self.variables for v in ['dens','velx','vely','velz']):
|
||||
self.variables['kinetic energy'] = 'ekin'
|
||||
if self.attributes['eqsys'] in [ 'mhd', 'srmhd' ] \
|
||||
and all(v in self.variables for v in ['magx','magy','magz']):
|
||||
self.variables['magnetic energy'] = 'emag'
|
||||
self.variables['magnetic pressure'] = 'emag'
|
||||
if all(v in self.variables for v in ['velx','vely','velz', 'magx','magy','magz']) in self.variables:
|
||||
self.variables['electric field'] = 'evec'
|
||||
self.variables['electric field magnitude'] = 'elec'
|
||||
|
||||
@ -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)
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out):: um, cm
|
||||
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,9 @@ 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 +171,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, cglm = 0.0d+00
|
||||
|
||||
! the lower limits for density and pressure to be treated as physical
|
||||
!
|
||||
@ -222,13 +230,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, cglm
|
||||
public :: nv, nf, ns
|
||||
public :: inx, iny, inz
|
||||
public :: idn, ivx, ivy, ivz, imx, imy, imz
|
||||
@ -361,11 +369,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 +410,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 +537,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 +585,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 +715,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 +897,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 +1895,58 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_hd_iso(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
integer :: i, j, k = 1
|
||||
real(kind=8) :: vl, vu
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
vm = 0.0d+00
|
||||
cm = 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))
|
||||
cm = max(cm, 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 +2321,60 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_hd_adi(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
integer :: i, j, k = 1
|
||||
real(kind=8) :: vl, vu, cc
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
vm = 0.0d+00
|
||||
cm = 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, 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:
|
||||
! -----------------------------
|
||||
!
|
||||
@ -2295,7 +2415,7 @@ module equations
|
||||
real(kind=8), dimension(5,5), save :: lvec, rvec
|
||||
!$omp threadprivate(first, lvec, rvec)
|
||||
|
||||
real(kind=8) :: vv, vh, c2, cc, vc, nd, fc, fh, f1, f2, fv, fx
|
||||
real(kind=8) :: vv, vh, c2, cc, vc, fc, fh, f1, f2, fv, fx
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -2682,6 +2802,65 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_mhd_iso(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
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
|
||||
|
||||
#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, 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:
|
||||
! ------------------------------
|
||||
!
|
||||
@ -2960,24 +3139,16 @@ module equations
|
||||
!
|
||||
subroutine prim2cons_mhd_adi(q, u, s)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8) :: ei, ek, em
|
||||
!
|
||||
real(kind=8) :: ei, ek, em, ep
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
|
||||
@ -2990,15 +3161,13 @@ module equations
|
||||
u(ibz,i) = q(ibz,i)
|
||||
u(ibp,i) = q(ibp,i)
|
||||
ei = gammam1i * q(ipr,i)
|
||||
ek = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i) &
|
||||
+ u(imz,i) * q(ivz,i))
|
||||
em = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
u(ien,i) = ei + ek + em
|
||||
ek = sum(u(imx:imz,i) * q(ivx:ivz,i))
|
||||
em = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
ep = q(ibp,i) * q(ibp,i)
|
||||
u(ien,i) = ei + 5.0d-01 * (ek + em + ep)
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
@ -3036,7 +3205,7 @@ module equations
|
||||
integer , intent(out) :: s
|
||||
|
||||
integer :: i, p
|
||||
real(kind=8) :: ei, ek, em
|
||||
real(kind=8) :: ei, ek, em, ep
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3053,9 +3222,10 @@ module equations
|
||||
q(iby,i) = u(iby,i)
|
||||
q(ibz,i) = u(ibz,i)
|
||||
q(ibp,i) = u(ibp,i)
|
||||
ek = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
|
||||
em = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
ei = u(ien,i) - (ek + em)
|
||||
ek = sum(u(imx:imz,i) * q(ivx:ivz,i))
|
||||
em = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
ep = q(ibp,i) * q(ibp,i)
|
||||
ei = u(ien,i) - 5.0d-01 * (ek + em + ep)
|
||||
if (ei > 0.0d+00) then
|
||||
q(ipr,i) = gammam1 * ei
|
||||
else
|
||||
@ -3245,6 +3415,66 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_mhd_adi(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
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
|
||||
|
||||
#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, 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:
|
||||
! ------------------------------
|
||||
!
|
||||
@ -3869,6 +4099,65 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_srhd_adi(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
integer :: i, j, k = 1
|
||||
real(kind=8) :: vl, vu, vv, ww, aa, cc, ss, fc
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
vm = 0.0d+00
|
||||
cm = 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, 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:
|
||||
! ----------------------------------
|
||||
!
|
||||
@ -5073,6 +5362,57 @@ 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 eigenvalue;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine get_maximum_speeds_srmhd_adi(qq, vm, cm)
|
||||
|
||||
use coordinates, only : nb, ne
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||||
real(kind=8) , intent(out) :: vm, cm
|
||||
|
||||
integer :: i, j, k = 1
|
||||
real(kind=8) :: vl, vu
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
vm = 0.0d+00
|
||||
cm = 1.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))
|
||||
end do
|
||||
end do
|
||||
#if NDIMS == 3
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine get_maximum_speeds_srmhd_adi
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine NR_VELOCITY_SRMHD_ADI_1D:
|
||||
! -----------------------------------
|
||||
!
|
||||
|
||||
@ -1170,7 +1170,7 @@ module evolution
|
||||
#if NDIMS == 3
|
||||
use coordinates, only : adz
|
||||
#endif /* NDIMS == 3 */
|
||||
use equations , only : nf, maxspeed, cmax, cmax2
|
||||
use equations , only : nf, get_maximum_speeds, cmax2, cmax, umax, cglm
|
||||
use helpers , only : print_message
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_maximum, reduce_sum
|
||||
@ -1183,7 +1183,7 @@ module evolution
|
||||
|
||||
logical :: flag
|
||||
integer :: l, m, n, status
|
||||
real(kind=8) :: cm, dx_min, fnorm, h0, h1
|
||||
real(kind=8) :: um, cm, dx_min, fnorm, h0, h1
|
||||
real(kind=8), dimension(3) :: d
|
||||
|
||||
real(kind=8), dimension(:,:), allocatable :: s
|
||||
@ -1200,6 +1200,7 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
umax = 0.0d+00
|
||||
cmax = eps
|
||||
m = 1
|
||||
n = get_dblocks()
|
||||
@ -1208,7 +1209,9 @@ module evolution
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
|
||||
cm = maxspeed(pdata%q)
|
||||
call get_maximum_speeds(pdata%q, um, cm)
|
||||
!$omp atomic update
|
||||
umax = max(umax, um)
|
||||
!$omp atomic update
|
||||
cmax = max(cmax, cm)
|
||||
!$omp atomic update
|
||||
@ -1217,10 +1220,12 @@ module evolution
|
||||
!$omp end parallel do
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(cmax)
|
||||
call reduce_maximum(m)
|
||||
#endif /* MPI */
|
||||
|
||||
cglm = cmax - umax
|
||||
cmax2 = cmax * cmax
|
||||
|
||||
#if NDIMS == 2
|
||||
@ -1398,7 +1403,7 @@ module evolution
|
||||
#if NDIMS == 3
|
||||
use coordinates , only : adz
|
||||
#endif /* NDIMS == 3 */
|
||||
use equations , only : maxspeed, cmax, cmax2
|
||||
use equations , only : get_maximum_speeds, cmax, cmax2, umax, cglm
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_maximum
|
||||
#endif /* MPI */
|
||||
@ -1409,12 +1414,13 @@ module evolution
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
integer :: l, n, m
|
||||
real(kind=8) :: cm, dx_min
|
||||
real(kind=8) :: um, cm, dx_min
|
||||
|
||||
real(kind=8), parameter :: eps = tiny(cmax)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
umax = 0.0d+00
|
||||
cmax = eps
|
||||
m = 1
|
||||
n = get_dblocks()
|
||||
@ -1423,7 +1429,9 @@ module evolution
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
|
||||
cm = maxspeed(pdata%q)
|
||||
call get_maximum_speeds(pdata%q, um, cm)
|
||||
!$omp atomic update
|
||||
umax = max(umax, um)
|
||||
!$omp atomic update
|
||||
cmax = max(cmax, cm)
|
||||
!$omp atomic update
|
||||
@ -1432,10 +1440,12 @@ module evolution
|
||||
!$omp end parallel do
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(cmax)
|
||||
call reduce_maximum(m)
|
||||
#endif /* MPI */
|
||||
|
||||
cglm = cmax - umax
|
||||
cmax2 = cmax * cmax
|
||||
|
||||
#if NDIMS == 2
|
||||
|
||||
@ -2060,11 +2060,12 @@ module schemes
|
||||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||||
|
||||
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
|
||||
dvl = smr - sml
|
||||
f(:,i) = sml * wcr(:) / dvl - smr * wcl(:) / dvl
|
||||
|
||||
else ! Bx = 0 -> Sₘ = 0
|
||||
|
||||
f(:,i) = (sl * wr(:) - sr *wl(:)) / srml
|
||||
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
|
||||
|
||||
end if
|
||||
|
||||
@ -2085,11 +2086,12 @@ module schemes
|
||||
|
||||
wcl(:) = (sml - sl) * uil(:) + wl(:)
|
||||
|
||||
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
|
||||
dvl = sr - sml
|
||||
f(:,i) = sml * wr(:) / dvl - sr * wcl(:) / dvl
|
||||
|
||||
else ! Bx = 0 -> Sₘ = 0
|
||||
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
|
||||
|
||||
end if ! sl* < 0
|
||||
|
||||
@ -2108,17 +2110,18 @@ module schemes
|
||||
|
||||
wcr(:) = (smr - sr) * uir(:) + wr(:)
|
||||
|
||||
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
|
||||
dvr = smr - sl
|
||||
f(:,i) = sl * wcr(:) / dvr - smr * wl(:) / dvr
|
||||
|
||||
else ! Bx = 0 -> Sₘ = 0
|
||||
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
|
||||
|
||||
end if ! sr* > 0
|
||||
|
||||
else ! sl* < sl & sr* > sr
|
||||
|
||||
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
|
||||
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
|
||||
|
||||
end if
|
||||
|
||||
@ -3244,8 +3247,8 @@ module schemes
|
||||
!
|
||||
subroutine riemann_hd_iso_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
|
||||
|
||||
use equations, only : idn, ivx, ivy, ivz
|
||||
use equations, only : csnd, csnd2
|
||||
use equations, only : idn, ivx, ivy, ivz, imx, imy, imz
|
||||
use equations, only : nf, csnd, csnd2
|
||||
|
||||
implicit none
|
||||
|
||||
@ -3256,18 +3259,18 @@ module schemes
|
||||
real(kind=8), dimension(4,4), save :: rm
|
||||
!$omp threadprivate(first, rm)
|
||||
|
||||
integer :: i
|
||||
real(kind=8) :: dnl, vxa, vya, vza
|
||||
integer :: i, p
|
||||
real(kind=8) :: dnl, dna, vxa, vya, vza
|
||||
|
||||
real(kind=8), dimension(4) :: v, lm, tm
|
||||
real(kind=8), dimension(nf) :: v, lm, tm, lv
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
if (first) then
|
||||
rm(:,:) = 0.0d+00
|
||||
rm(1,:) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||||
rm(3,2) = 1.0d+00
|
||||
rm(4,3) = 1.0d+00
|
||||
rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ]
|
||||
rm(3,:) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ]
|
||||
rm(4,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ]
|
||||
|
||||
first = .false.
|
||||
end if
|
||||
@ -3277,16 +3280,25 @@ module schemes
|
||||
! the intermediate state averages
|
||||
!
|
||||
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))
|
||||
|
||||
! the difference in entropy variables between states
|
||||
!
|
||||
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(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) &
|
||||
- (vxa * v(ivx) + vya * v(ivy) + vza * v(ivz))
|
||||
|
||||
! the eigenvalue diagonal matrix
|
||||
!
|
||||
lm(1) = vxa - csnd
|
||||
lm(1) = vxa + csnd
|
||||
lm(2) = vxa
|
||||
lm(3) = vxa
|
||||
lm(4) = vxa + csnd
|
||||
lm(4) = vxa - csnd
|
||||
|
||||
! the scaling diagonal matrix
|
||||
!
|
||||
@ -3297,32 +3309,32 @@ module schemes
|
||||
|
||||
! the average of the right eigenvector matrix
|
||||
!
|
||||
rm(2,1) = vxa - csnd
|
||||
rm(2,4) = vxa + csnd
|
||||
rm(2,1) = vxa + csnd
|
||||
rm(3,1) = vya
|
||||
rm(3,4) = vya
|
||||
rm(4,1) = vza
|
||||
rm(2,4) = vxa - csnd
|
||||
rm(3,4) = vya
|
||||
rm(4,4) = vza
|
||||
|
||||
! the difference in entropy variables between states
|
||||
!
|
||||
v(1) = csnd2 * diff(log(ql(idn,i)) , log(qr(idn,i)) ) &
|
||||
- 5.0d-01 * diff(sum(ql(ivx:ivz,i)**2), sum(qr(ivx:ivz,i)**2))
|
||||
v(2) = diff(ql(ivx,i), qr(ivx,i))
|
||||
v(3) = diff(ql(ivy,i), qr(ivy,i))
|
||||
v(4) = diff(ql(ivz,i), qr(ivz,i))
|
||||
|
||||
! KEPEC flux
|
||||
!
|
||||
f(idn,i) = dnl * vxa
|
||||
f(ivx,i) = f(idn,i) * vxa + csnd2 * dnl
|
||||
f(ivy,i) = f(idn,i) * vya
|
||||
f(ivz,i) = f(idn,i) * vza
|
||||
f(imx,i) = f(idn,i) * vxa + csnd2 * dna
|
||||
f(imy,i) = f(idn,i) * vya
|
||||
f(imz,i) = f(idn,i) * vza
|
||||
|
||||
! KEPES correction
|
||||
!
|
||||
f(:,i) = f(:,i) &
|
||||
- 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
|
||||
lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
|
||||
if (vxa >= 0.0d+00) then
|
||||
do p = 1, nf
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
else
|
||||
do p = nf, 1, -1
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
@ -3357,8 +3369,8 @@ module schemes
|
||||
!
|
||||
subroutine riemann_hd_adi_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
|
||||
|
||||
use equations, only : idn, ivx, ivy, ivz, ipr
|
||||
use equations, only : adiabatic_index
|
||||
use equations, only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien
|
||||
use equations, only : nf, adiabatic_index
|
||||
|
||||
implicit none
|
||||
|
||||
@ -3370,11 +3382,11 @@ module schemes
|
||||
real(kind=8), dimension(5,5), save :: rm
|
||||
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
|
||||
|
||||
integer :: i
|
||||
real(kind=8) :: btl, btr, bta, dnl, prl, pra, vxa, vya, vza, vva
|
||||
real(kind=8) :: csn, ent, sl, sr, tl, tr
|
||||
integer :: i, p
|
||||
real(kind=8) :: bl, br, btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra
|
||||
real(kind=8) :: csn, ent
|
||||
|
||||
real(kind=8), dimension(5) :: v, lm, tm
|
||||
real(kind=8), dimension(nf) :: v, lm, tm, lv
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3382,39 +3394,58 @@ module schemes
|
||||
adi_m1 = adiabatic_index - 1.0d+00
|
||||
adi_m1x = adi_m1 / adiabatic_index
|
||||
|
||||
rm(:,:) = 0.0d+00
|
||||
rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ]
|
||||
rm(3,3) = 1.0d+00
|
||||
rm(4,4) = 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, size(ql, 2)
|
||||
|
||||
! various parameters
|
||||
!
|
||||
bl = ql(idn,i) / ql(ipr,i)
|
||||
br = qr(idn,i) / qr(ipr,i)
|
||||
|
||||
! the intermediate state averages
|
||||
!
|
||||
btl = ql(idn,i) / ql(ipr,i)
|
||||
btr = qr(idn,i) / qr(ipr,i)
|
||||
bta = lmean(btl, btr)
|
||||
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 = dnl / amean(btl, btr)
|
||||
prl = dnl / bta
|
||||
pra = dna / bta
|
||||
prl = dnl / btl
|
||||
vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i))
|
||||
ent = 1.0d+00 / adi_m1x / btl + vva
|
||||
|
||||
! the eigenvalue diagonal matrix
|
||||
! the difference in entropy variables between states
|
||||
!
|
||||
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
|
||||
+ (log(br) - log(bl)) / adi_m1 &
|
||||
- 5.0d-01 * (br * sum(qr(ivx:ivz,i)**2) &
|
||||
- bl * sum(ql(ivx:ivz,i)**2))
|
||||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||||
v(ipr) = bl - br
|
||||
|
||||
! the eigenvalues
|
||||
!
|
||||
csn = sqrt(adiabatic_index * pra / dnl)
|
||||
ent = 1.0d+00 / adi_m1x / bta + vva
|
||||
|
||||
lm(1) = vxa - csn
|
||||
! the diagonal matrix of eigenvalues
|
||||
!
|
||||
lm(1) = vxa + csn
|
||||
lm(2) = vxa
|
||||
lm(3) = vxa
|
||||
lm(4) = vxa
|
||||
lm(5) = vxa + csn
|
||||
lm(5) = vxa - csn
|
||||
|
||||
! the scaling diagonal matrix
|
||||
!
|
||||
@ -3426,43 +3457,37 @@ module schemes
|
||||
|
||||
! the average of the right eigenvector matrix
|
||||
!
|
||||
rm(2,1) = vxa - csn
|
||||
rm(2,2) = vxa
|
||||
rm(2,5) = vxa + csn
|
||||
rm(2,1) = vxa + csn
|
||||
rm(3,1) = vya
|
||||
rm(3,2) = vya
|
||||
rm(3,5) = vya
|
||||
rm(4,1) = vza
|
||||
rm(2,2) = vxa
|
||||
rm(3,2) = vya
|
||||
rm(4,2) = vza
|
||||
rm(2,5) = vxa - csn
|
||||
rm(3,5) = vya
|
||||
rm(4,5) = vza
|
||||
rm(5,:) = [ ent - vxa * csn, vva, vya, vza, ent + vxa * csn ]
|
||||
|
||||
! the difference in entropy variables between states
|
||||
!
|
||||
sl = - adi_m1 * log(ql(idn,i)) - log(btl)
|
||||
sr = - adi_m1 * log(qr(idn,i)) - log(btr)
|
||||
tl = (adiabatic_index - sl) / adi_m1 &
|
||||
- 5.0d-01 * btl * sum(ql(ivx:ivz,i)**2)
|
||||
tr = (adiabatic_index - sr) / adi_m1 &
|
||||
- 5.0d-01 * btr * sum(qr(ivx:ivz,i)**2)
|
||||
v(1) = diff(tl, tr)
|
||||
v(2) = diff(btl * ql(ivx,i), btr * qr(ivx,i))
|
||||
v(3) = diff(btl * ql(ivy,i), btr * qr(ivy,i))
|
||||
v(4) = diff(btl * ql(ivz,i), btr * qr(ivz,i))
|
||||
v(5) = - diff(btl, btr)
|
||||
rm(5,:) = [ ent + vxa * csn, vva, vya, vza, ent - vxa * csn ]
|
||||
|
||||
! KEPEC flux
|
||||
!
|
||||
f(idn,i) = dnl * vxa
|
||||
f(ivx,i) = f(idn,i) * vxa + pra
|
||||
f(ivy,i) = f(idn,i) * vya
|
||||
f(ivz,i) = f(idn,i) * vza
|
||||
f(ipr,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
|
||||
f(imx,i) = f(idn,i) * vxa + pra
|
||||
f(imy,i) = f(idn,i) * vya
|
||||
f(imz,i) = f(idn,i) * vza
|
||||
f(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa
|
||||
|
||||
! KEPES correction
|
||||
!
|
||||
f(:,i) = f(:,i) &
|
||||
- 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
|
||||
lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
|
||||
if (vxa >= 0.0d+00) then
|
||||
do p = 1, nf
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
else
|
||||
do p = nf, 1, -1
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
@ -3500,7 +3525,7 @@ module schemes
|
||||
subroutine riemann_mhd_iso_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
|
||||
|
||||
use equations, only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp
|
||||
use equations, only : csnd, csnd2, cmax
|
||||
use equations, only : nf, csnd, csnd2, ch => cglm
|
||||
|
||||
implicit none
|
||||
|
||||
@ -3511,15 +3536,15 @@ module schemes
|
||||
real(kind=8), dimension(8,8), save :: rm
|
||||
!$omp threadprivate(first, rm)
|
||||
|
||||
integer :: i
|
||||
integer :: i, p
|
||||
real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa
|
||||
real(kind=8) :: dnl, bp
|
||||
real(kind=8) :: dnl, sqd, bp
|
||||
real(kind=8) :: bb2, bp2, b1, b2, b3, bb
|
||||
real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3
|
||||
real(kind=8) :: alf2, als2, alf, als
|
||||
real(kind=8) :: f1, f2, f3, sgnb1, v2l, v2r, b2l, b2r
|
||||
|
||||
real(kind=8), dimension(8) :: v, lm, tm
|
||||
real(kind=8), dimension(nf) :: v, lm, tm, lv
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3535,24 +3560,23 @@ module schemes
|
||||
|
||||
do i = 1, size(ql, 2)
|
||||
|
||||
! the difference in entropy variables between the states
|
||||
! various parameters
|
||||
!
|
||||
v2l = sum(ql(ivx:ivz,i)**2)
|
||||
v2r = sum(qr(ivx:ivz,i)**2)
|
||||
b2l = sum(ql(ibx:ibz,i)**2)
|
||||
b2r = sum(qr(ibx:ibz,i)**2)
|
||||
|
||||
f1 = csnd2 * (log(ql(idn,i)) + 1.0d+00) - 5.0d-01 * v2l
|
||||
f2 = csnd2 * (log(qr(idn,i)) + 1.0d+00) - 5.0d-01 * v2r
|
||||
|
||||
v(idn) = diff(f1, f2)
|
||||
v(ivx) = diff(ql(ivx,i), qr(ivx,i))
|
||||
v(ivy) = diff(ql(ivy,i), qr(ivy,i))
|
||||
v(ivz) = diff(ql(ivz,i), qr(ivz,i))
|
||||
v(ibx) = diff(ql(ibx,i), qr(ibx,i))
|
||||
v(iby) = diff(ql(iby,i), qr(iby,i))
|
||||
v(ibz) = diff(ql(ibz,i), qr(ibz,i))
|
||||
v(ibp) = diff(ql(ibp,i), qr(ibp,i))
|
||||
! the difference in entropy variables between the states
|
||||
!
|
||||
v(idn) = csnd2 * (log(qr(idn,i)) - log(ql(idn,i))) - 5.0d-01 * (v2r - v2l)
|
||||
v(ivx) = qr(ivx,i) - ql(ivx,i)
|
||||
v(ivy) = qr(ivy,i) - ql(ivy,i)
|
||||
v(ivz) = qr(ivz,i) - ql(ivz,i)
|
||||
v(ibx) = qr(ibx,i) - ql(ibx,i)
|
||||
v(iby) = qr(iby,i) - ql(iby,i)
|
||||
v(ibz) = qr(ibz,i) - ql(ibz,i)
|
||||
v(ibp) = qr(ibp,i) - ql(ibp,i)
|
||||
|
||||
! the intermediate state averages
|
||||
!
|
||||
@ -3569,9 +3593,10 @@ module schemes
|
||||
|
||||
! the eigenvalues and related parameters
|
||||
!
|
||||
b1 = bxa / sqrt(dnl)
|
||||
b2 = bya / sqrt(dnl)
|
||||
b3 = bza / sqrt(dnl)
|
||||
sqd = sqrt(dnl)
|
||||
b1 = bxa / sqd
|
||||
b2 = bya / sqd
|
||||
b3 = bza / sqd
|
||||
bp2 = b2**2 + b3**2
|
||||
bp = sqrt(bp2)
|
||||
bb2 = b1**2 + bp2
|
||||
@ -3624,7 +3649,7 @@ module schemes
|
||||
!
|
||||
f1 = dnl * alf
|
||||
f2 = dnl * als * cs * sgnb1
|
||||
f3 = als * csnd * sqrt(dnl)
|
||||
f3 = als * csnd * sqd
|
||||
rm(1,1) = f1
|
||||
rm(2,1) = f1 * (vxa + cf)
|
||||
rm(3,1) = f1 * vya - f2 * x2
|
||||
@ -3656,7 +3681,7 @@ module schemes
|
||||
!
|
||||
f1 = dnl * als
|
||||
f2 = dnl * alf * cf * sgnb1
|
||||
f3 = - alf * csnd * sqrt(dnl)
|
||||
f3 = - alf * csnd * sqd
|
||||
rm(1,3) = f1
|
||||
rm(2,3) = f1 * (vxa + cs)
|
||||
rm(3,3) = f1 * vya + f2 * x2
|
||||
@ -3676,8 +3701,8 @@ module schemes
|
||||
lm(1) = vxa + cf
|
||||
lm(2) = vxa + ca
|
||||
lm(3) = vxa + cs
|
||||
lm(4) = vxa + cmax
|
||||
lm(5) = vxa - cmax
|
||||
lm(4) = vxa + ch
|
||||
lm(5) = vxa - ch
|
||||
lm(6) = vxa - cs
|
||||
lm(7) = vxa - ca
|
||||
lm(8) = vxa - cf
|
||||
@ -3688,15 +3713,23 @@ module schemes
|
||||
f(ivx,i) = f(idn,i) * vxa - bxa * bxa + pta
|
||||
f(ivy,i) = f(idn,i) * vya - bxa * bya
|
||||
f(ivz,i) = f(idn,i) * vza - bxa * bza
|
||||
f(ibx,i) = cmax * bpa
|
||||
f(ibx,i) = ch * bpa
|
||||
f(iby,i) = vxa * bya - bxa * vya
|
||||
f(ibz,i) = vxa * bza - bxa * vza
|
||||
f(ibp,i) = cmax * bxa
|
||||
f(ibp,i) = ch * bxa
|
||||
|
||||
! KEPES correction
|
||||
!
|
||||
f(:,i) = f(:,i) &
|
||||
- 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
|
||||
lv = 5.0d-01 * (abs(lm) * tm) * matmul(v, rm)
|
||||
if (vxa >= 0.0d+00) then
|
||||
do p = 1, nf
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
else
|
||||
do p = nf, 1, -1
|
||||
f(:,i) = f(:,i) - lv(p) * rm(:,p)
|
||||
end do
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
@ -3707,7 +3740,7 @@ module schemes
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RIEMANN_MHD_ADI_KEPES:
|
||||
! -------------------------------
|
||||
! --------------------------------
|
||||
!
|
||||
! Subroutine solves one dimensional adiabatic magnetohydrodynamic
|
||||
! Riemann problem using the entropy stable KEPES method.
|
||||
@ -3732,7 +3765,8 @@ module schemes
|
||||
subroutine riemann_mhd_adi_kepes(ql, qr, ul, ur, fl, fr, cl, cr, f)
|
||||
|
||||
use equations, only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
use equations, only : adiabatic_index, cmax
|
||||
use equations, only : imx, imy, imz, ien
|
||||
use equations, only : nf, adiabatic_index, ch => cglm
|
||||
|
||||
implicit none
|
||||
|
||||
@ -3744,7 +3778,7 @@ module schemes
|
||||
real(kind=8), dimension(9,9), save :: rm
|
||||
!$omp threadprivate(first, adi_m1, adi_m1x, rm)
|
||||
|
||||
integer :: i
|
||||
integer :: i, p
|
||||
real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa
|
||||
real(kind=8) :: dnl, prl, pta, vva, br, bl, bp
|
||||
real(kind=8) :: btl, bta, eka, ema, ub2, uba
|
||||
@ -3753,7 +3787,7 @@ module schemes
|
||||
real(kind=8) :: alf2, als2, alf, als, wps, wms, wpf, wmf
|
||||
real(kind=8) :: f1, f2, f3, f4, sgnb1, v2l, v2r, b2l, b2r
|
||||
|
||||
real(kind=8), dimension(9) :: v, lm, tm
|
||||
real(kind=8), dimension(nf) :: v, lm, tm
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3775,27 +3809,25 @@ module schemes
|
||||
|
||||
! the difference in entropy variables between the states
|
||||
!
|
||||
v2l = sum(ql(ivx:ivz,i)**2)
|
||||
v2r = sum(qr(ivx:ivz,i)**2)
|
||||
b2l = sum(ql(ibx:ibz,i)**2)
|
||||
b2r = sum(qr(ibx:ibz,i)**2)
|
||||
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)
|
||||
f1 = - adi_m1 * log(ql(idn,i)) - log(bl)
|
||||
f2 = - adi_m1 * log(qr(idn,i)) - log(br)
|
||||
f1 = (adiabatic_index - f1) / adi_m1 - 5.0d-01 * bl * v2l
|
||||
f2 = (adiabatic_index - f2) / adi_m1 - 5.0d-01 * br * v2r
|
||||
|
||||
v(idn) = diff(f1, f2)
|
||||
v(ivx) = diff(bl * ql(ivx,i), br * qr(ivx,i))
|
||||
v(ivy) = diff(bl * ql(ivy,i), br * qr(ivy,i))
|
||||
v(ivz) = diff(bl * ql(ivz,i), br * qr(ivz,i))
|
||||
v(ipr) = -diff(bl, br)
|
||||
v(ibx) = diff(bl * ql(ibx,i), br * qr(ibx,i))
|
||||
v(iby) = diff(bl * ql(iby,i), br * qr(iby,i))
|
||||
v(ibz) = diff(bl * ql(ibz,i), br * qr(ibz,i))
|
||||
v(ibp) = diff(bl * ql(ibp,i), br * qr(ibp,i))
|
||||
v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) &
|
||||
+ (log(br) - log(bl)) / adi_m1 &
|
||||
- 5.0d-01 * (br * v2r - bl * v2l)
|
||||
v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i)
|
||||
v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i)
|
||||
v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i)
|
||||
v(ipr) = bl - br
|
||||
v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i)
|
||||
v(iby) = br * qr(iby,i) - bl * ql(iby,i)
|
||||
v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i)
|
||||
v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i)
|
||||
|
||||
! the intermediate state averages
|
||||
!
|
||||
@ -3907,8 +3939,8 @@ module schemes
|
||||
rm(7,1) = f3 * x2
|
||||
rm(8,1) = f3 * x3
|
||||
|
||||
rm(1,9) = rm(1,1)
|
||||
rm(2,9) = rm(1,9) * (vxa - cf)
|
||||
rm(1,9) = f1
|
||||
rm(2,9) = f1 * (vxa - cf)
|
||||
rm(3,9) = f1 * vya + f2 * x2
|
||||
rm(4,9) = f1 * vza + f2 * x3
|
||||
rm(5,9) = wmf
|
||||
@ -3968,9 +4000,9 @@ module schemes
|
||||
lm(1) = vxa + cf
|
||||
lm(2) = vxa + ca
|
||||
lm(3) = vxa + cs
|
||||
lm(4) = vxa + cmax
|
||||
lm(4) = vxa + ch
|
||||
lm(5) = vxa
|
||||
lm(6) = vxa - cmax
|
||||
lm(6) = vxa - ch
|
||||
lm(7) = vxa - cs
|
||||
lm(8) = vxa - ca
|
||||
lm(9) = vxa - cf
|
||||
@ -3978,23 +4010,22 @@ module schemes
|
||||
! KEPEC flux
|
||||
!
|
||||
f(idn,i) = dnl * vxa
|
||||
f(ivx,i) = f(idn,i) * vxa - bxa * bxa + pta
|
||||
f(ivy,i) = f(idn,i) * vya - bxa * bya
|
||||
f(ivz,i) = f(idn,i) * vza - bxa * bza
|
||||
f(ibx,i) = cmax * bpa
|
||||
f(imx,i) = f(idn,i) * vxa - bxa * bxa + pta
|
||||
f(imy,i) = f(idn,i) * vya - bxa * bya
|
||||
f(imz,i) = f(idn,i) * vza - bxa * bza
|
||||
f(ibx,i) = ch * bpa
|
||||
f(iby,i) = vxa * bya - bxa * vya
|
||||
f(ibz,i) = vxa * bza - bxa * vza
|
||||
f(ibp,i) = cmax * bxa
|
||||
f(ipr,i) = f(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
|
||||
+ f(ivx,i) * vxa + f(ivy,i) * vya + f(ivz,i) * vza &
|
||||
+ f(ibx,i) * bxa + f(iby,i) * bya + f(ibz,i) * bza &
|
||||
+ f(ibp,i) * bpa - ub2 + bxa * uba &
|
||||
- cmax * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
|
||||
f(ibp,i) = ch * bxa
|
||||
f(ien,i) = f(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) &
|
||||
+ (f(imx,i) * vxa + f(imy,i) * vya + f(imz,i) * vza) &
|
||||
+ (f(ibx,i) * bxa + f(iby,i) * bya + f(ibz,i) * bza) &
|
||||
+ f(ibp,i) * bpa - ub2 + bxa * uba &
|
||||
- ch * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i))
|
||||
|
||||
! KEPES correction
|
||||
!
|
||||
f(:,i) = f(:,i) &
|
||||
- 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(transpose(rm), v))
|
||||
f(:,i) = f(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm))
|
||||
|
||||
end do
|
||||
|
||||
@ -4107,19 +4138,25 @@ module schemes
|
||||
|
||||
real(kind=8), intent(in) :: l, r
|
||||
|
||||
real(kind=8) :: u, rr
|
||||
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
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
rr = r * r
|
||||
u = (l * (l - 2.0d+00 * r) + rr) / (l * (l + 2.0d+00 * r) + rr)
|
||||
d = r - l
|
||||
s = r + l
|
||||
u = (d / s)**2
|
||||
if (u < 1.0d-04) then
|
||||
lmean = (l + r) / (c1 + u * (c2 + u * (c3 + u * c4)))
|
||||
lmean = s / (c1 + u * (c2 + u * (c3 + u * c4)))
|
||||
else
|
||||
lmean = (r - l) / log(r / l)
|
||||
if (d >= 0.0d+00) then
|
||||
s = log(r / l)
|
||||
else
|
||||
s = log(l / r)
|
||||
end if
|
||||
lmean = abs(d) / s
|
||||
end if
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -4152,33 +4189,6 @@ module schemes
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end function amean
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! function DIFF:
|
||||
! -------------
|
||||
!
|
||||
! Function calculate the difference.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! l, r - the values to calculate the difference of;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
real(kind=8) function diff(l, r)
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), intent(in) :: l, r
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
diff = r - l
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end function diff
|
||||
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
@ -121,6 +121,9 @@ module sources
|
||||
case("heglm", "HEGLM")
|
||||
glm_type = 2
|
||||
glm_name = "HEGLM"
|
||||
case("kepes", "KEPES")
|
||||
glm_type = 3
|
||||
glm_name = "KEPES"
|
||||
case default
|
||||
glm_type = 0
|
||||
glm_name = "none"
|
||||
@ -556,7 +559,7 @@ module sources
|
||||
dh(2) = ady(pdata%meta%level)
|
||||
dh(3) = adz(pdata%meta%level)
|
||||
|
||||
! add the EGLM-MHD source terms
|
||||
! add the EGLM-MHD and KEPES source terms
|
||||
!
|
||||
if (glm_type > 0) then
|
||||
|
||||
@ -571,23 +574,7 @@ module sources
|
||||
du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
|
||||
du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
|
||||
|
||||
! update the energy equation
|
||||
!
|
||||
if (ien > 0 .and. ibp > 0) then
|
||||
|
||||
! calculate the gradient of divergence potential
|
||||
!
|
||||
call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
|
||||
|
||||
! add the divergence potential source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = - B.(∇ψ)
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
|
||||
|
||||
end if ! ien > 0
|
||||
|
||||
! add the HEGLM-MHD source terms
|
||||
! add the HEGLM-MHD and KEPES source terms
|
||||
!
|
||||
if (glm_type > 1) then
|
||||
|
||||
@ -616,7 +603,41 @@ module sources
|
||||
|
||||
end if ! glm_type > 1
|
||||
|
||||
end if ! glmtype > 0
|
||||
if (ibp > 0 .and. (ien > 0 .or. glm_type == 3)) then
|
||||
|
||||
! calculate the gradient of divergence correcting field
|
||||
!
|
||||
call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
|
||||
|
||||
db(:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
|
||||
|
||||
! update the divergence correcting field
|
||||
! d/dt ψ+ ∇.F = - (v.∇)ψ
|
||||
!
|
||||
du(ibp,:,:,:) = du(ibp,:,:,:) - db(:,:,:)
|
||||
|
||||
! update the energy equation
|
||||
!
|
||||
if (ien > 0) then
|
||||
if (glm_type == 3) then
|
||||
|
||||
! add the divergence potential source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = - ψ(v.∇)ψ
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) - pdata%q(ibp,:,:,:) * db(:,:,:)
|
||||
|
||||
else
|
||||
|
||||
! add the divergence potential source term to the energy equation, i.e.
|
||||
! d/dt E + ∇.F = - B.(∇ψ)
|
||||
!
|
||||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||||
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
|
||||
|
||||
end if
|
||||
end if ! ien > 0
|
||||
end if ! glm_type == 3
|
||||
end if ! glm_type > 0
|
||||
|
||||
! if anomalous resistivity is enabled
|
||||
!
|
||||
|
||||
@ -38,12 +38,14 @@ module statistics
|
||||
! munit - a file handler for the mesh statistics file;
|
||||
! cunit - a file handler for the conserved variable integrals file;
|
||||
! sunit - a file handler for the variable statistics file;
|
||||
! funit - a file handler for the forcing statistics file;
|
||||
! eunit - a file handler for the integration errors file;
|
||||
!
|
||||
integer(kind=4), save :: munit = 10
|
||||
integer(kind=4), save :: cunit = 11
|
||||
integer(kind=4), save :: sunit = 12
|
||||
integer(kind=4), save :: eunit = 13
|
||||
integer(kind=4), save :: funit = 13
|
||||
integer(kind=4), save :: eunit = 14
|
||||
|
||||
! iintd - the number of steps between subsequent intervals;
|
||||
!
|
||||
@ -96,6 +98,7 @@ module statistics
|
||||
use coordinates, only : toplev, ncells, bcells, nghosts, domain_base_dims
|
||||
use equations , only : pvars, nf
|
||||
use evolution , only : error_control
|
||||
use forcing , only : forcing_enabled
|
||||
use mpitools , only : master
|
||||
#ifdef MPI
|
||||
use mpitools , only : nprocs
|
||||
@ -234,10 +237,9 @@ module statistics
|
||||
#endif /* __INTEL_COMPILER */
|
||||
end if
|
||||
|
||||
write(cunit,"('#',a8,13(1x,a18))") 'step', 'time', 'dt' &
|
||||
write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt' &
|
||||
, 'mass', 'momx', 'momy', 'momz' &
|
||||
, 'ener', 'ekin', 'emag', 'eint' &
|
||||
, 'einj', 'erat', 'arms'
|
||||
, 'ener', 'ekin', 'emag', 'eint', 'entr'
|
||||
write(cunit,"('#')")
|
||||
|
||||
! prepare the variable statistics file
|
||||
@ -279,6 +281,40 @@ module statistics
|
||||
, 'mean(Malf)', 'min(Malf)', 'max(Malf)'
|
||||
write(sunit,"('#')")
|
||||
|
||||
! prepare the forcing statistics file
|
||||
!
|
||||
if (forcing_enabled) then
|
||||
if (append) then
|
||||
write(fname, "('forcing-statistics.dat')")
|
||||
inquire(file=fname, exist=flag)
|
||||
flag = flag .and. nrun > 1
|
||||
else
|
||||
write(fname, "('forcing-statistics_',i2.2,'.dat')") nrun
|
||||
flag = .false.
|
||||
end if
|
||||
|
||||
if (flag) then
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=funit, file=fname, form='formatted', status='old', &
|
||||
position='append', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=funit, file=fname, form='formatted', status='old', &
|
||||
position='append')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
write(cunit,"('#')")
|
||||
else
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=funit, file=fname, form='formatted', &
|
||||
status='replace', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=funit, file=fname, form='formatted', status='replace')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
end if
|
||||
|
||||
write(funit,"('#',a18,4(1x,a18))") 'time', 'einj', 'erat', 'arms'
|
||||
write(funit,"('#')")
|
||||
end if
|
||||
|
||||
! prepare the integration errors file
|
||||
!
|
||||
if (error_control) then
|
||||
@ -344,6 +380,7 @@ module statistics
|
||||
subroutine finalize_statistics(status)
|
||||
|
||||
use evolution, only : error_control
|
||||
use forcing , only : forcing_enabled
|
||||
use mpitools , only : master
|
||||
|
||||
implicit none
|
||||
@ -358,6 +395,7 @@ module statistics
|
||||
close(munit)
|
||||
close(cunit)
|
||||
close(sunit)
|
||||
if (forcing_enabled) close(funit)
|
||||
if (error_control) close(eunit)
|
||||
end if
|
||||
|
||||
@ -386,11 +424,11 @@ module statistics
|
||||
use coordinates, only : toplev
|
||||
use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp
|
||||
use equations , only : ien, imx, imy, imz
|
||||
use equations , only : magnetized, adiabatic_index, csnd
|
||||
use equations , only : magnetized, adiabatic_index, csnd, csnd2
|
||||
use equations , only : errors
|
||||
use evolution , only : error_control
|
||||
use evolution , only : step, time, dt, dth, dte
|
||||
use forcing , only : einj, rinj, arms
|
||||
use forcing , only : forcing_enabled, einj, rinj, arms
|
||||
use helpers , only : print_message, flush_and_sync
|
||||
use mesh , only : changed
|
||||
use mpitools , only : master
|
||||
@ -566,7 +604,7 @@ module statistics
|
||||
inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! sum up total energy
|
||||
! sum up total energy and entropy
|
||||
!
|
||||
if (ien > 0) then
|
||||
#if NDIMS == 3
|
||||
@ -574,7 +612,26 @@ module statistics
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
#if NDIMS == 3
|
||||
tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne,nb:ne) &
|
||||
/ pdata%q(idn,nb:ne,nb:ne,nb:ne)**adiabatic_index) &
|
||||
* pdata%q(idn,nb:ne,nb:ne,nb:ne) / (adiabatic_index - 1.0d+00)
|
||||
#else /* NDIMS == 3 */
|
||||
tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne, : ) &
|
||||
/ pdata%q(idn,nb:ne,nb:ne, : )**adiabatic_index) &
|
||||
* pdata%q(idn,nb:ne,nb:ne, : ) / (adiabatic_index - 1.0d+00)
|
||||
#endif /* NDIMS == 3 */
|
||||
else
|
||||
#if NDIMS == 3
|
||||
tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
|
||||
* pdata%q(idn,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne, : )) &
|
||||
* pdata%q(idn,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
inarr(9) = inarr(9) + sum(tmp(:,:,:)) * dvol
|
||||
|
||||
! sum up kinetic energy
|
||||
!
|
||||
@ -602,12 +659,21 @@ module statistics
|
||||
+ pdata%u(iby,nb:ne,nb:ne, : )**2 &
|
||||
+ pdata%u(ibz,nb:ne,nb:ne, : )**2) * dvolh
|
||||
#endif /* NDIMS == 3 */
|
||||
if (ien <= 0) then
|
||||
#if NDIMS == 3
|
||||
inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne,nb:ne)**2) * dvolh
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne, : )**2) * dvolh
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
end if
|
||||
|
||||
! sum up the injected energy and injection rate
|
||||
!
|
||||
inarr( 9) = einj
|
||||
inarr(10) = rinj
|
||||
if (forcing_enabled) then
|
||||
inarr(11) = einj
|
||||
inarr(12) = rinj
|
||||
end if
|
||||
|
||||
! get average, minimum and maximum values of density
|
||||
!
|
||||
@ -720,9 +786,13 @@ module statistics
|
||||
call reduce_maximum(mxarr(1:narr))
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate the internal energy
|
||||
! calculate the internal energy, or update the entropy
|
||||
!
|
||||
if (ien > 0) inarr(8) = inarr(5) - inarr(6) - inarr(7)
|
||||
if (ien > 0) then
|
||||
inarr(8) = inarr(5) - inarr(6) - inarr(7)
|
||||
else
|
||||
inarr(9) = inarr(9) + inarr(6) + inarr(7)
|
||||
end if
|
||||
|
||||
! normalize the averages by the volume of domain
|
||||
!
|
||||
@ -731,7 +801,7 @@ module statistics
|
||||
! write down the integrals and statistics to appropriate files
|
||||
!
|
||||
if (master) then
|
||||
write(cunit,"(i9,13(1x,1es18.8e3))") step, time, dt, inarr(1:10), arms
|
||||
write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9)
|
||||
write(sunit,"(i9,23(1x,1es18.8e3))") step, time &
|
||||
, avarr(1), mnarr(1), mxarr(1) &
|
||||
, avarr(2), mnarr(2), mxarr(2) &
|
||||
@ -744,6 +814,10 @@ module statistics
|
||||
call flush_and_sync(cunit)
|
||||
call flush_and_sync(sunit)
|
||||
|
||||
if (forcing_enabled) then
|
||||
write(funit,"(4(1x,1es18.8e3))") time, inarr(11:12), arms
|
||||
call flush_and_sync(funit)
|
||||
end if
|
||||
if (error_control) then
|
||||
write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:)
|
||||
call flush_and_sync(eunit)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user