diff --git a/sources/equations.F90 b/sources/equations.F90 index 0fec063..c1f05d0 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -71,20 +71,10 @@ module equations real(kind=8), dimension(:,:) , intent(out) :: f real(kind=8), dimension(:,:), optional, intent(out) :: c end subroutine - function maxspeed_iface(qq) result(maxspeed) - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - real(kind=8) :: maxspeed - end function subroutine get_maximum_speeds_iface(qq, um, cm) real(kind=8), dimension(:,:,:,:), intent(in) :: qq real(kind=8) , intent(out):: um, cm end subroutine - subroutine esystem_roe_iface(x, y, q, c, r, l) - real(kind=8) , intent(in) :: x, y - real(kind=8), dimension(:) , intent(in) :: q - real(kind=8), dimension(:) , intent(inout) :: c - real(kind=8), dimension(:,:), intent(inout) :: l, r - end subroutine subroutine nr_iterate_iface(mm, bb, mb, en, dn, w, vv, info) real(kind=8), intent(in) :: mm, bb, mb, en, dn real(kind=8), intent(inout) :: w, vv @@ -101,17 +91,9 @@ module equations ! procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null() -! pointer to the maxspeed procedure -! - procedure(maxspeed_iface) , pointer, save :: maxspeed => null() - procedure(get_maximum_speeds_iface), pointer, save :: & get_maximum_speeds => null() -! pointer to the Roe eigensystem procedure -! - procedure(esystem_roe_iface), pointer, save :: eigensystem_roe => null() - ! pointer to the variable conversion method ! procedure(nr_iterate_iface) , pointer, save :: nr_iterate => null() @@ -230,8 +212,7 @@ module equations public :: prim2cons_mhd_adi, fluxspeed_mhd_adi public :: prim2cons_srhd_adi, fluxspeed_srhd_adi public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi - public :: maxspeed, reset_maxspeed, get_maximum_speeds - public :: eigensystem_roe + public :: get_maximum_speeds public :: update_primitive_variables public :: fix_unphysical_cells, correct_unphysical_states public :: adiabatic_index, relativistic, magnetized @@ -372,8 +353,6 @@ module equations prim2cons => prim2cons_hd_iso cons2prim => cons2prim_hd_iso fluxspeed => fluxspeed_hd_iso - maxspeed => maxspeed_hd_iso - eigensystem_roe => esystem_roe_hd_iso get_maximum_speeds => get_maximum_speeds_hd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -413,8 +392,6 @@ module equations prim2cons => prim2cons_hd_adi cons2prim => cons2prim_hd_adi fluxspeed => fluxspeed_hd_adi - maxspeed => maxspeed_hd_adi - eigensystem_roe => esystem_roe_hd_adi get_maximum_speeds => get_maximum_speeds_hd_adi ! warn about the unimplemented equation of state @@ -540,8 +517,6 @@ module equations prim2cons => prim2cons_mhd_iso cons2prim => cons2prim_mhd_iso fluxspeed => fluxspeed_mhd_iso - maxspeed => maxspeed_mhd_iso - eigensystem_roe => esystem_roe_mhd_iso get_maximum_speeds => get_maximum_speeds_mhd_iso case("adi", "ADI", "adiabatic", "ADIABATIC") @@ -588,8 +563,6 @@ module equations prim2cons => prim2cons_mhd_adi cons2prim => cons2prim_mhd_adi fluxspeed => fluxspeed_mhd_adi - maxspeed => maxspeed_mhd_adi - eigensystem_roe => esystem_roe_mhd_adi get_maximum_speeds => get_maximum_speeds_mhd_adi case default @@ -718,7 +691,6 @@ module equations prim2cons => prim2cons_srhd_adi cons2prim => cons2prim_srhd_adi fluxspeed => fluxspeed_srhd_adi - maxspeed => maxspeed_srhd_adi get_maximum_speeds => get_maximum_speeds_srhd_adi ! warn about the unimplemented equation of state @@ -900,7 +872,6 @@ module equations prim2cons => prim2cons_srmhd_adi cons2prim => cons2prim_srmhd_adi fluxspeed => fluxspeed_srmhd_adi - maxspeed => maxspeed_srmhd_adi get_maximum_speeds => get_maximum_speeds_srmhd_adi ! warn about the unimplemented equation of state @@ -1200,8 +1171,6 @@ module equations nullify(prim2cons) nullify(cons2prim) nullify(fluxspeed) - nullify(maxspeed) - nullify(eigensystem_roe) ! deallocate variable indices ! @@ -1299,62 +1268,6 @@ module equations ! !=============================================================================== ! -! subroutine RESET_MAXSPEED: -! ------------------------- -! -! Subroutine resets the maximum speed in the domain to zero. -! -! -!=============================================================================== -! - subroutine reset_maxspeed() - -! local variables are not implicit by default -! - implicit none -! -!------------------------------------------------------------------------------- -! -! reset the maximum speed -! - cmax = 0.0d+00 - -!------------------------------------------------------------------------------- -! - end subroutine reset_maxspeed -! -!=============================================================================== -! -! function GET_MAXSPEED: -! ----------------- -! -! Function returns the maximum speed in the domain. -! -! -!=============================================================================== -! - real(kind=8) function get_maxspeed() - -! local variables are not implicit by default -! - implicit none -! -!------------------------------------------------------------------------------- -! -! return the maximum speed -! - get_maxspeed = cmax - -! return the value -! - return - -!------------------------------------------------------------------------------- -! - end function get_maxspeed -! -!=============================================================================== -! ! subroutine UPDATE_PRIMITIVE_VARIABLES: ! ------------------------------------- ! @@ -1841,67 +1754,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_HD_ISO: -! ------------------------ -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! q - the array of primitive variables; -! -! -!=============================================================================== -! - function maxspeed_hd_iso(qq) result(maxspeed) - - use coordinates, only : nb, ne - - implicit none - - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - - real(kind=8) :: maxspeed - - integer :: i, j, k - real(kind=8) :: vv, v - -!------------------------------------------------------------------------------- -! - maxspeed = 0.0d+00 - -#if NDIMS == 2 - k = 1 -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - do k = nb, ne -#endif /* NDIMS == 3 */ - do j = nb, ne - do i = nb, ne - -! calculate the velocity amplitude -! - vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k)) - v = sqrt(vv) - -! calculate the maximum speed -! - maxspeed = max(maxspeed, v + csnd) - - end do ! i = nb, ne - end do ! j = nb, ne -#if NDIMS == 3 - end do ! k = nb, ne -#endif /* NDIMS == 3 */ - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_hd_iso -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_HD_ISO: ! ------------------------------------ ! @@ -1955,103 +1807,6 @@ module equations ! end subroutine get_maximum_speeds_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) - - implicit none - - real(kind=8) , intent(in) :: x, y - real(kind=8), dimension(:) , intent(in) :: q - real(kind=8), dimension(:) , intent(inout) :: c - real(kind=8), dimension(:,:), intent(inout) :: l, r - - logical , save :: first = .true. - real(kind=8), dimension(4,4), save :: lvec, rvec -!$omp threadprivate(first, lvec, rvec) - - real(kind=8) :: vc - -!------------------------------------------------------------------------------- -! - if (first) then - lvec(:,:) = 0.0d+00 - rvec(:,:) = 0.0d+00 - - lvec(ivx,1) = - 5.0d-01 / csnd - lvec(ivy,2) = 1.0d+00 - lvec(ivz,3) = 1.0d+00 - lvec(ivx,4) = - lvec(ivx,1) - - rvec(1,idn) = 1.0d+00 - rvec(2,ivy) = 1.0d+00 - rvec(3,ivz) = 1.0d+00 - rvec(4,idn) = 1.0d+00 - - first = .false. - end if - -! 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 left eigenvectors matrix -! - vc = 5.0d-01 * q(ivx) / csnd - lvec(idn,1) = 5.0d-01 + vc - lvec(idn,2) = - q(ivy) - lvec(idn,3) = - q(ivz) - lvec(idn,4) = 5.0d-01 - vc - -! update the varying elements of the right eigenvectors matrix -! - rvec(1,ivx) = c(1) - rvec(1,ivy) = q(ivy) - rvec(1,ivz) = q(ivz) - - rvec(4,ivx) = c(4) - rvec(4,ivy) = q(ivy) - rvec(4,ivz) = q(ivz) - -! copy matrices of eigenvectors -! - l(:,:) = lvec(:,:) - r(:,:) = rvec(:,:) - -!------------------------------------------------------------------------------- -! - end subroutine esystem_roe_hd_iso -! !******************************************************************************* ! ! ADIABATIC HYDRODYNAMIC EQUATIONS @@ -2257,71 +2012,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_HD_ADI: -! ------------------------ -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! q - the array of primitive variables; -! -! -!=============================================================================== -! - function maxspeed_hd_adi(qq) result(maxspeed) - - use coordinates, only : nb, ne - - implicit none - - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - - real(kind=8) :: maxspeed - - integer :: i, j, k - real(kind=8) :: vv, v, c - -!------------------------------------------------------------------------------- -! - maxspeed = 0.0d+00 - -#if NDIMS == 2 - k = 1 -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - do k = nb, ne -#endif /* NDIMS == 3 */ - do j = nb, ne - do i = nb, ne - -! calculate the velocity amplitude -! - vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k)) - v = sqrt(vv) - -! calculate the adiabatic speed of sound -! - c = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)) - -! calculate the maximum speed -! - maxspeed = max(maxspeed, v + c) - - end do ! i = nb, ne - end do ! j = nb, ne -#if NDIMS == 3 - end do ! k = nb, ne -#endif /* NDIMS == 3 */ - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_hd_adi -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_HD_ADI: ! ------------------------------------ ! @@ -2377,144 +2067,6 @@ module equations ! end subroutine get_maximum_speeds_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) - - implicit none - - real(kind=8) , intent(in) :: x, y - real(kind=8), dimension(:) , intent(in) :: q - real(kind=8), dimension(:) , intent(inout) :: c - real(kind=8), dimension(:,:), intent(inout) :: l, r - - logical , save :: first = .true. - real(kind=8), dimension(5,5), save :: lvec, rvec -!$omp threadprivate(first, lvec, rvec) - - real(kind=8) :: vv, vh, c2, cc, vc, fc, fh, f1, f2, fv, fx - -!------------------------------------------------------------------------------- -! - if (first) then - lvec(:,:) = 0.0d+00 - rvec(:,:) = 0.0d+00 - - lvec(ivy,2) = 1.0d+00 - lvec(ivz,3) = 1.0d+00 - - rvec(1,idn) = 1.0d+00 - rvec(2,ivy) = 1.0d+00 - rvec(3,ivz) = 1.0d+00 - rvec(4,idn) = 1.0d+00 - rvec(5,idn) = 1.0d+00 - - first = .false. - end if - -! calculate characteristic speeds and useful variables -! - vv = sum(q(ivx:ivz)**2) - vh = 5.0d-01 * vv - c2 = gammam1 * (q(ien) - vh) - cc = sqrt(c2) - vc = q(ivx) * cc - fc = gammam1 / c2 - fh = 5.0d-01 * fc - f1 = 5.0d-01 * vc / c2 - f2 = 5.0d-01 * cc / c2 - fv = fh * vh - fx = fh * q(ivx) - -! 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 left eigenvectors matrix -! - lvec(idn,1) = fv + f1 - lvec(ivx,1) = - fx - f2 - lvec(ivy,1) = - fh * q(ivy) - lvec(ivz,1) = - fh * q(ivz) - lvec(ien,1) = fh - - lvec(idn,2) = - q(ivy) - - lvec(idn,3) = - q(ivz) - - lvec(idn,4) = 1.0d+00 - fh * vv - lvec(ivx,4) = fc * q(ivx) - lvec(ivy,4) = fc * q(ivy) - lvec(ivz,4) = fc * q(ivz) - lvec(ien,4) = - fc - - lvec(idn,5) = fv - f1 - lvec(ivx,5) = - fx + f2 - lvec(ivy,5) = lvec(ivy,1) - lvec(ivz,5) = lvec(ivz,1) - lvec(ien,5) = fh - -! update the varying elements of the right eigenvectors matrix -! - rvec(1,ivx) = q(ivx) - cc - rvec(1,ivy) = q(ivy) - rvec(1,ivz) = q(ivz) - rvec(1,ien) = q(ien) - vc - - rvec(2,ien) = q(ivy) - - rvec(3,ien) = q(ivz) - - rvec(4,ivx) = q(ivx) - rvec(4,ivy) = q(ivy) - rvec(4,ivz) = q(ivz) - rvec(4,ien) = vh - - rvec(5,ivx) = q(ivx) + cc - rvec(5,ivy) = q(ivy) - rvec(5,ivz) = q(ivz) - rvec(5,ien) = q(ien) + vc - -! copy matrices of eigenvectors -! - l(:,:) = lvec(:,:) - r(:,:) = rvec(:,:) - -!------------------------------------------------------------------------------- -! - end subroutine esystem_roe_hd_adi -! !******************************************************************************* ! ! ISOTHERMAL MAGNETOHYDRODYNAMIC EQUATIONS @@ -2732,71 +2284,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_MHD_ISO: -! ------------------------- -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! q - the array of primitive variables; -! -!=============================================================================== -! - function maxspeed_mhd_iso(qq) result(maxspeed) - - use coordinates, only : nb, ne - - implicit none - - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - - real(kind=8) :: maxspeed - - integer :: i, j, k - real(kind=8) :: vv, bb, v, c - -!------------------------------------------------------------------------------- -! - maxspeed = 0.0d+00 - -#if NDIMS == 2 - k = 1 -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - do k = nb, ne -#endif /* NDIMS == 3 */ - do j = nb, ne - do i = nb, ne - -! calculate the velocity amplitude -! - vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k)) - v = sqrt(vv) - bb = sum(qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k)) - -! calculate the fast magnetosonic speed -! - c = sqrt(csnd2 + bb / qq(idn,i,j,k)) - -! calculate the maximum of speed -! - maxspeed = max(maxspeed, v + c) - - end do ! i = nb, ne - end do ! j = nb, ne -#if NDIMS == 3 - end do ! k = nb, ne -#endif /* NDIMS == 3 */ - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_mhd_iso -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_MHD_ISO: ! ------------------------------------- ! @@ -2857,261 +2344,6 @@ module equations ! end subroutine get_maximum_speeds_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) - - implicit none - - real(kind=8) , intent(in) :: x, y - real(kind=8), dimension(:) , intent(in) :: q - real(kind=8), dimension(:) , intent(inout) :: c - real(kind=8), dimension(:,:), intent(inout) :: l, r - - logical , save :: first = .true. - real(kind=8), dimension(8,8), save :: lvec, rvec -!$omp threadprivate(first, lvec, rvec) - - real(kind=8) :: ca, ca2, ct2, cf, cf2, cs, cs2, cf2_cs2 - real(kind=8) :: br, br2, brs, br2s, tsum, tdif, twid_c2, twid_c - real(kind=8) :: bty, btz, btys, btzs, bt2s, alf, als, norm - real(kind=8) :: sqrtd, sgn, qf, qs, af_prime, as_prime - real(kind=8) :: cff, css, af, as, afpb, aspb, q2s, q3s, vqstr - -!------------------------------------------------------------------------------- -! - if (first) then - lvec(:,:) = 0.0d+00 - rvec(:,:) = 0.0d+00 - - first = .false. - end if - -! coefficients for eigenvalues -! - ca2 = q(ibx) * q(ibx) / q(idn) - ca = sqrt(ca2) - br2 = sum(q(iby:ibz)**2) - br2s = br2 * y - ct2 = br2s / q(idn) - - twid_c2 = csnd2 + x - tsum = ca2 + ct2 + twid_c2 - tdif = ca2 + ct2 - twid_c2 - cf2_cs2 = sqrt(tdif * tdif + 4.0d+00 * twid_c2 * ct2) - cf2 = 0.5d+00 * (tsum + cf2_cs2) - cf = sqrt(cf2) - cs2 = twid_c2 * ca2 / cf2 - cs = sqrt(cs2) - -! 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 .or. c(7) <= 0.0d+00) return - -! remaining coefficients -! - br = sqrt(br2) - brs = sqrt(br2s) - if (abs(br) > 0.0d+00) then - bty = q(iby) / br - btz = q(ibz) / br - else - bty = 1.0d+00 - btz = 0.0d+00 - end if - btys = bty / sqrt(y) - btzs = btz / sqrt(y) - bt2s = btys * btys + btzs * btzs - - if (.not. abs(cf2 - cs2) > 0.0d+00) then - alf = 1.0d+00 - als = 0.0d+00 - else if ((twid_c2 - cs2) <= 0.0d+00) then - alf = 0.0d+00 - als = 1.0d+00 - else if ((cf2 - twid_c2) <= 0.0d+00) then - alf = 1.0d+00 - als = 0.0d+00 - else - alf = sqrt((twid_c2 - cs2) / (cf2 - cs2)) - als = sqrt((cf2 - twid_c2) / (cf2 - cs2)) - end if - - sqrtd = sqrt(q(idn)) - sgn = sign(1.0d+00, q(ibx)) - twid_c = sqrt(twid_c2) - qf = cf * alf * sgn - qs = cs * als * sgn - af_prime = twid_c * alf / sqrtd - as_prime = twid_c * als / sqrtd - -! === update the varying elements of the right eigenvectors matrix -! -! left-going fast wave -! - rvec(1,idn) = alf - rvec(1,ivx) = alf * c(1) - rvec(1,ivy) = alf * q(ivy) + qs * btys - rvec(1,ivz) = alf * q(ivz) + qs * btzs - rvec(1,iby) = as_prime * btys - rvec(1,ibz) = as_prime * btzs - -! left-going Alfvèn wave -! - rvec(2,ivy) = - btz - rvec(2,ivz) = bty - rvec(2,iby) = - btz * sgn / sqrtd - rvec(2,ibz) = bty * sgn / sqrtd - -! left-going slow wave -! - rvec(3,idn) = als - rvec(3,ivx) = als * c(3) - rvec(3,ivy) = als * q(ivy) - qf * btys - rvec(3,ivz) = als * q(ivz) - qf * btzs - rvec(3,iby) = - af_prime * btys - rvec(3,ibz) = - af_prime * btzs - -! right-going slow wave -! - rvec(5,idn) = als - rvec(5,ivx) = als * c(5) - rvec(5,ivy) = als * q(ivy) + qf * btys - rvec(5,ivz) = als * q(ivz) + qf * btzs - rvec(5,iby) = rvec(3,iby) - rvec(5,ibz) = rvec(3,ibz) - -! right-going Alfvèn wave -! - rvec(6,ivy) = btz - rvec(6,ivz) = - bty - rvec(6,iby) = rvec(2,iby) - rvec(6,ibz) = rvec(2,ibz) - -! right-going fast wave -! - rvec(7,idn) = alf - rvec(7,ivx) = alf * c(7) - rvec(7,ivy) = alf * q(ivy) - qs * btys - rvec(7,ivz) = alf * q(ivz) - qs * btzs - rvec(7,iby) = rvec(1,iby) - rvec(7,ibz) = rvec(1,ibz) - -! === update the varying elements of the left eigenvectors matrix -! - norm = 2.0d+00 * twid_c2 - cff = alf * cf / norm - css = als * cs / norm - qf = qf / norm - qs = qs / norm - af = af_prime * q(idn) / norm - as = as_prime * q(idn) / norm - afpb = af_prime * brs / norm - aspb = as_prime * brs / norm - - q2s = btys / bt2s - q3s = btzs / bt2s - vqstr = q(ivy) * q2s + q(ivz) * q3s - -! left-going fast wave -! - lvec(idn,1) = cff * c(7) - qs * vqstr - aspb - lvec(ivx,1) = - cff - lvec(ivy,1) = qs * q2s - lvec(ivz,1) = qs * q3s - lvec(iby,1) = as * q2s - lvec(ibz,1) = as * q3s - -! left-going Alfvèn wave -! - lvec(idn,2) = 0.5d+00 * (q(ivy) * btz - q(ivz) * bty) - lvec(ivy,2) = - 0.5d+00 * btz - lvec(ivz,2) = 0.5d+00 * bty - lvec(iby,2) = - 0.5d+00 * sqrtd * btz * sgn - lvec(ibz,2) = 0.5d+00 * sqrtd * bty * sgn - -! left-going slow wave -! - lvec(idn,3) = css * c(5) + qf * vqstr + afpb - lvec(ivx,3) = - css - lvec(ivy,3) = - qf * q2s - lvec(ivz,3) = - qf * q3s - lvec(iby,3) = - af * q2s - lvec(ibz,3) = - af * q3s - -! right-going slow wave -! - lvec(idn,5) = - css * c(3) - qf * vqstr + afpb - lvec(ivx,5) = css - lvec(ivy,5) = - lvec(ivy,3) - lvec(ivz,5) = - lvec(ivz,3) - lvec(iby,5) = lvec(iby,3) - lvec(ibz,5) = lvec(ibz,3) - -! right-going Alfvèn wave -! - lvec(idn,6) = - lvec(idn,2) - lvec(ivy,6) = - lvec(ivy,2) - lvec(ivz,6) = - lvec(ivz,2) - lvec(iby,6) = lvec(iby,2) - lvec(ibz,6) = lvec(ibz,2) - -! right-going fast wave -! - lvec(idn,7) = - cff * c(1) + qs * vqstr - aspb - lvec(ivx,7) = cff - lvec(ivy,7) = - lvec(ivy,1) - lvec(ivz,7) = - lvec(ivz,1) - lvec(iby,7) = lvec(iby,1) - lvec(ibz,7) = lvec(ibz,1) - -! copy matrices of eigenvectors -! - l(:,:) = lvec(:,:) - r(:,:) = rvec(:,:) - -!------------------------------------------------------------------------------- -! - end subroutine esystem_roe_mhd_iso -! !******************************************************************************* ! ! ADIABATIC MAGNETOHYDRODYNAMIC EQUATIONS @@ -3339,71 +2571,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_MHD_ADI: -! ------------------------- -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! q - the array of primitive variables; -! -!=============================================================================== -! - function maxspeed_mhd_adi(qq) result(maxspeed) - - use coordinates, only : nb, ne - - implicit none - - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - - real(kind=8) :: maxspeed - - integer :: i, j, k - real(kind=8) :: vv, bb, v, c - -!------------------------------------------------------------------------------- -! - maxspeed = 0.0d+00 - -#if NDIMS == 2 - k = 1 -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - do k = nb, ne -#endif /* NDIMS == 3 */ - do j = nb, ne - do i = nb, ne - -! calculate the velocity amplitude -! - vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k)) - v = sqrt(vv) - bb = sum(qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k)) - -! calculate the fast magnetosonic speed -! - c = sqrt((adiabatic_index * qq(ipr,i,j,k) + bb) / qq(idn,i,j,k)) - -! calculate the maximum of speed -! - maxspeed = max(maxspeed, v + c) - - end do ! i = nb, ne - end do ! j = nb, ne -#if NDIMS == 3 - end do ! k = nb, ne -#endif /* NDIMS == 3 */ - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_mhd_adi -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_MHD_ADI: ! ------------------------------------- ! @@ -3465,301 +2632,6 @@ module equations ! end subroutine get_maximum_speeds_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) - - implicit none - - real(kind=8) , intent(in) :: x, y - real(kind=8), dimension(:) , intent(in) :: q - real(kind=8), dimension(:) , intent(inout) :: c - real(kind=8), dimension(:,:), intent(inout) :: l, r - - logical , save :: first = .true. - real(kind=8) , save :: gammam2 - real(kind=8), dimension(9,9), save :: lvec, rvec -!$omp threadprivate(first, gammam2, lvec, rvec) - - real(kind=8) :: ca, ca2, cf, cf2, cs, cs2, cf2_cs2, ct2 - real(kind=8) :: v2, br2, br, br2s, brs, hp - real(kind=8) :: tsum, tdif, twid_a2, twid_a - real(kind=8) :: bty, btz, btys, btzs, bt2s, vbet - real(kind=8) :: alf, als, af_prime, as_prime - real(kind=8) :: sqrtd, sgn, qf, qs, afpbb, aspbb - real(kind=8) :: qa, qb, qc, qd, q2s, q3s - real(kind=8) :: norm, cff, css, af, as, afpb, aspb, vqstr - - real(kind=8), parameter :: eps = epsilon(1.0d+00) - -!------------------------------------------------------------------------------- -! - if (first) then - gammam2 = adiabatic_index - 2.0d+00 - - lvec(:, : ) = 0.0d+00 - rvec(:, : ) = 0.0d+00 - rvec(4,idn) = 1.0d+00 - - first = .false. - end if - -! coefficients -! - ca2 = q(ibx) * q(ibx) / q(idn) - ca = sqrt(ca2) - v2 = sum(q(ivx:ivz)**2) - br2 = sum(q(iby:ibz)**2) - br2s = (gammam1 - gammam2 * y) * br2 - hp = q(ien) - (ca2 + br2 / q(idn)) - twid_a2 = max(eps, (gammam1 * (hp - 0.5d+00 * v2) - gammam2 * x)) - ct2 = br2s / q(idn) - tsum = ca2 + ct2 + twid_a2 - tdif = ca2 + ct2 - twid_a2 - cf2_cs2 = sqrt((tdif * tdif + 4.0d+00 * twid_a2 * ct2)) - cf2 = 0.5d+00 * (tsum + cf2_cs2) - cf = sqrt(cf2) - cs2 = twid_a2 * ca2 / cf2 - cs = sqrt(cs2) - -! 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) - -! eigenvectors only for the case of waves propagating in both direction -! - if (c(1) >= 0.0d+00 .or. c(8) <= 0.0d+00) return - -! remaining coefficients -! - br = sqrt(br2) - brs = sqrt(br2s) - if (abs(br) > 0.0d+00) then - bty = q(iby) / br - btz = q(ibz) / br - else - bty = 1.0d+00 - btz = 0.0d+00 - end if - btys = bty / sqrt(gammam1 - gammam2 * y) - btzs = btz / sqrt(gammam1 - gammam2 * y) - bt2s = btys * btys + btzs * btzs - vbet = q(ivy) * btys + q(ivz) * btzs - - if ( .not. abs(cf2 - cs2) > 0.0d+00 ) then - alf = 1.0d+00 - als = 0.0d+00 - else if ( (twid_a2 - cs2) <= 0.0d+00 ) then - alf = 0.0d+00 - als = 1.0d+00 - else if ( (cf2 - twid_a2) <= 0.0d+00 ) then - alf = 1.0d+00 - als = 0.0d+00 - else - alf = sqrt((twid_a2 - cs2) / (cf2 - cs2)) - als = sqrt((cf2 - twid_a2) / (cf2 - cs2)) - end if - - sqrtd = sqrt(q(idn)) - sgn = sign(1.0d+00, q(ibx)) - twid_a = sqrt(twid_a2) - qf = cf * alf * sgn - qs = cs * als * sgn - af_prime = twid_a * alf / sqrtd - as_prime = twid_a * als / sqrtd - afpbb = af_prime * brs * bt2s - aspbb = as_prime * brs * bt2s - -! === update the varying elements of the right eigenvectors matrix -! - rvec(1,idn) = alf - rvec(3,idn) = als - rvec(6,idn) = als - rvec(8,idn) = alf - - rvec(1,ivx) = alf * c(1) - rvec(3,ivx) = als * c(3) - rvec(4,ivx) = q(ivx) - rvec(6,ivx) = als * c(6) - rvec(8,ivx) = alf * c(8) - - qa = alf * q(ivy) - qb = als * q(ivy) - qc = qs * btys - qd = qf * btys - - rvec(1,ivy) = qa + qc - rvec(2,ivy) = - btz - rvec(3,ivy) = qb - qd - rvec(4,ivy) = q(ivy) - rvec(6,ivy) = qb + qd - rvec(7,ivy) = btz - rvec(8,ivy) = qa - qc - - qa = alf * q(ivz) - qb = als * q(ivz) - qc = qs * btzs - qd = qf * btzs - - rvec(1,ivz) = qa + qc - rvec(2,ivz) = bty - rvec(3,ivz) = qb - qd - rvec(4,ivz) = q(ivz) - rvec(6,ivz) = qb + qd - rvec(7,ivz) = - bty - rvec(8,ivz) = qa - qc - - rvec(1,ipr) = alf * (hp - q(ivx) * cf) + qs * vbet + aspbb - rvec(2,ipr) = -(q(ivy) * btz - q(ivz) * bty) - rvec(3,ipr) = als * (hp - q(ivx) * cs) - qf * vbet - afpbb - rvec(4,ipr) = 0.5d+00 * v2 + gammam2 * x / gammam1 - rvec(6,ipr) = als * (hp + q(ivx) * cs) + qf * vbet - afpbb - rvec(7,ipr) = - rvec(2,ipr) - rvec(8,ipr) = alf * (hp + q(ivx) * cf) - qs * vbet + aspbb - - rvec(1,iby) = as_prime * btys - rvec(2,iby) = - btz * sgn / sqrtd - rvec(3,iby) = - af_prime * btys - rvec(6,iby) = rvec(3,iby) - rvec(7,iby) = rvec(2,iby) - rvec(8,iby) = rvec(1,iby) - - rvec(1,ibz) = as_prime * btzs - rvec(2,ibz) = bty * sgn / sqrtd - rvec(3,ibz) = - af_prime * btzs - rvec(6,ibz) = rvec(3,ibz) - rvec(7,ibz) = rvec(2,ibz) - rvec(8,ibz) = rvec(1,ibz) - -! === update the varying elements of the left eigenvectors matrix -! - norm = 2.0d+00 * twid_a2 - cff = alf * cf / norm - css = als * cs / norm - qf = qf / norm - qs = qs / norm - af = af_prime * q(idn) / norm - as = as_prime * q(idn) / norm - afpb = af_prime * brs / norm - aspb = as_prime * brs / norm - - norm = norm / gammam1 - alf = alf / norm - als = als / norm - q2s = btys / bt2s - q3s = btzs / bt2s - vqstr = (q(ivy) * q2s + q(ivz) * q3s) - -! left-going fast wave -! - lvec(idn,1) = alf * (v2 - hp) + cff * (cf + q(ivx)) - qs * vqstr - aspb - lvec(ipr,1) = alf - lvec(ivx,1) = - alf * q(ivx) - cff - lvec(ivy,1) = - alf * q(ivy) + qs * q2s - lvec(ivz,1) = - alf * q(ivz) + qs * q3s - lvec(iby,1) = as * q2s - alf * q(iby) - lvec(ibz,1) = as * q3s - alf * q(ibz) - -! left-going Alfvèn wave -! - lvec(idn,2) = 5.0d-01 * (q(ivy) * btz - q(ivz) * bty) - lvec(ivy,2) = - 5.0d-01 * btz - lvec(ivz,2) = 5.0d-01 * bty - lvec(iby,2) = - 5.0d-01 * sqrtd * btz * sgn - lvec(ibz,2) = 5.0d-01 * sqrtd * bty * sgn - -! left-going slow wave -! - lvec(idn,3) = als * (v2 - hp) + css * (cs + q(ivx)) + qf * vqstr + afpb - lvec(ipr,3) = als - lvec(ivx,3) = - als * q(ivx) - css - lvec(ivy,3) = - als * q(ivy) - qf * q2s - lvec(ivz,3) = - als * q(ivz) - qf * q3s - lvec(iby,3) = - af * q2s - als * q(iby) - lvec(ibz,3) = - af * q3s - als * q(ibz) - -! entropy wave -! - lvec(idn,4) = 1.0d+00 - (5.0d-01 * v2 - gammam2 * x / gammam1) / twid_a2 - lvec(ipr,4) = - 1.0d+00 / twid_a2 - lvec(ivx,4) = q(ivx) / twid_a2 - lvec(ivy,4) = q(ivy) / twid_a2 - lvec(ivz,4) = q(ivz) / twid_a2 - lvec(iby,4) = q(iby) / twid_a2 - lvec(ibz,4) = q(ibz) / twid_a2 - -! right-going slow wave -! - lvec(idn,6) = als * (v2 - hp) + css * (cs - q(ivx)) - qf * vqstr + afpb - lvec(ipr,6) = als - lvec(ivx,6) = - als * q(ivx) + css - lvec(ivy,6) = - als * q(ivy) + qf * q2s - lvec(ivz,6) = - als * q(ivz) + qf * q3s - lvec(iby,6) = lvec(iby,3) - lvec(ibz,6) = lvec(ibz,3) - -! right-going Alfvèn wave -! - lvec(idn,7) = - lvec(idn,2) - lvec(ivy,7) = - lvec(ivy,2) - lvec(ivz,7) = - lvec(ivz,2) - lvec(iby,7) = lvec(iby,2) - lvec(ibz,7) = lvec(ibz,2) - -! right-going fast wave -! - lvec(idn,8) = alf * (v2 - hp) + cff * (cf - q(ivx)) + qs * vqstr - aspb - lvec(ipr,8) = alf - lvec(ivx,8) = - alf * q(ivx) + cff - lvec(ivy,8) = - alf * q(ivy) - qs * q2s - lvec(ivz,8) = - alf * q(ivz) - qs * q3s - lvec(iby,8) = lvec(iby,1) - lvec(ibz,8) = lvec(ibz,1) - -! copy matrices of eigenvectors -! - l(:,:) = lvec(:,:) - r(:,:) = rvec(:,:) - -!------------------------------------------------------------------------------- -! - end subroutine esystem_roe_mhd_adi -! !******************************************************************************* ! ! ADIABATIC SPECIAL RELATIVITY HYDRODYNAMIC EQUATIONS @@ -4014,74 +2886,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_SRHD_ADI: -! -------------------------- -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! q - the array of primitive variables; -! -!=============================================================================== -! - function maxspeed_srhd_adi(qq) result(maxspeed) - - use coordinates, only : nb, ne - - implicit none - - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - - real(kind=8) :: maxspeed - - integer :: i, j, k - real(kind=8) :: vv, v, cc, ww, c2, ss, fc - -!------------------------------------------------------------------------------- -! - maxspeed = 0.0d+00 - -#if NDIMS == 2 - k = 1 -#endif /* NDIMS == 2 */ -#if NDIMS == 3 - do k = nb, ne -#endif /* NDIMS == 3 */ - do j = nb, ne - do i = nb, ne - -! calculate the velocity amplitude -! - vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k)) - v = sqrt(vv) - -! calculate the square of the sound speed -! - ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi - c2 = adiabatic_index * qq(ipr,i,j,k) / ww - ss = c2 * (1.0d+00 - vv) / (1.0d+00 - c2) - fc = 1.0d+00 + ss - cc = sqrt(ss * (fc - vv)) - -! calculate the maximum of speed -! - maxspeed = max(maxspeed, (v + cc) / fc) - - end do ! i = nb, ne - end do ! j = nb, ne -#if NDIMS == 3 - end do ! k = nb, ne -#endif /* NDIMS == 3 */ - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_srhd_adi -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_SRHD_ADI: ! -------------------------------------- ! @@ -5311,43 +4115,6 @@ module equations ! !=============================================================================== ! -! function MAXSPEED_SRMHD_ADI: -! --------------------------- -! -! Function scans the variable array and returns the maximum speed in within. -! -! Arguments: -! -! qq - the array of primitive variables; -! -!=============================================================================== -! - function maxspeed_srmhd_adi(qq) result(maxspeed) - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - real(kind=8), dimension(:,:,:,:), intent(in) :: qq - -! return value -! - real(kind=8) :: maxspeed -! -!------------------------------------------------------------------------------- -! - maxspeed = 1.0d+00 - - return - -!------------------------------------------------------------------------------- -! - end function maxspeed_srmhd_adi -! -!=============================================================================== -! ! subroutine GET_MAXIMUM_SPEEDS_SRMHD_ADI: ! --------------------------------------- ! diff --git a/sources/schemes.F90 b/sources/schemes.F90 index 62c9b70..6a80da6 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -2040,117 +2040,219 @@ module schemes subroutine riemann_mhd_iso_roe(ql, qr, fi) use coordinates, only : nn => bcells - use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ibx, iby, ibz, ibp - use equations , only : prim2cons, fluxspeed, eigensystem_roe + use equations , only : nf, csnd2 + use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & + ibx, iby, ibz, ibp + use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi - integer :: i, p - real(kind=8) :: sdl, sdr, sds - real(kind=8) :: xx, yy + integer :: i + real(kind=8) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2 + real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als + real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff + real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb + real(kind=8) :: f1, f2, f3, f4 real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr - real(kind=8), dimension( 2,nn) :: cl, cr - real(kind=8), dimension(nf ) :: qi, ci, al - real(kind=8), dimension(nf,nf) :: li, ri + real(kind=8), dimension(nf ) :: qi + real(kind=8), dimension(6 ) :: lm, al, df + + logical , save :: first = .true. + integer , dimension(6) , save :: ivs + real(kind=8), dimension(6,6), save :: rvec, lvec +!$omp threadprivate(first, ivs, rvec, lvec) !------------------------------------------------------------------------------- ! + if (first) then + + rvec(:,:) = 0.0d+00 + lvec(:,:) = 0.0d+00 + + ivs = [ idn, imx, imy, imz, iby, ibz ] + + first = .false. + end if + call prim2cons(ql, ul) call prim2cons(qr, ur) - call fluxspeed(ql, ul, fl, cl) - call fluxspeed(qr, ur, fr, cr) + call fluxspeed(ql, ul, fl) + call fluxspeed(qr, ur, fr) do i = 1, nn -! 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds + qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds + qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds qi(ibx) = ql(ibx,i) - qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby: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 - - fi(:,i) = fl(:,i) - - else if (ci(nf) <= 0.0d+00) then - - fi(:,i) = fr(:,i) + f1 = ql(iby,i) - qr(iby,i) + f2 = ql(ibz,i) - qr(ibz,i) + xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds) + yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn) + br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz) + br = sqrt(br2) + if (br2 > 0.0d+00) then + bty = qi(iby) / br + btz = qi(ibz) / br else + bty = 0.0d+00 + btz = 0.0d+00 + end if -! prepare fluxes which do not change across the states -! - fi(ibx,i) = fl(ibx,i) - fi(ibp,i) = fl(ibp,i) + cc2 = csnd2 + xx + ca2 = qi(ibx) * qi(ibx) / qi(idn) + f1 = br2 / qi(idn) + f2 = ca2 + f1 + f3 = f2 - cc2 + f4 = sqrt(f3 * f3 + 4.0d+00 * cc2 * f1) + cf2 = 5.0d-01 * (f2 + cc2 + f4) + cs2 = cc2 * ca2 / cf2 -! calculate wave amplitudes α = L.ΔU -! - al(:) = 0.0d+00 - do p = 1, nf - al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i)) - end do + if (cs2 < cc2 .and. cc2 < cf2) then + f1 = (cc2 - cs2) / (cf2 - cs2) + alf = sqrt(f1) + als = sqrt(1.0d+00 - f1) + else if (cc2 >= cf2) then + alf = 1.0d+00 + als = 0.0d+00 + else + alf = 0.0d+00 + als = 1.0d+00 + end if -! calculate the flux average -! - fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) - fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) - fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) - fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) - fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) - fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + cc = sqrt(cc2) + ca = sqrt(ca2) + cf = sign(sqrt(cf2), qi(ivx)) + cs = sign(sqrt(cs2), qi(ivx)) -! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) -! - if (qi(ivx) >= 0.0d+00) then - do p = 1, nf - xx = - 0.5d+00 * al(p) * abs(ci(p)) + lm(1) = qi(ivx) + cf + lm(2) = qi(ivx) + ca + lm(3) = qi(ivx) + cs + lm(4) = qi(ivx) - cs + lm(5) = qi(ivx) - ca + lm(6) = qi(ivx) - cf - fi(idn,i) = fi(idn,i) + xx * ri(p,idn) - fi(imx,i) = fi(imx,i) + xx * ri(p,imx) - fi(imy,i) = fi(imy,i) + xx * ri(p,imy) - fi(imz,i) = fi(imz,i) + xx * ri(p,imz) - fi(iby,i) = fi(iby,i) + xx * ri(p,iby) - fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz) - end do - else - do p = nf, 1, -1 - xx = - 0.5d+00 * al(p) * abs(ci(p)) + sqrtd = sqrt(qi(idn)) + f1 = sqrt(yy) + qf = cf * sign(alf, qi(ibx)) / f1 + qs = cs * sign(als, qi(ibx)) / f1 + f1 = cc / (f1 * sqrtd) + af_prime = alf * f1 + as_prime = als * f1 - fi(idn,i) = fi(idn,i) + xx * ri(p,idn) - fi(imx,i) = fi(imx,i) + xx * ri(p,imx) - fi(imy,i) = fi(imy,i) + xx * ri(p,imy) - fi(imz,i) = fi(imz,i) + xx * ri(p,imz) - fi(iby,i) = fi(iby,i) + xx * ri(p,iby) - fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz) - end do - end if + rvec(1,1) = alf + rvec(1,2) = alf * lm(1) + rvec(1,3) = alf * qi(ivy) - qs * bty + rvec(1,4) = alf * qi(ivz) - qs * btz + rvec(1,5) = as_prime * bty + rvec(1,6) = as_prime * btz - end if ! sl < 0 < sr + rvec(6,1) = alf + rvec(6,2) = alf * lm(6) + rvec(6,3) = alf * qi(ivy) + qs * bty + rvec(6,4) = alf * qi(ivz) + qs * btz + rvec(6,5) = rvec(1,5) + rvec(6,6) = rvec(1,6) + + rvec(3,1) = als + rvec(3,2) = als * lm(3) + rvec(3,3) = als * qi(ivy) + qf * bty + rvec(3,4) = als * qi(ivz) + qf * btz + rvec(3,5) = - af_prime * bty + rvec(3,6) = - af_prime * btz + + rvec(4,1) = als + rvec(4,2) = als * lm(4) + rvec(4,3) = als * qi(ivy) - qf * bty + rvec(4,4) = als * qi(ivz) - qf * btz + rvec(4,5) = rvec(3,5) + rvec(4,6) = rvec(3,6) + + f1 = sign(1.0d+00 / sqrtd, qi(ibx)) + rvec(2,3) = btz + rvec(2,4) = - bty + rvec(2,5) = - btz * f1 + rvec(2,6) = bty * f1 + + rvec(5,3) = - rvec(2,3) + rvec(5,4) = - rvec(2,4) + rvec(5,5) = rvec(2,5) + rvec(5,6) = rvec(2,6) + + norm = 2.0d+00 * cc2 + cff = alf * cf / norm + css = als * cs / norm + norm = norm / yy + qf = qf / norm + qs = qs / norm + f1 = qi(idn) / norm + af = af_prime * f1 + as = as_prime * f1 + f1 = br / norm + afpb = af_prime * f1 + aspb = as_prime * f1 + vqstr = (qi(ivy) * bty + qi(ivz) * btz) + + f1 = qs * vqstr + lvec(1,1) = - cff * lm(6) + f1 - aspb + lvec(1,2) = cff + lvec(1,3) = - qs * bty + lvec(1,4) = - qs * btz + lvec(1,5) = as * bty + lvec(1,6) = as * btz + + lvec(6,1) = cff * lm(1) - f1 - aspb + lvec(6,2) = - lvec(1,2) + lvec(6,3) = - lvec(1,3) + lvec(6,4) = - lvec(1,4) + lvec(6,5) = lvec(1,5) + lvec(6,6) = lvec(1,6) + + f1 = qf * vqstr + lvec(3,1) = - css * lm(4) - f1 + afpb + lvec(3,2) = css + lvec(3,3) = qf * bty + lvec(3,4) = qf * btz + lvec(3,5) = - af * bty + lvec(3,6) = - af * btz + + lvec(4,1) = css * lm(3) + f1 + afpb + lvec(4,2) = - lvec(3,2) + lvec(4,3) = - lvec(3,3) + lvec(4,4) = - lvec(3,4) + lvec(4,5) = lvec(3,5) + lvec(4,6) = lvec(3,6) + + f1 = sign(5.0d-01 * sqrtd, qi(ibx)) + lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty) + lvec(2,3) = 5.0d-01 * btz + lvec(2,4) = - 5.0d-01 * bty + lvec(2,5) = - f1 * btz + lvec(2,6) = f1 * bty + + lvec(5,1) = - lvec(2,1) + lvec(5,3) = - lvec(2,3) + lvec(5,4) = - lvec(2,4) + lvec(5,5) = lvec(2,5) + lvec(5,6) = lvec(2,6) + + al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i)) + df = matmul(al, rvec) + fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:)) + fi(ibx,i) = fl(ibx,i) + fi(ibp,i) = fl(ibp,i) end do @@ -2204,7 +2306,7 @@ module schemes real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3 - real(kind=8) :: alf, als, f1, f2, f3 + real(kind=8) :: alf, als, f1, f2, f3, f4 real(kind=8), dimension(nf) :: v, lm, tm @@ -2266,34 +2368,23 @@ module schemes x3 = 0.0d+00 end if - ca2 = bxa * bxa - f1 = (ca2 + bp2) / dnl + csnd2 - ca2 = ca2 / dnl - f2 = f1 * f1 - f3 = 4.0d+00 * csnd2 * ca2 - if (f2 > f3) then - f2 = sqrt(f2 - f3) - cf2 = 5.0d-01 * (f1 + f2) - cs2 = 5.0d-01 * (f1 - f2) - else - cf2 = 5.0d-01 * f1 - cs2 = cf2 - end if + ca2 = bxa * bxa / dnl + f1 = bp2 / dnl + f2 = ca2 + f1 + f3 = f2 - csnd2 + f4 = sqrt(f3 * f3 + 4.0d+00 * csnd2 * f1) + cf2 = 5.0d-01 * (f2 + csnd2 + f4) + cs2 = csnd2 * ca2 / cf2 - if (cf2 > cs2) then - f2 = cf2 - cs2 - if (csnd2 > cs2) then - alf = sqrt((csnd2 - cs2) / f2) - else - alf = 0.0d+00 - end if - if (cf2 > csnd2) then - als = sqrt((cf2 - csnd2) / f2) - else - als = 0.0d+00 - end if - else + if (cs2 < csnd2 .and. csnd2 < cf2) then + f1 = (csnd2 - cs2) / (cf2 - cs2) + alf = sqrt(f1) + als = sqrt(1.0d+00 - f1) + else if (csnd2 >= cf2) then alf = 1.0d+00 + als = 0.0d+00 + else + alf = 0.0d+00 als = 1.0d+00 end if @@ -3137,129 +3228,287 @@ module schemes subroutine riemann_mhd_adi_roe(ql, qr, fi) use coordinates, only : nn => bcells - use equations , only : nf, idn, ivx, ivz, imx, imy, imz, ipr, ien, & + use equations , only : nf, adiabatic_index + use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien, & ibx, iby, ibz, ibp - use equations , only : prim2cons, fluxspeed, eigensystem_roe + use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi - integer :: i, p + integer :: i real(kind=8) :: sdl, sdr, sds - real(kind=8) :: xx, yy - real(kind=8) :: pml, pmr + real(kind=8) :: xx, yy, pml, pmr + + real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs + real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty + real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm + real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb + real(kind=8) :: f1, f2, f3, f4, f5 real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr - real(kind=8), dimension( 2,nn) :: cl, cr - real(kind=8), dimension(nf ) :: qi, ci, al - real(kind=8), dimension(nf,nf) :: li, ri + real(kind=8), dimension(nf ) :: qi + real(kind=8), dimension(7 ) :: lm, al, df + + logical , save :: first = .true. + integer , dimension(7) , save :: ivs + real(kind=8) , save :: adi_m1, adi_m2, adi_m2d1 + real(kind=8), dimension(7,7), save :: rvec, lvec +!$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec) !------------------------------------------------------------------------------- ! + if (first) then + adi_m1 = adiabatic_index - 1.0d+00 + adi_m2 = adiabatic_index - 2.0d+00 + adi_m2d1 = adi_m2 / adi_m1 + + rvec(:,:) = 0.0d+00 + lvec(:,:) = 0.0d+00 + rvec(4,1) = 1.0d+00 + + ivs = [ idn, imx, imy, imz, ipr, iby, ibz ] + + first = .false. + end if + call prim2cons(ql, ul) call prim2cons(qr, ur) - call fluxspeed(ql, ul, fl, cl) - call fluxspeed(qr, ur, fr, cr) + call fluxspeed(ql, ul, fl) + call fluxspeed(qr, ur, fr) do i = 1, nn -! calculate variables for the Roe intermediate state -! + pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i)) + pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i)) 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:ivz) = sdl * ql(ivx:ivz,i) / sds + sdr * qr(ivx:ivz,i) / sds - qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl & + qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds + qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl & + (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds + qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds qi(ibx) = ql(ibx,i) - qi(iby:ibz) = sdr * ql(iby:ibz,i) / sds + sdl * qr(iby: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 - - fi(:,i) = fl(:,i) - - else if (ci(nf) <= 0.0d+00) then - - fi(:,i) = fr(:,i) + f1 = qr(iby,i) - ql(iby,i) + f2 = qr(ibz,i) - ql(ibz,i) + xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds) + yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn) + br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz) + if (br2 > 0.0d+00) then + br = sqrt(br2) + bty = qi(iby) / br + btz = qi(ibz) / br else + br = 0.0d+00 + bty = 0.0d+00 + btz = 0.0d+00 + end if -! prepare fluxes which do not change across the states -! - fi(ibx,i) = fl(ibx,i) - fi(ibp,i) = fl(ibp,i) + v2 = sum(qi(ivx:ivz) * qi(ivx:ivz)) + v2h = 5.0d-01 * v2 + ayy = adi_m1 - adi_m2 * yy + sqty = sqrt(ayy) + ca2 = qi(ibx) * qi(ibx) + hp = qi(ien) - (ca2 + br2) / qi(idn) + ca2 = ca2 / qi(idn) + f1 = adi_m1 * (hp - v2h) + f2 = adi_m2 * xx + if (f1 > f2) then + cc2 = f1 - f2 + cc = sqrt(cc2) + f1 = ayy * br2 / qi(idn) + f2 = ca2 + f1 + f3 = f2 + cc2 + f4 = f2 - cc2 + f5 = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1) + cf2 = 5.0d-01 * (f3 + f5) + cs2 = cc2 * ca2 / cf2 + else + cf2 = ca2 + ayy * br2 / qi(idn) + cc2 = 0.0d+00 + cc = 0.0d+00 + cs2 = 0.0d+00 + end if -! calculate wave amplitudes α = L.ΔU -! - al(:) = 0.0d+00 - do p = 1, nf - al(:) = al(:) + li(p,:) * (ur(p,i) - ul(p,i)) - end do + if (cs2 < cc2 .and. cc2 < cf2) then + f1 = (cc2 - cs2) / (cf2 - cs2) + alf = sqrt(f1) + als = sqrt(1.0d+00 - f1) + else if (cc2 >= cf2) then + alf = 1.0d+00 + als = 0.0d+00 + else + alf = 0.0d+00 + als = 1.0d+00 + end if -! calculate the flux average -! - fi(idn,i) = 0.5d+00 * (fl(idn,i) + fr(idn,i)) - fi(imx,i) = 0.5d+00 * (fl(imx,i) + fr(imx,i)) - fi(imy,i) = 0.5d+00 * (fl(imy,i) + fr(imy,i)) - fi(imz,i) = 0.5d+00 * (fl(imz,i) + fr(imz,i)) - fi(ien,i) = 0.5d+00 * (fl(ien,i) + fr(ien,i)) - fi(iby,i) = 0.5d+00 * (fl(iby,i) + fr(iby,i)) - fi(ibz,i) = 0.5d+00 * (fl(ibz,i) + fr(ibz,i)) + cf = sign(sqrt(cf2), qi(ivx)) + ca = sign(sqrt(ca2), qi(ivx)) + cs = sign(sqrt(cs2), qi(ivx)) -! correct the flux by adding the characteristic wave contribution: ∑(α|λ|K) -! - if (qi(ivx) >= 0.0d+00) then - do p = 1, nf - xx = - 0.5d+00 * al(p) * abs(ci(p)) + lm(1) = qi(ivx) + cf + lm(2) = qi(ivx) + ca + lm(3) = qi(ivx) + cs + lm(4) = qi(ivx) + lm(5) = qi(ivx) - cs + lm(6) = qi(ivx) - ca + lm(7) = qi(ivx) - cf - fi(idn,i) = fi(idn,i) + xx * ri(p,idn) - fi(imx,i) = fi(imx,i) + xx * ri(p,imx) - fi(imy,i) = fi(imy,i) + xx * ri(p,imy) - fi(imz,i) = fi(imz,i) + xx * ri(p,imz) - fi(ien,i) = fi(ien,i) + xx * ri(p,ien) - fi(iby,i) = fi(iby,i) + xx * ri(p,iby) - fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz) - end do - else - do p = nf, 1, -1 - xx = - 0.5d+00 * al(p) * abs(ci(p)) + vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty - fi(idn,i) = fi(idn,i) + xx * ri(p,idn) - fi(imx,i) = fi(imx,i) + xx * ri(p,imx) - fi(imy,i) = fi(imy,i) + xx * ri(p,imy) - fi(imz,i) = fi(imz,i) + xx * ri(p,imz) - fi(ien,i) = fi(ien,i) + xx * ri(p,ien) - fi(iby,i) = fi(iby,i) + xx * ri(p,iby) - fi(ibz,i) = fi(ibz,i) + xx * ri(p,ibz) - end do - end if + sqrtd = sqrt(qi(idn)) + qf = cf * sign(alf, qi(ibx)) + qs = cs * sign(als, qi(ibx)) + f1 = cc / sqrtd + af_prime = f1 * alf + as_prime = f1 * als + f1 = br / sqty + afpbb = af_prime * f1 + aspbb = as_prime * f1 - end if ! sl < 0 < sr + f1 = qs / sqty + f2 = qs * vbet + f3 = as_prime / sqty + rvec(1,1) = alf + rvec(1,2) = alf * lm(1) + rvec(1,3) = alf * qi(ivy) - f1 * bty + rvec(1,4) = alf * qi(ivz) - f1 * btz + rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb + rvec(1,6) = f3 * bty + rvec(1,7) = f3 * btz + + rvec(7,1) = alf + rvec(7,2) = alf * lm(7) + rvec(7,3) = alf * qi(ivy) + f1 * bty + rvec(7,4) = alf * qi(ivz) + f1 * btz + rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb + rvec(7,6) = rvec(1,6) + rvec(7,7) = rvec(1,7) + + f1 = qf / sqty + f2 = qf * vbet + f3 = - af_prime / sqty + rvec(3,1) = als + rvec(3,2) = als * lm(3) + rvec(3,3) = als * qi(ivy) + f1 * bty + rvec(3,4) = als * qi(ivz) + f1 * btz + rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb + rvec(3,6) = f3 * bty + rvec(3,7) = f3 * btz + + rvec(5,1) = als + rvec(5,2) = als * lm(5) + rvec(5,3) = als * qi(ivy) - f1 * bty + rvec(5,4) = als * qi(ivz) - f1 * btz + rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb + rvec(5,6) = rvec(3,6) + rvec(5,7) = rvec(3,7) + + f1 = sign(1.0d+00, qi(ivx)) + f2 = sign(1.0d+00 / sqrtd, qi(ibx)) + rvec(2,3) = f1 * btz + rvec(2,4) = - f1 * bty + rvec(2,5) = f1 * (qi(ivy) * btz - qi(ivz) * bty) + rvec(2,6) = - f2 * btz + rvec(2,7) = f2 * bty + + rvec(6,3) = - rvec(2,3) + rvec(6,4) = - rvec(2,4) + rvec(6,5) = - rvec(2,5) + rvec(6,6) = rvec(2,6) + rvec(6,7) = rvec(2,7) + + rvec(4,2) = qi(ivx) + rvec(4,3) = qi(ivy) + rvec(4,4) = qi(ivz) + rvec(4,5) = v2h + adi_m2d1 * xx + + norm = 2.0d+00 * cc2 + f1 = sqty * br / norm + afpb = af_prime * f1 + aspb = as_prime * f1 + sqty = sqty / norm + vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty + + f1 = adi_m1 * alf / norm + f2 = qs * sqty + f3 = as_prime * qi(idn) * sqty + f4 = alf * cf / norm + f5 = qs * vqstr + lvec(1,1) = f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb + lvec(1,5) = f1 + lvec(1,2) = - f1 * qi(ivx) + f4 + lvec(1,3) = - f1 * qi(ivy) - f2 * bty + lvec(1,4) = - f1 * qi(ivz) - f2 * btz + lvec(1,6) = f3 * bty - f1 * qi(iby) + lvec(1,7) = f3 * btz - f1 * qi(ibz) + + lvec(7,1) = f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb + lvec(7,5) = f1 + lvec(7,2) = - f1 * qi(ivx) - f4 + lvec(7,3) = - f1 * qi(ivy) + f2 * bty + lvec(7,4) = - f1 * qi(ivz) + f2 * btz + lvec(7,6) = lvec(1,6) + lvec(7,7) = lvec(1,7) + + f1 = adi_m1 * als / norm + f2 = qf * sqty + f3 = af_prime * qi(idn) * sqty + f4 = als * cs / norm + f5 = qf * vqstr + lvec(3,1) = f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb + lvec(3,5) = f1 + lvec(3,2) = - f1 * qi(ivx) + f4 + lvec(3,3) = - f1 * qi(ivy) + f2 * bty + lvec(3,4) = - f1 * qi(ivz) + f2 * btz + lvec(3,6) = - f3 * bty - f1 * qi(iby) + lvec(3,6) = - f3 * btz - f1 * qi(ibz) + + lvec(5,1) = f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb + lvec(5,5) = f1 + lvec(5,2) = - f1 * qi(ivx) - f4 + lvec(5,3) = - f1 * qi(ivy) - f2 * bty + lvec(5,4) = - f1 * qi(ivz) - f2 * btz + lvec(5,6) = lvec(3,6) + lvec(5,7) = lvec(3,7) + + f1 = sign(5.0d-01, qi(ivx)) * bty + f2 = sign(5.0d-01, qi(ivx)) * btz + f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx)) + lvec(2,1) = qi(ivz) * f1 - qi(ivy) * f2 + lvec(2,3) = f2 + lvec(2,4) = - f1 + lvec(2,6) = - f2 * f3 + lvec(2,7) = f1 * f3 + + lvec(6,1) = - lvec(2,1) + lvec(6,3) = - lvec(2,3) + lvec(6,4) = - lvec(2,4) + lvec(6,6) = lvec(2,6) + lvec(6,7) = lvec(2,7) + + f1 = 1.0d+00 / cc2 + lvec(4,1) = 1.0d+00 - f1 * (v2h - adi_m2d1 * xx) + lvec(4,5) = - f1 + lvec(4,2) = f1 * qi(ivx) + lvec(4,3) = f1 * qi(ivy) + lvec(4,4) = f1 * qi(ivz) + lvec(4,6) = f1 * qi(iby) + lvec(4,7) = f1 * qi(ibz) + + al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i)) + df = matmul(al, rvec) + fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:)) + fi(ibx,i) = fl(ibx,i) + fi(ibp,i) = fl(ibp,i) end do @@ -3293,7 +3542,7 @@ module schemes subroutine riemann_mhd_adi_kepes(ql, qr, fi) use coordinates, only : nn => bcells - use equations , only : nf, adiabatic_index, ch => cglm + use equations , only : nf, adiabatic_index, cglm use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & ibx, iby, ibz, ibp, ipr, ien @@ -3302,53 +3551,76 @@ module schemes real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi - logical , save :: first = .true. - real(kind=8) , save :: adi_m1, adi_m1x - real(kind=8), dimension(9,9), save :: rm -!$omp threadprivate(first, adi_m1, adi_m1x, rm) - integer :: i - real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa - real(kind=8) :: dnl, prl, pta, vva, br, bl, bp - real(kind=8) :: btl, bta, eka, ema, ub2, uba - real(kind=8) :: aa2, al2, ab2, bb2, bp2, aa, al, ab, b1, b2, b3, bb - real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3 - real(kind=8) :: alf2, als2, alf, als, wps, wms, wpf, wmf - real(kind=8) :: f1, f2, f3, f4, sgnb1, v2l, v2r, b2l, b2r + real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa, dnl, prl, pta + real(kind=8) :: v2l, v2r, b2l, b2r, bl, br, bta, btl, bp2, eka, ema, vva + real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, aa2, x2, x3, ub2, uba + real(kind=8) :: alf, als, das, daf, csq, wps, wms, wpf, wmf + real(kind=8) :: f1, f2, f3, f4, f5 real(kind=8), dimension(nf) :: v, lm, tm + logical , save :: first = .true. + real(kind=8) , save :: adi_m1, adi_m1x, adi_m1xi + real(kind=8), dimension(9,9), save :: rm +!$omp threadprivate(first, adi_m1, adi_m1x, adi_m1xi, rm) + !------------------------------------------------------------------------------- ! if (first) then - adi_m1 = adiabatic_index - 1.0d+00 - adi_m1x = adi_m1 / adiabatic_index + adi_m1 = adiabatic_index - 1.0d+00 + adi_m1x = adi_m1 / adiabatic_index + adi_m1xi = adiabatic_index / adi_m1 rm(:,:) = 0.0d+00 - rm(6,4) = 1.0d+00 - rm(9,4) = 1.0d+00 rm(1,5) = 1.0d+00 + rm(6,4) = 1.0d+00 rm(6,6) = 1.0d+00 - rm(9,6) = - 1.0d+00 first = .false. end if do i = 1, nn -! the difference in entropy variables between the states -! v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i)) v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i)) b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i)) b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i)) - bl = ql(idn,i) / ql(ipr,i) br = qr(idn,i) / qr(ipr,i) - v(idn) = (log(qr(idn,i)) - log(ql(idn,i))) & - + (log(br) - log(bl)) / adi_m1 & - - 5.0d-01 * (br * v2r - bl * v2l) + btl = lmean(bl, br) + bta = amean(bl, br) + dnl = lmean(ql(idn,i), qr(idn,i)) + dna = amean(ql(idn,i), qr(idn,i)) + vxa = amean(ql(ivx,i), qr(ivx,i)) + vya = amean(ql(ivy,i), qr(ivy,i)) + vza = amean(ql(ivz,i), qr(ivz,i)) + bxa = amean(ql(ibx,i), qr(ibx,i)) + bya = amean(ql(iby,i), qr(iby,i)) + bza = amean(ql(ibz,i), qr(ibz,i)) + bpa = amean(ql(ibp,i), qr(ibp,i)) + prl = lmean(ql(ipr,i), qr(ipr,i)) + pra = dna / bta + eka = 5.0d-01 * amean(v2l, v2r) + ema = 5.0d-01 * amean(b2l, b2r) + pta = pra + ema + vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i)) + ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r) + uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), & + sum(qr(ivx:ivz,i) * qr(ibx:ibz,i))) + + v(idn) = 5.0d-01 * (bl * v2l - br * v2r) + if (qr(idn,i) > ql(idn,i)) then + v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i)) + else if (ql(idn,i) > qr(idn,i)) then + v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i)) + end if + if (br > bl) then + v(idn) = v(idn) + log(br / bl) / adi_m1 + else if (bl > br) then + v(idn) = v(idn) - log(bl / br) / adi_m1 + end if v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i) v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i) v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i) @@ -3358,92 +3630,55 @@ module schemes v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i) v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i) -! the intermediate state averages -! - btl = lmean(bl, br) - bta = amean(bl, br) - - dna = amean(ql(idn,i), qr(idn,i)) - vxa = amean(ql(ivx,i), qr(ivx,i)) - vya = amean(ql(ivy,i), qr(ivy,i)) - vza = amean(ql(ivz,i), qr(ivz,i)) - bxa = amean(ql(ibx,i), qr(ibx,i)) - bya = amean(ql(iby,i), qr(iby,i)) - bza = amean(ql(ibz,i), qr(ibz,i)) - bpa = amean(ql(ibp,i), qr(ibp,i)) - pra = dna / bta - dnl = lmean(ql(idn,i), qr(idn,i)) - prl = lmean(ql(ipr,i), qr(ipr,i)) - eka = 5.0d-01 * amean(v2l, v2r) - ema = 5.0d-01 * amean(b2l, b2r) - pta = pra + ema - vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i)) - ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r) - uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), & - sum(qr(ivx:ivz,i) * qr(ibx:ibz,i))) - -! the eigenvalues and related parameters -! - aa2 = adiabatic_index * pra / dnl - aa = sqrt(aa2) - al2 = adiabatic_index * prl / dnl - al = sqrt(al2) - ab2 = adiabatic_index / bta - ab = sqrt(ab2) - - b1 = bxa / sqrt(dnl) - b2 = bya / sqrt(dnl) - b3 = bza / sqrt(dnl) - bp2 = b2**2 + b3**2 - bp = sqrt(bp2) - bb2 = b1**2 + bp2 - bb = sqrt(bb2) - sgnb1 = sign(1.0d+00, b1) - - if (bp > 0.0d+00) then - x2 = b2 / bp - x3 = b3 / bp + bp2 = bya * bya + bza * bza + if (bp2 > 0.0d+00) then + f1 = sqrt(bp2) + x2 = bya / f1 + x3 = bza / f1 else x2 = 0.0d+00 x3 = 0.0d+00 end if - ca2 = b1**2 - ca = sqrt(ca2) - f1 = aa2 + bb2 - f3 = max(0.0d+00, f1**2 - 4.0d+00 * aa2 * ca2) - f2 = sqrt(f3) - cf2 = 5.0d-01 * (f1 + f2) - cf = sqrt(cf2) - cs2 = 5.0d-01 * (f1 - f2) - cs = sqrt(cs2) + aa2 = adiabatic_index * pra / dnl + ca2 = bxa * bxa / dnl + f1 = bp2 / dnl + f2 = ca2 + f1 + f3 = f2 - aa2 + f4 = sqrt(f3 * f3 + 4.0d+00 * aa2 * f1) + cf2 = 5.0d-01 * (f2 + aa2 + f4) + cs2 = aa2 * ca2 / cf2 - f2 = cf2 - cs2 - if (f2 > 0.0d+00) then - alf2 = max(0.0d+00, aa2 - cs2) / f2 - als2 = max(0.0d+00, cf2 - aa2) / f2 + if (cs2 < aa2 .and. aa2 < cf2) then + f1 = (aa2 - cs2) / (cf2 - cs2) + alf = sqrt(f1) + als = sqrt(1.0d+00 - f1) + else if (aa2 >= cf2) then + alf = 1.0d+00 + als = 0.0d+00 else - alf2 = 0.0d+00 - als2 = 0.0d+00 + alf = 0.0d+00 + als = 1.0d+00 end if - alf = sqrt(alf2) - als = sqrt(als2) - f3 = vva + al2 / adi_m1 - f4 = sgnb1 * (vya * x2 + vza * x3) - f1 = dnl * (als * f3 - alf * ab * bp) - f2 = dnl * (als * cs * vxa + alf * cf * f4) - wps = f1 + f2 - wms = f1 - f2 - f1 = dnl * (alf * f3 + als * ab * bp) - f2 = dnl * (alf * cf * vxa - als * cs * f4) - wpf = f1 + f2 - wmf = f1 - f2 + cf = sign(sqrt(cf2), vxa) + cs = sign(sqrt(cs2), vxa) + ca = sign(sqrt(ca2), vxa) + ch = sign( cglm, vxa) + + lm(1) = vxa + cf + lm(2) = vxa + ca + lm(3) = vxa + cs + lm(4) = vxa + ch + lm(5) = vxa + lm(6) = vxa - ch + lm(7) = vxa - cs + lm(8) = vxa - ca + lm(9) = vxa - cf -! the scaling diagonal matrix -! f1 = 5.0d-01 / bta - tm(1) = 5.0d-01 / (adiabatic_index * dnl) + f2 = adiabatic_index * dnl + tm(1) = 5.0d-01 / f2 tm(2) = f1 / dnl**2 tm(3) = tm(1) tm(4) = f1 @@ -3453,15 +3688,27 @@ module schemes tm(8) = tm(2) tm(9) = tm(1) -! === the average of the right eigenvector matrix === -! -! right (1) and left (9) fast magnetosonic eigenvectors -! - f1 = dnl * alf - f2 = dnl * als * cs * sgnb1 - f3 = als * ab * sqrt(dnl) + das = dnl * als + daf = dnl * alf + csq = sqrt(f2 / bta) + + f3 = vva + adi_m1xi * prl / dnl + f4 = vya * x2 + vza * x3 + f5 = sqrt(adiabatic_index * bp2 / (dnl * bta)) + f1 = dnl * (als * f3 - alf * f5) + f2 = dnl * (als * cs * vxa + sign(alf, bxa) * cf * f4) + wps = f1 + f2 + wms = f1 - f2 + f1 = dnl * (alf * f3 + als * f5) + f2 = dnl * (alf * cf * vxa - sign(als, bxa) * cs * f4) + wpf = f1 + f2 + wmf = f1 - f2 + + f1 = daf + f2 = sign(das, bxa) * cs + f3 = als * csq rm(1,1) = f1 - rm(2,1) = f1 * (vxa + cf) + rm(2,1) = f1 * lm(1) rm(3,1) = f1 * vya - f2 * x2 rm(4,1) = f1 * vza - f2 * x3 rm(5,1) = wpf @@ -3469,16 +3716,14 @@ module schemes rm(8,1) = f3 * x3 rm(1,9) = f1 - rm(2,9) = f1 * (vxa - cf) + rm(2,9) = f1 * lm(9) rm(3,9) = f1 * vya + f2 * x2 rm(4,9) = f1 * vza + f2 * x3 rm(5,9) = wmf rm(7,9) = rm(7,1) rm(8,9) = rm(8,1) -! right (2) and left (8) Alfvénic eigenvectors -! - f1 = dnl * sqrt(dna) + f1 = sign(dnl * sqrt(dna), vxa) rm(3,2) = f1 * x3 rm(4,2) = - f1 * x2 rm(5,2) = - f1 * (x2 * vza - x3 * vya) @@ -3491,13 +3736,11 @@ module schemes rm(7,8) = rm(7,2) rm(8,8) = rm(8,2) -! right (3) and left (7) slow magnetosonic eigenvectors -! - f1 = dnl * als - f2 = dnl * alf * cf * sgnb1 - f3 = - alf * ab * sqrt(dnl) + f1 = das + f2 = sign(daf, bxa) * cf + f3 = - alf * csq rm(1,3) = f1 - rm(2,3) = f1 * (vxa + cs) + rm(2,3) = f1 * lm(3) rm(3,3) = f1 * vya + f2 * x2 rm(4,3) = f1 * vza + f2 * x3 rm(5,3) = wps @@ -3505,55 +3748,38 @@ module schemes rm(8,3) = f3 * x3 rm(1,7) = f1 - rm(2,7) = f1 * (vxa - cs) + rm(2,7) = f1 * lm(7) rm(3,7) = f1 * vya - f2 * x2 rm(4,7) = f1 * vza - f2 * x3 rm(5,7) = wms rm(7,7) = rm(7,3) rm(8,7) = rm(8,3) -! right (4) and left (6) GLM eigenvectors -! - rm(5,4) = bxa + bpa - rm(5,6) = bxa - bpa + f1 = sign(1.0d+00, vxa) + rm(5,4) = bxa + bpa * f1 + rm(5,6) = bxa - bpa * f1 + rm(9,4) = f1 + rm(9,6) = - f1 -! entropy (5) eigenvectors -! rm(2,5) = vxa rm(3,5) = vya rm(4,5) = vza rm(5,5) = vva -! the eigenvalue diagonal matrix -! - lm(1) = vxa + cf - lm(2) = vxa + ca - lm(3) = vxa + cs - lm(4) = vxa + ch - lm(5) = vxa - lm(6) = vxa - ch - lm(7) = vxa - cs - lm(8) = vxa - ca - lm(9) = vxa - cf - -! KEPEC flux -! fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta fi(imy,i) = fi(idn,i) * vya - bxa * bya fi(imz,i) = fi(idn,i) * vza - bxa * bza - fi(ibx,i) = ch * bpa + fi(ibx,i) = cglm * bpa fi(iby,i) = vxa * bya - bxa * vya fi(ibz,i) = vxa * bza - bxa * vza - fi(ibp,i) = ch * bxa + fi(ibp,i) = cglm * bxa fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) & + (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) & + (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) & + fi(ibp,i) * bpa - ub2 + bxa * uba & - - ch * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i)) + - cglm * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i)) -! KEPES correction -! fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) end do