diff --git a/src/equations.F90 b/src/equations.F90 index 108b75b..bde35e8 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -63,16 +63,20 @@ 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() + 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 @@ -237,10 +246,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_iso - cons2prim => cons2prim_hd_iso - fluxspeed => fluxspeed_hd_iso - maxspeed => maxspeed_hd_iso + prim2cons => prim2cons_hd_iso + cons2prim => cons2prim_hd_iso + fluxspeed => fluxspeed_hd_iso + maxspeed => maxspeed_hd_iso + eigensystem_roe => esystem_roe_hd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -259,10 +269,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_adi - cons2prim => cons2prim_hd_adi - fluxspeed => fluxspeed_hd_adi - maxspeed => maxspeed_hd_adi + prim2cons => prim2cons_hd_adi + 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 ! @@ -337,10 +348,11 @@ module equations ! set pointers to the subroutines ! - prim2cons => prim2cons_mhd_iso - cons2prim => cons2prim_mhd_iso - fluxspeed => fluxspeed_mhd_iso - maxspeed => maxspeed_mhd_iso + prim2cons => prim2cons_mhd_iso + cons2prim => cons2prim_mhd_iso + fluxspeed => fluxspeed_mhd_iso + maxspeed => maxspeed_mhd_iso + eigensystem_roe => esystem_roe_mhd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -359,10 +371,11 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_mhd_adi - cons2prim => cons2prim_mhd_adi - fluxspeed => fluxspeed_mhd_adi - maxspeed => maxspeed_mhd_adi + prim2cons => prim2cons_mhd_adi + 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 !=============================================================================== ! diff --git a/src/io.F90 b/src/io.F90 index 7c95c6c..cecf5a3 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -61,17 +61,17 @@ module io ! (in hours); ! hsnap - the problem time interval for regular snapshot storing; ! - character(len=255), save :: respath = "./" - character , save :: ftype = "p" - integer , save :: nrest = -1 - integer(kind=4) , save :: irest = 1 - integer(kind=4) , save :: isnap = -1 - real , save :: hrest = 6.0d+00 - real , save :: hsnap = 1.0d+00 + character(len=255), save :: respath = "./" + character , save :: ftype = "p" + integer , save :: nrest = -1 + integer(kind=4) , save :: irest = 1 + integer(kind=4) , save :: isnap = -1 + real , save :: hrest = 6.0d+00 + real , save :: hsnap = 1.0d+00 ! 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) diff --git a/src/schemes.F90 b/src/schemes.F90 index 47904f3..d1d53a5 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -46,13 +46,17 @@ module schemes #ifdef PROFILE ! timer indices ! - integer , save :: imi, imu, imf, imr + integer , save :: imi, imu, imf, ims, imr #endif /* PROFILE */ ! pointer to the flux update procedure ! procedure(update_flux_hd_iso), pointer, save :: update_flux => null() +! pointer to the state reconstruction +! + procedure(states_hd_iso) , pointer, save :: states => null() + ! pointer to the Riemann solver ! procedure(riemann_hd_iso_hll), pointer, save :: riemann => null() @@ -93,8 +97,8 @@ module schemes ! include external procedures and variables ! - use equations , only : eqsys, eos - use parameters, only : get_parameter_string + use equations , only : eqsys, eos + use parameters , only : get_parameter_string ! local variables are not implicit by default ! @@ -118,6 +122,7 @@ module schemes call set_timer('schemes:: initialization' , imi) call set_timer('schemes:: increment update', imu) call set_timer('schemes:: flux update' , imf) + call set_timer('schemes:: Riemann states' , ims) call set_timer('schemes:: Riemann solver' , imr) ! start accounting time for module initialization/finalization @@ -146,11 +151,21 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_hd_iso + states => states_hd_iso ! select the Riemann solver ! select case(trim(solver)) + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_hd_iso_roe ! in the case of unknown Riemann solver, revert to HLL ! @@ -171,6 +186,7 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_hd_adi + states => states_hd_adi ! select the Riemann solver ! @@ -186,6 +202,16 @@ module schemes ! riemann => riemann_hd_adi_hllc + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_hd_adi_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -215,6 +241,7 @@ module schemes ! set pointers to the subroutines ! update_flux => update_flux_mhd_iso + states => states_mhd_iso ! select the Riemann solver ! @@ -240,6 +267,16 @@ module schemes ! riemann => riemann_mhd_iso_hlldm + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_mhd_iso_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -259,6 +296,7 @@ module schemes ! set pointers to subroutines ! update_flux => update_flux_mhd_adi + states => states_mhd_adi ! select the Riemann solver ! @@ -284,6 +322,16 @@ module schemes ! riemann => riemann_mhd_adi_hlld + case("roe", "ROE") + +! set the solver name +! + name_sol = "ROE" + +! set pointers to subroutines +! + riemann => riemann_mhd_adi_roe + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -354,6 +402,7 @@ module schemes ! nullify procedure pointers ! nullify(update_flux) + nullify(states) nullify(riemann) #ifdef PROFILE @@ -386,8 +435,8 @@ module schemes ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv ! local variables are not implicit by default ! @@ -451,6 +500,8 @@ module schemes !! !=============================================================================== ! +!***** ISOTHERMAL HYDRODYNAMICS ***** +! !=============================================================================== ! ! subroutine UPDATE_FLUX_HD_ISO: @@ -467,16 +518,15 @@ module schemes ! q - the array of primitive variables; ! f - the array of numerical fluxes; ! -! !=============================================================================== ! subroutine update_flux_hd_iso(idir, dx, q, f) ! include external variables ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz ! local variables are not implicit by default ! @@ -491,14 +541,14 @@ module schemes ! local variables ! - integer :: i, j, k + integer :: i, j, k ! local temporary arrays ! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy #if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- @@ -530,9 +580,13 @@ module schemes qx(ivy,1:im) = q(ivy,1:im,j,k) qx(ivz,1:im) = q(ivz,1:im,j,k) +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) ! update the array of fluxes ! @@ -558,9 +612,13 @@ module schemes qy(ivy,1:jm) = q(ivz,i,1:jm,k) qy(ivz,1:jm) = q(ivx,i,1:jm,k) +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) ! update the array of fluxes ! @@ -587,9 +645,13 @@ module schemes qz(ivy,1:km) = q(ivx,i,j,1:km) qz(ivz,1:km) = q(ivy,i,j,1:km) +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + ! call one dimensional Riemann solver in order to obtain numerical fluxes ! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) ! update the array of fluxes ! @@ -616,554 +678,71 @@ module schemes ! !=============================================================================== ! -! subroutine UPDATE_FLUX_HD_ADI: -! ----------------------------- +! subroutine STATES_HD_ISO: +! ------------------------ ! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. +! Subroutine reconstructs the Riemann states. ! ! Arguments: ! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; ! !=============================================================================== ! - subroutine update_flux_hd_adi(idir, dx, q, f) + subroutine states_hd_iso(n, h, q, ql, qr) -! include external variables +! include external procedures ! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + use equations , only : nv + use equations , only : idn + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! implicit none -! input arguments +! subroutine arguments ! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + integer , intent(in) :: n + real(kind=8) , intent(in) :: h + real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr ! local variables ! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ + integer :: p ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for flux update +! start accounting time for the state reconstruction ! - call start_timer(imf) + call start_timer(ims) #endif /* PROFILE */ -! reset the flux array +! reconstruct the left and right states of primitive variables ! - f(:,:,:,:) = 0.0d+00 + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do -! select the directional flux to compute +! check if the reconstruction gives negative values of density, +! if so, correct the states ! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ipr,1:im) = q(ipr,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ien,1:im,j,k) = fx(ien,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ipr,1:jm) = q(ipr,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ien,i,1:jm,k) = fy(ien,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ipr,1:km) = q(ipr,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ien,i,j,1:km) = fz(ien,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) #ifdef PROFILE -! stop accounting time for flux update +! stop accounting time for the state reconstruction ! - call stop_timer(imf) + call stop_timer(ims) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! - end subroutine update_flux_hd_adi -! -!=============================================================================== -! -! subroutine UPDATE_FLUX_MHD_ADI: -! ------------------------------ -! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. -! -! Arguments: -! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! -! -!=============================================================================== -! - subroutine update_flux_mhd_iso(idir, dx, q, f) - -! include external variables -! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz - use equations , only : ibx, iby, ibz, ibp - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f - -! local variables -! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for flux update -! - call start_timer(imf) -#endif /* PROFILE */ - -! reset the flux array -! - f(:,:,:,:) = 0.0d+00 - -! select the directional flux to compute -! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ibx,1:im) = q(ibx,1:im,j,k) - qx(iby,1:im) = q(iby,1:im,j,k) - qx(ibz,1:im) = q(ibz,1:im,j,k) - qx(ibp,1:im) = q(ibp,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ibx,1:im,j,k) = fx(ibx,1:im) - f(iby,1:im,j,k) = fx(iby,1:im) - f(ibz,1:im,j,k) = fx(ibz,1:im) - f(ibp,1:im,j,k) = fx(ibp,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ibx,1:jm) = q(iby,i,1:jm,k) - qy(iby,1:jm) = q(ibz,i,1:jm,k) - qy(ibz,1:jm) = q(ibx,i,1:jm,k) - qy(ibp,1:jm) = q(ibp,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ibx,i,1:jm,k) = fy(ibz,1:jm) - f(iby,i,1:jm,k) = fy(ibx,1:jm) - f(ibz,i,1:jm,k) = fy(iby,1:jm) - f(ibp,i,1:jm,k) = fy(ibp,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ibx,1:km) = q(ibz,i,j,1:km) - qz(iby,1:km) = q(ibx,i,j,1:km) - qz(ibz,1:km) = q(iby,i,j,1:km) - qz(ibp,1:km) = q(ibp,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ibx,i,j,1:km) = fz(iby,1:km) - f(iby,i,j,1:km) = fz(ibz,1:km) - f(ibz,i,j,1:km) = fz(ibx,1:km) - f(ibp,i,j,1:km) = fz(ibp,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select - -#ifdef PROFILE -! stop accounting time for flux update -! - call stop_timer(imf) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_flux_mhd_iso -! -!=============================================================================== -! -! subroutine UPDATE_FLUX_MHD_ADI: -! ------------------------------ -! -! Subroutine solves the Riemann problem along each direction and calculates -! the numerical fluxes, which are used later to calculate the conserved -! variable increment. -! -! Arguments: -! -! idir - direction along which the flux is calculated; -! dx - the spatial step; -! q - the array of primitive variables; -! f - the array of numerical fluxes; -! -! -!=============================================================================== -! - subroutine update_flux_mhd_adi(idir, dx, q, f) - -! include external variables -! - use coordinates, only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien - use equations , only : ibx, iby, ibz, ibp - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - integer , intent(in) :: idir - real(kind=8) , intent(in) :: dx - real(kind=8), dimension(nv,im,jm,km), intent(in) :: q - real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f - -! local variables -! - integer :: i, j, k - -! local temporary arrays -! - real(kind=8), dimension(nv,im) :: qx, fx - real(kind=8), dimension(nv,jm) :: qy, fy -#if NDIMS == 3 - real(kind=8), dimension(nv,km) :: qz, fz -#endif /* NDIMS == 3 */ -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for flux update -! - call start_timer(imf) -#endif /* PROFILE */ - -! reset the flux array -! - f(:,:,:,:) = 0.0d+00 - -! select the directional flux to compute -! - select case(idir) - case(1) - -! calculate the flux along the X-direction -! - do k = kbl, keu - do j = jbl, jeu - -! copy directional variable vectors to pass to the one dimensional solver -! - qx(idn,1:im) = q(idn,1:im,j,k) - qx(ivx,1:im) = q(ivx,1:im,j,k) - qx(ivy,1:im) = q(ivy,1:im,j,k) - qx(ivz,1:im) = q(ivz,1:im,j,k) - qx(ibx,1:im) = q(ibx,1:im,j,k) - qx(iby,1:im) = q(iby,1:im,j,k) - qx(ibz,1:im) = q(ibz,1:im,j,k) - qx(ibp,1:im) = q(ibp,1:im,j,k) - qx(ipr,1:im) = q(ipr,1:im,j,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im)) - -! update the array of fluxes -! - f(idn,1:im,j,k) = fx(idn,1:im) - f(imx,1:im,j,k) = fx(imx,1:im) - f(imy,1:im,j,k) = fx(imy,1:im) - f(imz,1:im,j,k) = fx(imz,1:im) - f(ibx,1:im,j,k) = fx(ibx,1:im) - f(iby,1:im,j,k) = fx(iby,1:im) - f(ibz,1:im,j,k) = fx(ibz,1:im) - f(ibp,1:im,j,k) = fx(ibp,1:im) - f(ien,1:im,j,k) = fx(ien,1:im) - - end do ! j = jbl, jeu - end do ! k = kbl, keu - - case(2) - -! calculate the flux along the Y direction -! - do k = kbl, keu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qy(idn,1:jm) = q(idn,i,1:jm,k) - qy(ivx,1:jm) = q(ivy,i,1:jm,k) - qy(ivy,1:jm) = q(ivz,i,1:jm,k) - qy(ivz,1:jm) = q(ivx,i,1:jm,k) - qy(ibx,1:jm) = q(iby,i,1:jm,k) - qy(iby,1:jm) = q(ibz,i,1:jm,k) - qy(ibz,1:jm) = q(ibx,i,1:jm,k) - qy(ibp,1:jm) = q(ibp,i,1:jm,k) - qy(ipr,1:jm) = q(ipr,i,1:jm,k) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm)) - -! update the array of fluxes -! - f(idn,i,1:jm,k) = fy(idn,1:jm) - f(imx,i,1:jm,k) = fy(imz,1:jm) - f(imy,i,1:jm,k) = fy(imx,1:jm) - f(imz,i,1:jm,k) = fy(imy,1:jm) - f(ibx,i,1:jm,k) = fy(ibz,1:jm) - f(iby,i,1:jm,k) = fy(ibx,1:jm) - f(ibz,i,1:jm,k) = fy(iby,1:jm) - f(ibp,i,1:jm,k) = fy(ibp,1:jm) - f(ien,i,1:jm,k) = fy(ien,1:jm) - - end do ! i = ibl, ieu - end do ! k = kbl, keu - -#if NDIMS == 3 - case(3) - -! calculate the flux along the Z direction -! - do j = jbl, jeu - do i = ibl, ieu - -! copy directional variable vectors to pass to the one dimensional solver -! - qz(idn,1:km) = q(idn,i,j,1:km) - qz(ivx,1:km) = q(ivz,i,j,1:km) - qz(ivy,1:km) = q(ivx,i,j,1:km) - qz(ivz,1:km) = q(ivy,i,j,1:km) - qz(ibx,1:km) = q(ibz,i,j,1:km) - qz(iby,1:km) = q(ibx,i,j,1:km) - qz(ibz,1:km) = q(iby,i,j,1:km) - qz(ibp,1:km) = q(ibp,i,j,1:km) - qz(ipr,1:km) = q(ipr,i,j,1:km) - -! call one dimensional Riemann solver in order to obtain numerical fluxes -! - call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km)) - -! update the array of fluxes -! - f(idn,i,j,1:km) = fz(idn,1:km) - f(imx,i,j,1:km) = fz(imy,1:km) - f(imy,i,j,1:km) = fz(imz,1:km) - f(imz,i,j,1:km) = fz(imx,1:km) - f(ibx,i,j,1:km) = fz(iby,1:km) - f(iby,i,j,1:km) = fz(ibz,1:km) - f(ibz,i,j,1:km) = fz(ibx,1:km) - f(ibp,i,j,1:km) = fz(ibp,1:km) - f(ien,i,j,1:km) = fz(ien,1:km) - - end do ! i = ibl, ieu - end do ! j = jbl, jeu -#endif /* NDIMS == 3 */ - - end select - -#ifdef PROFILE -! stop accounting time for flux update -! - call stop_timer(imf) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_flux_mhd_adi + end subroutine states_hd_iso ! !=============================================================================== ! @@ -1175,10 +754,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1189,14 +767,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_hd_iso_hll(n, h, q, f) + subroutine riemann_hd_iso_hll(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -1205,18 +782,17 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -1228,17 +804,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1299,36 +864,358 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_HD_ADI_HLL: +! subroutine RIEMANN_HD_ISO_ROE: ! ----------------------------- ! ! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! the Roe's method. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! - subroutine riemann_hd_adi_hll(n, h, q, f) + subroutine riemann_hd_iso_roe(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ipr, ivx - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + else + do p = nv, 1, - 1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_hd_iso_roe +! +!=============================================================================== +! +!***** ADIABATIC HYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_HD_ADI: +! ----------------------------- +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_hd_adi(idir, dx, q, f) + +! include external variables +! + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ipr,1:im) = q(ipr,1:im,j,k) + +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ien,1:im,j,k) = fx(ien,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ipr,1:jm) = q(ipr,i,1:jm,k) + +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ien,i,1:jm,k) = fy(ien,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ipr,1:km) = q(ipr,i,j,1:km) + +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ien,i,j,1:km) = fz(ien,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_hd_adi +! +!=============================================================================== +! +! subroutine STATES_HD_ADI: +! ------------------------ +! +! Subroutine reconstructs the Riemann states. +! +! Arguments: +! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; +! +!=============================================================================== +! + subroutine states_hd_adi(n, h, q, ql, qr) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ipr + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -1339,16 +1226,91 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr + +! local variables +! + integer :: p +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the state reconstruction +! + call start_timer(ims) +#endif /* PROFILE */ + +! reconstruct the left and right states of primitive variables +! + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do + +! check if the reconstruction gives negative values of density or pressure, +! if so, correct the states +! + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) + call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) + +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_hd_adi +! +!=============================================================================== +! +! subroutine RIEMANN_HD_ADI_HLL: +! ----------------------------- +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_hd_adi_hll(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -1360,18 +1322,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1441,10 +1391,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1454,14 +1403,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_hd_adi_hllc(n, h, q, f) + subroutine riemann_hd_adi_hllc(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -1470,20 +1418,19 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm real(kind=8) :: srml, slmm, srmm real(kind=8) :: dn, pr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -1495,18 +1442,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1623,37 +1558,380 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_MHD_ISO_HLL: +! subroutine RIEMANN_HD_ADI_ROE: ! ----------------------------- ! ! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! the Roe's method. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! [1] Roe, P. L. +! "Approximate Riemann Solvers, Parameter Vectors, and Difference +! Schemes", +! Journal of Computational Physics, 1981, 43, pp. 357-372 +! [2] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! - subroutine riemann_mhd_iso_hll(n, h, q, f) + subroutine riemann_hd_adi_roe(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ibx, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ien + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ipr) = ((ul(ien,i) + ql(ipr,i)) / sdl & + + (ur(ien,i) + qr(ipr,i)) / sdr) / sds + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(0.0d+00, 0.0d+00, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(1:nv,i) = 0.5d+00 * (fl(1:nv,i) + fr(1:nv,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + else + do p = nv, 1, - 1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + f(1:nv,i) = f(1:nv,i) + xx * ri(p,1:nv) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_hd_adi_roe +! +!=============================================================================== +! +!***** ISOTHERMAL MAGNETOHYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_MHD_ISO: +! ------------------------------ +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_mhd_iso(idir, dx, q, f) + +! include external variables +! + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz + use equations , only : ibx, iby, ibz, ibp + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ibx,1:im) = q(ibx,1:im,j,k) + qx(iby,1:im) = q(iby,1:im,j,k) + qx(ibz,1:im) = q(ibz,1:im,j,k) + qx(ibp,1:im) = q(ibp,1:im,j,k) + +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ibx,1:im,j,k) = fx(ibx,1:im) + f(iby,1:im,j,k) = fx(iby,1:im) + f(ibz,1:im,j,k) = fx(ibz,1:im) + f(ibp,1:im,j,k) = fx(ibp,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ibx,1:jm) = q(iby,i,1:jm,k) + qy(iby,1:jm) = q(ibz,i,1:jm,k) + qy(ibz,1:jm) = q(ibx,i,1:jm,k) + qy(ibp,1:jm) = q(ibp,i,1:jm,k) + +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ibx,i,1:jm,k) = fy(ibz,1:jm) + f(iby,i,1:jm,k) = fy(ibx,1:jm) + f(ibz,i,1:jm,k) = fy(iby,1:jm) + f(ibp,i,1:jm,k) = fy(ibp,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ibx,1:km) = q(ibz,i,j,1:km) + qz(iby,1:km) = q(ibx,i,j,1:km) + qz(ibz,1:km) = q(iby,i,j,1:km) + qz(ibp,1:km) = q(ibp,i,j,1:km) + +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ibx,i,j,1:km) = fz(iby,1:km) + f(iby,i,j,1:km) = fz(ibz,1:km) + f(ibz,i,j,1:km) = fz(ibx,1:km) + f(ibp,i,j,1:km) = fz(ibp,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_mhd_iso +! +!=============================================================================== +! +! subroutine STATES_MHD_ISO: +! ------------------------- +! +! Subroutine reconstructs the Riemann states. +! +! Arguments: +! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; +! +!=============================================================================== +! + subroutine states_mhd_iso(n, h, q, ql, qr) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ibx, ibp + use equations , only : cmax + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -1664,25 +1942,19 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q - real(kind=8), dimension(nv,n), intent(out) :: f + real(kind=8), dimension(nv,n), intent(out) :: ql, qr ! local variables ! - integer :: p, i - real(kind=8) :: sl, sr, srml - -! local arrays to store the states -! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr - real(kind=8), dimension(nv) :: wl, wr - real(kind=8), dimension(n) :: cl, cr + integer :: i, p + real(kind=8) :: bx, bp ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for Riemann solver +! start accounting time for the state reconstruction ! - call start_timer(imr) + call start_timer(ims) #endif /* PROFILE */ ! reconstruct the left and right states of primitive variables @@ -1693,18 +1965,95 @@ module schemes ! obtain the state values for Bx and Psi for the GLM-MHD equations ! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) + do i = 1, n -! check if the reconstruction doesn't give the negative density or pressure, + bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) & + - (qr(ibp,i) - ql(ibp,i)) / cmax) + bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) & + - (qr(ibx,i) - ql(ibx,i)) * cmax) + + ql(ibx,i) = bx + qr(ibx,i) = bx + ql(ibp,i) = bp + qr(ibp,i) = bp + + end do ! i = 1, n + +! check if the reconstruction gives negative values of density, ! if so, correct the states ! call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_mhd_iso +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ISO_HLL: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_mhd_iso_hll(n, ql, qr, f) + +! include external variables and procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: i + real(kind=8) :: sl, sr, srml + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(nv) :: wl, wr + real(kind=8), dimension(n) :: cl, cr +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -1773,10 +2122,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -1787,15 +2135,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_iso_hlld(n, h, q, f) + subroutine riemann_mhd_iso_hlld(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -1804,19 +2150,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -1828,26 +2173,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -2166,10 +2491,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -2179,15 +2503,13 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_iso_hlldm(n, h, q, f) + subroutine riemann_mhd_iso_hlldm(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, imx, imy, imz, ibx, iby, ibz, ibp + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -2196,19 +2518,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm real(kind=8) :: bx, b2, dn, dnl, dnr, dvl, dvr, ca ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -2220,26 +2541,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -2548,37 +2849,416 @@ module schemes ! !=============================================================================== ! -! subroutine RIEMANN_MHD_ADI_HLL: +! subroutine RIEMANN_MHD_ISO_ROE: ! ------------------------------ ! ! Subroutine solves one dimensional Riemann problem using -! the Harten-Lax-van Leer (HLL) method. +! the Roe's method. ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! -! [1] Harten, A., Lax, P. D. & Van Leer, B. -! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic -! Conservation Laws", -! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! [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] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! - subroutine riemann_mhd_adi_hll(n, h, q, f) + subroutine riemann_mhd_iso_roe(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ipr, ivx, ibx, ibp - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp + use equations , only : imx, imy, imz + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: xx, yy + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ibx) = ql(ibx,i) + qi(iby) = (sdr * ql(iby,i) + sdl * qr(iby,i)) / sds + qi(ibz) = (sdr * ql(ibz,i) + sdl * qr(ibz,i)) / sds + qi(ibp) = ql(ibp,i) + +! prepare coefficients +! + xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 & + + (ql(ibz,i) - qr(ibz,i))**2) / sds**2 + yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn) + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! prepare fluxes which do not change across the states +! + f(ibx,i) = fl(ibx,i) + f(ibp,i) = fl(ibp,i) + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) + f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) + f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) + f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) + f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) + f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + else + do p = nv, 1, -1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_iso_roe +! +!=============================================================================== +! +!***** ADIABATIC MAGNETOHYDRODYNAMICS ***** +! +!=============================================================================== +! +! subroutine UPDATE_FLUX_MHD_ADI: +! ------------------------------ +! +! Subroutine solves the Riemann problem along each direction and calculates +! the numerical fluxes, which are used later to calculate the conserved +! variable increment. +! +! Arguments: +! +! idir - direction along which the flux is calculated; +! dx - the spatial step; +! q - the array of primitive variables; +! f - the array of numerical fluxes; +! +!=============================================================================== +! + subroutine update_flux_mhd_adi(idir, dx, q, f) + +! include external variables +! + use coordinates , only : im, jm, km, ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien + use equations , only : ibx, iby, ibz, ibp + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + integer , intent(in) :: idir + real(kind=8) , intent(in) :: dx + real(kind=8), dimension(nv,im,jm,km), intent(in) :: q + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: f + +! local variables +! + integer :: i, j, k + +! local temporary arrays +! + real(kind=8), dimension(nv,im) :: qx, qxl, qxr, fx + real(kind=8), dimension(nv,jm) :: qy, qyl, qyr, fy +#if NDIMS == 3 + real(kind=8), dimension(nv,km) :: qz, qzl, qzr, fz +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for flux update +! + call start_timer(imf) +#endif /* PROFILE */ + +! reset the flux array +! + f(:,:,:,:) = 0.0d+00 + +! select the directional flux to compute +! + select case(idir) + case(1) + +! calculate the flux along the X-direction +! + do k = kbl, keu + do j = jbl, jeu + +! copy directional variable vectors to pass to the one dimensional solver +! + qx(idn,1:im) = q(idn,1:im,j,k) + qx(ivx,1:im) = q(ivx,1:im,j,k) + qx(ivy,1:im) = q(ivy,1:im,j,k) + qx(ivz,1:im) = q(ivz,1:im,j,k) + qx(ibx,1:im) = q(ibx,1:im,j,k) + qx(iby,1:im) = q(iby,1:im,j,k) + qx(ibz,1:im) = q(ibz,1:im,j,k) + qx(ibp,1:im) = q(ibp,1:im,j,k) + qx(ipr,1:im) = q(ipr,1:im,j,k) + +! reconstruct Riemann states +! + call states(im, dx, qx(1:nv,1:im), qxl(1:nv,1:im), qxr(1:nv,1:im)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(im, qxl(1:nv,1:im), qxr(1:nv,1:im), fx(1:nv,1:im)) + +! update the array of fluxes +! + f(idn,1:im,j,k) = fx(idn,1:im) + f(imx,1:im,j,k) = fx(imx,1:im) + f(imy,1:im,j,k) = fx(imy,1:im) + f(imz,1:im,j,k) = fx(imz,1:im) + f(ibx,1:im,j,k) = fx(ibx,1:im) + f(iby,1:im,j,k) = fx(iby,1:im) + f(ibz,1:im,j,k) = fx(ibz,1:im) + f(ibp,1:im,j,k) = fx(ibp,1:im) + f(ien,1:im,j,k) = fx(ien,1:im) + + end do ! j = jbl, jeu + end do ! k = kbl, keu + + case(2) + +! calculate the flux along the Y direction +! + do k = kbl, keu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qy(idn,1:jm) = q(idn,i,1:jm,k) + qy(ivx,1:jm) = q(ivy,i,1:jm,k) + qy(ivy,1:jm) = q(ivz,i,1:jm,k) + qy(ivz,1:jm) = q(ivx,i,1:jm,k) + qy(ibx,1:jm) = q(iby,i,1:jm,k) + qy(iby,1:jm) = q(ibz,i,1:jm,k) + qy(ibz,1:jm) = q(ibx,i,1:jm,k) + qy(ibp,1:jm) = q(ibp,i,1:jm,k) + qy(ipr,1:jm) = q(ipr,i,1:jm,k) + +! reconstruct Riemann states +! + call states(jm, dx, qy(1:nv,1:jm), qyl(1:nv,1:jm), qyr(1:nv,1:jm)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(jm, qyl(1:nv,1:jm), qyr(1:nv,1:jm), fy(1:nv,1:jm)) + +! update the array of fluxes +! + f(idn,i,1:jm,k) = fy(idn,1:jm) + f(imx,i,1:jm,k) = fy(imz,1:jm) + f(imy,i,1:jm,k) = fy(imx,1:jm) + f(imz,i,1:jm,k) = fy(imy,1:jm) + f(ibx,i,1:jm,k) = fy(ibz,1:jm) + f(iby,i,1:jm,k) = fy(ibx,1:jm) + f(ibz,i,1:jm,k) = fy(iby,1:jm) + f(ibp,i,1:jm,k) = fy(ibp,1:jm) + f(ien,i,1:jm,k) = fy(ien,1:jm) + + end do ! i = ibl, ieu + end do ! k = kbl, keu + +#if NDIMS == 3 + case(3) + +! calculate the flux along the Z direction +! + do j = jbl, jeu + do i = ibl, ieu + +! copy directional variable vectors to pass to the one dimensional solver +! + qz(idn,1:km) = q(idn,i,j,1:km) + qz(ivx,1:km) = q(ivz,i,j,1:km) + qz(ivy,1:km) = q(ivx,i,j,1:km) + qz(ivz,1:km) = q(ivy,i,j,1:km) + qz(ibx,1:km) = q(ibz,i,j,1:km) + qz(iby,1:km) = q(ibx,i,j,1:km) + qz(ibz,1:km) = q(iby,i,j,1:km) + qz(ibp,1:km) = q(ibp,i,j,1:km) + qz(ipr,1:km) = q(ipr,i,j,1:km) + +! reconstruct Riemann states +! + call states(km, dx, qz(1:nv,1:km), qzl(1:nv,1:km), qzr(1:nv,1:km)) + +! call one dimensional Riemann solver in order to obtain numerical fluxes +! + call riemann(km, qzl(1:nv,1:km), qzr(1:nv,1:km), fz(1:nv,1:km)) + +! update the array of fluxes +! + f(idn,i,j,1:km) = fz(idn,1:km) + f(imx,i,j,1:km) = fz(imy,1:km) + f(imy,i,j,1:km) = fz(imz,1:km) + f(imz,i,j,1:km) = fz(imx,1:km) + f(ibx,i,j,1:km) = fz(iby,1:km) + f(iby,i,j,1:km) = fz(ibz,1:km) + f(ibz,i,j,1:km) = fz(ibx,1:km) + f(ibp,i,j,1:km) = fz(ibp,1:km) + f(ien,i,j,1:km) = fz(ien,1:km) + + end do ! i = ibl, ieu + end do ! j = jbl, jeu +#endif /* NDIMS == 3 */ + + end select + +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_flux_mhd_adi +! +!=============================================================================== +! +! subroutine STATES_MHD_ADI: +! ------------------------- +! +! Subroutine reconstructs the Riemann states. +! +! Arguments: +! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! ql, qr - the reconstructed Riemann states; +! +!=============================================================================== +! + subroutine states_mhd_adi(n, h, q, ql, qr) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ipr, ibx, ibp + use equations , only : cmax + use interpolations , only : reconstruct, fix_positivity ! local variables are not implicit by default ! @@ -2589,16 +3269,108 @@ module schemes integer , intent(in) :: n real(kind=8) , intent(in) :: h real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(out) :: ql, qr + +! local variables +! + integer :: i, p + real(kind=8) :: bx, bp +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the state reconstruction +! + call start_timer(ims) +#endif /* PROFILE */ + +! reconstruct the left and right states of primitive variables +! + do p = 1, nv + call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) + end do ! p = 1, nv + +! obtain the state values for Bx and Psi for the GLM-MHD equations +! + do i = 1, n + + bx = 0.5d+00 * ((qr(ibx,i) + ql(ibx,i)) & + - (qr(ibp,i) - ql(ibp,i)) / cmax) + bp = 0.5d+00 * ((qr(ibp,i) + ql(ibp,i)) & + - (qr(ibx,i) - ql(ibx,i)) * cmax) + + ql(ibx,i) = bx + qr(ibx,i) = bx + ql(ibp,i) = bp + qr(ibp,i) = bp + + end do ! i = 1, n + +! check if the reconstruction gives negative values of density or density, +! if so, correct the states +! + call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) + call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) + +#ifdef PROFILE +! stop accounting time for the state reconstruction +! + call stop_timer(ims) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine states_mhd_adi +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ADI_HLL: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Harten-Lax-van Leer (HLL) method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! References: +! +! [1] Harten, A., Lax, P. D. & Van Leer, B. +! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic +! Conservation Laws", +! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61 +! +!=============================================================================== +! + subroutine riemann_mhd_adi_hll(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : ivx + use equations , only : prim2cons, fluxspeed + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, srml ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr real(kind=8), dimension(n) :: cl, cr ! @@ -2610,27 +3382,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -2700,10 +3451,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -2719,16 +3469,14 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_adi_hllc(n, h, q, f) + subroutine riemann_mhd_adi_hllc(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr - use equations , only : imx, imy, imz, ien - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -2737,19 +3485,18 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -2761,27 +3508,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -3004,10 +3730,9 @@ module schemes ! ! Arguments: ! -! n - the length of input vectors; -! h - the spatial step; -! q - the input array of primitive variables; -! f - the output array of fluxes; +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; ! ! References: ! @@ -3018,16 +3743,14 @@ module schemes ! !=============================================================================== ! - subroutine riemann_mhd_adi_hlld(n, h, q, f) + subroutine riemann_mhd_adi_hlld(n, ql, qr, f) ! include external procedures ! - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr - use equations , only : imx, imy, imz, ien - use equations , only : cmax - use equations , only : prim2cons, fluxspeed - use interpolations, only : reconstruct, fix_positivity + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed ! local variables are not implicit by default ! @@ -3036,20 +3759,19 @@ module schemes ! subroutine arguments ! integer , intent(in) :: n - real(kind=8) , intent(in) :: h - real(kind=8), dimension(nv,n), intent(in) :: q + real(kind=8), dimension(nv,n), intent(in) :: ql, qr real(kind=8), dimension(nv,n), intent(out) :: f ! local variables ! - integer :: p, i + integer :: i real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb, dv real(kind=8) :: dnl, dnr, cal, car, sml, smr ! local arrays to store the states ! - real(kind=8), dimension(nv,n) :: ql, qr, ul, ur, fl, fr + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr real(kind=8), dimension(nv) :: wl, wr, wcl, wcr, ui real(kind=8), dimension(n) :: cl, cr ! @@ -3061,27 +3783,6 @@ module schemes call start_timer(imr) #endif /* PROFILE */ -! reconstruct the left and right states of primitive variables -! - do p = 1, nv - call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) - end do - -! obtain the state values for Bx and Psi for the GLM-MHD equations -! - cl(:) = 0.5d+00 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax) - cr(:) = 0.5d+00 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) - ql(ibx,:) = cl(:) - qr(ibx,:) = cl(:) - ql(ibp,:) = cr(:) - qr(ibp,:) = cr(:) - -! check if the reconstruction doesn't give the negative density or pressure, -! if so, correct the states -! - call fix_positivity(n, q(idn,:), ql(idn,:), qr(idn,:)) - call fix_positivity(n, q(ipr,:), ql(ipr,:), qr(ipr,:)) - ! calculate corresponding conserved variables of the left and right states ! call prim2cons(n, ql(:,:), ul(:,:)) @@ -3592,6 +4293,196 @@ module schemes !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_adi_hlld +! +!=============================================================================== +! +! subroutine RIEMANN_MHD_ADI_ROE: +! ------------------------------ +! +! Subroutine solves one dimensional Riemann problem using +! the Roe's method. +! +! Arguments: +! +! n - the length of input vectors; +! ql, qr - the array of primitive variables at the Riemann states; +! f - the output array of fluxes; +! +! 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] Toro, E. F., +! "Riemann Solvers and Numerical Methods for Fluid dynamics", +! Springer-Verlag Berlin Heidelberg, 2009 +! +!=============================================================================== +! + subroutine riemann_mhd_adi_roe(n, ql, qr, f) + +! include external procedures +! + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use equations , only : imx, imy, imz, ien + use equations , only : prim2cons, fluxspeed, eigensystem_roe + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer , intent(in) :: n + real(kind=8), dimension(nv,n), intent(in) :: ql, qr + real(kind=8), dimension(nv,n), intent(out) :: f + +! local variables +! + integer :: p, i + real(kind=8) :: sdl, sdr, sds + real(kind=8) :: pml, pmr + real(kind=8) :: xx, yy + +! local arrays to store the states +! + real(kind=8), dimension(nv,n) :: ul, ur, fl, fr + real(kind=8), dimension(n) :: cl, cr + real(kind=8), dimension(nv) :: qi, ci, al + real(kind=8), dimension(nv,nv) :: li, ri +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for Riemann solver +! + call start_timer(imr) +#endif /* PROFILE */ + +! calculate corresponding conserved variables of the left and right states +! + call prim2cons(n, ql(:,:), ul(:,:)) + call prim2cons(n, qr(:,:), ur(:,:)) + +! calculate the physical fluxes and speeds at the states +! + call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) + call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) + +! iterate over all points +! + do i = 1, n + +! calculate variables for the Roe intermediate state +! + sdl = sqrt(ql(idn,i)) + sdr = sqrt(qr(idn,i)) + sds = sdl + sdr + +! prepare magnetic pressures +! + pml = 0.5d+00 * sum(ql(ibx:ibz,i)**2) + pmr = 0.5d+00 * sum(qr(ibx:ibz,i)**2) + +! prepare the Roe intermediate state vector (eq. 11.60 in [2]) +! + qi(idn) = sdl * sdr + qi(ivx) = (sdl * ql(ivx,i) + sdr * qr(ivx,i)) / sds + qi(ivy) = (sdl * ql(ivy,i) + sdr * qr(ivy,i)) / sds + qi(ivz) = (sdl * ql(ivz,i) + sdr * qr(ivz,i)) / sds + qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl & + + (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds + qi(ibx) = ql(ibx,i) + qi(iby) = (sdr * ql(iby,i) + sdl * qr(iby,i)) / sds + qi(ibz) = (sdr * ql(ibz,i) + sdl * qr(ibz,i)) / sds + qi(ibp) = ql(ibp,i) + +! prepare coefficients +! + xx = 0.5d+00 * ((ql(iby,i) - qr(iby,i))**2 & + + (ql(ibz,i) - qr(ibz,i))**2) / sds**2 + yy = 0.5d+00 * (ql(idn,i) + qr(idn,i)) / qi(idn) + +! obtain eigenvalues and eigenvectors +! + call eigensystem_roe(xx, yy, qi(:), ci(:), ri(:,:), li(:,:)) + +! return upwind fluxes +! + if (ci(1) >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (ci(nv) <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else + +! prepare fluxes which do not change across the states +! + f(ibx,i) = fl(ibx,i) + f(ibp,i) = fl(ibp,i) + +! calculate wave amplitudes α = L.ΔU +! + al(1:nv) = 0.0d+00 + do p = 1, nv + al(1:nv) = al(1:nv) + li(p,1:nv) * (ur(p,i) - ul(p,i)) + end do + +! calculate the flux average +! + f(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) + f(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) + f(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) + f(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) + f(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i)) + f(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) + f(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + +! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) +! + if (qi(ivx) >= 0.0d+00) then + do p = 1, nv + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(ien,i) = f(ien,i) + xx * ri(p,ien) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + else + do p = nv, 1, -1 + xx = - 0.5d+00 * al(p) * abs(ci(p)) + + f(idn,i) = f(idn,i) + xx * ri(p,idn) + f(imx,i) = f(imx,i) + xx * ri(p,imx) + f(imy,i) = f(imy,i) + xx * ri(p,imy) + f(imz,i) = f(imz,i) + xx * ri(p,imz) + f(ien,i) = f(ien,i) + xx * ri(p,ien) + f(iby,i) = f(iby,i) + xx * ri(p,iby) + f(ibz,i) = f(ibz,i) + xx * ri(p,ibz) + end do + end if + + end if ! sl < 0 < sr + + end do ! i = 1, n + +#ifdef PROFILE +! stop accounting time for Riemann solver +! + call stop_timer(imr) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_adi_roe !=============================================================================== !