gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EQUATIONS: Take into account ½ψ² in total energy for adiabatic MHD.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-01-06 16:32:58 -03:00
parent 4f3a05f9ba
commit 25ea3c3daf
1 changed files with 11 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
sources/equations.F90
View File

@ -2960,24 +2960,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 +2982,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 +3026,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 +3043,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