Merge branch 'master' into reconnection
This commit is contained in:
commit
58075613e6
@ -236,11 +236,11 @@ module blocks
|
||||
type block_data
|
||||
! pointers to the previous and next data blocks
|
||||
!
|
||||
type(block_data), pointer :: prev, next
|
||||
type(block_data) , pointer :: prev, next
|
||||
|
||||
! a pointer to the associated meta block
|
||||
!
|
||||
type(block_meta), pointer :: meta
|
||||
type(block_meta) , pointer :: meta
|
||||
|
||||
! a pointer to the current conserved variable
|
||||
! array
|
||||
@ -263,6 +263,11 @@ module blocks
|
||||
!
|
||||
real(kind=8), dimension(:,:,:,:) , allocatable :: fx, fy, fz
|
||||
|
||||
! flag indicating if all block cells
|
||||
! are physical
|
||||
!
|
||||
logical :: physical
|
||||
|
||||
end type block_data
|
||||
|
||||
! define the INFO block structure
|
||||
|
||||
@ -1295,7 +1295,7 @@ module equations
|
||||
logical , intent(in) :: ghosts
|
||||
integer , intent(out) :: status
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: i, j, k, s
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1309,11 +1309,8 @@ module equations
|
||||
do k = 1, nn
|
||||
#endif /* NDIMS == 3 */
|
||||
do j = 1, nn
|
||||
|
||||
call cons2prim(uu(:,:,j,k), qq(:,:,j,k), status)
|
||||
|
||||
if (status /= 0) go to 100
|
||||
|
||||
call cons2prim(uu(:,:,j,k), qq(:,:,j,k), s)
|
||||
status = max(status, s)
|
||||
end do
|
||||
#if NDIMS == 3
|
||||
end do
|
||||
@ -1326,11 +1323,8 @@ module equations
|
||||
do k = nb, ne
|
||||
#endif /* NDIMS == 3 */
|
||||
do j = nb, ne
|
||||
|
||||
call cons2prim(uu(:,nb:ne,j,k), qq(:,nb:ne,j,k), status)
|
||||
|
||||
if (status /= 0) go to 100
|
||||
|
||||
call cons2prim(uu(:,nb:ne,j,k), qq(:,nb:ne,j,k), s)
|
||||
status = max(status, s)
|
||||
end do ! j = nb, ne
|
||||
#if NDIMS == 3
|
||||
end do ! k = nb, ne
|
||||
@ -1371,8 +1365,6 @@ module equations
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_primitive_variables
|
||||
@ -1388,14 +1380,15 @@ module equations
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! it - the time step number;
|
||||
! id - the block id where the states are being checked;
|
||||
! qq - the output array of primitive variables;
|
||||
! uu - the input array of conservative variables;
|
||||
! it - the time step number;
|
||||
! id - the block id where the states are being checked;
|
||||
! qq - the output array of primitive variables;
|
||||
! uu - the input array of conservative variables;
|
||||
! status - the call status;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine correct_unphysical_states(it, id, qq, uu)
|
||||
subroutine correct_unphysical_states(it, id, qq, uu, status)
|
||||
|
||||
use coordinates, only : nn => bcells
|
||||
use helpers , only : print_message
|
||||
@ -1405,6 +1398,7 @@ module equations
|
||||
integer(kind=4) , intent(in) :: it, id
|
||||
real(kind=8), dimension(:,:,:,:), intent(inout) :: qq
|
||||
real(kind=8), dimension(:,:,:,:), intent(inout) :: uu
|
||||
integer , intent(out) :: status
|
||||
|
||||
character(len=255) :: msg, sfmt
|
||||
character(len=16) :: sit, sid, snc
|
||||
@ -1425,7 +1419,8 @@ module equations
|
||||
integer , dimension(:,:), allocatable :: idx
|
||||
real(kind=8), dimension(:,:), allocatable :: q, u
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::correct_unphysical_states()'
|
||||
character(len=*), parameter :: loc = &
|
||||
'EQUATIONS::correct_unphysical_states()'
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1442,37 +1437,31 @@ module equations
|
||||
ku = 1
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! search for negative density or pressure
|
||||
!
|
||||
physical(:,:,:) = qq(idn,:,:,:) > 0.0d+00
|
||||
if (ipr > 0) then
|
||||
physical(:,:,:) = physical(:,:,:) .and. qq(ipr,:,:,:) > 0.0d+00
|
||||
physical(:,:,:) = qq(idn,:,:,:) > 0.0d+00 .and. qq(ipr,:,:,:) > 0.0d+00
|
||||
else
|
||||
physical(:,:,:) = qq(idn,:,:,:) > 0.0d+00
|
||||
end if
|
||||
|
||||
! apply averaging for unphysical states
|
||||
!
|
||||
if (.not. all(physical)) then
|
||||
nc = count(.not. physical)
|
||||
|
||||
! count unphysical cells
|
||||
!
|
||||
nc = count(.not. physical)
|
||||
if (nc > 0) then
|
||||
|
||||
! inform about the encountered unphysical states
|
||||
!
|
||||
write(sit,'(i15)') it
|
||||
write(sid,'(i15)') id
|
||||
write(snc,'(i15)') nc
|
||||
sfmt = '("Unphysical cells in block ID:",a," (",a," counted)' &
|
||||
sfmt = '("Unphysical cells in block ID:",a," (",a," counted)' &
|
||||
// ' at time step ",a,".")'
|
||||
write(msg,sfmt) trim(adjustl(sid)), trim(adjustl(snc)), trim(adjustl(sit))
|
||||
call print_message(loc, msg)
|
||||
|
||||
! allocate temporary vectors for unphysical states
|
||||
!
|
||||
allocate(q(nv,nc), u(nv,nc), idx(3,nc))
|
||||
allocate(q(nv,nc), u(nv,nc), idx(3,nc), stat=status)
|
||||
if (status /= 0) then
|
||||
call print_message(loc, &
|
||||
"Could not allocate vectors for unphysical cells!")
|
||||
return
|
||||
end if
|
||||
|
||||
! iterate over block cells
|
||||
!
|
||||
n = 0
|
||||
#if NDIMS == 3
|
||||
do k = 1, nn
|
||||
@ -1480,13 +1469,11 @@ module equations
|
||||
do j = 1, nn
|
||||
do i = 1, nn
|
||||
|
||||
! fix unphysical states
|
||||
!
|
||||
if (.not. physical(i,j,k)) then
|
||||
|
||||
n = n + 1
|
||||
|
||||
idx(:,n) = (/ i, j, k /)
|
||||
idx(:,n) = [ i, j, k ]
|
||||
|
||||
! increase the region until we find at least three physical cells, but no more
|
||||
! than 4 cells away
|
||||
@ -1542,26 +1529,22 @@ module equations
|
||||
msg = "Applying lower bounds for positive variables."
|
||||
call print_message(loc, msg)
|
||||
|
||||
q(:,n) = qq(:,i,j,k)
|
||||
q(idn ,n) = max(dmin, qq(idn,i,j,k))
|
||||
q( : ,n) = qq( : ,i,j,k)
|
||||
q(idn,n) = max(dmin, qq(idn,i,j,k))
|
||||
if (ipr > 0) q(ipr,n) = max(pmin, qq(ipr,i,j,k))
|
||||
|
||||
end if ! not sufficient number of physical cells for averaging
|
||||
|
||||
end if ! not physical
|
||||
end if
|
||||
|
||||
end do ! i = 1, nn
|
||||
end do ! j = 1, nn
|
||||
end do
|
||||
end do
|
||||
#if NDIMS == 3
|
||||
end do ! k = 1, nn
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! convert the vector of primitive variables to conservative ones
|
||||
!
|
||||
call prim2cons(q(:,1:nc), u(:,1:nc), .true.)
|
||||
|
||||
! update block variables
|
||||
!
|
||||
do n = 1, nc
|
||||
i = idx(1,n)
|
||||
j = idx(2,n)
|
||||
@ -1571,11 +1554,14 @@ module equations
|
||||
uu(:,i,j,k) = u(:,n)
|
||||
end do
|
||||
|
||||
! deallocate temporary vectors
|
||||
!
|
||||
deallocate(q, u, idx)
|
||||
deallocate(q, u, idx, stat=status)
|
||||
if (status /= 0) then
|
||||
call print_message(loc, &
|
||||
"Could not deallocate vectors of unphysical cells!")
|
||||
return
|
||||
end if
|
||||
|
||||
end if ! there are unphysical cells
|
||||
end if
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1671,40 +1657,48 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_hd_iso(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: u
|
||||
real(kind=8), dimension(:,:), intent(out) :: q
|
||||
integer , intent(out) :: s
|
||||
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
|
||||
#ifdef DEBUG
|
||||
character(len=256) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_hd_iso()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
s = 0
|
||||
|
||||
do i = 1, size(u,2)
|
||||
|
||||
if (u(idn,i) > 0.0d+00) then
|
||||
q(idn,i) = u(idn,i)
|
||||
q(ivx,i) = u(imx,i) / u(idn,i)
|
||||
q(ivy,i) = u(imy,i) / u(idn,i)
|
||||
q(ivz,i) = u(imz,i) / u(idn,i)
|
||||
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,4(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_hd_iso
|
||||
@ -1917,21 +1911,30 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_hd_adi(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: u
|
||||
real(kind=8), dimension(:,:), intent(out) :: q
|
||||
integer , intent(out) :: s
|
||||
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
real(kind=8) :: ek, ei
|
||||
|
||||
#ifdef DEBUG
|
||||
character(len=256) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_hd_adi()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
s = 0
|
||||
|
||||
do i = 1, size(u,2)
|
||||
|
||||
if (u(idn,i) > 0.0d+00) then
|
||||
q(idn,i) = u(idn,i)
|
||||
q(ivx,i) = u(imx,i) / u(idn,i)
|
||||
@ -1942,23 +1945,29 @@ module equations
|
||||
if (ei > 0.0d+00) then
|
||||
q(ipr,i) = gammam1 * ei
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(ipr,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, "Conversion to physical primitive state " // &
|
||||
"resulted in negative pressure!")
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_hd_adi
|
||||
@ -2176,20 +2185,29 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_mhd_iso(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: u
|
||||
real(kind=8), dimension(:,:), intent(out) :: q
|
||||
integer , intent(out) :: s
|
||||
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
|
||||
#ifdef DEBUG
|
||||
character(len=256) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_mhd_iso()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
s = 0
|
||||
|
||||
do i = 1, size(u,2)
|
||||
|
||||
if (u(idn,i) > 0.0d+00) then
|
||||
q(idn,i) = u(idn,i)
|
||||
q(ivx,i) = u(imx,i) / u(idn,i)
|
||||
@ -2199,21 +2217,20 @@ module equations
|
||||
q(iby,i) = u(iby,i)
|
||||
q(ibz,i) = u(ibz,i)
|
||||
q(ibp,i) = u(ibp,i)
|
||||
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,8(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_mhd_iso
|
||||
@ -2449,21 +2466,30 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_mhd_adi(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: u
|
||||
real(kind=8), dimension(:,:), intent(out) :: q
|
||||
integer , intent(out) :: s
|
||||
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
real(kind=8) :: ei, ek, em, ep
|
||||
|
||||
#ifdef DEBUG
|
||||
character(len=256) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_mhd_adi()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
s = 0
|
||||
|
||||
do i = 1, size(u,2)
|
||||
|
||||
if (u(idn,i) > 0.0d+00) then
|
||||
q(idn,i) = u(idn,i)
|
||||
q(ivx,i) = u(imx,i) / u(idn,i)
|
||||
@ -2480,24 +2506,29 @@ module equations
|
||||
if (ei > 0.0d+00) then
|
||||
q(ipr,i) = gammam1 * ei
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(ipr,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, "Conversion to physical primitive state " // &
|
||||
"resulted in negative pressure!")
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
go to 100
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_mhd_adi
|
||||
@ -2748,7 +2779,9 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_srhd_adi(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
@ -2757,13 +2790,15 @@ module equations
|
||||
integer , intent(out) :: s
|
||||
|
||||
logical :: info
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
real(kind=8) :: mm, bb, mb, en, dn
|
||||
real(kind=8) :: w , vv, vm, vs
|
||||
|
||||
character(len=80) :: msg
|
||||
#ifdef DEBUG
|
||||
character(len=256) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -2771,58 +2806,51 @@ module equations
|
||||
|
||||
do i = 1, size(u,2)
|
||||
|
||||
! prepare variables which do not change during the Newton-Ralphson iterations
|
||||
!
|
||||
mm = sum(u(imx:imz,i) * u(imx:imz,i))
|
||||
en = u(ien,i) + u(idn,i)
|
||||
dn = u(idn,i)
|
||||
|
||||
! find the exact W using the Newton-Ralphson interative method
|
||||
!
|
||||
call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
|
||||
|
||||
! if info is .true., the solution was found
|
||||
!
|
||||
if (info) then
|
||||
|
||||
! prepare coefficients
|
||||
!
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
|
||||
! calculate the primitive variables
|
||||
!
|
||||
q(idn,i) = dn * vs
|
||||
q(ivx,i) = u(imx,i) / w
|
||||
q(ivy,i) = u(imy,i) / w
|
||||
q(ivz,i) = u(imz,i) / w
|
||||
q(ipr,i) = w - en
|
||||
|
||||
else ! cannot find physical solution
|
||||
if (q(ipr,i) <= 0.0d+00) then
|
||||
s = 1
|
||||
q(ipr,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, "Conversion to physical primitive state " // &
|
||||
"resulted in negative pressure!")
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
write(msg,"(a,3(1x,1es24.16e3))") "D, |m|², E = ", dn, mm, en
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
call print_message(loc, "Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(:,i)
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
write(msg,"(a,3(1x,1es24.16e3))") "D, |m|², E = ", dn, mm, en
|
||||
call print_message(loc, msg)
|
||||
|
||||
s = 1
|
||||
go to 100
|
||||
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_srhd_adi
|
||||
@ -3758,7 +3786,9 @@ module equations
|
||||
!
|
||||
subroutine cons2prim_srmhd_adi(u, q, s)
|
||||
|
||||
#ifdef DEBUG
|
||||
use helpers, only : print_message
|
||||
#endif /* DEBUG */
|
||||
|
||||
implicit none
|
||||
|
||||
@ -3767,48 +3797,37 @@ module equations
|
||||
integer , intent(out) :: s
|
||||
|
||||
logical :: info
|
||||
integer :: i, p
|
||||
integer :: i
|
||||
real(kind=8) :: mm, mb, bb, en, dn
|
||||
real(kind=8) :: w, wt, vv, vm, vs, fc
|
||||
|
||||
#ifdef DEBUG
|
||||
character(len=80) :: msg
|
||||
|
||||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srmhd_adi()'
|
||||
#endif /* DEBUG */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
s = 0
|
||||
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(u,2)
|
||||
|
||||
! prepare variables which do not change during the Newton-Ralphson iterations
|
||||
! (|B|², |M|² and B.M and their multiplications)
|
||||
!
|
||||
mm = sum(u(imx:imz,i) * u(imx:imz,i))
|
||||
mb = sum(u(imx:imz,i) * u(ibx:ibz,i))
|
||||
bb = sum(u(ibx:ibz,i) * u(ibx:ibz,i))
|
||||
en = u(ien,i) + u(idn,i)
|
||||
dn = u(idn,i)
|
||||
|
||||
! find the exact W using an Newton-Ralphson interative method
|
||||
!
|
||||
call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
|
||||
|
||||
! if info is .true., the solution was found
|
||||
!
|
||||
if (info) then
|
||||
|
||||
! prepare coefficients
|
||||
!
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
wt = w + bb
|
||||
fc = mb / w
|
||||
|
||||
! calculate the primitive variables
|
||||
!
|
||||
q(idn,i) = dn * vs
|
||||
q(ivx,i) = (u(imx,i) + fc * u(ibx,i)) / wt
|
||||
q(ivy,i) = (u(imy,i) + fc * u(iby,i)) / wt
|
||||
@ -3819,48 +3838,36 @@ module equations
|
||||
q(ibp,i) = u(ibp,i)
|
||||
q(ipr,i) = w - en + 0.5d+00 * (bb + (bb * mm - mb * mb) / wt**2)
|
||||
|
||||
! check if the pressure is positive, if not, print a warning and replace it
|
||||
! with the minimum allowed value pmin
|
||||
!
|
||||
if (q(ipr,i) <= 0.0d+00) then
|
||||
|
||||
s = 1
|
||||
q(ipr,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, "Conversion to physical primitive state " // &
|
||||
"resulted in negative pressure!")
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(:,i)
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
write(msg,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = ", &
|
||||
write(msg,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = ", &
|
||||
dn, mm, mb, bb, en, w
|
||||
call print_message(loc, msg)
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
|
||||
s = 1
|
||||
go to 100
|
||||
|
||||
end if ! p <= 0
|
||||
|
||||
else ! unphysical state
|
||||
|
||||
call print_message(loc, "Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(:,i)
|
||||
if (ns > 0) q(isl:isu,i) = u(isl:isu,i) / u(idn,i)
|
||||
else
|
||||
s = 1
|
||||
q(idn,i) = 0.0d+00
|
||||
#ifdef DEBUG
|
||||
call print_message(loc, &
|
||||
"Conversion to physical primitive state failed!")
|
||||
write(msg,"(a,9(1x,1es24.16e3))") "U = ", u(1:nf,i)
|
||||
call print_message(loc, msg)
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ", &
|
||||
write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ", &
|
||||
dn, mm, mb, bb, en
|
||||
call print_message(loc, msg)
|
||||
|
||||
s = 1
|
||||
go to 100
|
||||
|
||||
end if ! unphysical state
|
||||
|
||||
#endif /* DEBUG */
|
||||
end if
|
||||
end do
|
||||
|
||||
if (ns > 0) then
|
||||
do p = isl, isu
|
||||
q(p,:) = u(p,:) / u(idn,:)
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine cons2prim_srmhd_adi
|
||||
|
||||
@ -4077,8 +4077,7 @@ module evolution
|
||||
!
|
||||
subroutine update_variables(tm, dtm, status)
|
||||
|
||||
use blocks , only : block_meta, list_meta
|
||||
use blocks , only : block_data, list_data, data_blocks
|
||||
use blocks , only : block_meta, block_data, data_blocks, list_meta
|
||||
use blocks , only : set_neighbors_update, get_dblocks
|
||||
use boundaries, only : boundary_variables
|
||||
use equations , only : update_primitive_variables
|
||||
@ -4105,72 +4104,71 @@ module evolution
|
||||
n = get_dblocks()
|
||||
|
||||
call boundary_variables(tm, dtm, status)
|
||||
if (status /= 0) go to 100
|
||||
|
||||
!$omp parallel do default(shared) private(pdata,s)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
|
||||
if (pdata%meta%update .or. pdata%meta%boundary) then
|
||||
call update_primitive_variables(pdata%u, pdata%q, .true., s)
|
||||
!$omp critical
|
||||
if (s /= 0) status = 1
|
||||
!$omp end critical
|
||||
end if
|
||||
|
||||
end do
|
||||
!$omp end parallel do
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(status)
|
||||
#endif /* MPI */
|
||||
if (status /= 0) go to 100
|
||||
|
||||
!$omp parallel do default(shared) private(pdata)
|
||||
!$omp parallel do default(shared) private(pdata,pmeta,s)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (pdata%meta%update) &
|
||||
call update_shapes(pdata, tm, dtm)
|
||||
if (pmeta%update .or. pmeta%boundary) then
|
||||
call update_primitive_variables(pdata%u, pdata%q, .true., s)
|
||||
|
||||
pdata%physical = s == 0
|
||||
else
|
||||
pdata%physical = .true.
|
||||
end if
|
||||
end do
|
||||
!$omp end parallel do
|
||||
|
||||
if (fix_unphysical_cells) then
|
||||
!$omp parallel do default(shared) private(pdata,pmeta,s)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (.not. pdata%physical) then
|
||||
call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u, s)
|
||||
!$omp critical
|
||||
if (s /= 0) status = 1
|
||||
!$omp end critical
|
||||
end if
|
||||
|
||||
end do
|
||||
!$omp end parallel do
|
||||
else
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
if (.not. pdata%physical) status = 1
|
||||
end do
|
||||
end if
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(status)
|
||||
#endif /* MPI */
|
||||
if (status /= 0) go to 100
|
||||
|
||||
!$omp parallel do default(shared) private(pdata,pmeta)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (pmeta%update .or. pmeta%boundary) call update_shapes(pdata, tm, dtm)
|
||||
end do
|
||||
!$omp end parallel do
|
||||
|
||||
100 continue
|
||||
|
||||
pmeta => list_meta
|
||||
do while (associated(pmeta))
|
||||
pmeta%boundary = .false.
|
||||
pmeta => pmeta%next
|
||||
end do
|
||||
|
||||
if (fix_unphysical_cells) then
|
||||
|
||||
! if an unphysical cell appeared in a block while updating its primitive
|
||||
! variables it could be propagated to its neighbors through boundary update;
|
||||
! mark all neighbors of such a block to be verified and corrected for
|
||||
! unphysical cells too
|
||||
!
|
||||
pmeta => list_meta
|
||||
do while (associated(pmeta))
|
||||
|
||||
if (pmeta%update) call set_neighbors_update(pmeta)
|
||||
|
||||
pmeta => pmeta%next
|
||||
end do
|
||||
|
||||
! verify and correct, if necessary, unphysical cells in recently updated blocks
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (pmeta%update) &
|
||||
call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
end if
|
||||
|
||||
100 continue
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_variables
|
||||
@ -4334,7 +4332,7 @@ module evolution
|
||||
use equations , only : nv, pvars, cvars
|
||||
use ieee_arithmetic, only : ieee_is_nan
|
||||
|
||||
use blocks , only : block_meta, list_meta
|
||||
use blocks , only : block_meta
|
||||
use blocks , only : block_data, list_data
|
||||
|
||||
implicit none
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user