gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into binaries

Browse Source
This commit is contained in:
Grzegorz Kowal 2017-04-10 08:06:22 -03:00
commit 7bc5807c80
6 changed files with 441 additions and 151 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/gravity.F90
View File

@ -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
!

211
src/interpolations.F90
View File

@ -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
@ -868,7 +889,7 @@ module interpolations
! -------------------------
!
! Subroutine reconstructs both side interfaces of variable using
! multidimentional Gaussian Process method.
! multidimensional Gaussian Process method.
!
! Arguments:
!
@ -1028,12 +1049,194 @@ 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) = 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) = 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) = limiter_minmod(0.5d+00, 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
! update the interpolated states
!
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
#if NDIMS == 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
!-------------------------------------------------------------------------------
!
end subroutine mlp_limiting
!
!===============================================================================
!
! subroutine RECONSTRUCT:
! ----------------------
!

40
src/io.F90
View File

@ -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,24 @@ module io
call start_timer(ios)
#endif /* PROFILE */
! reset the return flag
!
iret = 0
! 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
@ -411,6 +423,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
@ -421,12 +437,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 +500,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 +514,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

319
src/schemes.F90
View File

@ -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,12 +2263,16 @@ module schemes
end if ! one degeneracy
else ! Bₓ² = 0
else if (b2 > 0.0d+00) then ! 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
! 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
!
dn = dn / srml
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
else ! Bₓ² = 0
! take the right flux depending on the sign of the advection speed
!
@ -2277,14 +2280,14 @@ module schemes
! 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
ui(iby) = wl(iby) / slmm
ui(ibz) = wl(ibz) / slmm
ui(ibp) = ql(ibp,i)
ui(ibp) = ul(ibp,i)
! the left intermediate flux
!
@ -2294,14 +2297,14 @@ 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
ui(iby) = wr(iby) / srmm
ui(ibz) = wr(ibz) / srmm
ui(ibp) = qr(ibp,i)
ui(ibp) = ur(ibp,i)
! the right intermediate flux
!
@ -2310,16 +2313,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
@ -2384,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
!
@ -2475,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
@ -2508,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
!
@ -2525,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
!
@ -2542,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
!
@ -2557,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
!
@ -2584,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
!
@ -2617,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
!
@ -2649,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
@ -2667,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
!
@ -2684,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
!
@ -2692,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
@ -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)
@ -3729,7 +3729,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
@ -3755,7 +3755,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
@ -3781,7 +3781,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
@ -3805,7 +3805,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
@ -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
@ -3836,8 +3836,8 @@ module schemes
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = ql(ibp,i)
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / cal
ui(ibp) = ul(ibp,i)
ui(ien) = (wcl(ien) + sm * pt - bx * vb) / (sml - sm)
! the inmost left intermediate flux
!
@ -3854,8 +3854,8 @@ module schemes
ui(ibx) = bx
ui(iby) = by
ui(ibz) = bz
ui(ibp) = qr(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car
ui(ibp) = ur(ibp,i)
ui(ien) = (wcr(ien) + sm * pt - bx * vb) / (smr - sm)
! the inmost right intermediate flux
!
@ -3895,7 +3895,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
@ -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
@ -3922,23 +3922,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) / (sml - sm)
! 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(:)
@ -3947,11 +3947,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
@ -3976,7 +3976,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
@ -4003,23 +4003,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) / (smr - sm)
! 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(:)
@ -4028,11 +4028,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
@ -4048,30 +4048,83 @@ 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(imy) = ui(idn) * ql(ivy,i)
ui(imz) = ui(idn) * ql(ivz,i)
ui(idn) = dnl
ui(imx) = dnl * 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)
ui(ibp) = ul(ibp,i)
ui(ien) = (wl(ien) + sm * pt) / slmm
! the left intermediate flux
@ -4080,20 +4133,16 @@ 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(imy) = ui(idn) * qr(ivy,i)
ui(imz) = ui(idn) * qr(ivz,i)
ui(idn) = dnr
ui(imx) = dnr * 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)
ui(ibp) = ur(ibp,i)
ui(ien) = (wr(ien) + sm * pt) / srmm
! the right intermediate flux
@ -4102,18 +4151,10 @@ module schemes
else ! sm = 0
! intermediate flux is constant across the contact discontinuity and all except
! the parallel momentum flux are zero
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
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

8
src/sources.F90
View File

@ -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

4
src/user_problem.F90
View File

@ -831,6 +831,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
!