Merge branch 'master' into gem-reconnection-challenge
This commit is contained in:
commit
a436a03e9a
@ -2,7 +2,7 @@
|
||||
# **The AMUN Code**
|
||||
## Copyright (C) 2008-2023 Grzegorz Kowal
|
||||
|
||||
[](https://ampere-orbis.nsupdate.info/gkowal/amun-code)
|
||||
[](https://ampere.clarokowal.com/gkowal/amun-code)
|
||||
|
||||
AMUN is a parallel code to perform numerical simulations in fluid approximation
|
||||
on uniform or non-uniform (adaptive) meshes. The goal in developing this code is
|
||||
|
||||
@ -53,13 +53,19 @@ module equations
|
||||
integer(kind=4), save :: nf = 0
|
||||
integer(kind=4), save :: ns = 0
|
||||
|
||||
! subroutine interfaces
|
||||
!
|
||||
interface prim2cons
|
||||
module procedure prim2cons_cell
|
||||
module procedure prim2cons_stencil
|
||||
end interface
|
||||
|
||||
! interfaces for procedure pointers
|
||||
!
|
||||
interface
|
||||
subroutine prim2cons_iface(q, u, s)
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
subroutine prim2cons_iface(q, u)
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
end subroutine
|
||||
subroutine cons2prim_iface(u, q, s)
|
||||
real(kind=8), dimension(:,:), intent(in) :: u
|
||||
@ -84,12 +90,12 @@ module equations
|
||||
|
||||
! pointers to the conversion procedures
|
||||
!
|
||||
procedure(prim2cons_iface) , pointer, save :: prim2cons => null()
|
||||
procedure(cons2prim_iface) , pointer, save :: cons2prim => null()
|
||||
procedure(prim2cons_iface) , pointer, save :: prim2cons_ptr => null()
|
||||
procedure(cons2prim_iface) , pointer, save :: cons2prim => null()
|
||||
|
||||
! pointer to the flux procedure
|
||||
!
|
||||
procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null()
|
||||
procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null()
|
||||
|
||||
procedure(get_maximum_speeds_iface), pointer, save :: &
|
||||
get_maximum_speeds => null()
|
||||
@ -205,13 +211,7 @@ module equations
|
||||
! declare public variables and subroutines
|
||||
!
|
||||
public :: initialize_equations, finalize_equations, print_equations
|
||||
public :: cons2prim, prim2cons, fluxspeed
|
||||
public :: prim2cons_hd_iso, fluxspeed_hd_iso
|
||||
public :: prim2cons_hd_adi, fluxspeed_hd_adi
|
||||
public :: prim2cons_mhd_iso, fluxspeed_mhd_iso
|
||||
public :: prim2cons_mhd_adi, fluxspeed_mhd_adi
|
||||
public :: prim2cons_srhd_adi, fluxspeed_srhd_adi
|
||||
public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi
|
||||
public :: prim2cons, cons2prim, fluxspeed
|
||||
public :: get_maximum_speeds
|
||||
public :: update_primitive_variables
|
||||
public :: fix_unphysical_cells, correct_unphysical_states
|
||||
@ -350,7 +350,7 @@ module equations
|
||||
|
||||
! set pointers to subroutines
|
||||
!
|
||||
prim2cons => prim2cons_hd_iso
|
||||
prim2cons_ptr => prim2cons_hd_iso
|
||||
cons2prim => cons2prim_hd_iso
|
||||
fluxspeed => fluxspeed_hd_iso
|
||||
get_maximum_speeds => get_maximum_speeds_hd_iso
|
||||
@ -389,7 +389,7 @@ module equations
|
||||
|
||||
! set pointers to subroutines
|
||||
!
|
||||
prim2cons => prim2cons_hd_adi
|
||||
prim2cons_ptr => prim2cons_hd_adi
|
||||
cons2prim => cons2prim_hd_adi
|
||||
fluxspeed => fluxspeed_hd_adi
|
||||
get_maximum_speeds => get_maximum_speeds_hd_adi
|
||||
@ -514,7 +514,7 @@ module equations
|
||||
|
||||
! set pointers to the subroutines
|
||||
!
|
||||
prim2cons => prim2cons_mhd_iso
|
||||
prim2cons_ptr => prim2cons_mhd_iso
|
||||
cons2prim => cons2prim_mhd_iso
|
||||
fluxspeed => fluxspeed_mhd_iso
|
||||
get_maximum_speeds => get_maximum_speeds_mhd_iso
|
||||
@ -560,7 +560,7 @@ module equations
|
||||
|
||||
! set pointers to subroutines
|
||||
!
|
||||
prim2cons => prim2cons_mhd_adi
|
||||
prim2cons_ptr => prim2cons_mhd_adi
|
||||
cons2prim => cons2prim_mhd_adi
|
||||
fluxspeed => fluxspeed_mhd_adi
|
||||
get_maximum_speeds => get_maximum_speeds_mhd_adi
|
||||
@ -688,7 +688,7 @@ module equations
|
||||
|
||||
! set pointers to subroutines
|
||||
!
|
||||
prim2cons => prim2cons_srhd_adi
|
||||
prim2cons_ptr => prim2cons_srhd_adi
|
||||
cons2prim => cons2prim_srhd_adi
|
||||
fluxspeed => fluxspeed_srhd_adi
|
||||
get_maximum_speeds => get_maximum_speeds_srhd_adi
|
||||
@ -869,7 +869,7 @@ module equations
|
||||
|
||||
! set pointers to subroutines
|
||||
!
|
||||
prim2cons => prim2cons_srmhd_adi
|
||||
prim2cons_ptr => prim2cons_srmhd_adi
|
||||
cons2prim => cons2prim_srmhd_adi
|
||||
fluxspeed => fluxspeed_srmhd_adi
|
||||
get_maximum_speeds => get_maximum_speeds_srmhd_adi
|
||||
@ -1166,9 +1166,7 @@ module equations
|
||||
!
|
||||
status = 0
|
||||
|
||||
! release the procedure pointers
|
||||
!
|
||||
nullify(prim2cons)
|
||||
nullify(prim2cons_ptr)
|
||||
nullify(cons2prim)
|
||||
nullify(fluxspeed)
|
||||
|
||||
@ -1483,6 +1481,7 @@ module equations
|
||||
|
||||
n = n + 1
|
||||
|
||||
#ifdef DEBUG
|
||||
#if NDIMS == 3
|
||||
sfmt = '("n: ",i0,", [i,j,k]: ",3(1x,i0),", [x,y,z]: ",3es12.4)'
|
||||
write(msg,sfmt) n, i, j, k, x(i), y(j), z(k)
|
||||
@ -1491,6 +1490,7 @@ module equations
|
||||
write(msg,sfmt) n, i, j, x(i), y(j)
|
||||
#endif /* NDIMS == 3 */
|
||||
call print_message(msg)
|
||||
#endif /* DEBUG */
|
||||
|
||||
idx(:,n) = [ i, j, k ]
|
||||
|
||||
@ -1560,7 +1560,7 @@ module equations
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
call prim2cons(q(:,1:nc), u(:,1:nc), .true.)
|
||||
call prim2cons_stencil(q(:,1:nc), u(:,1:nc), .true.)
|
||||
|
||||
do n = 1, nc
|
||||
i = idx(1,n)
|
||||
@ -1590,6 +1590,79 @@ module equations
|
||||
!!
|
||||
!===============================================================================
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine PRIM2CONS_CELL:
|
||||
! -------------------------
|
||||
!
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_cell(q, u)
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
call prim2cons_ptr(q, u)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine prim2cons_cell
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine PRIM2CONS_STENCIL:
|
||||
! ----------------------------
|
||||
!
|
||||
! Subroutine converts primitive to conservative variables along a stencil.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_stencil(q, u, s)
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
|
||||
integer :: i, p
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
call prim2cons_ptr(q(:,i), u(:,i))
|
||||
end do
|
||||
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine prim2cons_stencil
|
||||
!
|
||||
!*******************************************************************************
|
||||
!
|
||||
! ISOTHERMAL HYDRODYNAMIC EQUATIONS
|
||||
@ -1601,56 +1674,28 @@ module equations
|
||||
! subroutine PRIM2CONS_HD_ISO:
|
||||
! ---------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_hd_iso(q, u, s)
|
||||
subroutine prim2cons_hd_iso(q, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
|
||||
u(idn,i) = q(idn,i)
|
||||
u(imx,i) = q(idn,i) * q(ivx,i)
|
||||
u(imy,i) = q(idn,i) * q(ivy,i)
|
||||
u(imz,i) = q(idn,i) * q(ivz,i)
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn)
|
||||
u(imx) = q(idn) * q(ivx)
|
||||
u(imy) = q(idn) * q(ivy)
|
||||
u(imz) = q(idn) * q(ivz)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1850,61 +1895,33 @@ module equations
|
||||
! subroutine PRIM2CONS_HD_ADI:
|
||||
! ---------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_hd_adi(q, u, s)
|
||||
subroutine prim2cons_hd_adi(q, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8) :: ek, ei
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
|
||||
u(idn,i) = q(idn,i)
|
||||
u(imx,i) = q(idn,i) * q(ivx,i)
|
||||
u(imy,i) = q(idn,i) * q(ivy,i)
|
||||
u(imz,i) = q(idn,i) * q(ivz,i)
|
||||
ek = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i) &
|
||||
+ u(imz,i) * q(ivz,i))
|
||||
ei = gammam1i * q(ipr,i)
|
||||
u(ien,i) = ei + ek
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn)
|
||||
u(imx) = q(idn) * q(ivx)
|
||||
u(imy) = q(idn) * q(ivy)
|
||||
u(imz) = q(idn) * q(ivz)
|
||||
ek = 5.0d-01 * (u(imx) * q(ivx) + u(imy) * q(ivy) + u(imz) * q(ivz))
|
||||
ei = gammam1i * q(ipr)
|
||||
u(ien) = ei + ek
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -2125,60 +2142,32 @@ module equations
|
||||
! subroutine PRIM2CONS_MHD_ISO:
|
||||
! ----------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_mhd_iso(q, u, s)
|
||||
subroutine prim2cons_mhd_iso(q, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
|
||||
u(idn,i) = q(idn,i)
|
||||
u(imx,i) = q(idn,i) * q(ivx,i)
|
||||
u(imy,i) = q(idn,i) * q(ivy,i)
|
||||
u(imz,i) = q(idn,i) * q(ivz,i)
|
||||
u(ibx,i) = q(ibx,i)
|
||||
u(iby,i) = q(iby,i)
|
||||
u(ibz,i) = q(ibz,i)
|
||||
u(ibp,i) = q(ibp,i)
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn)
|
||||
u(imx) = q(idn) * q(ivx)
|
||||
u(imy) = q(idn) * q(ivy)
|
||||
u(imz) = q(idn) * q(ivz)
|
||||
u(ibx) = q(ibx)
|
||||
u(iby) = q(iby)
|
||||
u(ibz) = q(ibz)
|
||||
u(ibp) = q(ibp)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -2410,56 +2399,39 @@ module equations
|
||||
! subroutine PRIM2CONS_MHD_ADI:
|
||||
! ----------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_mhd_adi(q, u, s)
|
||||
subroutine prim2cons_mhd_adi(q, u)
|
||||
|
||||
implicit none
|
||||
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
|
||||
integer :: i, p
|
||||
real(kind=8) :: ei, ek, em, ep
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
|
||||
u(idn,i) = q(idn,i)
|
||||
u(imx,i) = q(idn,i) * q(ivx,i)
|
||||
u(imy,i) = q(idn,i) * q(ivy,i)
|
||||
u(imz,i) = q(idn,i) * q(ivz,i)
|
||||
u(ibx,i) = q(ibx,i)
|
||||
u(iby,i) = q(iby,i)
|
||||
u(ibz,i) = q(ibz,i)
|
||||
u(ibp,i) = q(ibp,i)
|
||||
ei = gammam1i * q(ipr,i)
|
||||
ek = sum(u(imx:imz,i) * q(ivx:ivz,i))
|
||||
em = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
ep = q(ibp,i) * q(ibp,i)
|
||||
u(ien,i) = ei + 5.0d-01 * (ek + em + ep)
|
||||
|
||||
end do
|
||||
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn)
|
||||
u(imx) = q(idn) * q(ivx)
|
||||
u(imy) = q(idn) * q(ivy)
|
||||
u(imz) = q(idn) * q(ivz)
|
||||
u(ibx) = q(ibx)
|
||||
u(iby) = q(iby)
|
||||
u(ibz) = q(ibz)
|
||||
u(ibp) = q(ibp)
|
||||
ei = gammam1i * q(ipr)
|
||||
ek = sum(u(imx:imz) * q(ivx:ivz))
|
||||
em = sum(q(ibx:ibz) * q(ibx:ibz))
|
||||
ep = q(ibp) * q(ibp)
|
||||
u(ien) = ei + 5.0d-01 * (ek + em + ep)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -2712,67 +2684,36 @@ module equations
|
||||
! subroutine PRIM2CONS_SRHD_ADI:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_srhd_adi(q, u, s)
|
||||
subroutine prim2cons_srhd_adi(q, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8) :: vv, vm, vs, ww
|
||||
!
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
vv = sum(q(ivx:ivz) * q(ivx:ivz))
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
ww = (q(idn) + q(ipr) / gammaxi) / vm
|
||||
|
||||
! calculate the square of velocity, the Lorentz factor and specific enthalpy
|
||||
!
|
||||
vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
ww = (q(idn,i) + q(ipr,i) / gammaxi) / vm
|
||||
|
||||
! calculate conservative variables
|
||||
!
|
||||
u(idn,i) = q(idn,i) / vs
|
||||
u(imx,i) = ww * q(ivx,i)
|
||||
u(imy,i) = ww * q(ivy,i)
|
||||
u(imz,i) = ww * q(ivz,i)
|
||||
u(ien,i) = ww - q(ipr,i) - u(idn,i)
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn) / vs
|
||||
u(imx) = ww * q(ivx)
|
||||
u(imy) = ww * q(ivy)
|
||||
u(imz) = ww * q(ivz)
|
||||
u(ien) = ww - q(ipr) - u(idn)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -3710,76 +3651,43 @@ module equations
|
||||
! subroutine PRIM2CONS_SRMHD_ADI:
|
||||
! ------------------------------
|
||||
!
|
||||
! Subroutine converts primitive variables to their corresponding
|
||||
! conservative representation.
|
||||
! Subroutine converts primitive to conservative variables.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! q - the input array of primitive variables;
|
||||
! u - the output array of conservative variables;
|
||||
! s - an optional flag indicating that passive scalars have
|
||||
! to be calculated too;
|
||||
! q - primitive variables (input);
|
||||
! u - conservative variables (output);
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine prim2cons_srmhd_adi(q, u, s)
|
||||
subroutine prim2cons_srmhd_adi(q, u)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input/output arguments
|
||||
!
|
||||
real(kind=8), dimension(:,:), intent(in) :: q
|
||||
real(kind=8), dimension(:,:), intent(out) :: u
|
||||
logical , optional , intent(in) :: s
|
||||
real(kind=8), dimension(:), intent(in) :: q
|
||||
real(kind=8), dimension(:), intent(out) :: u
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, p
|
||||
real(kind=8) :: vv, bb, vb, vm, vs, ww, wt
|
||||
!
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! iterate over all positions
|
||||
!
|
||||
do i = 1, size(q,2)
|
||||
vv = sum(q(ivx:ivz) * q(ivx:ivz))
|
||||
bb = sum(q(ibx:ibz) * q(ibx:ibz))
|
||||
vb = sum(q(ivx:ivz) * q(ibx:ibz))
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
ww = (q(idn) + q(ipr) / gammaxi) / vm
|
||||
wt = ww + bb
|
||||
|
||||
! calculate the square of velocity, the quare of magnetic field, the scalar
|
||||
! product of velocity and magnetic field, the Lorentz factor, specific and
|
||||
! total enthalpies
|
||||
!
|
||||
vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
|
||||
bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||||
vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
|
||||
vm = 1.0d+00 - vv
|
||||
vs = sqrt(vm)
|
||||
ww = (q(idn,i) + q(ipr,i) / gammaxi) / vm
|
||||
wt = ww + bb
|
||||
|
||||
! calculate conservative variables
|
||||
!
|
||||
u(idn,i) = q(idn,i) / vs
|
||||
u(imx,i) = wt * q(ivx,i) - vb * q(ibx,i)
|
||||
u(imy,i) = wt * q(ivy,i) - vb * q(iby,i)
|
||||
u(imz,i) = wt * q(ivz,i) - vb * q(ibz,i)
|
||||
u(ibx,i) = q(ibx,i)
|
||||
u(iby,i) = q(iby,i)
|
||||
u(ibz,i) = q(ibz,i)
|
||||
u(ibp,i) = q(ibp,i)
|
||||
u(ien,i) = wt - q(ipr,i) - u(idn,i) - 0.5d+00 * (vm * bb + vb * vb)
|
||||
|
||||
end do
|
||||
|
||||
! update primitive passive scalars
|
||||
!
|
||||
if (ns > 0 .and. present(s)) then
|
||||
if (s) then
|
||||
do p = isl, isu
|
||||
u(p,:) = q(p,:) * u(idn,:)
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
u(idn) = q(idn) / vs
|
||||
u(imx) = wt * q(ivx) - vb * q(ibx)
|
||||
u(imy) = wt * q(ivy) - vb * q(iby)
|
||||
u(imz) = wt * q(ivz) - vb * q(ibz)
|
||||
u(ibx) = q(ibx)
|
||||
u(iby) = q(iby)
|
||||
u(ibz) = q(ibz)
|
||||
u(ibp) = q(ibp)
|
||||
u(ien) = wt - q(ipr) - u(idn) - 5.0d-01 * (vm * bb + vb * vb)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
|
||||
@ -75,6 +75,7 @@ module evolution
|
||||
real(kind=8) , save :: dth = 1.0d+00
|
||||
real(kind=8) , save :: dte = 0.0d+00
|
||||
real(kind=8) , save :: dtp = huge(dt) ! dt for precise snapshots
|
||||
real(kind=8) , save :: dtsafe = 1.0d-02
|
||||
|
||||
! the absolute and relative tolerances, limiting factors, the maximum error,
|
||||
! the maximum number of passes for the adaptive step,
|
||||
@ -195,6 +196,7 @@ module evolution
|
||||
call get_parameter("factor" , fac )
|
||||
call get_parameter("minimum_factor" , facmin)
|
||||
call get_parameter("maximum_factor" , facmax)
|
||||
call get_parameter("dt_safe_factor" , dtsafe)
|
||||
|
||||
! select the integration method, check the correctness of the integration
|
||||
! parameters and adjust the CFL coefficient if necessary
|
||||
@ -2947,7 +2949,7 @@ module evolution
|
||||
|
||||
else
|
||||
if (status == 0) then
|
||||
dt = dte
|
||||
dt = max(dte, dtsafe * dth)
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
else
|
||||
dt = 2.5d-01 * dt
|
||||
@ -3202,7 +3204,7 @@ module evolution
|
||||
|
||||
else
|
||||
if (status == 0) then
|
||||
dt = dte
|
||||
dt = max(dte, dtsafe * dth)
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
else
|
||||
dt = 2.5d-01 * dt
|
||||
@ -3525,7 +3527,7 @@ module evolution
|
||||
errs(2) = errs(1)
|
||||
else
|
||||
if (status == 0) then
|
||||
dt = dte
|
||||
dt = max(dte, dtsafe * dth)
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
else
|
||||
dt = 2.5d-01 * dt
|
||||
@ -3786,7 +3788,7 @@ module evolution
|
||||
if (fsal) fsal_update = .true.
|
||||
else
|
||||
if (status == 0) then
|
||||
dt = dte
|
||||
dt = max(dte, dtsafe * dth)
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
else
|
||||
dt = 2.5d-01 * dt
|
||||
@ -4103,6 +4105,15 @@ module evolution
|
||||
|
||||
n = get_dblocks()
|
||||
|
||||
!$omp parallel do default(shared) private(pdata,pmeta)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (pmeta%update) call update_shapes(pdata, tm, dtm)
|
||||
end do
|
||||
!$omp end parallel do
|
||||
|
||||
call boundary_variables(tm, dtm, status)
|
||||
|
||||
#ifdef MPI
|
||||
@ -4110,6 +4121,15 @@ module evolution
|
||||
#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%boundary) call update_shapes(pdata, tm, dtm)
|
||||
end do
|
||||
!$omp end parallel do
|
||||
|
||||
!$omp parallel do default(shared) private(pdata,pmeta,s)
|
||||
do l = 1, n
|
||||
pdata => data_blocks(l)%ptr
|
||||
@ -4136,6 +4156,8 @@ module evolution
|
||||
!$omp critical
|
||||
if (s /= 0) status = 1
|
||||
!$omp end critical
|
||||
|
||||
call update_shapes(pdata, tm, dtm)
|
||||
end if
|
||||
|
||||
end do
|
||||
@ -4152,15 +4174,6 @@ module evolution
|
||||
#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
|
||||
|
||||
@ -5709,8 +5709,8 @@ module interpolations
|
||||
|
||||
end if
|
||||
|
||||
qmx = max(min(u(i), u(ip1), qmd), min(u(i), qul, qlc))
|
||||
qmn = min(max(u(i), u(ip1), qmd), max(u(i), qul, qlc))
|
||||
qmn = max(min(u(i), u(ip1), qmd), min(u(i), qul, qlc))
|
||||
qmx = min(max(u(i), u(ip1), qmd), max(u(i), qul, qlc))
|
||||
|
||||
if (p .and. qmn <= 0.0d+00) then
|
||||
if (d(i) <= 0.0d+00 .and. d(ip1) >= 0.0d+00) then
|
||||
|
||||
@ -1752,14 +1752,6 @@ module mesh
|
||||
end if
|
||||
|
||||
end if ! pmeta belongs to the current process
|
||||
|
||||
! redistribute blocks equally among all processors
|
||||
!
|
||||
call redistribute_blocks(status)
|
||||
if (status /= 0) then
|
||||
call print_message(loc, "Cannot redistribute blocks!")
|
||||
go to 100
|
||||
end if
|
||||
#endif /* MPI */
|
||||
|
||||
end if ! leaf at current level selected for refinement
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user