Merge branch 'master' into reconnection
This commit is contained in:
commit
6f5c27293c
@ -63,17 +63,21 @@ module equations
|
||||
|
||||
! pointers to the conversion procedures
|
||||
!
|
||||
procedure(prim2cons_hd_iso), pointer, save :: prim2cons => null()
|
||||
procedure(cons2prim_hd_iso), pointer, save :: cons2prim => null()
|
||||
procedure(prim2cons_hd_iso) , pointer, save :: prim2cons => null()
|
||||
procedure(cons2prim_hd_iso) , pointer, save :: cons2prim => null()
|
||||
|
||||
! pointer to the flux procedure
|
||||
!
|
||||
procedure(fluxspeed_hd_iso), pointer, save :: fluxspeed => null()
|
||||
procedure(fluxspeed_hd_iso) , pointer, save :: fluxspeed => null()
|
||||
|
||||
! pointer to the maxspeed procedure
|
||||
!
|
||||
procedure(maxspeed_hd_iso) , pointer, save :: maxspeed => null()
|
||||
|
||||
! pointer to the Roe eigensystem procedure
|
||||
!
|
||||
procedure(esystem_roe_hd_iso), pointer, save :: eigensystem_roe => null()
|
||||
|
||||
|
||||
! the system of equations and the equation of state
|
||||
!
|
||||
@ -97,6 +101,10 @@ module equations
|
||||
!
|
||||
character(len=4), dimension(:), allocatable, save :: pvars, cvars
|
||||
|
||||
! eigenvectors
|
||||
!
|
||||
real(kind=8), dimension(:,:,:), allocatable, save :: evroe
|
||||
|
||||
! adiabatic heat ratio
|
||||
!
|
||||
real(kind=8) , save :: gamma = 5.0d+00 / 3.0d+00
|
||||
@ -123,6 +131,7 @@ module equations
|
||||
public :: prim2cons, cons2prim
|
||||
public :: fluxspeed
|
||||
public :: maxspeed, reset_maxspeed, get_maxspeed
|
||||
public :: eigensystem_roe
|
||||
public :: update_primitive_variables
|
||||
public :: gamma
|
||||
public :: csnd
|
||||
@ -241,6 +250,7 @@ module equations
|
||||
cons2prim => cons2prim_hd_iso
|
||||
fluxspeed => fluxspeed_hd_iso
|
||||
maxspeed => maxspeed_hd_iso
|
||||
eigensystem_roe => esystem_roe_hd_iso
|
||||
|
||||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||||
|
||||
@ -263,6 +273,7 @@ module equations
|
||||
cons2prim => cons2prim_hd_adi
|
||||
fluxspeed => fluxspeed_hd_adi
|
||||
maxspeed => maxspeed_hd_adi
|
||||
eigensystem_roe => esystem_roe_hd_adi
|
||||
|
||||
! warn about the unimplemented equation of state
|
||||
!
|
||||
@ -341,6 +352,7 @@ module equations
|
||||
cons2prim => cons2prim_mhd_iso
|
||||
fluxspeed => fluxspeed_mhd_iso
|
||||
maxspeed => maxspeed_mhd_iso
|
||||
eigensystem_roe => esystem_roe_mhd_iso
|
||||
|
||||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||||
|
||||
@ -363,6 +375,7 @@ module equations
|
||||
cons2prim => cons2prim_mhd_adi
|
||||
fluxspeed => fluxspeed_mhd_adi
|
||||
maxspeed => maxspeed_mhd_adi
|
||||
eigensystem_roe => esystem_roe_mhd_adi
|
||||
|
||||
case default
|
||||
|
||||
@ -435,6 +448,10 @@ module equations
|
||||
!
|
||||
csnd2 = csnd * csnd
|
||||
|
||||
! allocate space for Roe eigenvectors
|
||||
!
|
||||
allocate(evroe(2,nv,nv))
|
||||
|
||||
! print information about the equation module
|
||||
!
|
||||
if (verbose) then
|
||||
@ -490,12 +507,17 @@ module equations
|
||||
if (allocated(pvars)) deallocate(pvars)
|
||||
if (allocated(cvars)) deallocate(cvars)
|
||||
|
||||
! deallocate Roe eigenvectors
|
||||
!
|
||||
if (allocated(evroe)) deallocate(evroe)
|
||||
|
||||
! release the procedure pointers
|
||||
!
|
||||
nullify(prim2cons)
|
||||
nullify(cons2prim)
|
||||
nullify(fluxspeed)
|
||||
nullify(maxspeed )
|
||||
nullify(maxspeed)
|
||||
nullify(eigensystem_roe)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
@ -911,6 +933,120 @@ module equations
|
||||
!
|
||||
end function maxspeed_hd_iso
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine ESYSTEM_ROE_HD_ISO:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine computes eigenvalues and eigenvectors for a given set of
|
||||
! equations and input variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! x - ratio of the perpendicular magnetic field component difference
|
||||
! y - ratio of the density
|
||||
! q - the intermediate Roe state vector;
|
||||
! c - the vector of eigenvalues;
|
||||
! r - the matrix of right eigenvectors;
|
||||
! l - the matrix of left eigenvectors;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Roe, P. L.
|
||||
! "Approximate Riemann Solvers, Parameter Vectors, and Difference
|
||||
! Schemes",
|
||||
! Journal of Computational Physics, 1981, 43, pp. 357-372
|
||||
! [2] Stone, J. M. & Gardiner, T. A.,
|
||||
! "ATHENA: A New Code for Astrophysical MHD",
|
||||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine esystem_roe_hd_iso(x, y, q, c, r, l)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8) , intent(in) :: x, y
|
||||
real(kind=8), dimension(nv) , intent(in) :: q
|
||||
real(kind=8), dimension(nv) , intent(inout) :: c
|
||||
real(kind=8), dimension(nv,nv), intent(inout) :: l, r
|
||||
|
||||
! local variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: ch
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! prepare the internal arrays at the first run
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! prepare constants
|
||||
!
|
||||
ch = 0.5d+00 / csnd
|
||||
|
||||
! reset all elements
|
||||
!
|
||||
evroe(:, : ,:) = 0.0d+00
|
||||
|
||||
! initiate the matrix of left eigenvectors
|
||||
!
|
||||
evroe(1,ivx,1) = - ch
|
||||
evroe(1,ivy,2) = 1.0d+00
|
||||
evroe(1,ivz,3) = 1.0d+00
|
||||
evroe(1,ivx,4) = ch
|
||||
|
||||
! initiate the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,1,idn) = 1.0d+00
|
||||
evroe(2,2,ivy) = 1.0d+00
|
||||
evroe(2,3,ivz) = 1.0d+00
|
||||
evroe(2,4,idn) = 1.0d+00
|
||||
|
||||
! unset the first execution flag
|
||||
!
|
||||
first = .false.
|
||||
|
||||
end if ! first execution
|
||||
|
||||
! prepare eigenvalues
|
||||
!
|
||||
c(1) = q(ivx) - csnd
|
||||
c(2) = q(ivx)
|
||||
c(3) = q(ivx)
|
||||
c(4) = q(ivx) + csnd
|
||||
|
||||
! update the varying elements of the matrix of left eigenvectors
|
||||
!
|
||||
evroe(1,idn,1) = ch * c(4)
|
||||
evroe(1,idn,2) = - q(ivy)
|
||||
evroe(1,idn,3) = - q(ivz)
|
||||
evroe(1,idn,4) = - ch * c(1)
|
||||
|
||||
! update the varying elements of the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,1,ivx) = c(1)
|
||||
evroe(2,1,ivy) = q(ivy)
|
||||
evroe(2,1,ivz) = q(ivz)
|
||||
|
||||
evroe(2,4,ivx) = c(4)
|
||||
evroe(2,4,ivy) = q(ivy)
|
||||
evroe(2,4,ivz) = q(ivz)
|
||||
|
||||
! copy matrices of eigenvectors
|
||||
!
|
||||
l(1:nv,1:nv) = evroe(1,1:nv,1:nv)
|
||||
r(1:nv,1:nv) = evroe(2,1:nv,1:nv)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine esystem_roe_hd_iso
|
||||
!
|
||||
!*******************************************************************************
|
||||
!
|
||||
! ADIABATIC HYDRODYNAMIC EQUATIONS
|
||||
@ -1205,6 +1341,158 @@ module equations
|
||||
!
|
||||
end function maxspeed_hd_adi
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine ESYSTEM_ROE_HD_ADI:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine computes eigenvalues and eigenvectors for a given set of
|
||||
! equations and input variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! x - ratio of the perpendicular magnetic field component difference
|
||||
! y - ratio of the density
|
||||
! q - the intermediate Roe state vector;
|
||||
! c - the vector of eigenvalues;
|
||||
! r - the matrix of right eigenvectors;
|
||||
! l - the matrix of left eigenvectors;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Roe, P. L.
|
||||
! "Approximate Riemann Solvers, Parameter Vectors, and Difference
|
||||
! Schemes",
|
||||
! Journal of Computational Physics, 1981, 43, pp. 357-372
|
||||
! [2] Stone, J. M. & Gardiner, T. A.,
|
||||
! "ATHENA: A New Code for Astrophysical MHD",
|
||||
! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine esystem_roe_hd_adi(x, y, q, c, r, l)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8) , intent(in) :: x, y
|
||||
real(kind=8), dimension(nv) , intent(in) :: q
|
||||
real(kind=8), dimension(nv) , intent(inout) :: c
|
||||
real(kind=8), dimension(nv,nv), intent(inout) :: l, r
|
||||
|
||||
! local variables
|
||||
!
|
||||
logical, save :: first = .true.
|
||||
real(kind=8) :: vv, vh, c2, na, cc, vc, ng, nd, nw, nh, nc
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! prepare the internal arrays at the first run
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! reset all elements
|
||||
!
|
||||
evroe(:, : ,:) = 0.0d+00
|
||||
|
||||
! initiate the matrix of left eigenvectors
|
||||
!
|
||||
evroe(1,ivy,2) = 1.0d+00
|
||||
evroe(1,ivz,3) = 1.0d+00
|
||||
|
||||
! initiate the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,1,idn) = 1.0d+00
|
||||
evroe(2,2,ivy) = 1.0d+00
|
||||
evroe(2,3,ivz) = 1.0d+00
|
||||
evroe(2,4,idn) = 1.0d+00
|
||||
evroe(2,5,idn) = 1.0d+00
|
||||
|
||||
! unset the first execution flag
|
||||
!
|
||||
first = .false.
|
||||
|
||||
end if ! first execution
|
||||
|
||||
! calculate characteristic speeds and useful variables
|
||||
!
|
||||
vv = sum(q(ivx:ivz)**2)
|
||||
vh = 0.5d+00 * vv
|
||||
c2 = gammam1 * (q(ien) - vh)
|
||||
na = 0.5d+00 / c2
|
||||
cc = sqrt(c2)
|
||||
vc = q(ivx) * cc
|
||||
ng = na * gammam1
|
||||
nd = 2.0d+00 * ng
|
||||
nw = na * vc
|
||||
nh = na * gammam1 * vh
|
||||
nc = na * cc
|
||||
|
||||
! prepare eigenvalues
|
||||
!
|
||||
c(1) = q(ivx) - cc
|
||||
c(2) = q(ivx)
|
||||
c(3) = q(ivx)
|
||||
c(4) = q(ivx)
|
||||
c(5) = q(ivx) + cc
|
||||
|
||||
! update the varying elements of the matrix of left eigenvectors
|
||||
!
|
||||
evroe(1,idn,1) = nh + nw
|
||||
evroe(1,ivx,1) = - ng * q(ivx) - nc
|
||||
evroe(1,ivy,1) = - ng * q(ivy)
|
||||
evroe(1,ivz,1) = - ng * q(ivz)
|
||||
evroe(1,ien,1) = ng
|
||||
|
||||
evroe(1,idn,2) = - q(ivy)
|
||||
|
||||
evroe(1,idn,3) = - q(ivz)
|
||||
|
||||
evroe(1,idn,4) = 1.0d+00 - ng * vv
|
||||
evroe(1,ivx,4) = nd * q(ivx)
|
||||
evroe(1,ivy,4) = nd * q(ivy)
|
||||
evroe(1,ivz,4) = nd * q(ivz)
|
||||
evroe(1,ien,4) = - nd
|
||||
|
||||
evroe(1,idn,5) = nh - nw
|
||||
evroe(1,ivx,5) = - ng * q(ivx) + nc
|
||||
evroe(1,ivy,5) = - ng * q(ivy)
|
||||
evroe(1,ivz,5) = - ng * q(ivz)
|
||||
evroe(1,ien,5) = ng
|
||||
|
||||
! update the varying elements of the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,1,ivx) = q(ivx) - cc
|
||||
evroe(2,1,ivy) = q(ivy)
|
||||
evroe(2,1,ivz) = q(ivz)
|
||||
evroe(2,1,ien) = q(ien) - vc
|
||||
|
||||
evroe(2,2,ien) = q(ivy)
|
||||
|
||||
evroe(2,3,ien) = q(ivz)
|
||||
|
||||
evroe(2,4,ivx) = q(ivx)
|
||||
evroe(2,4,ivy) = q(ivy)
|
||||
evroe(2,4,ivz) = q(ivz)
|
||||
evroe(2,4,ien) = vh
|
||||
|
||||
evroe(2,5,ivx) = q(ivx) + cc
|
||||
evroe(2,5,ivy) = q(ivy)
|
||||
evroe(2,5,ivz) = q(ivz)
|
||||
evroe(2,5,ien) = q(ien) + vc
|
||||
|
||||
! copy matrices of eigenvectors
|
||||
!
|
||||
l(1:nv,1:nv) = evroe(1,1:nv,1:nv)
|
||||
r(1:nv,1:nv) = evroe(2,1:nv,1:nv)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine esystem_roe_hd_adi
|
||||
!
|
||||
!*******************************************************************************
|
||||
!
|
||||
! ISOTHERMAL MAGNETOHYDRODYNAMIC EQUATIONS
|
||||
@ -1520,6 +1808,276 @@ module equations
|
||||
!
|
||||
end function maxspeed_mhd_iso
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine ESYSTEM_ROE_MHD_ISO:
|
||||
! ------------------------------
|
||||
!
|
||||
! Subroutine computes eigenvalues and eigenvectors for a given set of
|
||||
! equations and input variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! x - ratio of the perpendicular magnetic field component difference
|
||||
! y - ratio of the density
|
||||
! q - the intermediate Roe state vector;
|
||||
! c - the vector of eigenvalues;
|
||||
! r - the matrix of right eigenvectors;
|
||||
! l - the matrix of left eigenvectors;
|
||||
!
|
||||
! 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] Balsara, D. S.
|
||||
! "Linearized Formulation of the Riemann Problem for Adiabatic and
|
||||
! Isothermal Magnetohydrodynamics",
|
||||
! The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine esystem_roe_mhd_iso(x, y, q, c, r, l)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8) , intent(in) :: x, y
|
||||
real(kind=8), dimension(nv) , intent(in) :: q
|
||||
real(kind=8), dimension(nv) , intent(inout) :: c
|
||||
real(kind=8), dimension(nv,nv), intent(inout) :: l, r
|
||||
|
||||
! saved variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
|
||||
! local variables
|
||||
!
|
||||
real(kind=8) :: di, btsq, bt_starsq, casq, twid_csq
|
||||
real(kind=8) :: ct2, tsum, tdif, cf2_cs2, cfsq, cf, cssq, cs, ca
|
||||
real(kind=8) :: bt, bt_star, bet2, bet3, bet2_star, bet3_star, bet_starsq
|
||||
real(kind=8) :: alpha_f, alpha_s
|
||||
real(kind=8) :: sqrtd, s, twid_c, qf, qs, af_prime, as_prime
|
||||
real(kind=8) :: norm, cff, css, af, as, afpb, aspb, q2_star, q3_star, vqstr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! prepare the internal arrays at the first run
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! reset all elements
|
||||
!
|
||||
evroe(:, : ,:) = 0.0d+00
|
||||
|
||||
! unset the first execution flag
|
||||
!
|
||||
first = .false.
|
||||
|
||||
end if ! first execution
|
||||
|
||||
! prepare coefficients for eigenvalues
|
||||
!
|
||||
di = 1.0d+00 / q(idn)
|
||||
casq = q(ibx) * q(ibx) * di
|
||||
ca = sqrt(casq)
|
||||
btsq = q(iby) * q(iby) + q(ibz) * q(ibz)
|
||||
bt_starsq = btsq * y
|
||||
twid_csq = csnd2 + x
|
||||
ct2 = bt_starsq * di
|
||||
tsum = casq + ct2 + twid_csq
|
||||
tdif = casq + ct2 - twid_csq
|
||||
cf2_cs2 = sqrt(tdif * tdif + 4.0d+00 * twid_csq * ct2)
|
||||
cfsq = 0.5d+00 * (tsum + cf2_cs2)
|
||||
cf = sqrt(cfsq)
|
||||
cssq = twid_csq * casq / cfsq
|
||||
cs = sqrt(cssq)
|
||||
|
||||
! prepare eigenvalues
|
||||
!
|
||||
c(1) = q(ivx) - cf
|
||||
c(2) = q(ivx) - ca
|
||||
c(3) = q(ivx) - cs
|
||||
c(4) = q(ivx)
|
||||
c(5) = q(ivx) + cs
|
||||
c(6) = q(ivx) + ca
|
||||
c(7) = q(ivx) + cf
|
||||
c(8) = c(7)
|
||||
|
||||
! calculate the eigenvectors only if the waves propagate in both direction
|
||||
!
|
||||
if (c(1) >= 0.0d+00) return
|
||||
if (c(7) <= 0.0d+00) return
|
||||
|
||||
! prepare remaining coefficients for eigenvectors
|
||||
!
|
||||
bt = sqrt(btsq)
|
||||
bt_star = sqrt(bt_starsq)
|
||||
if (bt == 0.0d+00) then
|
||||
bet2 = 1.0d+00
|
||||
bet3 = 0.0d+00
|
||||
else
|
||||
bet2 = q(iby) / bt
|
||||
bet3 = q(ibz) / bt
|
||||
end if
|
||||
bet2_star = bet2 / sqrt(y)
|
||||
bet3_star = bet3 / sqrt(y)
|
||||
bet_starsq = bet2_star * bet2_star + bet3_star * bet3_star
|
||||
|
||||
if ((cfsq - cssq) == 0.0d+00) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
else if ((twid_csq - cssq) <= 0.0d+00) then
|
||||
alpha_f = 0.0d+00
|
||||
alpha_s = 1.0d+00
|
||||
else if ((cfsq - twid_csq) <= 0.0d+00) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
else
|
||||
alpha_f = sqrt((twid_csq - cssq) / (cfsq - cssq))
|
||||
alpha_s = sqrt((cfsq - twid_csq) / (cfsq - cssq))
|
||||
end if
|
||||
|
||||
sqrtd = sqrt(q(idn))
|
||||
s = sign(1.0d+00, q(ibx))
|
||||
twid_c = sqrt(twid_csq)
|
||||
qf = cf * alpha_f * s
|
||||
qs = cs * alpha_s * s
|
||||
af_prime = twid_c * alpha_f / sqrtd
|
||||
as_prime = twid_c * alpha_s / sqrtd
|
||||
|
||||
! update the varying elements of the matrix of right eigenvectors
|
||||
!
|
||||
! left-going fast wave
|
||||
!
|
||||
evroe(2,1,idn) = alpha_f
|
||||
evroe(2,1,ivx) = alpha_f * c(1)
|
||||
evroe(2,1,ivy) = alpha_f * q(ivy) + qs * bet2_star
|
||||
evroe(2,1,ivz) = alpha_f * q(ivz) + qs * bet3_star
|
||||
evroe(2,1,iby) = as_prime * bet2_star
|
||||
evroe(2,1,ibz) = as_prime * bet3_star
|
||||
|
||||
! left-going Alfvèn wave
|
||||
!
|
||||
evroe(2,2,ivy) = - bet3
|
||||
evroe(2,2,ivz) = bet2
|
||||
evroe(2,2,iby) = - bet3 * s / sqrtd
|
||||
evroe(2,2,ibz) = bet2 * s / sqrtd
|
||||
|
||||
! left-going slow wave
|
||||
!
|
||||
evroe(2,3,idn) = alpha_s
|
||||
evroe(2,3,ivx) = alpha_s * c(3)
|
||||
evroe(2,3,ivy) = alpha_s * q(ivy) - qf * bet2_star
|
||||
evroe(2,3,ivz) = alpha_s * q(ivz) - qf * bet3_star
|
||||
evroe(2,3,iby) = - af_prime * bet2_star
|
||||
evroe(2,3,ibz) = - af_prime * bet3_star
|
||||
|
||||
! right-going slow wave
|
||||
!
|
||||
evroe(2,5,idn) = alpha_s
|
||||
evroe(2,5,ivx) = alpha_s * c(5)
|
||||
evroe(2,5,ivy) = alpha_s * q(ivy) + qf * bet2_star
|
||||
evroe(2,5,ivz) = alpha_s * q(ivz) + qf * bet3_star
|
||||
evroe(2,5,iby) = evroe(2,3,iby)
|
||||
evroe(2,5,ibz) = evroe(2,3,ibz)
|
||||
|
||||
! right-going Alfvèn wave
|
||||
!
|
||||
evroe(2,6,ivy) = bet3
|
||||
evroe(2,6,ivz) = - bet2
|
||||
evroe(2,6,iby) = evroe(2,2,iby)
|
||||
evroe(2,6,ibz) = evroe(2,2,ibz)
|
||||
|
||||
! right-going fast wave
|
||||
!
|
||||
evroe(2,7,idn) = alpha_f
|
||||
evroe(2,7,ivx) = alpha_f * c(7)
|
||||
evroe(2,7,ivy) = alpha_f * q(ivy) - qs * bet2_star
|
||||
evroe(2,7,ivz) = alpha_f * q(ivz) - qs * bet3_star
|
||||
evroe(2,7,iby) = evroe(2,1,iby)
|
||||
evroe(2,7,ibz) = evroe(2,1,ibz)
|
||||
|
||||
! update the varying elements of the matrix of left eigenvectors
|
||||
!
|
||||
norm = 0.5d+00 / twid_csq
|
||||
cff = norm * alpha_f * cf
|
||||
css = norm * alpha_s * cs
|
||||
qf = qf * norm
|
||||
qs = qs * norm
|
||||
af = norm * af_prime * q(idn)
|
||||
as = norm * as_prime * q(idn)
|
||||
afpb = norm * af_prime * bt_star
|
||||
aspb = norm * as_prime * bt_star
|
||||
|
||||
q2_star = bet2_star / bet_starsq
|
||||
q3_star = bet3_star / bet_starsq
|
||||
vqstr = q(ivy) * q2_star + q(ivz) * q3_star
|
||||
|
||||
! left-going fast wave
|
||||
!
|
||||
evroe(1,idn,1) = cff * c(7) - qs * vqstr - aspb
|
||||
evroe(1,ivx,1) = - cff
|
||||
evroe(1,ivy,1) = qs * q2_star
|
||||
evroe(1,ivz,1) = qs * q3_star
|
||||
evroe(1,iby,1) = as * q2_star
|
||||
evroe(1,ibz,1) = as * q3_star
|
||||
|
||||
! left-going Alfvèn wave
|
||||
!
|
||||
evroe(1,idn,2) = 0.5d+00 * (q(ivy) * bet3 - q(ivz) * bet2)
|
||||
evroe(1,ivy,2) = - 0.5d+00 * bet3
|
||||
evroe(1,ivz,2) = 0.5d+00 * bet2
|
||||
evroe(1,iby,2) = - 0.5d+00 * sqrtd * bet3 * s
|
||||
evroe(1,ibz,2) = 0.5d+00 * sqrtd * bet2 * s
|
||||
|
||||
! left-going slow wave
|
||||
!
|
||||
evroe(1,idn,3) = css * c(5) + qf * vqstr + afpb
|
||||
evroe(1,ivx,3) = - css
|
||||
evroe(1,ivy,3) = - qf * q2_star
|
||||
evroe(1,ivz,3) = - qf * q3_star
|
||||
evroe(1,iby,3) = - af * q2_star
|
||||
evroe(1,ibz,3) = - af * q3_star
|
||||
|
||||
! right-going slow wave
|
||||
!
|
||||
evroe(1,idn,5) = - css * c(3) - qf * vqstr + afpb
|
||||
evroe(1,ivx,5) = css
|
||||
evroe(1,ivy,5) = - evroe(1,ivy,3)
|
||||
evroe(1,ivz,5) = - evroe(1,ivz,3)
|
||||
evroe(1,iby,5) = evroe(1,iby,3)
|
||||
evroe(1,ibz,5) = evroe(1,ibz,3)
|
||||
|
||||
! right-going Alfvèn wave
|
||||
!
|
||||
evroe(1,idn,6) = - evroe(1,idn,2)
|
||||
evroe(1,ivy,6) = - evroe(1,ivy,2)
|
||||
evroe(1,ivz,6) = - evroe(1,ivz,2)
|
||||
evroe(1,iby,6) = evroe(1,iby,2)
|
||||
evroe(1,ibz,6) = evroe(1,ibz,2)
|
||||
|
||||
! right-going fast wave
|
||||
!
|
||||
evroe(1,idn,7) = - cff * c(1) + qs * vqstr - aspb
|
||||
evroe(1,ivx,7) = cff
|
||||
evroe(1,ivy,7) = - evroe(1,ivy,1)
|
||||
evroe(1,ivz,7) = - evroe(1,ivz,1)
|
||||
evroe(1,iby,7) = evroe(1,iby,1)
|
||||
evroe(1,ibz,7) = evroe(1,ibz,1)
|
||||
|
||||
! copy matrices of eigenvectors
|
||||
!
|
||||
l(1:nv,1:nv) = evroe(1,1:nv,1:nv)
|
||||
r(1:nv,1:nv) = evroe(2,1:nv,1:nv)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine esystem_roe_mhd_iso
|
||||
!
|
||||
!*******************************************************************************
|
||||
!
|
||||
! ADIABATIC MAGNETOHYDRODYNAMIC EQUATIONS
|
||||
@ -1847,6 +2405,329 @@ module equations
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end function maxspeed_mhd_adi
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine ESYSTEM_ROE_MHD_ADI:
|
||||
! ------------------------------
|
||||
!
|
||||
! Subroutine computes eigenvalues and eigenvectors for a given set of
|
||||
! equations and input variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! x - ratio of the perpendicular magnetic field component difference
|
||||
! y - ratio of the density
|
||||
! q - the intermediate Roe state vector;
|
||||
! c - the vector of eigenvalues;
|
||||
! r - the matrix of right eigenvectors;
|
||||
! l - the matrix of left eigenvectors;
|
||||
!
|
||||
! 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] Balsara, D. S.
|
||||
! "Linearized Formulation of the Riemann Problem for Adiabatic and
|
||||
! Isothermal Magnetohydrodynamics",
|
||||
! The Astrophysical Journal Suplement Series, 1998, 116, pp. 119-131
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine esystem_roe_mhd_adi(x, y, q, c, r, l)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
real(kind=8) , intent(in) :: x, y
|
||||
real(kind=8), dimension(nv) , intent(in) :: q
|
||||
real(kind=8), dimension(nv) , intent(inout) :: c
|
||||
real(kind=8), dimension(nv,nv), intent(inout) :: l, r
|
||||
|
||||
! saved variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: gammam2
|
||||
|
||||
! local variables
|
||||
!
|
||||
real(kind=8) :: di, vsq, btsq, bt_starsq, casq, hp, twid_asq
|
||||
real(kind=8) :: ct2, tsum, tdif, cf2_cs2, cfsq, cf, cssq, cs, ca, bt
|
||||
real(kind=8) :: bt_star, bet2, bet3, bet2_star, bet3_star, bet_starsq, vbet
|
||||
real(kind=8) :: alpha_f, alpha_s, af_prime, as_prime
|
||||
real(kind=8) :: sqrtd, isqrtd, s, twid_a, qf, qs, afpbb, aspbb
|
||||
real(kind=8) :: qa, qb, qc, qd
|
||||
real(kind=8) :: norm, cff, css, af, as, afpb, aspb
|
||||
real(kind=8) :: q2_star, q3_star, vqstr
|
||||
|
||||
! local parameters
|
||||
!
|
||||
real(kind=8), parameter :: eps = epsilon(di)
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! prepare the internal arrays at the first run
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! prepare coefficients
|
||||
!
|
||||
gammam2 = gamma - 2.0d+00
|
||||
|
||||
! reset all elements
|
||||
!
|
||||
evroe(:, : ,:) = 0.0d+00
|
||||
|
||||
! initiate the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,4,idn) = 1.0d+00
|
||||
|
||||
! unset the first execution flag
|
||||
!
|
||||
first = .false.
|
||||
|
||||
end if ! first execution
|
||||
|
||||
! prepare coefficients for eigenvalues
|
||||
!
|
||||
di = 1.0d+00 / q(idn)
|
||||
casq = q(ibx) * q(ibx) * di
|
||||
ca = sqrt(casq)
|
||||
vsq = q(ivx) * q(ivx) + q(ivy) * q(ivy) + q(ivz) * q(ivz)
|
||||
btsq = q(iby) * q(iby) + q(ibz) * q(ibz)
|
||||
bt_starsq = (gammam1 - gammam2 * y) * btsq
|
||||
hp = q(ien) - (casq + btsq * di)
|
||||
twid_asq = max(eps, (gammam1 * (hp - 0.5d+00 * vsq) - gammam2 * x))
|
||||
ct2 = bt_starsq * di
|
||||
tsum = casq + ct2 + twid_asq
|
||||
tdif = casq + ct2 - twid_asq
|
||||
cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_asq * ct2))
|
||||
cfsq = 0.5d+00 * (tsum + cf2_cs2)
|
||||
cf = sqrt(cfsq)
|
||||
cssq = twid_asq * casq / cfsq
|
||||
cs = sqrt(cssq)
|
||||
|
||||
! prepare eigenvalues
|
||||
!
|
||||
c(1) = q(ivx) - cf
|
||||
c(2) = q(ivx) - ca
|
||||
c(3) = q(ivx) - cs
|
||||
c(4) = q(ivx)
|
||||
c(5) = q(ivx)
|
||||
c(6) = q(ivx) + cs
|
||||
c(7) = q(ivx) + ca
|
||||
c(8) = q(ivx) + cf
|
||||
c(9) = c(8)
|
||||
|
||||
! calculate the eigenvectors only if the waves propagate in both direction
|
||||
!
|
||||
if (c(1) >= 0.0d+00) return
|
||||
if (c(8) <= 0.0d+00) return
|
||||
|
||||
! prepare remaining coefficients for eigenvectors
|
||||
!
|
||||
bt = sqrt(btsq)
|
||||
bt_star = sqrt(bt_starsq)
|
||||
if (bt == 0.0d+00) then
|
||||
bet2 = 1.0d+00
|
||||
bet3 = 0.0d+00
|
||||
else
|
||||
bet2 = q(iby) / bt
|
||||
bet3 = q(ibz) / bt
|
||||
end if
|
||||
bet2_star = bet2 / sqrt(gammam1 - gammam2 * y)
|
||||
bet3_star = bet3 / sqrt(gammam1 - gammam2 * y)
|
||||
bet_starsq = bet2_star * bet2_star + bet3_star * bet3_star
|
||||
vbet = q(ivy) * bet2_star + q(ivz) * bet3_star
|
||||
|
||||
if ( (cfsq - cssq) == 0.0d+00 ) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
else if ( (twid_asq - cssq) <= 0.0d+00 ) then
|
||||
alpha_f = 0.0d+00
|
||||
alpha_s = 1.0d+00
|
||||
else if ( (cfsq - twid_asq) <= 0.0d+00 ) then
|
||||
alpha_f = 1.0d+00
|
||||
alpha_s = 0.0d+00
|
||||
else
|
||||
alpha_f = sqrt((twid_asq - cssq) / (cfsq - cssq))
|
||||
alpha_s = sqrt((cfsq - twid_asq) / (cfsq - cssq))
|
||||
end if
|
||||
|
||||
sqrtd = sqrt(q(idn))
|
||||
isqrtd = 1.0d+00 / sqrtd
|
||||
s = sign(1.0d+00, q(ibx))
|
||||
twid_a = sqrt(twid_asq)
|
||||
qf = cf * alpha_f * s
|
||||
qs = cs * alpha_s * s
|
||||
af_prime = twid_a * alpha_f * isqrtd
|
||||
as_prime = twid_a * alpha_s * isqrtd
|
||||
afpbb = af_prime * bt_star * bet_starsq
|
||||
aspbb = as_prime * bt_star * bet_starsq
|
||||
|
||||
! update the varying elements of the matrix of right eigenvectors
|
||||
!
|
||||
evroe(2,1,idn) = alpha_f
|
||||
evroe(2,3,idn) = alpha_s
|
||||
evroe(2,6,idn) = alpha_s
|
||||
evroe(2,8,idn) = alpha_f
|
||||
|
||||
evroe(2,1,ivx) = alpha_f * c(1)
|
||||
evroe(2,3,ivx) = alpha_s * c(3)
|
||||
evroe(2,4,ivx) = q(ivx)
|
||||
evroe(2,6,ivx) = alpha_s * c(6)
|
||||
evroe(2,8,ivx) = alpha_f * c(8)
|
||||
|
||||
qa = alpha_f * q(ivy)
|
||||
qb = alpha_s * q(ivy)
|
||||
qc = qs * bet2_star
|
||||
qd = qf * bet2_star
|
||||
|
||||
evroe(2,1,ivy) = qa + qc
|
||||
evroe(2,2,ivy) = - bet3
|
||||
evroe(2,3,ivy) = qb - qd
|
||||
evroe(2,4,ivy) = q(ivy)
|
||||
evroe(2,6,ivy) = qb + qd
|
||||
evroe(2,7,ivy) = bet3
|
||||
evroe(2,8,ivy) = qa - qc
|
||||
|
||||
qa = alpha_f * q(ivz)
|
||||
qb = alpha_s * q(ivz)
|
||||
qc = qs * bet3_star
|
||||
qd = qf * bet3_star
|
||||
|
||||
evroe(2,1,ivz) = qa + qc
|
||||
evroe(2,2,ivz) = bet2
|
||||
evroe(2,3,ivz) = qb - qd
|
||||
evroe(2,4,ivz) = q(ivz)
|
||||
evroe(2,6,ivz) = qb + qd
|
||||
evroe(2,7,ivz) = - bet2
|
||||
evroe(2,8,ivz) = qa - qc
|
||||
|
||||
evroe(2,1,ipr) = alpha_f * (hp - q(ivx) * cf) + qs * vbet + aspbb
|
||||
evroe(2,2,ipr) = -(q(ivy) * bet3 - q(ivz) * bet2)
|
||||
evroe(2,3,ipr) = alpha_s * (hp - q(ivx) * cs) - qf * vbet - afpbb
|
||||
evroe(2,4,ipr) = 0.5d+00 * vsq + gammam2 * x / gammam1
|
||||
evroe(2,6,ipr) = alpha_s * (hp + q(ivx) * cs) + qf * vbet - afpbb
|
||||
evroe(2,7,ipr) = - evroe(2,2,ipr)
|
||||
evroe(2,8,ipr) = alpha_f * (hp + q(ivx) * cf) - qs * vbet + aspbb
|
||||
|
||||
evroe(2,1,iby) = as_prime * bet2_star
|
||||
evroe(2,2,iby) = - bet3 * s * isqrtd
|
||||
evroe(2,3,iby) = - af_prime * bet2_star
|
||||
evroe(2,6,iby) = evroe(2,3,iby)
|
||||
evroe(2,7,iby) = evroe(2,2,iby)
|
||||
evroe(2,8,iby) = evroe(2,1,iby)
|
||||
|
||||
evroe(2,1,ibz) = as_prime * bet3_star
|
||||
evroe(2,2,ibz) = bet2 * s * isqrtd
|
||||
evroe(2,3,ibz) = - af_prime * bet3_star
|
||||
evroe(2,6,ibz) = evroe(2,3,ibz)
|
||||
evroe(2,7,ibz) = evroe(2,2,ibz)
|
||||
evroe(2,8,ibz) = evroe(2,1,ibz)
|
||||
|
||||
! update the varying elements of the matrix of left eigenvectors
|
||||
!
|
||||
norm = 0.5d+00 / twid_asq
|
||||
cff = norm * alpha_f * cf
|
||||
css = norm * alpha_s * cs
|
||||
qf = qf * norm
|
||||
qs = qs * norm
|
||||
af = norm * af_prime * q(idn)
|
||||
as = norm * as_prime * q(idn)
|
||||
afpb = norm * af_prime * bt_star
|
||||
aspb = norm * as_prime * bt_star
|
||||
|
||||
norm = norm * gammam1
|
||||
alpha_f = alpha_f * norm
|
||||
alpha_s = alpha_s * norm
|
||||
q2_star = bet2_star / bet_starsq
|
||||
q3_star = bet3_star / bet_starsq
|
||||
vqstr = (q(ivy) * q2_star + q(ivz) * q3_star)
|
||||
norm = 2.0d+00 * norm
|
||||
|
||||
! left-going fast wave
|
||||
!
|
||||
evroe(1,idn,1) = alpha_f * (vsq - hp) + cff * (cf + q(ivx)) &
|
||||
- qs * vqstr - aspb
|
||||
evroe(1,ivx,1) = - alpha_f * q(ivx) - cff
|
||||
evroe(1,ivy,1) = - alpha_f * q(ivy) + qs * q2_star
|
||||
evroe(1,ivz,1) = - alpha_f * q(ivz) + qs * q3_star
|
||||
evroe(1,ipr,1) = alpha_f
|
||||
evroe(1,iby,1) = as * q2_star - alpha_f * q(iby)
|
||||
evroe(1,ibz,1) = as * q3_star - alpha_f * q(ibz)
|
||||
|
||||
! left-going Alfvèn wave
|
||||
!
|
||||
evroe(1,idn,2) = 0.5d+00 * (q(ivy) * bet3 - q(ivz) * bet2)
|
||||
evroe(1,ivy,2) = - 0.5d+00 * bet3
|
||||
evroe(1,ivz,2) = 0.5d+00 * bet2
|
||||
evroe(1,iby,2) = - 0.5d+00 * sqrtd * bet3 * s
|
||||
evroe(1,ibz,2) = 0.5d+00 * sqrtd * bet2 * s
|
||||
|
||||
! left-going slow wave
|
||||
!
|
||||
evroe(1,idn,3) = alpha_s * (vsq - hp) + css * (cs + q(ivx)) &
|
||||
+ qf * vqstr + afpb
|
||||
evroe(1,ivx,3) = - alpha_s * q(ivx) - css
|
||||
evroe(1,ivy,3) = - alpha_s * q(ivy) - qf * q2_star
|
||||
evroe(1,ivz,3) = - alpha_s * q(ivz) - qf * q3_star
|
||||
evroe(1,ipr,3) = alpha_s
|
||||
evroe(1,iby,3) = - af * q2_star - alpha_s * q(iby)
|
||||
evroe(1,ibz,3) = - af * q3_star - alpha_s * q(ibz)
|
||||
|
||||
! entropy wave
|
||||
!
|
||||
evroe(1,idn,4) = 1.0d+00 - norm * (0.5d+00 * vsq - gammam2 * x / gammam1)
|
||||
evroe(1,ivx,4) = norm * q(ivx)
|
||||
evroe(1,ivy,4) = norm * q(ivy)
|
||||
evroe(1,ivz,4) = norm * q(ivz)
|
||||
evroe(1,ipr,4) = - norm
|
||||
evroe(1,iby,4) = norm * q(iby)
|
||||
evroe(1,ibz,4) = norm * q(ibz)
|
||||
|
||||
! right-going slow wave
|
||||
!
|
||||
evroe(1,idn,6) = alpha_s * (vsq - hp) + css * (cs - q(ivx)) &
|
||||
- qf * vqstr + afpb
|
||||
evroe(1,ivx,6) = - alpha_s * q(ivx) + css
|
||||
evroe(1,ivy,6) = - alpha_s * q(ivy) + qf * q2_star
|
||||
evroe(1,ivz,6) = - alpha_s * q(ivz) + qf * q3_star
|
||||
evroe(1,ipr,6) = alpha_s
|
||||
evroe(1,iby,6) = evroe(1,iby,3)
|
||||
evroe(1,ibz,6) = evroe(1,ibz,3)
|
||||
|
||||
! right-going Alfvèn wave
|
||||
!
|
||||
evroe(1,idn,7) = - evroe(1,idn,2)
|
||||
evroe(1,ivy,7) = - evroe(1,ivy,2)
|
||||
evroe(1,ivz,7) = - evroe(1,ivz,2)
|
||||
evroe(1,iby,7) = evroe(1,iby,2)
|
||||
evroe(1,ibz,7) = evroe(1,ibz,2)
|
||||
|
||||
! right-going fast wave
|
||||
!
|
||||
evroe(1,idn,8) = alpha_f * (vsq - hp) + cff * (cf - q(ivx)) &
|
||||
+ qs * vqstr - aspb
|
||||
evroe(1,ivx,8) = - alpha_f * q(ivx) + cff
|
||||
evroe(1,ivy,8) = - alpha_f * q(ivy) - qs * q2_star
|
||||
evroe(1,ivz,8) = - alpha_f * q(ivz) - qs * q3_star
|
||||
evroe(1,ipr,8) = alpha_f
|
||||
evroe(1,iby,8) = evroe(1,iby,1)
|
||||
evroe(1,ibz,8) = evroe(1,ibz,1)
|
||||
|
||||
! copy matrices of eigenvectors
|
||||
!
|
||||
l(1:nv,1:nv) = evroe(1,1:nv,1:nv)
|
||||
r(1:nv,1:nv) = evroe(2,1:nv,1:nv)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine esystem_roe_mhd_adi
|
||||
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
17
src/io.F90
17
src/io.F90
@ -71,7 +71,7 @@ module io
|
||||
|
||||
! flags to determine the way of data writing
|
||||
!
|
||||
logical , save :: with_ghosts = .false.
|
||||
logical , save :: with_ghosts = .true.
|
||||
|
||||
! local variables to store the number of processors and maximum level read from
|
||||
! the restart file
|
||||
@ -141,7 +141,7 @@ module io
|
||||
|
||||
! local variables
|
||||
!
|
||||
character(len=255) :: ghosts = "off"
|
||||
character(len=255) :: ghosts = "on"
|
||||
integer :: dd, hh, mm, ss
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -171,15 +171,15 @@ module io
|
||||
|
||||
! get the flag determining if the ghost cells are stored
|
||||
!
|
||||
call get_parameter_string ("include_ghosts", ghosts )
|
||||
call get_parameter_string ("include_ghosts" , ghosts )
|
||||
|
||||
! check ghost cell storing flag
|
||||
!
|
||||
select case(trim(ghosts))
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
with_ghosts = .true.
|
||||
case default
|
||||
case ("off", "OFF", "n", "N", "false", "FALSE", "no", "NO")
|
||||
with_ghosts = .false.
|
||||
case default
|
||||
with_ghosts = .true.
|
||||
end select
|
||||
|
||||
! return the run number
|
||||
@ -193,6 +193,11 @@ module io
|
||||
"snapshot type", "primitive variables"
|
||||
if (ftype == 'c') write (*,"(4x,a13,10x,'=',1x,a)") &
|
||||
"snapshot type", "conservative variables"
|
||||
if (with_ghosts) then
|
||||
write (*,"(4x,a21,2x,'=',1x,a)") "with ghosts cells ", "on"
|
||||
else
|
||||
write (*,"(4x,a21,2x,'=',1x,a)") "with ghosts cells ", "off"
|
||||
end if
|
||||
write (*,"(4x,a21,2x,'=',1x,e9.2)") "snapshot interval ", hsnap
|
||||
if (hrest > 0.0d+00) then
|
||||
dd = int(hrest / 2.4d+01)
|
||||
|
||||
2421
src/schemes.F90
2421
src/schemes.F90
File diff suppressed because it is too large
Load Diff
Loading…
x
Reference in New Issue
Block a user