From 3735112d9abc04a52a66d554d569b4f765489071 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 20 Dec 2013 13:53:45 -0200 Subject: [PATCH 1/2] SCHEMES: Implement isothermal MHD HLLD Riemann solver. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 391 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 391 insertions(+) diff --git a/src/schemes.F90 b/src/schemes.F90 index 126efdd..f65d044 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -195,6 +195,16 @@ module schemes ! select case(trim(solver)) + case ("hlld", "HLLD") + +! set the solver name +! + name_sol = "HLLD" + +! set the solver pointer +! + riemann => riemann_mhd_iso_hlld + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -1593,6 +1603,387 @@ module schemes ! !=============================================================================== ! +! subroutine RIEMANN_MHD_ISO_HLLD: +! ------------------------------- +! +! Subroutine solves one dimensional Riemann problem using the isothermal HLLD +! method with correct state averaging and degeneracies treatement. +! +! Arguments: +! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! f - the output array of fluxes; +! +! References: +! +! [1] Kowal, G., +! "An isothermal MHD Riemann solver with correct state averaging and +! degeneracies treatement", +! Journal of Computational Physics, in preparation +! +!=============================================================================== +! + subroutine riemann_mhd_iso_hlld(n, h, q, 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 + +! local variables are not implicit by default +! + implicit none + +! 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(out) :: f + +! local variables +! + integer :: p, 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) :: wl, wr, wcl, wcr, ui + real(kind=8), dimension(n) :: cl, cr +! +!------------------------------------------------------------------------------- +! +! 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(:,:)) + 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 + +! estimate the minimum and maximum speeds +! + sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) + sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) + +! calculate the HLL flux +! + if (sl >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (sr <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else ! sl < 0 < sr + +! calculate speed difference +! + srml = sr - sl + +! calculate vectors of the left and right-going waves +! + wl(:) = sl * ul(:,i) - fl(:,i) + wr(:) = sr * ur(:,i) - fr(:,i) + +! the advection speed in the intermediate states +! + dn = wr(idn) - wl(idn) + sm = (wr(imx) - wl(imx)) / dn + +! square of Bₓ, i.e. Bₓ² +! + bx = ql(ibx,i) + b2 = ql(ibx,i) * qr(ibx,i) + +! speed differences +! + slmm = sl - sm + srmm = sr - sm + +! left and right state densities +! + dnl = wl(idn) / slmm + dnr = wr(idn) / srmm + +! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to +! the HLL one +! + if (b2 > 0.0d+00) then ! Bₓ² > 0 + +! left and right Alfvén speeds +! + sml = sm - sqrt(b2 / dnl) + smr = sm + sqrt(b2 / dnr) + +! calculate denominators in order to check degeneracy +! + dvl = slmm * wl(idn) - b2 + dvr = srmm * wr(idn) - b2 + +! check degeneracy Sl* -> Sl or Sr* -> Sr +! + if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy + +! choose the correct state depending on the speed signs +! + if (sml >= 0.0d+00) then ! sl* ≥ 0 + +! conservative variables for the outer left intermediate state +! + ui(idn) = dnl + ui(imx) = dnl * sm + ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! the outer left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else if (smr <= 0.0d+00) then ! sr* ≤ 0 + +! conservative variables for the outer right intermediate state +! + ui(idn) = dnr + ui(imx) = dnr * sm + ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! the outer right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sl* < 0 < sr* + +! conservative variables for the outer left intermediate state +! + ui(idn) = dnl + ui(imx) = dnl * sm + ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! vector of the left-going Alfvén wave +! + wcl(:) = (sml - sl) * ui(:) + wl(:) + +! conservative variables for the outer right intermediate state +! + ui(idn) = dnr + ui(imx) = dnr * sm + ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! vector of the right-going Alfvén wave +! + wcr(:) = (smr - sr) * ui(:) + wr(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml) + + end if ! sl* < 0 < sr* + + else ! one degeneracy + +! separate degeneracies +! + if (dvl > 0.0d+00) then ! sr* > sr + +! conservative variables for the outer left intermediate state +! + ui(idn) = dnl + ui(imx) = dnl * sm + ui(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! choose the correct state depending on the speed signs +! + if (sml >= 0.0d+00) then ! sl* ≥ 0 + +! the outer left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else ! sl* < 0 + +! vector of the left-going Alfvén wave +! + wcl(:) = (sml - sl) * ui(:) + wl(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml) + + end if ! sl* < 0 + + else if (dvr > 0.0d+00) then ! sl* < sl + +! conservative variables for the outer right intermediate state +! + ui(idn) = dnr + ui(imx) = dnr * sm + ui(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! choose the correct state depending on the speed signs +! + if (smr <= 0.0d+00) then ! sr* ≤ 0 + +! the outer right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sr* > 0 + +! vector of the right-going Alfvén wave +! + wcr(:) = (smr - sr) * ui(:) + wr(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl) + + end if ! sr* > 0 + + else ! sl* < sl & sr* > sr + +! both outer states are degenerate so revert to the HLL flux +! + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml + + end if ! sl* < sl & sr* > sr + + end if ! one degeneracy + + else ! Bₓ² = 0 + +! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and +! still we can have a discontinuity in the perpendicular components +! + dn = dn / srml + +! take the right flux depending on the sign of the advection speed +! + if (sm > 0.0d+00) then ! sm > 0 + +! conservative variables for the left intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = wl(imy) / slmm + ui(imz) = wl(imz) / slmm + ui(ibx) = 0.0d+00 + ui(iby) = wl(iby) / slmm + ui(ibz) = wl(ibz) / slmm + ui(ibp) = ql(ibp,i) + +! the left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else if (sm < 0.0d+00) then ! sm < 0 + +! conservative variables for the right intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = wr(imy) / srmm + ui(imz) = wr(imz) / srmm + ui(ibx) = 0.0d+00 + ui(iby) = wr(iby) / srmm + ui(ibz) = wr(ibz) / srmm + ui(ibp) = qr(ibp,i) + +! the right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sm = 0 + +! the intermediate flux; since the advection speed is zero, perpendicular +! components do not change +! + f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml + f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml + f(imy,i) = 0.0d+00 + f(imz,i) = 0.0d+00 + f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml + f(iby,i) = 0.0d+00 + f(ibz,i) = 0.0d+00 + f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml + + end if ! sm = 0 + + end if ! Bₓ² = 0 + + end if ! sl < 0 < sr + + end do ! i = 1, n + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_iso_hlld +! +!=============================================================================== +! ! subroutine RIEMANN_MHD_ADI_HLL: ! ------------------------------ ! From 9517e6ce4c82b136ede3a75195bc0e870253b3eb Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 20 Dec 2013 14:01:33 -0200 Subject: [PATCH 2/2] SCHEMES: Implement isothermal MHD HLLD-M Riemann solver. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 388 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 388 insertions(+) diff --git a/src/schemes.F90 b/src/schemes.F90 index f65d044..bc4500f 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -205,6 +205,16 @@ module schemes ! riemann => riemann_mhd_iso_hlld + case ("hlldm", "hlld-m", "HLLDM", "HLLD-M") + +! set the solver name +! + name_sol = "HLLD-M" + +! set the solver pointer +! + riemann => riemann_mhd_iso_hlldm + ! in the case of unknown Riemann solver, revert to HLL ! case default @@ -1984,6 +1994,384 @@ module schemes ! !=============================================================================== ! +! subroutine RIEMANN_MHD_ISO_HLLDM: +! -------------------------------- +! +! Subroutine solves one dimensional Riemann problem using the isothermal HLLD +! method described by Mignone. +! +! Arguments: +! +! n - the length of input vectors; +! h - the spatial step; +! q - the input array of primitive variables; +! f - the output array of fluxes; +! +! References: +! +! [1] Mignone, A., +! "A simple and accurate Riemann solver for isothermal MHD", +! Journal of Computational Physics, 2007, 225, pp. 1427-1441 +! +!=============================================================================== +! + subroutine riemann_mhd_iso_hlldm(n, h, q, 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 + +! local variables are not implicit by default +! + implicit none + +! 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(out) :: f + +! local variables +! + integer :: p, 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) :: wl, wr, wcl, wcr, ui + real(kind=8), dimension(n) :: cl, cr +! +!------------------------------------------------------------------------------- +! +! 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(:,:)) + 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 + +! estimate the minimum and maximum speeds +! + sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) + sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) + +! calculate the HLL flux +! + if (sl >= 0.0d+00) then + + f(:,i) = fl(:,i) + + else if (sr <= 0.0d+00) then + + f(:,i) = fr(:,i) + + else ! sl < 0 < sr + +! calculate speed difference +! + srml = sr - sl + +! calculate vectors of the left and right-going waves +! + wl(:) = sl * ul(:,i) - fl(:,i) + wr(:) = sr * ur(:,i) - fr(:,i) + +! the advection speed in the intermediate states +! + dn = wr(idn) - wl(idn) + sm = (wr(imx) - wl(imx)) / dn + +! square of Bₓ, i.e. Bₓ² +! + bx = ql(ibx,i) + b2 = ql(ibx,i) * qr(ibx,i) + +! speed differences +! + slmm = sl - sm + srmm = sr - sm + +! calculate density of the intermediate states +! + dn = dn / srml + +! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to +! the HLL one +! + if (b2 > 0.0d+00) then ! Bₓ² > 0 + +! left and right Alfvén speeds +! + ca = sqrt(b2 / dn) + sml = sm - ca + smr = sm + ca + +! calculate denominators in order to check degeneracy +! + dvl = slmm * slmm * dn - b2 + dvr = srmm * srmm * dn - b2 + +! check degeneracy Sl* -> Sl or Sr* -> Sr +! + if (min(dvl, dvr) > 0.0d+00) then ! no degeneracy + +! choose the correct state depending on the speed signs +! + if (sml >= 0.0d+00) then ! sl* ≥ 0 + +! conservative variables for the outer left intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! the outer left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else if (smr <= 0.0d+00) then ! sr* ≤ 0 + +! conservative variables for the outer right intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! the outer right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sl* < 0 < sr* + +! conservative variables for the outer left intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! vector of the left-going Alfvén wave +! + wcl(:) = (sml - sl) * ui(:) + wl(:) + +! conservative variables for the outer right intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! vector of the right-going Alfvén wave +! + wcr(:) = (smr - sr) * ui(:) + wr(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml) + + end if ! sl* < 0 < sr* + + else ! one degeneracy + +! separate degeneracies +! + if (dvl > 0.0d+00) then ! sr* > sr + +! conservative variables for the outer left intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (slmm * wl(imy) - bx * wl(iby)) / dvl + ui(imz) = dn * (slmm * wl(imz) - bx * wl(ibz)) / dvl + ui(ibx) = bx + ui(iby) = (slmm * dn * wl(iby) - bx * wl(imy)) / dvl + ui(ibz) = (slmm * dn * wl(ibz) - bx * wl(imz)) / dvl + ui(ibp) = ql(ibp,i) + +! choose the correct state depending on the speed signs +! + if (sml >= 0.0d+00) then ! sl* ≥ 0 + +! the outer left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else ! sl* < 0 + +! vector of the left-going Alfvén wave +! + wcl(:) = (sml - sl) * ui(:) + wl(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml) + + end if ! sl* < 0 + + else if (dvr > 0.0d+00) then ! sl* < sl + +! conservative variables for the outer right intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = dn * (srmm * wr(imy) - bx * wr(iby)) / dvr + ui(imz) = dn * (srmm * wr(imz) - bx * wr(ibz)) / dvr + ui(ibx) = bx + ui(iby) = (srmm * dn * wr(iby) - bx * wr(imy)) / dvr + ui(ibz) = (srmm * dn * wr(ibz) - bx * wr(imz)) / dvr + ui(ibp) = qr(ibp,i) + +! choose the correct state depending on the speed signs +! + if (smr <= 0.0d+00) then ! sr* ≤ 0 + +! the outer right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sr* > 0 + +! vector of the right-going Alfvén wave +! + wcr(:) = (smr - sr) * ui(:) + wr(:) + +! the flux corresponding to the middle intermediate state +! + f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl) + + end if ! sr* > 0 + + else ! sl* < sl & sr* > sr + +! both outer states are degenerate so revert to the HLL flux +! + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml + + end if ! sl* < sl & sr* > sr + + end if ! one degeneracy + + else ! Bₓ² = 0 + +! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and +! still we can have a discontinuity in the perpendicular components +! +! take the right flux depending on the sign of the advection speed +! + if (sm > 0.0d+00) then ! sm > 0 + +! conservative variables for the left intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = wl(imy) / slmm + ui(imz) = wl(imz) / slmm + ui(ibx) = 0.0d+00 + ui(iby) = wl(iby) / slmm + ui(ibz) = wl(ibz) / slmm + ui(ibp) = ql(ibp,i) + +! the left intermediate flux +! + f(:,i) = sl * ui(:) - wl(:) + + else if (sm < 0.0d+00) then ! sm < 0 + +! conservative variables for the right intermediate state +! + ui(idn) = dn + ui(imx) = dn * sm + ui(imy) = wr(imy) / srmm + ui(imz) = wr(imz) / srmm + ui(ibx) = 0.0d+00 + ui(iby) = wr(iby) / srmm + ui(ibz) = wr(ibz) / srmm + ui(ibp) = qr(ibp,i) + +! the right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sm = 0 + +! the intermediate flux; since the advection speed is zero, perpendicular +! components do not change +! + f(idn,i) = (sl * wr(idn) - sr * wl(idn)) / srml + f(imx,i) = (sl * wr(imx) - sr * wl(imx)) / srml + f(imy,i) = 0.0d+00 + f(imz,i) = 0.0d+00 + f(ibx,i) = (sl * wr(ibx) - sr * wl(ibx)) / srml + f(iby,i) = 0.0d+00 + f(ibz,i) = 0.0d+00 + f(ibp,i) = (sl * wr(ibp) - sr * wl(ibp)) / srml + + end if ! sm = 0 + + end if ! Bₓ² = 0 + + end if ! sl < 0 < sr + + end do ! i = 1, n + +!------------------------------------------------------------------------------- +! + end subroutine riemann_mhd_iso_hlldm +! +!=============================================================================== +! ! subroutine RIEMANN_MHD_ADI_HLL: ! ------------------------------ !