From 322b8a141b3d411d21b8a9410a2a219a860d3e58 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 9 Mar 2017 16:26:28 -0300 Subject: [PATCH 01/13] USER_PROBLEM: Reset acc in gravitational_acceleration_user(). Signed-off-by: Grzegorz Kowal --- src/user_problem.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/user_problem.F90 b/src/user_problem.F90 index 13624a6..a466f41 100644 --- a/src/user_problem.F90 +++ b/src/user_problem.F90 @@ -407,6 +407,10 @@ module user_problem call start_timer(img) #endif /* PROFILE */ +! reset gravitational acceleration +! + acc(:) = 0.0d+00 + #ifdef PROFILE ! stop accounting time for the gravitational acceleration calculation ! From 31bf3006b7c507a448f3d3d9da1fbcb75cde4032 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 9 Mar 2017 16:28:30 -0300 Subject: [PATCH 02/13] GRAVITY: Reset acc in gacc_none(). Signed-off-by: Grzegorz Kowal --- src/gravity.F90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/gravity.F90 b/src/gravity.F90 index 200fb64..82824ec 100644 --- a/src/gravity.F90 +++ b/src/gravity.F90 @@ -216,11 +216,11 @@ module gravity ! ! t, dt - time and the time increment; ! x, y, z - rectangular coordinates; -! gacc - vector of the gravitational acceleration; +! acc - vector of the gravitational acceleration; ! !=============================================================================== ! - subroutine gacc_none(t, dt, x, y, z, gacc) + subroutine gacc_none(t, dt, x, y, z, acc) ! local variables are not implicit by default ! @@ -230,7 +230,7 @@ module gravity ! real(kind=8) , intent(in) :: t, dt real(kind=8) , intent(in) :: x, y, z - real(kind=8), dimension(3), intent(out) :: gacc + real(kind=8), dimension(3), intent(out) :: acc ! !------------------------------------------------------------------------------- ! @@ -240,6 +240,10 @@ module gravity call start_timer(imc) #endif /* PROFILE */ +! reset gravitational acceleration +! + acc(:) = 0.0d+00 + #ifdef PROFILE ! stop accounting time for the gravitational acceleration calculation ! From fe0e43b095160338631a565aded722b39fbdccfd Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 9 Mar 2017 16:31:16 -0300 Subject: [PATCH 03/13] IO: Reset iret in read_restart_snapshot() and write_restart_snapshot(). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/io.F90 b/src/io.F90 index 1381bb2..9dfbf8c 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -360,6 +360,10 @@ module io call start_timer(ios) #endif /* PROFILE */ +! reset the return flag +! + iret = 0 + #ifdef HDF5 ! read HDF5 restart file and rebuild the meta and data block structures ! @@ -411,6 +415,10 @@ module io ! !------------------------------------------------------------------------------- ! +! reset the return flag +! + iret = 0 + ! check if conditions for storing the restart snapshot have been met ! if (hrest < 5.0d-02 .or. thrs < irest * hrest) return From 5213004734de1264d74e31777bcb01fa3e7729a7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 10 Mar 2017 12:04:42 -0300 Subject: [PATCH 04/13] IO: Account time for disk operations (for snapshots). Signed-off-by: Grzegorz Kowal --- src/io.F90 | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/src/io.F90 b/src/io.F90 index 1381bb2..2a013b2 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -33,9 +33,7 @@ module io ! import external subroutines ! use blocks, only : pointer_meta -#ifdef PROFILE use timers, only : set_timer, start_timer, stop_timer -#endif /* PROFILE */ ! module variables are not implicit by default ! @@ -89,9 +87,10 @@ module io #endif /* HDF5 */ end interface -#ifdef PROFILE ! timer indices ! + integer , save :: iio +#ifdef PROFILE integer , save :: ioi, iow, ios #endif /* PROFILE */ @@ -204,9 +203,10 @@ module io ! !------------------------------------------------------------------------------- ! -#ifdef PROFILE ! set timer descriptions ! + call set_timer('SNAPSHOTS I/O' , iio) +#ifdef PROFILE call set_timer('io:: initialization' , ioi) call set_timer('io:: snapshot writing', iow) call set_timer('io:: snapshot reading', ios) @@ -360,12 +360,20 @@ module io call start_timer(ios) #endif /* PROFILE */ +! start accounting time for I/O +! + call start_timer(iio) + #ifdef HDF5 ! read HDF5 restart file and rebuild the meta and data block structures ! call read_restart_snapshot_h5(iret) #endif /* HDF5 */ +! stop accounting time for I/O +! + call stop_timer(iio) + ! calculate the shift of the snapshot counter, and the next snapshot time ! ishift = int(time / hsnap) - isnap + 1 @@ -421,12 +429,20 @@ module io call start_timer(iow) #endif /* PROFILE */ +! start accounting time for I/O +! + call start_timer(iio) + #ifdef HDF5 ! store restart file ! call write_restart_snapshot_h5(nrun, iret) #endif /* HDF5 */ +! stop accounting time for I/O +! + call stop_timer(iio) + ! increase the restart snapshot counter ! irest = irest + 1 @@ -476,6 +492,10 @@ module io call start_timer(iow) #endif /* PROFILE */ +! start accounting time for I/O +! + call start_timer(iio) + #ifdef HDF5 ! store variable snapshot file ! @@ -486,6 +506,10 @@ module io end if #endif /* HDF5 */ +! stop accounting time for I/O +! + call stop_timer(iio) + ! increase the snapshot counter and calculate the next snapshot time ! isnap = isnap + 1 From b092c80aaeb647f845ad8c152d8428532bd39ca7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 15 Mar 2017 10:18:39 -0300 Subject: [PATCH 05/13] INTERPOLATIONS: Fix misspelling. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index c04521d..cc54e52 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -868,7 +868,7 @@ module interpolations ! ------------------------- ! ! Subroutine reconstructs both side interfaces of variable using -! multidimentional Gaussian Process method. +! multidimensional Gaussian Process method. ! ! Arguments: ! From aad62c9c519b5f283be8f07571e9331637e1181f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 16 Mar 2017 09:50:10 -0300 Subject: [PATCH 06/13] INTERPOLATIONS: Implement Multi-dimensional Limiting Process. This additional multi-dimensional limiting can be applied to any reconstruction method. It is controlled by the parameter 'mlp_limiting' by setting 'on' or 'off'. The implementation is based on Gerlinger (2012). Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 206 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 203 insertions(+), 3 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index cc54e52..02adebb 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -86,6 +86,7 @@ module interpolations ! flags for reconstruction corrections ! + logical , save :: mlp = .false. logical , save :: positivity = .false. logical , save :: clip = .false. @@ -137,6 +138,7 @@ module interpolations character(len=255) :: tlimiter = "mm" character(len=255) :: plimiter = "mm" character(len=255) :: climiter = "mm" + character(len=255) :: mlp_limiting = "off" character(len=255) :: positivity_fix = "off" character(len=255) :: clip_extrema = "off" character(len=255) :: name_rec = "" @@ -169,6 +171,7 @@ module interpolations call get_parameter_string ("clip_extrema" , clip_extrema ) call get_parameter_string ("extrema_limiter" , climiter ) call get_parameter_string ("prolongation_limiter", plimiter ) + call get_parameter_string ("mlp_limiting" , mlp_limiting ) call get_parameter_integer("nghosts" , ng ) call get_parameter_integer("ngp" , ngp ) call get_parameter_real ("sgp" , sgp ) @@ -390,6 +393,12 @@ module interpolations ! check additional reconstruction limiting ! + select case(trim(mlp_limiting)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + mlp = .true. + case default + mlp = .false. + end select select case(trim(positivity_fix)) case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") positivity = .true. @@ -410,6 +419,7 @@ module interpolations write (*,"(4x,a14, 9x,'=',1x,a)") "reconstruction" , trim(name_rec) write (*,"(4x,a11,12x,'=',1x,a)") "TVD limiter" , trim(name_tlim) write (*,"(4x,a20, 3x,'=',1x,a)") "prolongation limiter", trim(name_plim) + write (*,"(4x,a12,11x,'=',1x,a)") "MLP limiting" , trim(mlp_limiting) write (*,"(4x,a14, 9x,'=',1x,a)") "fix positivity" , trim(positivity_fix) write (*,"(4x,a12,11x,'=',1x,a)") "clip extrema" , trim(clip_extrema) if (clip) then @@ -673,9 +683,12 @@ module interpolations ! local variables ! - integer :: i, im1, ip1 - integer :: j, jm1, jp1 - integer :: k, km1, kp1 + integer :: i, im1, ip1 + integer :: j, jm1, jp1 + integer :: k, km1, kp1 + +! local vectors +! real(kind=8), dimension(NDIMS) :: dql, dqr, dq ! !------------------------------------------------------------------------------- @@ -751,6 +764,10 @@ module interpolations end do ! j = jbl, jeu end do ! k = kbl, keu +! apply Multi-dimensional Limiting Process +! + if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:)) + !------------------------------------------------------------------------------- ! end subroutine interfaces_tvd @@ -858,6 +875,10 @@ module interpolations end if +! apply Multi-dimensional Limiting Process +! + if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:)) + !------------------------------------------------------------------------------- ! end subroutine interfaces_dir @@ -1028,12 +1049,191 @@ module interpolations end do ! j = jbl, jeu end do ! k = kbl, keu +! apply Multi-dimensional Limiting Process +! + if (mlp) call mlp_limiting(q(:,:,:), qi(:,:,:,:,:)) + !------------------------------------------------------------------------------- ! end subroutine interfaces_mgp ! !=============================================================================== ! +! subroutine MLP_LIMITING: +! ----------------------- +! +! Subroutine applies the multi-dimensional limiting process to +! the reconstructed states in order to control oscillations due to +! the Gibbs phenomena near discontinuities. +! +! Arguments: +! +! q - the variable array; +! qi - the array of reconstructed interfaces (2 in each direction); +! +! References: +! +! [1] Gerlinger, P., +! "Multi-dimensional limiting for high-order schemes including +! turbulence and combustion", +! Journal of Computational Physics, +! 2012, vol. 231, pp. 2199-2228, +! http://dx.doi.org/10.1016/j.jcp.2011.10.024 +! +!=============================================================================== +! + subroutine mlp_limiting(q, qi) + +! include external procedures +! + use coordinates , only : im , jm , km + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), dimension(im,jm,km) , intent(in) :: q + real(kind=8), dimension(im,jm,km,2,NDIMS), intent(inout) :: qi + +! local variables +! + integer :: i, im1, ip1 + integer :: j, jm1, jp1 + integer :: k, km1, kp1 + integer :: m +#if NDIMS == 3 + integer :: n, np1, np2 +#endif /* NDIMS == 3 */ + real(kind=8) :: qmn, qmx, dqc + real(kind=8) :: fc, fl + +! local vectors +! + real(kind=8), dimension(NDIMS) :: dql, dqr, dq + real(kind=8), dimension(NDIMS) :: dqm, ap +#if NDIMS == 3 + real(kind=8), dimension(NDIMS) :: hh, uu +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! + do k = kbl, keu +#if NDIMS == 3 + km1 = k - 1 + kp1 = k + 1 +#endif /* NDIMS == 3 */ + do j = jbl, jeu + jm1 = j - 1 + jp1 = j + 1 + do i = ibl, ieu + im1 = i - 1 + ip1 = i + 1 + +! calculate the minmod TVD derivatives +! + dql(1) = q(i ,j,k) - q(im1,j,k) + dqr(1) = q(ip1,j,k) - q(i ,j,k) + dq (1) = minmod(dql(1), dqr(1)) + + dql(2) = q(i,j ,k) - q(i,jm1,k) + dqr(2) = q(i,jp1,k) - q(i,j ,k) + dq (2) = minmod(dql(2), dqr(2)) + +#if NDIMS == 3 + dql(3) = q(i,j,k ) - q(i,j,km1) + dqr(3) = q(i,j,kp1) - q(i,j,k ) + dq (3) = minmod(dql(3), dqr(3)) +#endif /* NDIMS == 3 */ + +! calculate dqc +! +#if NDIMS == 3 + qmn = minval(q(im1:ip1,jm1:jp1,km1:kp1)) + qmx = maxval(q(im1:ip1,jm1:jp1,km1:kp1)) +#else /* NDIMS == 3 */ + qmn = minval(q(im1:ip1,jm1:jp1,k)) + qmx = maxval(q(im1:ip1,jm1:jp1,k)) +#endif /* NDIMS == 3 */ + dqc = min(qmx - q(i,j,k), q(i,j,k) - qmn) + +! if needed, apply the multi-dimensional limiting process +! + if (sum(abs(dq(1:NDIMS))) > dqc) then + + dqm(1) = abs(q(ip1,j,k) - q(im1,j,k)) + dqm(2) = abs(q(i,jp1,k) - q(i,jm1,k)) +#if NDIMS == 3 + dqm(3) = abs(q(i,j,kp1) - q(i,j,km1)) +#endif /* NDIMS == 3 */ + + fc = dqc / max(1.0d-16, sum(abs(dqm(1:NDIMS)))) + do m = 1, NDIMS + ap(m) = fc * abs(dqm(m)) + end do + +#if NDIMS == 3 + do n = 1, NDIMS + hh(n) = max(ap(n) - abs(dq(n)), 0.0d+00) + np1 = mod(n , NDIMS) + 1 + np2 = mod(n + 1, NDIMS) + 1 + uu(n) = ap(n) - hh(n) + uu(np1) = ap(np1) + 0.5d+00 * hh(n) + uu(np2) = ap(np2) + 0.5d+00 * hh(n) + fc = hh(n) / (hh(n) + 1.0d-16) + fl = fc * (max(ap(np1) - abs(dq(np1)), 0.0d+00) & + - max(ap(np2) - abs(dq(np2)), 0.0d+00)) + ap(n ) = uu(n) + ap(np1) = uu(np1) - fl + ap(np2) = uu(np2) + fl + end do +#else /* NDIMS == 3 */ + fl = max(ap(1) - abs(dq(1)), 0.0d+00) & + - max(ap(2) - abs(dq(2)), 0.0d+00) + ap(1) = ap(1) - fl + ap(2) = ap(2) + fl +#endif /* NDIMS == 3 */ + + do m = 1, NDIMS + dq(m) = sign(ap(m), dq(m)) + end do + end if + +! calculate the limited variable increments +! + dql(1) = minmod(dq(1), qi(i ,j,k,1,1) - q(i,j,k)) + dqr(1) = minmod(dq(1), - qi(im1,j,k,2,1) + q(i,j,k)) + dql(2) = minmod(dq(2), qi(i,j ,k,1,2) - q(i,j,k)) + dqr(2) = minmod(dq(2), - qi(i,jm1,k,2,2) + q(i,j,k)) +#if NDIMS == 3 + dql(3) = minmod(dq(3), qi(i,j,k ,1,3) - q(i,j,k)) + dqr(3) = minmod(dq(3), - qi(i,j,km1,2,3) + q(i,j,k)) +#endif /* NDIMS == 3 */ + +! update the interpolated states +! + qi(i ,j,k,1,1) = q(i,j,k) + dql(1) + qi(im1,j,k,2,1) = q(i,j,k) - dqr(1) + + qi(i,j ,k,1,2) = q(i,j,k) + dql(2) + qi(i,jm1,k,2,2) = q(i,j,k) - dqr(2) +#if NDIMS == 3 + qi(i,j,k ,1,3) = q(i,j,k) + dql(3) + qi(i,j,km1,2,3) = q(i,j,k) - dqr(3) +#endif /* NDIMS == 3 */ + + end do ! i = ibl, ieu + end do ! j = jbl, jeu + end do ! k = kbl, keu + +!------------------------------------------------------------------------------- +! + end subroutine mlp_limiting +! +!=============================================================================== +! ! subroutine RECONSTRUCT: ! ---------------------- ! From 3cd24e3b0f26fce53b5e775916f0e96cb23ad767 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 29 Mar 2017 10:40:31 -0300 Subject: [PATCH 07/13] =?UTF-8?q?SOURCES:=20Fix=20calculation=20of=20?= =?UTF-8?q?=CE=B7=20J=C2=B2=20term.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In the sum of the J² we used Fortran subroutine sum() and the sum was done along the 2nd index. However, the interpretation of this index is uncertain when the first array rank is 1. Therefore, instead of using subroutine sum(), the calculation is done directly. After this change there are no strange effects appearing on the boundaries of the blocks. Signed-off-by: Grzegorz Kowal --- src/sources.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/sources.F90 b/src/sources.F90 index 88d198f..1c13354 100644 --- a/src/sources.F90 +++ b/src/sources.F90 @@ -582,11 +582,17 @@ module sources call curl(dh(:), pdata%q(ibx:ibz,1:im,1:jm,1:km) & , tmp(inz,inx:inz,1:im,1:jm,1:km)) +! calculate J² +! + db(1:im,1:jm,1:km) = tmp(inz,inx,1:im,1:jm,1:km)**2 & + + tmp(inz,iny,1:im,1:jm,1:km)**2 & + + tmp(inz,inz,1:im,1:jm,1:km)**2 + ! add the second resistive source term to the energy equation, i.e. ! d/dt E + ∇.F = η J² ! du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) & - + resistivity * sum(tmp(inz,inx:inz,1:im,1:jm,1:km)**2, 2) + + resistivity * db(1:im,1:jm,1:km) end if ! energy equation present From b4b72a783d19dbe0bfd4a0515abef3aa184ce57b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 5 Apr 2017 07:23:38 -0300 Subject: [PATCH 08/13] SCHEMES: A small fix for isothermal HLLD when Bx = 0. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 24 ++++++------------------ 1 file changed, 6 insertions(+), 18 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index d07c90d..426940d 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -2266,19 +2266,14 @@ module schemes 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(idn) = dnl + ui(imx) = dnl * sm ui(imy) = wl(imy) / slmm ui(imz) = wl(imz) / slmm ui(ibx) = 0.0d+00 @@ -2294,8 +2289,8 @@ module schemes ! conservative variables for the right intermediate state ! - ui(idn) = dn - ui(imx) = dn * sm + ui(idn) = dnr + ui(imx) = dnr * sm ui(imy) = wr(imy) / srmm ui(imz) = wr(imz) / srmm ui(ibx) = 0.0d+00 @@ -2310,16 +2305,9 @@ module schemes else ! sm = 0 ! the intermediate flux; since the advection speed is zero, perpendicular -! components do not change +! components do not change, so revert it to HLL ! - 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 + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sm = 0 From b41b5a000b9833d389cccb57e11c490a456c1df2 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 5 Apr 2017 10:23:15 -0300 Subject: [PATCH 09/13] INTERPOLATIONS: Fix MLP limiting. We should take half of the TVD limited derivatives in order to compare properly with the high order interpolated derivative. Fix it. Also, TVD limit both states if the high order interpolation of any of them overshoot. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 47 ++++++++++++++++++++++-------------------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 02adebb..80fa8c1 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -1136,16 +1136,16 @@ module interpolations ! dql(1) = q(i ,j,k) - q(im1,j,k) dqr(1) = q(ip1,j,k) - q(i ,j,k) - dq (1) = minmod(dql(1), dqr(1)) + dq (1) = limiter_minmod(0.5d+00, dql(1), dqr(1)) dql(2) = q(i,j ,k) - q(i,jm1,k) dqr(2) = q(i,jp1,k) - q(i,j ,k) - dq (2) = minmod(dql(2), dqr(2)) + dq (2) = limiter_minmod(0.5d+00, dql(2), dqr(2)) #if NDIMS == 3 dql(3) = q(i,j,k ) - q(i,j,km1) dqr(3) = q(i,j,kp1) - q(i,j,k ) - dq (3) = minmod(dql(3), dqr(3)) + dq (3) = limiter_minmod(0.5d+00, dql(3), dqr(3)) #endif /* NDIMS == 3 */ ! calculate dqc @@ -1197,33 +1197,36 @@ module interpolations #endif /* NDIMS == 3 */ do m = 1, NDIMS - dq(m) = sign(ap(m), dq(m)) + dq(m) = sign(ap(m), dq(m)) end do - end if - -! calculate the limited variable increments -! - dql(1) = minmod(dq(1), qi(i ,j,k,1,1) - q(i,j,k)) - dqr(1) = minmod(dq(1), - qi(im1,j,k,2,1) + q(i,j,k)) - dql(2) = minmod(dq(2), qi(i,j ,k,1,2) - q(i,j,k)) - dqr(2) = minmod(dq(2), - qi(i,jm1,k,2,2) + q(i,j,k)) -#if NDIMS == 3 - dql(3) = minmod(dq(3), qi(i,j,k ,1,3) - q(i,j,k)) - dqr(3) = minmod(dq(3), - qi(i,j,km1,2,3) + q(i,j,k)) -#endif /* NDIMS == 3 */ ! update the interpolated states ! - qi(i ,j,k,1,1) = q(i,j,k) + dql(1) - qi(im1,j,k,2,1) = q(i,j,k) - dqr(1) + dql(1) = qi(i ,j,k,1,1) - q(i,j,k) + dqr(1) = qi(im1,j,k,2,1) - q(i,j,k) + if (max(abs(dql(1)), abs(dqr(1))) > abs(dq(1))) then + qi(i ,j,k,1,1) = q(i,j,k) + dq(1) + qi(im1,j,k,2,1) = q(i,j,k) - dq(1) + end if + + dql(2) = qi(i,j ,k,1,2) - q(i,j,k) + dqr(2) = qi(i,jm1,k,2,2) - q(i,j,k) + if (max(abs(dql(2)), abs(dqr(2))) > abs(dq(2))) then + qi(i,j ,k,1,2) = q(i,j,k) + dq(2) + qi(i,jm1,k,2,2) = q(i,j,k) - dq(2) + end if - qi(i,j ,k,1,2) = q(i,j,k) + dql(2) - qi(i,jm1,k,2,2) = q(i,j,k) - dqr(2) #if NDIMS == 3 - qi(i,j,k ,1,3) = q(i,j,k) + dql(3) - qi(i,j,km1,2,3) = q(i,j,k) - dqr(3) + dql(3) = qi(i,j,k ,1,3) - q(i,j,k)) + dqr(3) = qi(i,j,km1,2,3) - q(i,j,k)) + if (max(abs(dql(3)), abs(dqr(3))) > abs(dq(3))) then + qi(i,j,k ,1,3) = q(i,j,k) + dq(3) + qi(i,j,km1,2,3) = q(i,j,k) - dq(3) + end if #endif /* NDIMS == 3 */ + end if + end do ! i = ibl, ieu end do ! j = jbl, jeu end do ! k = kbl, keu From 3a6090499adda0969bb7bdeba09b6d23b9f09a5f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 7 Apr 2017 08:34:44 -0300 Subject: [PATCH 10/13] SCHEMES: Fix adiabatic HLLD solver for MHD. For one of the degenerate situations the state vector was not calculated. This introduced numerical instabilities. This change fixes it. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 98 +++++++++++++++++++++++-------------------------- 1 file changed, 45 insertions(+), 53 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index 426940d..71ff769 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -3717,7 +3717,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! the outer left intermediate flux @@ -3743,7 +3743,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! the outer right intermediate flux @@ -3769,7 +3769,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! vector of the left-going Alfvén wave @@ -3793,7 +3793,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! vector of the right-going Alfvén wave @@ -3824,7 +3824,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal ! the inmost left intermediate flux @@ -3842,7 +3842,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car ! the inmost right intermediate flux @@ -3883,7 +3883,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! choose the correct state depending on the speed signs @@ -3910,23 +3910,23 @@ module schemes bz = (wr(ibz) - wcl(ibz)) / dv vb = sm * bx + vy * by + vz * bz +! conservative variables for the inner left intermediate state +! + ui(idn) = dnl + ui(imx) = dnl * sm + ui(imy) = dnl * vy + ui(imz) = dnl * vz + ui(ibx) = bx + ui(iby) = by + ui(ibz) = bz + ui(ibp) = ul(ibp,i) + ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal + ! choose the correct state depending on the sign of contact discontinuity ! advection speed ! if (sm >= 0.0d+00) then ! sm ≥ 0 -! conservative variables for the inner left intermediate state -! - ui(idn) = dnl - ui(imx) = dnl * sm - ui(imy) = dnl * vy - ui(imz) = dnl * vz - ui(ibx) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = ql(ibp,i) - ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal - ! the inner left intermediate flux ! f(:,i) = sml * ui(:) - wcl(:) @@ -3935,11 +3935,11 @@ module schemes ! vector of the left-going Alfvén wave ! - wcl(:) = (sm - sml) * ui(:) + wcl(:) + wcr(:) = (sm - sml) * ui(:) + wcl(:) ! calculate the average flux over the right inner intermediate state ! - f(:,i) = (sm * wr(:) - sr * wcl(:)) / srmm + f(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm end if ! sm < 0 @@ -3964,7 +3964,7 @@ module schemes ui(ibx) = bx ui(iby) = by ui(ibz) = bz - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! choose the correct state depending on the speed signs @@ -3991,23 +3991,23 @@ module schemes bz = (wl(ibz) - wcr(ibz)) / dv vb = sm * bx + vy * by + vz * bz +! conservative variables for the inner left intermediate state +! + ui(idn) = dnr + ui(imx) = dnr * sm + ui(imy) = dnr * vy + ui(imz) = dnr * vz + ui(ibx) = bx + ui(iby) = by + ui(ibz) = bz + ui(ibp) = ur(ibp,i) + ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car + ! choose the correct state depending on the sign of contact discontinuity ! advection speed ! if (sm <= 0.0d+00) then ! sm ≤ 0 -! conservative variables for the inner left intermediate state -! - ui(idn) = dnr - ui(imx) = dnr * sm - ui(imy) = dnr * vy - ui(imz) = dnr * vz - ui(ibx) = bx - ui(iby) = by - ui(ibz) = bz - ui(ibp) = qr(ibp,i) - ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car - ! the inner right intermediate flux ! f(:,i) = smr * ui(:) - wcr(:) @@ -4016,11 +4016,11 @@ module schemes ! vector of the right-going Alfvén wave ! - wcr(:) = (sm - smr) * ui(:) + wcr(:) + wcl(:) = (sm - smr) * ui(:) + wcr(:) ! calculate the average flux over the left inner intermediate state ! - f(:,i) = (sm * wl(:) - sl * wcr(:)) / slmm + f(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm end if ! sm > 0 @@ -4054,12 +4054,12 @@ module schemes ! ui(idn) = wl(idn) / slmm ui(imx) = ui(idn) * sm - ui(imy) = ui(idn) * ql(ivy,i) - ui(imz) = ui(idn) * ql(ivz,i) + 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) + ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt) / slmm ! the left intermediate flux @@ -4076,12 +4076,12 @@ module schemes ! ui(idn) = wr(idn) / srmm ui(imx) = ui(idn) * sm - ui(imy) = ui(idn) * qr(ivy,i) - ui(imz) = ui(idn) * qr(ivz,i) + 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) + ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt) / srmm ! the right intermediate flux @@ -4090,18 +4090,10 @@ module schemes else ! sm = 0 -! intermediate flux is constant across the contact discontinuity and all except -! the parallel momentum flux are zero +! intermediate flux is constant across the contact discontinuity, so we end up +! having only one intermediate state, which is equivalent to the HLL solver ! - f(idn,i) = 0.0d+00 - f(imx,i) = - 0.5d+00 * (wl(imx) + wr(imx)) - f(imy,i) = 0.0d+00 - f(imz,i) = 0.0d+00 - f(ibx,i) = fl(ibx,i) - f(iby,i) = 0.0d+00 - f(ibz,i) = 0.0d+00 - f(ibp,i) = fl(ibp,i) - f(ien,i) = 0.0d+00 + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sm = 0 From a6490f31bbef1a6a24b9e613207523a03fe01452 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 10 Apr 2017 07:56:49 -0300 Subject: [PATCH 11/13] SCHEMES: Fix isothermal HLLD solver for weak Bx. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If the parallel component of magnetic field is too small, the intermediate states might produce numerical instabilities, since they are obtained by the division by a small factor. In order to remove this situation, we apply full HLLD solver for strong enough Alfvén wave. If it is weak, the HLL solver is applied. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index 71ff769..fc67d19 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -2091,10 +2091,9 @@ module schemes 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 +! apply full HLLD solver only if the Alfvén wave is strong enough ! - if (b2 > 0.0d+00) then ! Bₓ² > 0 + if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε ! left and right Alfvén speeds ! @@ -2123,7 +2122,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! the outer left intermediate flux ! @@ -2140,7 +2139,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! the outer right intermediate flux ! @@ -2157,7 +2156,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! vector of the left-going Alfvén wave ! @@ -2172,7 +2171,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! vector of the right-going Alfvén wave ! @@ -2199,7 +2198,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! choose the correct state depending on the speed signs ! @@ -2232,7 +2231,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! choose the correct state depending on the speed signs ! @@ -2264,6 +2263,15 @@ module schemes end if ! one degeneracy + else if (b2 > 0.0d+00) then ! Bₓ² > 0 + +! the Alfvén wave is too weak, so ignore it in order to not introduce any +! numerical instabilities; since Bₓ is not zero, the perpendicular components +! of velocity and magnetic field are continuous; we have only one intermediate +! state, therefore simply apply the HLL solver +! + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml + else ! Bₓ² = 0 ! take the right flux depending on the sign of the advection speed @@ -2279,7 +2287,7 @@ module schemes ui(ibx) = 0.0d+00 ui(iby) = wl(iby) / slmm ui(ibz) = wl(ibz) / slmm - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ! the left intermediate flux ! @@ -2296,7 +2304,7 @@ module schemes ui(ibx) = 0.0d+00 ui(iby) = wr(iby) / srmm ui(ibz) = wr(ibz) / srmm - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ! the right intermediate flux ! From 8dc157620c7596ed7d3317f3e761872507b12f81 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 10 Apr 2017 08:02:05 -0300 Subject: [PATCH 12/13] SCHEMES: Fix isothermal HLLD-M solver for weak Bx. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If the parallel component of magnetic field is too small, the intermediate states might produce numerical instabilities, since they are obtained by the division by a small factor. In order to remove this situation, we apply full HLLD solver for strong enough Alfvén wave. If it is weak, the HLL solver is applied. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 50 ++++++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index fc67d19..e88f153 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -2380,7 +2380,7 @@ module schemes ! integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm - real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr, ca + real(kind=8) :: bx, bp, b2, dn, dnl, dnr, dvl, dvr, ca, ca2 ! local arrays to store the states ! @@ -2471,14 +2471,17 @@ module schemes ! dn = dn / srml -! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to -! the HLL one +! get the square of the Alfvén speed ! - if (b2 > 0.0d+00) then ! Bₓ² > 0 + ca2 = b2 / dn + +! apply full HLLD solver only if the Alfvén wave is strong enough +! + if (ca2 > 1.0d-08 * max(sl**2, sr**2)) then ! Bₓ² > ε ! left and right Alfvén speeds ! - ca = sqrt(b2 / dn) + ca = sqrt(ca2) sml = sm - ca smr = sm + ca @@ -2504,7 +2507,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! the outer left intermediate flux ! @@ -2521,7 +2524,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! the outer right intermediate flux ! @@ -2538,7 +2541,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! vector of the left-going Alfvén wave ! @@ -2553,7 +2556,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! vector of the right-going Alfvén wave ! @@ -2580,7 +2583,7 @@ module schemes 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) + ui(ibp) = ul(ibp,i) ! choose the correct state depending on the speed signs ! @@ -2613,7 +2616,7 @@ module schemes 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) + ui(ibp) = ur(ibp,i) ! choose the correct state depending on the speed signs ! @@ -2645,6 +2648,15 @@ module schemes end if ! one degeneracy + else if (b2 > 0.0d+00) then ! Bₓ² > 0 + +! the Alfvén wave is very weak, so ignore it in order to not introduce any +! numerical instabilities; since Bₓ is not zero, the perpendicular components +! of velocity and magnetic field are continuous; we have only one intermediate +! state, therefore simply apply the HLL solver +! + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml + else ! Bₓ² = 0 ! in the vase of vanishing Bₓ there is no Alfvén wave, density is constant, and @@ -2663,7 +2675,7 @@ module schemes ui(ibx) = 0.0d+00 ui(iby) = wl(iby) / slmm ui(ibz) = wl(ibz) / slmm - ui(ibp) = ql(ibp,i) + ui(ibp) = ul(ibp,i) ! the left intermediate flux ! @@ -2680,7 +2692,7 @@ module schemes ui(ibx) = 0.0d+00 ui(iby) = wr(iby) / srmm ui(ibz) = wr(ibz) / srmm - ui(ibp) = qr(ibp,i) + ui(ibp) = ur(ibp,i) ! the right intermediate flux ! @@ -2688,17 +2700,9 @@ module schemes else ! sm = 0 -! the intermediate flux; since the advection speed is zero, perpendicular -! components do not change +! both states are equal so revert to the HLL flux ! - 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 + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sm = 0 From da2f2a1bf051a3ff31c8897179a6b46ce51e5e7a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 10 Apr 2017 08:03:50 -0300 Subject: [PATCH 13/13] SCHEMES: Fix adiabatic HLLD solver for weak Bx. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If the parallel component of magnetic field is too small, the intermediate states might produce numerical instabilities, since they are obtained by the division by a small factor. In order to remove this situation, we apply full HLLD solver for strong enough Alfvén wave. If it is weak, the HLLC solver is applied. Signed-off-by: Grzegorz Kowal --- src/schemes.F90 | 127 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 88 insertions(+), 39 deletions(-) diff --git a/src/schemes.F90 b/src/schemes.F90 index e88f153..95f4611 100644 --- a/src/schemes.F90 +++ b/src/schemes.F90 @@ -3668,35 +3668,35 @@ module schemes dn = wr(idn) - wl(idn) sm = (wr(imx) - wl(imx)) / dn +! speed differences +! + slmm = sl - sm + srmm = sr - sm + +! left and right state densities +! + dnl = wl(idn) / slmm + dnr = wr(idn) / srmm + ! square of Bₓ, i.e. Bₓ² ! bx = ql(ibx,i) b2 = ql(ibx,i) * qr(ibx,i) -! if there is an Alvén wave, apply the full HLLD solver, otherwise revert to -! the HLLC solver -! - if (b2 > 0.0d+00) then ! Bₓ² > 0 - ! the total pressure, constant across the contact discontinuity and Alfvén waves ! - pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2 + pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2 -! speed differences +! if Alfvén wave is strong enough, apply the full HLLD solver, otherwise +! revert to the HLLC one ! - slmm = sl - sm - srmm = sr - sm - -! left and right state densities -! - dnl = wl(idn) / slmm - dnr = wr(idn) / srmm + if (b2 > 1.0d-08 * max(dnl * sl**2, dnr * sr**2)) then ! Bₓ² > ε ! left and right Alfvén speeds ! - cal = - sqrt(b2 / dnl) - car = sqrt(b2 / dnr) - sml = sm + cal + cal = sqrt(b2 / dnl) + car = sqrt(b2 / dnr) + sml = sm - cal smr = sm + car ! calculate division factors (also used to determine degeneracies) @@ -3814,10 +3814,10 @@ module schemes ! prepare constant primitive variables of the intermediate states ! - dv = car * dnr - cal * dnl + dv = car * dnr + cal * dnl vy = (wcr(imy) - wcl(imy)) / dv vz = (wcr(imz) - wcl(imz)) / dv - dv = car - cal + dv = car + cal by = (wcr(iby) - wcl(iby)) / dv bz = (wcr(ibz) - wcl(ibz)) / dv vb = sm * bx + vy * by + vz * bz @@ -3837,7 +3837,7 @@ module schemes ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) - ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal + ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm) ! the inmost left intermediate flux ! @@ -3855,7 +3855,7 @@ module schemes ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) - ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car + ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm) ! the inmost right intermediate flux ! @@ -3914,7 +3914,7 @@ module schemes ! primitive variables for the inner left intermediate state ! - dv = srmm * dnr - cal * dnl + dv = srmm * dnr + cal * dnl vy = (wr(imy) - wcl(imy)) / dv vz = (wr(imz) - wcl(imz)) / dv dv = sr - sml @@ -3932,7 +3932,7 @@ module schemes ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) - ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal + ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm) ! choose the correct state depending on the sign of contact discontinuity ! advection speed @@ -4013,7 +4013,7 @@ module schemes ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) - ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car + ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm) ! choose the correct state depending on the sign of contact discontinuity ! advection speed @@ -4048,24 +4048,77 @@ module schemes end if ! one degeneracy - else ! Bₓ² = 0 + else if (b2 > 0.0d+00) then ! Bₓ² > 0 -! the total pressure, constant across the contact discontinuity +! constant intermediate state tangential components of velocity and +! magnetic field ! - pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + vy = (wr(imy) - wl(imy)) / dn + vz = (wr(imz) - wl(imz)) / dn + by = (wr(iby) - wl(iby)) / srml + bz = (wr(ibz) - wl(ibz)) / srml + +! scalar product of velocity and magnetic field for the intermediate states +! + vb = sm * bx + vy * by + vz * bz ! separate intermediate states depending on the sign of the advection speed ! if (sm > 0.0d+00) then ! sm > 0 -! the left speed difference +! conservative variables for the left intermediate state ! - slmm = sl - sm + ui(idn) = dnl + ui(imx) = dnl * sm + ui(imy) = dnl * vy + ui(imz) = dnl * vz + ui(ibx) = bx + ui(iby) = by + ui(ibz) = bz + ui(ibp) = ul(ibp,i) + ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm + +! 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) = dnr + ui(imx) = dnr * sm + ui(imy) = dnr * vy + ui(imz) = dnr * vz + ui(ibx) = bx + ui(iby) = by + ui(ibz) = bz + ui(ibp) = ur(ibp,i) + ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm + +! the right intermediate flux +! + f(:,i) = sr * ui(:) - wr(:) + + else ! sm = 0 + +! when Sₘ = 0 all variables are continuous, therefore the flux reduces +! to the HLL one +! + f(:,i) = (sl * wr(:) - sr * wl(:)) / srml + + end if ! sm = 0 + + else ! Bₓ² = 0 + +! separate intermediate states 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) = wl(idn) / slmm - ui(imx) = ui(idn) * sm + ui(idn) = dnl + ui(imx) = dnl * sm ui(imy) = wl(imy) / slmm ui(imz) = wl(imz) / slmm ui(ibx) = 0.0d+00 @@ -4080,14 +4133,10 @@ module schemes else if (sm < 0.0d+00) then ! sm < 0 -! the right speed difference -! - srmm = sr - sm - ! conservative variables for the right intermediate state ! - ui(idn) = wr(idn) / srmm - ui(imx) = ui(idn) * sm + ui(idn) = dnr + ui(imx) = dnr * sm ui(imy) = wr(imy) / srmm ui(imz) = wr(imz) / srmm ui(ibx) = 0.0d+00 @@ -4102,8 +4151,8 @@ module schemes else ! sm = 0 -! intermediate flux is constant across the contact discontinuity, so we end up -! having only one intermediate state, which is equivalent to the HLL solver +! when Sₘ = 0 all variables are continuous, therefore the flux reduces +! to the HLL one ! f(:,i) = (sl * wr(:) - sr * wl(:)) / srml