5136 lines
137 KiB
Fortran
5136 lines
137 KiB
Fortran
!!******************************************************************************
|
||
!!
|
||
!! This file is part of the AMUN source code, a program to perform
|
||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||
!! adaptive mesh.
|
||
!!
|
||
!! Copyright (C) 2008-2023 Grzegorz Kowal <grzegorz@amuncode.org>
|
||
!!
|
||
!! This program is free software: you can redistribute it and/or modify
|
||
!! it under the terms of the GNU General Public License as published by
|
||
!! the Free Software Foundation, either version 3 of the License, or
|
||
!! (at your option) any later version.
|
||
!!
|
||
!! This program is distributed in the hope that it will be useful,
|
||
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
!! GNU General Public License for more details.
|
||
!!
|
||
!! You should have received a copy of the GNU General Public License
|
||
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
!!
|
||
!!******************************************************************************
|
||
!!
|
||
!! module: EQUATIONS
|
||
!!
|
||
!! This module provides interface for the systems of equations. Any set of
|
||
!! equations gives us some basic informations, such as the number of variables,
|
||
!! the primitive and conservative variable definitions, the conversion between
|
||
!! those variables, the flux and characteristic speeds defined in terms of
|
||
!! primitive variables. All this information is provided by this module.
|
||
!!
|
||
!! In order to implement a new set of equations, we need to:
|
||
!!
|
||
!! 1) define the number of independent variables (or equations) nv;
|
||
!! 2) define the variable indices and names (both primitive and conservative);
|
||
!! 3) provide subroutines for primitive-conservative variable conversion and
|
||
!! point them to corresponding pointers;
|
||
!! 4) provide a subroutine to calculate physical fluxes and characteristic
|
||
!! speeds;
|
||
!! 5) provide a subroutine to calculate the maximum speed;
|
||
!! 6) optionally, define and read all physical constants related to a given
|
||
!! system;
|
||
!!
|
||
!!******************************************************************************
|
||
!
|
||
module equations
|
||
|
||
implicit none
|
||
|
||
! the number of variables, fluxes and passive scalars
|
||
!
|
||
integer(kind=4), save :: nv = 0
|
||
integer(kind=4), save :: nf = 0
|
||
integer(kind=4), save :: ns = 0
|
||
|
||
! 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
|
||
end subroutine
|
||
subroutine cons2prim_iface(u, q, s)
|
||
real(kind=8), dimension(:,:), intent(in) :: u
|
||
real(kind=8), dimension(:,:), intent(out) :: q
|
||
integer , intent(out) :: s
|
||
end subroutine
|
||
subroutine fluxspeed_iface(q, u, f, c)
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
end subroutine
|
||
subroutine get_maximum_speeds_iface(qq, um, cm)
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out):: um, cm
|
||
end subroutine
|
||
subroutine nr_iterate_iface(mm, bb, mb, en, dn, w, vv, info)
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
end subroutine
|
||
end interface
|
||
|
||
! pointers to the conversion procedures
|
||
!
|
||
procedure(prim2cons_iface) , pointer, save :: prim2cons => null()
|
||
procedure(cons2prim_iface) , pointer, save :: cons2prim => null()
|
||
|
||
! pointer to the flux procedure
|
||
!
|
||
procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null()
|
||
|
||
procedure(get_maximum_speeds_iface), pointer, save :: &
|
||
get_maximum_speeds => null()
|
||
|
||
! pointer to the variable conversion method
|
||
!
|
||
procedure(nr_iterate_iface) , pointer, save :: nr_iterate => null()
|
||
|
||
! the system of equations and the equation of state
|
||
!
|
||
character(len=32), save :: eqsys = "hydrodynamic"
|
||
character(len=32), save :: eos = "adiabatic"
|
||
|
||
! the flag indicating if the set of equations is relativistic or magnetized
|
||
!
|
||
logical , save :: relativistic = .false.
|
||
logical , save :: magnetized = .false.
|
||
|
||
! the variable conversion method
|
||
!
|
||
character(len=32), save :: c2p = "1Dw"
|
||
|
||
! the names of equations and methods
|
||
!
|
||
character(len=80), save :: name_eqsys = ""
|
||
character(len=80), save :: name_eos = ""
|
||
character(len=80), save :: name_c2p = ""
|
||
|
||
! direction indices
|
||
!
|
||
integer(kind=4) , save :: inx = 1, iny = 2, inz = 3
|
||
|
||
! variable indices
|
||
!
|
||
integer(kind=4) , save :: idn = -1
|
||
integer(kind=4) , save :: ivx = -1, ivy = -1, ivz = -1
|
||
integer(kind=4) , save :: imx = -1, imy = -1, imz = -1
|
||
integer(kind=4) , save :: ibx = -1, iby = -1, ibz = -1
|
||
integer(kind=4) , save :: ibp = -1
|
||
integer(kind=4) , save :: ipr = -1, ien = -1
|
||
integer(kind=4) , save :: isl = -1, isu = -1
|
||
|
||
! variable names
|
||
!
|
||
character(len=4), dimension(:), allocatable, save :: pvars, cvars
|
||
|
||
! variable boundary values
|
||
!
|
||
real(kind=8), dimension(:,:,:), allocatable, save :: qpbnd
|
||
|
||
! adiabatic heat ratio
|
||
!
|
||
real(kind=8), save :: adiabatic_index = 5.0d+00 / 3.0d+00
|
||
|
||
! additional adiabatic parameters
|
||
!
|
||
real(kind=8), save :: gammam1 = 2.0d+00 / 3.0d+00, gammam1i = 1.5d+00
|
||
real(kind=8), save :: gammaxi = 2.0d+00 / 5.0d+00
|
||
|
||
! isothermal speed of sound and its second power
|
||
!
|
||
real(kind=8), save :: csnd = 1.0d+00, csnd2 = 1.0d+00
|
||
|
||
! maximum speeds in the system
|
||
!
|
||
real(kind=8), save :: cmax = 0.0d+00, cmax2 = 0.0d+00
|
||
real(kind=8), save :: umax = 0.0d+00, cglm = 0.0d+00
|
||
|
||
! the lower limits for density and pressure to be treated as physical
|
||
!
|
||
real(kind=8), save :: dmin = 1.0d-16, pmin = 1.0d-16
|
||
|
||
! the upper limits for the Lorentz factor and corresponding |v|²
|
||
!
|
||
real(kind=8), save :: lmax = 1.0d+06
|
||
real(kind=8), save :: vmax = 0.999999999999d+00
|
||
|
||
! the upper bound for the sonic Mach number
|
||
!
|
||
real(kind=8), save :: msmax = 1.0d+03
|
||
real(kind=8), save :: msfac = 3.0d-06 / 5.0d+00
|
||
|
||
! the tolerance for Newton-Raphson interative method, the maximum number of
|
||
! iterations and the number of extra iterations for polishing
|
||
!
|
||
real(kind=8), save :: tol = 1.0d-10
|
||
integer , save :: nrmax = 100
|
||
integer , save :: nrext = 2
|
||
|
||
! flag for unphysical cells correction, the maximum distance of neighbors for
|
||
! averaging region, and the minimum number of cells for averaging
|
||
!
|
||
logical , save :: fix_unphysical_cells = .false.
|
||
integer , save :: ngavg = 2
|
||
integer , save :: npavg = 4
|
||
|
||
! the variable indices for update_flux() from module SCHEMES
|
||
!
|
||
integer, dimension(:,:), allocatable :: ivars
|
||
|
||
! the variable positivity indicator
|
||
!
|
||
logical, dimension(:), allocatable :: positive
|
||
|
||
! the array of maximum errors for conservative variables
|
||
!
|
||
real(kind=8), dimension(:), allocatable :: errors
|
||
|
||
! by default everything is private
|
||
!
|
||
private
|
||
|
||
! 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 :: get_maximum_speeds
|
||
public :: update_primitive_variables
|
||
public :: fix_unphysical_cells, correct_unphysical_states
|
||
public :: adiabatic_index, relativistic, magnetized
|
||
public :: csnd, csnd2
|
||
public :: cmax, cmax2, umax, cglm
|
||
public :: nv, nf, ns
|
||
public :: inx, iny, inz
|
||
public :: idn, ivx, ivy, ivz, imx, imy, imz
|
||
public :: ibx, iby, ibz, ibp, ipr, ien
|
||
public :: isl, isu
|
||
public :: eqsys, eos
|
||
public :: pvars, cvars
|
||
public :: qpbnd
|
||
public :: ivars, positive, errors
|
||
|
||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||
!
|
||
contains
|
||
!
|
||
!===============================================================================
|
||
!!
|
||
!!*** PUBLIC SUBROUTINES *****************************************************
|
||
!!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INITIALIZE_EQUATIONS:
|
||
! -------------------------------
|
||
!
|
||
! Subroutine initiate the module by setting module parameters and subroutine
|
||
! pointers.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! system - the equation system
|
||
! state - the equation of state
|
||
! verbose - a logical flag turning the information printing;
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine initialize_equations(system, state, verbose, status)
|
||
|
||
! include external procedures and variables
|
||
!
|
||
use parameters, only : get_parameter
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
character(len=32), intent(in) :: system, state
|
||
logical , intent(in) :: verbose
|
||
integer , intent(out) :: status
|
||
|
||
! local variables
|
||
!
|
||
integer :: p
|
||
character(len=32) :: unphysical_fix = "off"
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
! set the system of equations
|
||
!
|
||
eqsys = system
|
||
|
||
! set the equation of state
|
||
!
|
||
eos = state
|
||
|
||
! get the number of passive scalars
|
||
!
|
||
call get_parameter("nscalars" , ns )
|
||
ns = min(ns, 100)
|
||
|
||
! get the primitive variable solver
|
||
!
|
||
call get_parameter("primitive_solver", c2p)
|
||
|
||
! depending on the system of equations initialize the module variables
|
||
!
|
||
select case(trim(eqsys))
|
||
|
||
!--- HYDRODYNAMICS ---
|
||
!
|
||
case("hd", "HD", "hydro", "HYDRO", "hydrodynamic", "HYDRODYNAMIC")
|
||
|
||
! the name of equation system
|
||
!
|
||
name_eqsys = "HD"
|
||
eqsys = "hd"
|
||
|
||
! initialize the indices of density, and velocity and momenta components
|
||
!
|
||
idn = 1
|
||
ivx = 2
|
||
ivy = 3
|
||
ivz = 4
|
||
imx = 2
|
||
imy = 3
|
||
imz = 4
|
||
|
||
! depending on the equation of state complete the initialization
|
||
!
|
||
select case(trim(eos))
|
||
|
||
case("iso", "ISO", "isothermal", "ISOTHERMAL")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "isothermal"
|
||
eos = "iso"
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
!
|
||
nv = 4
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to subroutines
|
||
!
|
||
prim2cons => prim2cons_hd_iso
|
||
cons2prim => cons2prim_hd_iso
|
||
fluxspeed => fluxspeed_hd_iso
|
||
get_maximum_speeds => get_maximum_speeds_hd_iso
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "adiabatic"
|
||
eos = "adi"
|
||
|
||
! initialize the indices of pressure and total energy
|
||
!
|
||
ipr = 5
|
||
ien = 5
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
! - 1 component of pressure
|
||
!
|
||
nv = 5
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to subroutines
|
||
!
|
||
prim2cons => prim2cons_hd_adi
|
||
cons2prim => cons2prim_hd_adi
|
||
fluxspeed => fluxspeed_hd_adi
|
||
get_maximum_speeds => get_maximum_speeds_hd_adi
|
||
|
||
! warn about the unimplemented equation of state
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected equation of state is not " // &
|
||
"implemented for HD: '" // trim(eos) // "'."
|
||
write(*,"(1x,a)") "Available equations of state: 'iso', 'adi'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! proceed if everything is fine
|
||
!
|
||
if (status == 0) then
|
||
|
||
! initialize the number of fluxes
|
||
!
|
||
nf = nv
|
||
|
||
! set passive scalars index limits and update the number of variables
|
||
!
|
||
if (ns > 0) then
|
||
isl = nv + 1
|
||
isu = nv + ns
|
||
nv = isu
|
||
end if
|
||
|
||
! allocate arrays for variable names
|
||
!
|
||
allocate(pvars(nv), cvars(nv), stat = status)
|
||
|
||
if (status == 0) then
|
||
|
||
! fill in the primitive variable names
|
||
!
|
||
pvars(idn) = 'dens'
|
||
pvars(ivx) = 'velx'
|
||
pvars(ivy) = 'vely'
|
||
pvars(ivz) = 'velz'
|
||
if (ipr > 0) pvars(ipr) = 'pres'
|
||
|
||
! fill in the conservative variable names
|
||
!
|
||
cvars(idn) = 'dens'
|
||
cvars(imx) = 'momx'
|
||
cvars(imy) = 'momy'
|
||
cvars(imz) = 'momz'
|
||
if (ien > 0) cvars(ien) = 'ener'
|
||
|
||
end if ! status
|
||
end if ! status
|
||
|
||
!--- MAGNETOHYDRODYNAMICS ---
|
||
!
|
||
case("mhd", "MHD", "magnetohydrodynamic", "MAGNETOHYDRODYNAMIC")
|
||
|
||
! the name of equation system
|
||
!
|
||
name_eqsys = "MHD"
|
||
eqsys = "mhd"
|
||
|
||
! set magnetized flag
|
||
!
|
||
magnetized = .true.
|
||
|
||
! initialize the indices of density, and velocity and momenta components
|
||
!
|
||
idn = 1
|
||
ivx = 2
|
||
ivy = 3
|
||
ivz = 4
|
||
imx = 2
|
||
imy = 3
|
||
imz = 4
|
||
|
||
! depending on the equation of state complete the initialization
|
||
!
|
||
select case(trim(eos))
|
||
|
||
case("iso", "ISO", "isothermal", "ISOTHERMAL")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "isothermal"
|
||
eos = "iso"
|
||
|
||
! initialize the indices of magnetic field components and divergence potential
|
||
!
|
||
ibx = 5
|
||
iby = 6
|
||
ibz = 7
|
||
ibp = 8
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
! - 3 components of magnetic field
|
||
! - 1 component of divergence potential
|
||
!
|
||
nv = 8
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz, ibx, iby, ibz, ibp /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx, iby, ibz, ibx, ibp /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy, ibz, ibx, iby, ibp /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to the subroutines
|
||
!
|
||
prim2cons => prim2cons_mhd_iso
|
||
cons2prim => cons2prim_mhd_iso
|
||
fluxspeed => fluxspeed_mhd_iso
|
||
get_maximum_speeds => get_maximum_speeds_mhd_iso
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "adiabatic"
|
||
eos = "adi"
|
||
|
||
! initialize the indices of pressure, total energy, magnetic field components
|
||
! and divergence potential
|
||
!
|
||
ipr = 5
|
||
ien = 5
|
||
ibx = 6
|
||
iby = 7
|
||
ibz = 8
|
||
ibp = 9
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
! - 1 component of pressure
|
||
! - 3 components of magnetic field
|
||
! - 1 component of magnetic divergence potential
|
||
!
|
||
nv = 9
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr, iby, ibz, ibx, ibp /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr, ibz, ibx, iby, ibp /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to subroutines
|
||
!
|
||
prim2cons => prim2cons_mhd_adi
|
||
cons2prim => cons2prim_mhd_adi
|
||
fluxspeed => fluxspeed_mhd_adi
|
||
get_maximum_speeds => get_maximum_speeds_mhd_adi
|
||
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected equation of state is not " // &
|
||
"implemented for MHD: '" // trim(eos) // "'."
|
||
write(*,"(1x,a)") "Available equations of state: 'iso', 'adi'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! proceed if everything is fine
|
||
!
|
||
if (status == 0) then
|
||
|
||
! initialize the number of fluxes
|
||
!
|
||
nf = nv
|
||
|
||
! set passive scalars index limits and update the number of variables
|
||
!
|
||
if (ns > 0) then
|
||
isl = nv + 1
|
||
isu = nv + ns
|
||
nv = isu
|
||
end if
|
||
|
||
! allocate arrays for variable names
|
||
!
|
||
allocate(pvars(nv), cvars(nv), stat = status)
|
||
|
||
if (status == 0) then
|
||
|
||
! fill in the primitive variable names
|
||
!
|
||
pvars(idn) = 'dens'
|
||
pvars(ivx) = 'velx'
|
||
pvars(ivy) = 'vely'
|
||
pvars(ivz) = 'velz'
|
||
if (ipr > 0) pvars(ipr) = 'pres'
|
||
pvars(ibx) = 'magx'
|
||
pvars(iby) = 'magy'
|
||
pvars(ibz) = 'magz'
|
||
pvars(ibp) = 'bpsi'
|
||
|
||
! fill in the conservative variable names
|
||
!
|
||
cvars(idn) = 'dens'
|
||
cvars(imx) = 'momx'
|
||
cvars(imy) = 'momy'
|
||
cvars(imz) = 'momz'
|
||
if (ien > 0) cvars(ien) = 'ener'
|
||
cvars(ibx) = 'magx'
|
||
cvars(iby) = 'magy'
|
||
cvars(ibz) = 'magz'
|
||
cvars(ibp) = 'bpsi'
|
||
|
||
end if ! status
|
||
end if ! status
|
||
|
||
!--- SPECIAL RELATIVITY HYDRODYNAMICS ---
|
||
!
|
||
case("srhd", "SRHD")
|
||
|
||
! the name of equation system
|
||
!
|
||
name_eqsys = "Special Relativity HD"
|
||
eqsys = "srhd"
|
||
|
||
! set relativistic flag
|
||
!
|
||
relativistic = .true.
|
||
|
||
! initialize the indices of density, and velocity and momenta components
|
||
!
|
||
idn = 1
|
||
ivx = 2
|
||
ivy = 3
|
||
ivz = 4
|
||
imx = 2
|
||
imy = 3
|
||
imz = 4
|
||
|
||
! depending on the equation of state complete the initialization
|
||
!
|
||
select case(trim(eos))
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "adiabatic"
|
||
eos = "adi"
|
||
|
||
! initialize the indices of pressure and total energy
|
||
!
|
||
ipr = 5
|
||
ien = 5
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
! - 1 component of pressure
|
||
!
|
||
nv = 5
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to subroutines
|
||
!
|
||
prim2cons => prim2cons_srhd_adi
|
||
cons2prim => cons2prim_srhd_adi
|
||
fluxspeed => fluxspeed_srhd_adi
|
||
get_maximum_speeds => get_maximum_speeds_srhd_adi
|
||
|
||
! warn about the unimplemented equation of state
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected equation of state is not " // &
|
||
"implemented for SR-HD: '" // trim(eos) // "'."
|
||
write(*,"(1x,a)") "Available equations of state: 'adi'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! choose the conserved to primitive variable conversion method
|
||
!
|
||
select case(trim(c2p))
|
||
|
||
case("1Dw", "1dw", "1DW", "1D(w)", "1D(W)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "1D(W)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srhd_adi_1dw
|
||
|
||
case("2dwv", "2Dwv", "2D(w,v)", "2D(W,v)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "2D(W,v²)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srhd_adi_2dwv
|
||
|
||
case("2dwu", "2Dwu", "2D(w,u)", "2D(W,u)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "2D(W,u²)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srhd_adi_2dwu
|
||
|
||
! warn about the unimplemented method
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected conversion method is not " // &
|
||
"implemented " // trim(c2p) // " for SR-HD."
|
||
write(*,"(1x,a)") "Available conversions: '1Dw', '2Dwv', '2Dwu'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! proceed if everything is fine
|
||
!
|
||
if (status == 0) then
|
||
|
||
! initialize the number of fluxes
|
||
!
|
||
nf = nv
|
||
|
||
! set passive scalars index limits and update the number of variables
|
||
!
|
||
if (ns > 0) then
|
||
isl = nv + 1
|
||
isu = nv + ns
|
||
nv = isu
|
||
end if
|
||
|
||
! allocate arrays for variable names
|
||
!
|
||
allocate(pvars(nv), cvars(nv), stat = status)
|
||
|
||
if (status == 0) then
|
||
! fill in the primitive variable names
|
||
!
|
||
pvars(idn) = 'dens'
|
||
pvars(ivx) = 'velx'
|
||
pvars(ivy) = 'vely'
|
||
pvars(ivz) = 'velz'
|
||
if (ipr > 0) pvars(ipr) = 'pres'
|
||
|
||
! fill in the conservative variable names
|
||
!
|
||
cvars(idn) = 'dens'
|
||
cvars(imx) = 'momx'
|
||
cvars(imy) = 'momy'
|
||
cvars(imz) = 'momz'
|
||
if (ien > 0) cvars(ien) = 'ener'
|
||
|
||
end if ! status
|
||
end if ! status
|
||
|
||
!--- SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS ---
|
||
!
|
||
case("srmhd", "SRMHD")
|
||
|
||
! the name of equation system
|
||
!
|
||
name_eqsys = "Special Relativity MHD"
|
||
eqsys = "srmhd"
|
||
|
||
! set relativistic flag
|
||
!
|
||
relativistic = .true.
|
||
|
||
! set magnetized flag
|
||
!
|
||
magnetized = .true.
|
||
|
||
! initialize the indices of density, and velocity and momenta components
|
||
!
|
||
idn = 1
|
||
ivx = 2
|
||
ivy = 3
|
||
ivz = 4
|
||
imx = 2
|
||
imy = 3
|
||
imz = 4
|
||
|
||
! depending on the equation of state complete the initialization
|
||
!
|
||
select case(trim(eos))
|
||
|
||
case("adi", "ADI", "adiabatic", "ADIABATIC")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_eos = "adiabatic"
|
||
eos = "adi"
|
||
|
||
! initialize the indices of pressure, total energy, magnetic field components
|
||
! and divergence potential
|
||
!
|
||
ipr = 5
|
||
ien = 5
|
||
ibx = 6
|
||
iby = 7
|
||
ibz = 8
|
||
ibp = 9
|
||
|
||
! initialize the number of variables:
|
||
!
|
||
! - 1 component of density
|
||
! - 3 components of velocity
|
||
! - 1 component of pressure
|
||
! - 3 components of magnetic field
|
||
! - 1 component of magnetic divergence potential
|
||
!
|
||
nv = 9
|
||
|
||
! allocate the variable indices
|
||
!
|
||
allocate(ivars(NDIMS,nv), stat = status)
|
||
|
||
! fill up the variable indices
|
||
!
|
||
ivars(1,:) = (/ idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp /)
|
||
ivars(2,:) = (/ idn, ivy, ivz, ivx, ipr, iby, ibz, ibx, ibp /)
|
||
#if NDIMS == 3
|
||
ivars(3,:) = (/ idn, ivz, ivx, ivy, ipr, ibz, ibx, iby, ibp /)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! set pointers to subroutines
|
||
!
|
||
prim2cons => prim2cons_srmhd_adi
|
||
cons2prim => cons2prim_srmhd_adi
|
||
fluxspeed => fluxspeed_srmhd_adi
|
||
get_maximum_speeds => get_maximum_speeds_srmhd_adi
|
||
|
||
! warn about the unimplemented equation of state
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected equation of state is not " // &
|
||
"implemented for SR-MHD: '" // trim(eos) // "'."
|
||
write(*,"(1x,a)") "Available equations of state: 'adi'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! choose the conserved to primitive variable conversion method
|
||
!
|
||
select case(trim(c2p))
|
||
|
||
case("1Dw", "1dw", "1DW", "1D(w)", "1D(W)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "1D(W)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srmhd_adi_1dw
|
||
|
||
case("2dwv", "2Dwv", "2D(w,v)", "2D(W,v)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "2D(W,v²)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srmhd_adi_2dwv
|
||
|
||
case("2dwu", "2Dwu", "2D(w,u)", "2D(W,u)")
|
||
|
||
! the type of equation of state
|
||
!
|
||
name_c2p = "2D(W,u²)"
|
||
|
||
! set pointer to the conversion method
|
||
!
|
||
nr_iterate => nr_iterate_srmhd_adi_2dwu
|
||
|
||
! warn about the unimplemented method
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected conversion method is not " // &
|
||
"implemented " // trim(c2p) // " for SR-MHD."
|
||
write(*,"(1x,a)") "Available conversions: '1Dw', '2Dwv', '2Dwu'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! proceed if everything is fine
|
||
!
|
||
if (status == 0) then
|
||
|
||
! initialize the number of fluxes
|
||
!
|
||
nf = nv
|
||
|
||
! set passive scalars index limits and update the number of variables
|
||
!
|
||
if (ns > 0) then
|
||
isl = nv + 1
|
||
isu = nv + ns
|
||
nv = isu
|
||
end if
|
||
|
||
! allocate arrays for variable names
|
||
!
|
||
allocate(pvars(nv), cvars(nv), stat = status)
|
||
|
||
if (status == 0) then
|
||
|
||
! fill in the primitive variable names
|
||
!
|
||
pvars(idn) = 'dens'
|
||
pvars(ivx) = 'velx'
|
||
pvars(ivy) = 'vely'
|
||
pvars(ivz) = 'velz'
|
||
if (ipr > 0) pvars(ipr) = 'pres'
|
||
pvars(ibx) = 'magx'
|
||
pvars(iby) = 'magy'
|
||
pvars(ibz) = 'magz'
|
||
pvars(ibp) = 'bpsi'
|
||
|
||
! fill in the conservative variable names
|
||
!
|
||
cvars(idn) = 'dens'
|
||
cvars(imx) = 'momx'
|
||
cvars(imy) = 'momy'
|
||
cvars(imz) = 'momz'
|
||
if (ien > 0) cvars(ien) = 'ener'
|
||
cvars(ibx) = 'magx'
|
||
cvars(iby) = 'magy'
|
||
cvars(ibz) = 'magz'
|
||
cvars(ibp) = 'bpsi'
|
||
|
||
end if ! status
|
||
end if ! status
|
||
|
||
!--- EQUATION SYSTEM NOT IMPLEMENTED ---
|
||
!
|
||
case default
|
||
|
||
if (verbose) then
|
||
write(*,*)
|
||
write(*,"(1x,a)") "ERROR!"
|
||
write(*,"(1x,a)") "The selected equation system is not " // &
|
||
"implemented: " // trim(eqsys) // "."
|
||
write(*,"(1x,a)") "Available equation systems: 'hd' 'mhd'" // &
|
||
", 'srhd', 'srmhd'."
|
||
end if
|
||
status = 1
|
||
|
||
end select
|
||
|
||
! allocate positivity indicators
|
||
!
|
||
allocate(positive(nv), stat = status)
|
||
|
||
if (status == 0) then
|
||
|
||
positive( : ) = .false.
|
||
positive(idn) = .true.
|
||
if (ipr > 0) positive(ipr) = .true.
|
||
|
||
end if
|
||
|
||
! allocate an array for errors
|
||
!
|
||
allocate(errors(nf), stat = status)
|
||
|
||
if (status == 0) errors(:) = 0.0d+00
|
||
|
||
! proceed if everything is fine
|
||
!
|
||
if (status == 0) then
|
||
|
||
! generate the passive scalar names
|
||
!
|
||
if (ns > 0) then
|
||
do p = isl, isu
|
||
write(pvars(p), '("ps",i2.2)') p - isl
|
||
write(cvars(p), '("ps",i2.2)') p - isl
|
||
end do
|
||
end if
|
||
|
||
! obtain the adiabatic specific heat ratio
|
||
!
|
||
call get_parameter("adiabatic_index", adiabatic_index)
|
||
|
||
! calculate additional parameters
|
||
!
|
||
gammam1 = adiabatic_index - 1.0d+00
|
||
gammam1i = 1.0d+00 / gammam1
|
||
gammaxi = gammam1 / adiabatic_index
|
||
|
||
! obtain the isothermal sound speed
|
||
!
|
||
call get_parameter("sound_speed", csnd )
|
||
|
||
! calculate additional parameters
|
||
!
|
||
csnd2 = csnd * csnd
|
||
|
||
! allocate array for the boundary values
|
||
!
|
||
allocate(qpbnd(nv,3,2), stat = status)
|
||
|
||
if (status == 0) then
|
||
|
||
! set the boundary values
|
||
!
|
||
do p = 1, nv
|
||
|
||
! set the initial boundary values (1.0 for density and pressure, 0.0 otherwise)
|
||
!
|
||
if (pvars(p) == "dens" .or. pvars(p) == "pres") then
|
||
qpbnd(p,:,:) = 1.0d+00
|
||
else
|
||
qpbnd(p,:,:) = 0.0d+00
|
||
end if
|
||
|
||
! read the boundary values from the parameter file
|
||
!
|
||
call get_parameter(pvars(p) // "_bnd_xl", qpbnd(p,1,1))
|
||
call get_parameter(pvars(p) // "_bnd_xr", qpbnd(p,1,2))
|
||
call get_parameter(pvars(p) // "_bnd_yl", qpbnd(p,2,1))
|
||
call get_parameter(pvars(p) // "_bnd_yr", qpbnd(p,2,2))
|
||
call get_parameter(pvars(p) // "_bnd_zl", qpbnd(p,3,1))
|
||
call get_parameter(pvars(p) // "_bnd_zr", qpbnd(p,3,2))
|
||
|
||
end do ! over all variables
|
||
|
||
! get the minimum allowed density and pressure in the system, and the maximum
|
||
! Lorentz factor for special relativity
|
||
!
|
||
call get_parameter("dmin", dmin)
|
||
call get_parameter("pmin", pmin)
|
||
call get_parameter("lmax", lmax)
|
||
|
||
! calculate the maximum speed corresponding to the maximum Lorentz factor
|
||
!
|
||
vmax = 1.0d+00 - 1.0d+00 / lmax**2
|
||
|
||
! get the upper bound for the sonic Mach number
|
||
!
|
||
call get_parameter("msmax", msmax)
|
||
|
||
! calculate the sonic Mach number factor
|
||
!
|
||
msfac = 1.0d+00 / (adiabatic_index * msmax**2)
|
||
|
||
! get the tolerance
|
||
!
|
||
call get_parameter("tolerance" , tol )
|
||
|
||
! get the maximum number of Newton-Raphson method iterations
|
||
!
|
||
call get_parameter("nr_maxit" , nrmax)
|
||
call get_parameter("nr_extra" , nrext)
|
||
|
||
! get the state correction flags
|
||
!
|
||
call get_parameter("fix_unphysical_cells", unphysical_fix)
|
||
|
||
! check if the correction of unphysical cells is on
|
||
!
|
||
select case(trim(unphysical_fix))
|
||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||
fix_unphysical_cells = .true.
|
||
case default
|
||
fix_unphysical_cells = .false.
|
||
end select
|
||
|
||
! get parameters for unphysical cells correction
|
||
!
|
||
call get_parameter("ngavg", ngavg)
|
||
call get_parameter("npavg", npavg)
|
||
|
||
! correct the above parameters to reasonable values
|
||
!
|
||
ngavg = max(1, ngavg)
|
||
npavg = max(2, npavg)
|
||
|
||
end if ! status
|
||
end if ! status
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine initialize_equations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FINALIZE_EQUATIONS:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine releases memory used by the module.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! status - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine finalize_equations(status)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(out) :: status
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
! release the procedure pointers
|
||
!
|
||
nullify(prim2cons)
|
||
nullify(cons2prim)
|
||
nullify(fluxspeed)
|
||
|
||
! deallocate variable indices
|
||
!
|
||
if (allocated(ivars)) deallocate(ivars, stat = status)
|
||
|
||
! deallocate variable name arrays
|
||
!
|
||
if (allocated(pvars)) deallocate(pvars, stat = status)
|
||
if (allocated(cvars)) deallocate(cvars, stat = status)
|
||
|
||
! deallocate boundary values array
|
||
!
|
||
if (allocated(qpbnd)) deallocate(qpbnd, stat = status)
|
||
|
||
! deallocate positivity indicator
|
||
!
|
||
if (allocated(positive)) deallocate(positive, stat = status)
|
||
|
||
! deallocate the array for errors
|
||
!
|
||
if (allocated(errors)) deallocate(errors, stat = status)
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine finalize_equations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRINT_EQUATIONS:
|
||
! --------------------------
|
||
!
|
||
! Subroutine prints module parameters.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine print_equations(verbose)
|
||
|
||
! include external procedures and variables
|
||
!
|
||
use helpers, only : print_section, print_parameter
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
logical, intent(in) :: verbose
|
||
|
||
! local variables
|
||
!
|
||
character(len= 80) :: sfmt
|
||
character(len=255) :: msg
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
if (verbose) then
|
||
|
||
call print_section(verbose, "Physics")
|
||
call print_parameter(verbose, "equation system" , name_eqsys)
|
||
call print_parameter(verbose, "equation of state" , name_eos )
|
||
if (ipr > 0) then
|
||
call print_parameter(verbose, "adiabatic index", adiabatic_index)
|
||
else
|
||
call print_parameter(verbose, "sound speed" , csnd)
|
||
end if
|
||
call print_parameter(verbose, "number of variables" , nv )
|
||
call print_parameter(verbose, "number of fluxes" , nf )
|
||
call print_parameter(verbose, "number of passive scalars", ns )
|
||
write(sfmt,"(a,i0,a)") "(", nv, "(1x,a))"
|
||
write(msg,sfmt) cvars
|
||
call print_parameter(verbose, "conservative variables", msg )
|
||
write(msg,sfmt) pvars
|
||
call print_parameter(verbose, "primitive variables" , msg )
|
||
if (relativistic) then
|
||
call print_parameter(verbose, "variable conversion" , name_c2p )
|
||
end if
|
||
if (fix_unphysical_cells) then
|
||
call print_parameter(verbose, "fix unphysical cells" , "on" )
|
||
call print_parameter(verbose, "ngavg" , ngavg )
|
||
call print_parameter(verbose, "npavg" , npavg )
|
||
else
|
||
call print_parameter(verbose, "fix unphysical cells" , "off" )
|
||
end if
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine print_equations
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine UPDATE_PRIMITIVE_VARIABLES:
|
||
! -------------------------------------
|
||
!
|
||
! Subroutine updates primitive to conservative variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! uu - the input array of conservative variables;
|
||
! qq - the output array of primitive variables;
|
||
! ghosts - the flag indicating the conversion has to be done for
|
||
! the ghost zones too;
|
||
! status - the call status;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine update_primitive_variables(uu, qq, ghosts, status)
|
||
|
||
use coordinates, only : nn => bcells
|
||
use coordinates, only : nb, ne, nbl, neu
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(inout) :: uu
|
||
real(kind=8), dimension(:,:,:,:), intent(inout) :: qq
|
||
logical , intent(in) :: ghosts
|
||
integer , intent(out) :: status
|
||
|
||
integer :: i, j, k, s
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
status = 0
|
||
|
||
if (ghosts) then
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
call cons2prim(uu(:,:,j,k), qq(:,:,j,k), s)
|
||
status = max(status, s)
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
else
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
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
|
||
#endif /* NDIMS == 3 */
|
||
|
||
#if NDIMS == 3
|
||
do i = nbl, 1, -1
|
||
qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,nb,nb:ne,nb:ne)
|
||
end do
|
||
do i = neu, nn
|
||
qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,ne,nb:ne,nb:ne)
|
||
end do
|
||
do j = nbl, 1, -1
|
||
qq(1:nv, : , j,nb:ne) = qq(1:nv, : ,nb,nb:ne)
|
||
end do
|
||
do j = neu, nn
|
||
qq(1:nv, : , j,nb:ne) = qq(1:nv, : ,ne,nb:ne)
|
||
end do
|
||
do k = nbl, 1, -1
|
||
qq(1:nv, : , : , k) = qq(1:nv, : , : ,nb)
|
||
end do
|
||
do k = neu, nn
|
||
qq(1:nv, : , : , k) = qq(1:nv, : , : ,ne)
|
||
end do
|
||
#else /* NDIMS == 3 */
|
||
do i = nbl, 1, -1
|
||
qq(1:nv, i,nb:ne, : ) = qq(1:nv,nb,nb:ne, : )
|
||
end do
|
||
do i = neu, nn
|
||
qq(1:nv, i,nb:ne, : ) = qq(1:nv,ne,nb:ne, : )
|
||
end do
|
||
do j = nbl, 1, -1
|
||
qq(1:nv, : , j, : ) = qq(1:nv, : ,nb, : )
|
||
end do
|
||
do j = neu, nn
|
||
qq(1:nv, : , j, : ) = qq(1:nv, : ,ne, : )
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine update_primitive_variables
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CORRECT_UNPHYSICAL_STATES:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine seeks for unphysical states (cells with negative density or
|
||
! pressure) and try to fix them by averaging their values from physical
|
||
! neighbours.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! step - the time step number;
|
||
! pdata - the data block;
|
||
! status - the call status;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine correct_unphysical_states(step, pdata, status)
|
||
|
||
use blocks , only : block_data
|
||
use coordinates, only : nn => bcells
|
||
use helpers , only : print_message
|
||
#if NDIMS == 3
|
||
use coordinates, only : ax, ay, az
|
||
#else /* NDIMS == 3 */
|
||
use coordinates, only : ax, ay
|
||
#endif /* NDIMS == 3 */
|
||
|
||
implicit none
|
||
|
||
integer(kind=4) , intent(in) :: step
|
||
type(block_data), pointer, intent(inout) :: pdata
|
||
integer , intent(out) :: status
|
||
|
||
character(len=255) :: msg, sfmt
|
||
integer :: n, p, nc, np
|
||
integer :: i, il, iu
|
||
integer :: j, jl, ju
|
||
integer :: k
|
||
#if NDIMS == 3
|
||
integer :: kl, ku
|
||
#endif /* NDIMS == 3 */
|
||
|
||
#if NDIMS == 3
|
||
logical, dimension(nn,nn,nn) :: physical
|
||
#else /* NDIMS == 3 */
|
||
logical, dimension(nn,nn, 1) :: physical
|
||
#endif /* NDIMS == 3 */
|
||
|
||
integer , dimension(:,:), allocatable :: idx
|
||
real(kind=8), dimension(:,:), allocatable :: q, u
|
||
real(kind=8), dimension(nn) :: x, y
|
||
#if NDIMS == 3
|
||
real(kind=8), dimension(nn) :: z
|
||
#endif /* NDIMS == 3 */
|
||
|
||
character(len=*), parameter :: loc = &
|
||
'EQUATIONS::correct_unphysical_states()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
np = 0
|
||
i = 1
|
||
il = 1
|
||
iu = 1
|
||
j = 1
|
||
jl = 1
|
||
ju = 1
|
||
k = 1
|
||
#if NDIMS == 3
|
||
kl = 1
|
||
ku = 1
|
||
#endif /* NDIMS == 3 */
|
||
|
||
if (ipr > 0) then
|
||
physical(:,:,:) = pdata%q(idn,:,:,:) > 0.0d+00 .and. &
|
||
pdata%q(ipr,:,:,:) > 0.0d+00
|
||
else
|
||
physical(:,:,:) = pdata%q(idn,:,:,:) > 0.0d+00
|
||
end if
|
||
|
||
nc = count(.not. physical)
|
||
|
||
if (nc > 0) then
|
||
|
||
sfmt = '("Unphysical cells in block ID: ",i0," (",i0," counted) ' // &
|
||
'at time step ",i0,".")'
|
||
write(msg,sfmt) pdata%meta%id, nc, step
|
||
call print_message(loc, msg)
|
||
|
||
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
|
||
|
||
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
||
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
||
#if NDIMS == 3
|
||
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
n = 0
|
||
#if NDIMS == 3
|
||
do k = 1, nn
|
||
#endif /* NDIMS == 3 */
|
||
do j = 1, nn
|
||
do i = 1, nn
|
||
|
||
if (.not. physical(i,j,k)) then
|
||
|
||
n = n + 1
|
||
|
||
#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)
|
||
#else /* NDIMS == 3 */
|
||
sfmt = '("n: ",i0,", [i,j]: ",2(1x,i0),", [x,y]: ",2es12.4)'
|
||
write(msg,sfmt) n, i, j, x(i), y(j)
|
||
#endif /* NDIMS == 3 */
|
||
call print_message(msg)
|
||
|
||
idx(:,n) = [ i, j, k ]
|
||
|
||
! increase the region until we find at least three physical cells, but no more
|
||
! than 4 cells away
|
||
!
|
||
np = 0
|
||
p = 1
|
||
do while (np <= npavg .and. p <= ngavg)
|
||
il = max( 1, i - p)
|
||
iu = min(nn, i + p)
|
||
jl = max( 1, j - p)
|
||
ju = min(nn, j + p)
|
||
#if NDIMS == 3
|
||
kl = max( 1, k - p)
|
||
ku = min(nn, k + p)
|
||
|
||
np = count(physical(il:iu,jl:ju,kl:ku))
|
||
#else /* NDIMS == 3 */
|
||
|
||
np = count(physical(il:iu,jl:ju, 1 ))
|
||
#endif /* NDIMS == 3 */
|
||
|
||
p = p + 1
|
||
end do
|
||
|
||
! average primitive variables
|
||
!
|
||
if (np >= npavg) then
|
||
|
||
do p = 1, nv
|
||
#if NDIMS == 3
|
||
q(p,n) = sum(pdata%q(p,il:iu,jl:ju,kl:ku), &
|
||
physical(il:iu,jl:ju,kl:ku)) / np
|
||
#else /* NDIMS == 3 */
|
||
q(p,n) = sum(pdata%q(p,il:iu,jl:ju, 1 ), &
|
||
physical(il:iu,jl:ju, 1 )) / np
|
||
#endif /* NDIMS == 3 */
|
||
end do
|
||
|
||
else
|
||
|
||
! limit density or pressure to minimum value, since the averaging over
|
||
! neighbours failed
|
||
!
|
||
msg = "Not sufficient number of physical neighbors!"
|
||
call print_message(msg)
|
||
sfmt = '("Block ID:",i0,", cell position = ( ",3(i4," ")," ).")'
|
||
write(msg,sfmt) pdata%meta%id, i, j, k
|
||
call print_message(msg)
|
||
write(msg,"('Q = ',10(1x,1es24.16e3))") pdata%q(:,i,j,k)
|
||
call print_message(msg)
|
||
msg = "Applying lower bounds for positive variables."
|
||
call print_message(msg)
|
||
|
||
q( : ,n) = pdata%q( : ,i,j,k)
|
||
q(idn,n) = max(pdata%q(idn,i,j,k), dmin)
|
||
if (ipr > 0) q(ipr,n) = max(pdata%q(ipr,i,j,k), pmin)
|
||
|
||
end if ! not sufficient number of physical cells for averaging
|
||
|
||
end if
|
||
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
call prim2cons(q(:,1:nc), u(:,1:nc), .true.)
|
||
|
||
do n = 1, nc
|
||
i = idx(1,n)
|
||
j = idx(2,n)
|
||
k = idx(3,n)
|
||
|
||
pdata%q(:,i,j,k) = q(:,n)
|
||
pdata%u(:,i,j,k) = u(:,n)
|
||
end do
|
||
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine correct_unphysical_states
|
||
!
|
||
!===============================================================================
|
||
!!
|
||
!!*** PRIVATE SUBROUTINES ****************************************************
|
||
!!
|
||
!===============================================================================
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ISOTHERMAL HYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_HD_ISO:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_hd_iso(q, u, s)
|
||
|
||
! 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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! 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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_hd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_HD_ISO:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
|
||
#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
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_hd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_HD_ISO:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_hd_iso(q, u, f, c)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
! local variables
|
||
!
|
||
integer :: i
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate the hydrodynamic fluxes
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
f(idn,i) = u(imx,i)
|
||
f(imx,i) = q(ivx,i) * u(imx,i)
|
||
f(imy,i) = q(ivx,i) * u(imy,i)
|
||
f(imz,i) = q(ivx,i) * u(imz,i)
|
||
f(imx,i) = f(imx,i) + csnd2 * q(idn,i)
|
||
|
||
end do
|
||
|
||
! calculate the characteristic speeds
|
||
!
|
||
if (present(c)) then
|
||
|
||
do i = 1, size(q,2)
|
||
|
||
c(1,i) = q(ivx,i) - csnd
|
||
c(2,i) = q(ivx,i) + csnd
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_hd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_HD_ISO:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_hd_iso(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 0.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
cm = max(cm, abs(vl - csnd), abs(vu + csnd))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_hd_iso
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC HYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_HD_ADI:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_hd_adi(q, u, s)
|
||
|
||
! 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) :: 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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_hd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_HD_ADI:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
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)
|
||
q(ivy,i) = u(imy,i) / u(idn,i)
|
||
q(ivz,i) = u(imz,i) / u(idn,i)
|
||
ek = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
|
||
ei = u(ien,i) - ek
|
||
if (ei > 0.0d+00) then
|
||
q(ipr,i) = gammam1 * ei
|
||
else
|
||
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
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_hd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_HD_ADI:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_hd_adi(q, u, f, c)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
! local variables
|
||
!
|
||
integer :: i
|
||
real(kind=8) :: cs
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate the hydrodynamic fluxes
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
f(idn,i) = u(imx,i)
|
||
f(imx,i) = q(ivx,i) * u(imx,i)
|
||
f(imy,i) = q(ivx,i) * u(imy,i)
|
||
f(imz,i) = q(ivx,i) * u(imz,i)
|
||
f(imx,i) = f(imx,i) + q(ipr,i)
|
||
f(ien,i) = q(ivx,i) * (u(ien,i) + q(ipr,i))
|
||
|
||
end do
|
||
|
||
! calculate the characteristic speeds
|
||
!
|
||
if (present(c)) then
|
||
|
||
do i = 1, size(q,2)
|
||
|
||
cs = sqrt(adiabatic_index * q(ipr,i) / q(idn,i))
|
||
|
||
c(1,i) = q(ivx,i) - cs
|
||
c(2,i) = q(ivx,i) + cs
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_hd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_HD_ADI:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_hd_adi(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu, cc
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 0.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
|
||
cc = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k))
|
||
cm = max(cm, abs(vl - cc), abs(vu + cc))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_hd_adi
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ISOTHERMAL MAGNETOHYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_MHD_ISO:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_mhd_iso(q, u, s)
|
||
|
||
! 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
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! 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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_mhd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_MHD_ISO:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
|
||
#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)
|
||
q(ivy,i) = u(imy,i) / u(idn,i)
|
||
q(ivz,i) = u(imz,i) / u(idn,i)
|
||
q(ibx,i) = u(ibx,i)
|
||
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
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_mhd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_MHD_ISO:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_mhd_iso(q, u, f, c)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
! local variables
|
||
!
|
||
integer :: i
|
||
real(kind=8) :: by2, bz2, pt
|
||
real(kind=8) :: fa, fb, fc, cf
|
||
|
||
! local arrays
|
||
!
|
||
real(kind=8), dimension(size(q,2)) :: bx2, bb
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate the magnetohydrodynamic fluxes
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
bx2(i) = q(ibx,i) * q(ibx,i)
|
||
by2 = q(iby,i) * q(iby,i)
|
||
bz2 = q(ibz,i) * q(ibz,i)
|
||
bb(i) = bx2(i) + by2 + bz2
|
||
pt = csnd2 * q(idn,i) + 0.5d+00 * bb(i)
|
||
|
||
f(idn,i) = u(imx,i)
|
||
f(imx,i) = q(ivx,i) * u(imx,i) - bx2(i)
|
||
f(imy,i) = q(ivx,i) * u(imy,i) - q(ibx,i) * q(iby,i)
|
||
f(imz,i) = q(ivx,i) * u(imz,i) - q(ibx,i) * q(ibz,i)
|
||
f(imx,i) = f(imx,i) + pt
|
||
f(ibx,i) = q(ibp,i)
|
||
f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
|
||
f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
|
||
f(ibp,i) = cmax2 * q(ibx,i)
|
||
|
||
end do
|
||
|
||
! calculate the characteristic speeds
|
||
!
|
||
if (present(c)) then
|
||
|
||
do i = 1, size(q,2)
|
||
|
||
fa = csnd2 * q(idn,i)
|
||
fb = fa + bb(i)
|
||
fc = max(0.0d+00, fb * fb - 4.0d+00 * fa * bx2(i))
|
||
cf = sqrt(max(0.5d+00 * (fb + sqrt(fc)), bb(i)) / q(idn,i))
|
||
|
||
c(1,i) = q(ivx,i) - cf
|
||
c(2,i) = q(ivx,i) + cf
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_mhd_iso
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_MHD_ISO:
|
||
! -------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_mhd_iso(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu, cc, xx
|
||
|
||
real(kind=8), dimension(3) :: bb
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 0.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
|
||
bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
|
||
xx = csnd2 + bb(1) + bb(2) + bb(3)
|
||
cc = max(0.0d+00, xx**2 - 4.0d+00 * csnd2 * bb(1))
|
||
cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
|
||
cm = max(cm, abs(vl - cc), abs(vu + cc))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_mhd_iso
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC MAGNETOHYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_MHD_ADI:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_mhd_adi(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
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_mhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_MHD_ADI:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
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
|
||
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)
|
||
q(ivy,i) = u(imy,i) / u(idn,i)
|
||
q(ivz,i) = u(imz,i) / u(idn,i)
|
||
q(ibx,i) = u(ibx,i)
|
||
q(iby,i) = u(iby,i)
|
||
q(ibz,i) = u(ibz,i)
|
||
q(ibp,i) = u(ibp,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)
|
||
ei = u(ien,i) - 5.0d-01 * (ek + em + ep)
|
||
if (ei > 0.0d+00) then
|
||
q(ipr,i) = gammam1 * ei
|
||
else
|
||
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
|
||
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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_mhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_MHD_ADI:
|
||
! ----------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_mhd_adi(q, u, f, c)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
! local variables
|
||
!
|
||
integer :: i
|
||
real(kind=8) :: by2, bz2, pt
|
||
real(kind=8) :: vb
|
||
real(kind=8) :: fa, fb, fc, cf
|
||
|
||
! local arrays
|
||
!
|
||
real(kind=8), dimension(size(q,2)) :: bx2, bb
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate the magnetohydrodynamic fluxes
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
bx2(i) = q(ibx,i) * q(ibx,i)
|
||
by2 = q(iby,i) * q(iby,i)
|
||
bz2 = q(ibz,i) * q(ibz,i)
|
||
bb(i) = bx2(i) + by2 + bz2
|
||
vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
|
||
pt = q(ipr,i) + 0.5d+00 * bb(i)
|
||
|
||
f(idn,i) = u(imx,i)
|
||
f(imx,i) = q(ivx,i) * u(imx,i) - bx2(i)
|
||
f(imy,i) = q(ivx,i) * u(imy,i) - q(ibx,i) * q(iby,i)
|
||
f(imz,i) = q(ivx,i) * u(imz,i) - q(ibx,i) * q(ibz,i)
|
||
f(imx,i) = f(imx,i) + pt
|
||
f(ibx,i) = q(ibp,i)
|
||
f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
|
||
f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
|
||
f(ibp,i) = cmax2 * q(ibx,i)
|
||
f(ien,i) = q(ivx,i) * (u(ien,i) + pt) - q(ibx,i) * vb
|
||
|
||
end do
|
||
|
||
! calculate the characteristic speeds
|
||
!
|
||
if (present(c)) then
|
||
|
||
do i = 1, size(q,2)
|
||
|
||
fa = adiabatic_index * q(ipr,i)
|
||
fb = fa + bb(i)
|
||
fc = max(0.0d+00, fb * fb - 4.0d+00 * fa * bx2(i))
|
||
cf = sqrt(max(0.5d+00 * (fb + sqrt(fc)), bb(i)) / q(idn,i))
|
||
|
||
c(1,i) = q(ivx,i) - cf
|
||
c(2,i) = q(ivx,i) + cf
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_mhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_MHD_ADI:
|
||
! -------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_mhd_adi(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu, aa, cc, xx
|
||
|
||
real(kind=8), dimension(3) :: bb
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 0.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
|
||
aa = adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)
|
||
bb = qq(ibx:ibz,i,j,k) * qq(ibx:ibz,i,j,k) / qq(idn,i,j,k)
|
||
xx = aa + bb(1) + bb(2) + bb(3)
|
||
cc = max(0.0d+00, xx**2 - 4.0d+00 * aa * bb(1))
|
||
cc = sqrt(5.0d-01 * (xx + sqrt(cc)))
|
||
cm = max(cm, abs(vl - cc), abs(vu + cc))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_mhd_adi
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC SPECIAL RELATIVITY HYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_SRHD_ADI:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_srhd_adi(q, u, s)
|
||
|
||
! 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) :: vv, vm, vs, ww
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! iterate over all positions
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
! 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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_srhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_SRHD_ADI:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation using an interative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine cons2prim_srhd_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
|
||
|
||
logical :: info
|
||
integer :: i
|
||
real(kind=8) :: mm, bb, mb, en, dn
|
||
real(kind=8) :: w , vv, vm, vs
|
||
|
||
#ifdef DEBUG
|
||
character(len=256) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()'
|
||
#endif /* DEBUG */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
s = 0
|
||
|
||
do i = 1, size(u,2)
|
||
|
||
mm = sum(u(imx:imz,i) * u(imx:imz,i))
|
||
en = u(ien,i) + u(idn,i)
|
||
dn = u(idn,i)
|
||
|
||
call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
if (info) then
|
||
|
||
vm = 1.0d+00 - vv
|
||
vs = sqrt(vm)
|
||
|
||
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
|
||
|
||
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
|
||
|
||
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)
|
||
#endif /* DEBUG */
|
||
end if
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_srhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_SRHD_ADI:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A., Bodo, G.,
|
||
! "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society, 2005, 364, 126-136
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_srhd_adi(q, u, f, c)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
! local variables
|
||
!
|
||
integer :: i
|
||
real(kind=8) :: vv, ww, c2, ss, cc, fc
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate the relativistic hydrodynamic fluxes (eq. 2 in [1])
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
f(idn,i) = u(idn,i) * q(ivx,i)
|
||
f(imx,i) = u(imx,i) * q(ivx,i) + q(ipr,i)
|
||
f(imy,i) = u(imy,i) * q(ivx,i)
|
||
f(imz,i) = u(imz,i) * q(ivx,i)
|
||
f(ien,i) = u(imx,i) - f(idn,i)
|
||
|
||
end do
|
||
|
||
! calculate characteristic speeds (eq. 23 in [1])
|
||
!
|
||
if (present(c)) then
|
||
|
||
do i = 1, size(q,2)
|
||
|
||
ww = q(idn,i) + q(ipr,i) / gammaxi
|
||
c2 = adiabatic_index * q(ipr,i) / ww
|
||
vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
|
||
ss = c2 * (1.0d+00 - vv) / (1.0d+00 - c2)
|
||
fc = 1.0d+00 + ss
|
||
cc = sqrt(ss * (fc - q(ivx,i)**2))
|
||
|
||
c(1,i) = (q(ivx,i) - cc) / fc
|
||
c(2,i) = (q(ivx,i) + cc) / fc
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_srhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_SRHD_ADI:
|
||
! --------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_srhd_adi(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu, vv, ww, aa, cc, ss, fc
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 0.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
vv = sum(qq(ivx:ivz,i,j,k) * qq(ivx:ivz,i,j,k))
|
||
|
||
ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi
|
||
aa = adiabatic_index * qq(ipr,i,j,k) / ww
|
||
ss = aa * (1.0d+00 - vv) / (1.0d+00 - aa)
|
||
fc = 1.0d+00 + ss
|
||
cc = sqrt(ss * (fc - vv))
|
||
cm = max(cm, abs(vl - cc), abs(vu + cc))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_srhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_FUNCTION_SRHD_ADI_1D:
|
||
! ----------------------------------
|
||
!
|
||
! Subroutine calculate the value of function
|
||
!
|
||
! F(W) = W - P(W) - E
|
||
!
|
||
! for a given enthalpy W. It is used to estimate the initial guess.
|
||
!
|
||
! The pressure is
|
||
!
|
||
! P(W) = (γ - 1)/γ (W - D / sqrt(1 - |v|²(W))) (1 - |v|²(W))
|
||
!
|
||
! and the squared velocity is
|
||
!
|
||
! |v|²(W) = |m|² / W²
|
||
!
|
||
! Optional derivative is returned
|
||
!
|
||
! dF(W)/dW = 1 - dP(W)/dW
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en, dn, w - input coefficients for |M|² E, D, and W, respectively;
|
||
! f - the value of function F(W);
|
||
! df - optional derivative F'(W);
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8) , intent(in) :: mm, en, dn, w
|
||
real(kind=8) , intent(out) :: f
|
||
real(kind=8), optional, intent(out) :: df
|
||
|
||
! local variables
|
||
!
|
||
real(kind=8) :: vv, vm, vs
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vv = mm / (w * w)
|
||
vm = 1.0d+00 - vv
|
||
vs = sqrt(vm)
|
||
f = (1.0d+00 - gammaxi * vm) * w + gammaxi * dn * vs - en
|
||
if (present(df)) then
|
||
df = 1.0d+00 - gammaxi * (1.0d+00 + (1.0d+00 - dn / (vs * w)) * vv)
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_function_srhd_adi_1d
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRHD_ADI_1DW:
|
||
! ----------------------------------
|
||
!
|
||
! Subroutine finds a root W of equation
|
||
!
|
||
! F(W) = W - P(W) - E = 0
|
||
!
|
||
! using the Newton-Raphson 1Dw iterative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w, vv - input/output coefficients W and |v|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
! References:
|
||
!
|
||
! Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
|
||
! "Primitive Variable Solvers for Conservative General Relativistic
|
||
! Magnetohydrodynamics",
|
||
! The Astrophysical Journal, 2006, vol. 641, pp. 626-637
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
use helpers, only : print_message
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: f, df, dw
|
||
real(kind=8) :: err
|
||
|
||
character(len=80) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srhd_adi_1dw()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! check if the state is physical; this can save some time if unphysical state
|
||
! is considered
|
||
!
|
||
wl = sqrt(mm + dn * dn)
|
||
if (en > wl) then
|
||
|
||
! prepare the initial brackets
|
||
!
|
||
wl = wl + gammaxi * pmin
|
||
wu = en + pmin
|
||
|
||
! calculate the value of function for the lower bracket
|
||
!
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||
|
||
! the lower bracket gives negative function, so there is chance it bounds
|
||
! the root
|
||
!
|
||
if (fl < 0.0d+00) then
|
||
|
||
! make sure that the upper bracket is larger than the lower one and
|
||
! the function has positive value
|
||
!
|
||
if (wu <= wl) wu = 2.0d+00 * wl
|
||
|
||
! check if the brackets bound the root region, if not proceed until
|
||
! opposite function signs are found for the brackets
|
||
!
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
it = nrmax
|
||
keep = fl * fu > 0.0d+00
|
||
do while (keep)
|
||
it = it - 1
|
||
wl = wu
|
||
fl = fu
|
||
wu = 2.0d+00 * wu
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||
end do
|
||
|
||
! the upper bracket was found, so proceed with determining the root
|
||
!
|
||
if (it > 0 .and. fu >= 0.0d+00) then
|
||
|
||
! estimate the value of enthalpy close to the root and corresponding v²
|
||
!
|
||
w = wl - fl * (wu - wl) / (fu - fl)
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! iterate using the Newton-Raphson method in order to find a root w of the
|
||
! function
|
||
!
|
||
do while(keep)
|
||
|
||
! calculate F(W) and dF(W)/dW
|
||
!
|
||
call nr_function_srhd_adi_1d(mm, en, dn, w, f, df)
|
||
|
||
! update brackets
|
||
!
|
||
if (f > fl .and. f < 0.0d+00) then
|
||
wl = w
|
||
fl = f
|
||
end if
|
||
if (f < fu .and. f > 0.0d+00) then
|
||
wu = w
|
||
fu = f
|
||
end if
|
||
|
||
! calculate the increment dW, update the solution, and estimate the error
|
||
!
|
||
dw = f / df
|
||
w = w - dw
|
||
err = abs(dw / w)
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! if new W leaves the brackets, use the bisection method to estimate
|
||
! the new guess
|
||
!
|
||
if (w < wl .or. w > wu) then
|
||
w = 0.5d+00 * (wl + wu)
|
||
end if
|
||
|
||
it = it - 1
|
||
end do ! NR iterations
|
||
|
||
! print information about failed convergence or unphysical variables
|
||
!
|
||
if (err >= tol) then
|
||
call print_message(loc, "Convergence not reached!")
|
||
write(msg,"(a,1x,1es24.16e3)") "Error: ", err
|
||
call print_message(loc, msg)
|
||
end if
|
||
|
||
! calculate |V|² from W
|
||
!
|
||
vv = mm / (w * w)
|
||
|
||
else ! the upper brack not found
|
||
call print_message(loc, "Could not find the upper bracket!")
|
||
info = .false.
|
||
|
||
end if
|
||
|
||
else if (fl > 0.0d+00) then ! the root cannot be found, since it is below the lower bracket
|
||
call print_message(loc, "Positive function for lower bracket!")
|
||
info = .false.
|
||
else ! the lower bracket is a root, so return it
|
||
w = wl
|
||
info = .true.
|
||
end if
|
||
|
||
else ! the state is unphysical
|
||
call print_message(loc, "The state is not physical!")
|
||
info = .false.
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srhd_adi_1dw
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRHD_ADI_2DWV:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine finds a root (W,v²) of 2D equations
|
||
!
|
||
! F(W,v²) = W - E - P(W,v²) = 0
|
||
! G(W,v²) = W² v² - m² = 0
|
||
!
|
||
! using the Newton-Raphson 2D iterative method.
|
||
!
|
||
! All evaluated equations incorporate already the pressure of the form
|
||
!
|
||
! P(W,|v|²) = (γ - 1)/γ (W - Γ D) (1 - |v|²)
|
||
!
|
||
! in order to optimize calculations.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w, vv - input/output coefficients W and |v|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
! References:
|
||
!
|
||
! Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
|
||
! "Primitive Variable Solvers for Conservative General Relativistic
|
||
! Magnetohydrodynamics",
|
||
! The Astrophysical Journal, 2006, vol. 641, pp. 626-637
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srhd_adi_2dwv(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
use helpers, only : print_message
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: ww, vm, vs, gd, gv
|
||
real(kind=8) :: f, dfw, dfv
|
||
real(kind=8) :: g, dgw, dgv
|
||
real(kind=8) :: det, jfw, jfv, jgw, jgv
|
||
real(kind=8) :: dw, dv
|
||
real(kind=8) :: err
|
||
|
||
character(len=80) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srhd_adi_2dwv()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! prepare the initial brackets
|
||
!
|
||
wl = sqrt(mm + dn * dn) + gammaxi * pmin
|
||
wu = en + pmin
|
||
|
||
! make sure that the upper bracket is larger than the lower one
|
||
!
|
||
keep = wl >= wu
|
||
it = nrmax
|
||
do while(keep)
|
||
wu = 2.0d+00 * wu
|
||
it = it - 1
|
||
keep = (wl >= wu) .and. it > 0
|
||
end do
|
||
if (it <= 0) then
|
||
call print_message(loc, "Could not find the upper limit for enthalpy!")
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! check if the brackets bound the root region, if not proceed until
|
||
! opposite function signs are found for the brackets
|
||
!
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
|
||
keep = (fl * fu > 0.0d+00)
|
||
it = nrmax
|
||
|
||
do while (keep)
|
||
|
||
wl = wu
|
||
fl = fu
|
||
wu = 2.0d+00 * wu
|
||
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
|
||
it = it - 1
|
||
|
||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||
|
||
end do
|
||
if (it <= 0) then
|
||
call print_message(loc, "No initial brackets found!")
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! estimate the value of enthalpy close to the root and corresponding v²
|
||
!
|
||
w = wl - fl * (wu - wl) / (fu - fl)
|
||
vv = mm / (w * w)
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! iterate using the Newton-Raphson method in order to find the roots W and |v|²
|
||
! of functions
|
||
!
|
||
do while(keep)
|
||
|
||
! calculate W², (1 - |v|²), and the Lorentz factor
|
||
!
|
||
ww = w * w
|
||
vm = 1.0d+00 - vv
|
||
vs = sqrt(vm)
|
||
gd = gammaxi * dn
|
||
gv = 1.0d+00 - gammaxi * vm
|
||
|
||
! calculate F(W,|v|²) and G(W,|v|²)
|
||
!
|
||
f = gv * w - en + gd * vs
|
||
g = vv * ww - mm
|
||
|
||
! calculate dF(W,|v|²)/dW and dF(W,|v|²)/d|v|²
|
||
!
|
||
dfw = gv
|
||
dfv = gammaxi * w - 0.5d+00 * gd / vs
|
||
|
||
! calculate dG(W,|v|²)/dW and dG(W,|v|²)/d|v|²
|
||
!
|
||
dgw = 2.0d+00 * vv * w
|
||
dgv = ww
|
||
|
||
! invert the Jacobian J = | dF/dW, dF/d|v|² |
|
||
! | dG/dW, dG/d|v|² |
|
||
!
|
||
det = dfw * dgv - dfv * dgw
|
||
|
||
jfw = dgv / det
|
||
jgw = - dfv / det
|
||
jfv = - dgw / det
|
||
jgv = dfw / det
|
||
|
||
! calculate increments dW and d|v|²
|
||
!
|
||
dw = f * jfw + g * jgw
|
||
dv = f * jfv + g * jgv
|
||
|
||
! correct W and |v|²
|
||
!
|
||
w = w - dw
|
||
vv = vv - dv
|
||
|
||
! check if the new enthalpy and velocity are physical
|
||
!
|
||
if (w < wl) then
|
||
write(msg,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", &
|
||
w, wl
|
||
call print_message(loc, msg)
|
||
info = .false.
|
||
return
|
||
end if
|
||
if (vv < 0.0d+00 .or. vv >= 1.0d+00) then
|
||
write(msg,"(a,1x,1es24.16e3)") "Unphysical speed |v|²: ", vv
|
||
call print_message(loc, msg)
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! calculate the error
|
||
!
|
||
err = max(abs(dw / w), abs(dv))
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! decrease the number of remaining iterations
|
||
!
|
||
it = it - 1
|
||
|
||
end do ! NR iterations
|
||
|
||
! print information about failed convergence or unphysical variables
|
||
!
|
||
if (err >= tol) then
|
||
write(msg,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
|
||
call print_message(loc, msg)
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srhd_adi_2dwv
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRHD_ADI_2DWU:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine finds a root (W,u²) of 2D equations
|
||
!
|
||
! F(W,u²) = (W - E - P(W,u²)) (u² + 1) = 0
|
||
! G(W,u²) = W² u² - (u² + 1) m² = 0
|
||
!
|
||
! using the Newton-Raphson 2D iterative method.
|
||
!
|
||
! All evaluated equations incorporate already the pressure of the form
|
||
!
|
||
! P(W,|u|²) = (γ - 1)/γ (W - Γ D) / (1 + |u|²)
|
||
!
|
||
! in order to optimize calculations.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w, vv - input/output coefficients W and |v|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srhd_adi_2dwu(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
! local variables
|
||
!
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: ww, uu, up, gm, gd
|
||
real(kind=8) :: f, dfw, dfu
|
||
real(kind=8) :: g, dgw, dgu
|
||
real(kind=8) :: det, jfw, jfu, jgw, jgu
|
||
real(kind=8) :: dw, du
|
||
real(kind=8) :: err
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! prepare the initial brackets
|
||
!
|
||
wl = sqrt(mm + dn * dn) + gammaxi * pmin
|
||
wu = en + pmin
|
||
|
||
! make sure that the upper bracket is larger than the lower one
|
||
!
|
||
keep = wl >= wu
|
||
it = nrmax
|
||
do while(keep)
|
||
wu = 2.0d+00 * wu
|
||
it = it - 1
|
||
keep = (wl >= wu) .and. it > 0
|
||
end do
|
||
if (it <= 0) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)") "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srhd_adi_1dw()"
|
||
write(*,"(a)" ) "Could not find the upper limit for enthalpy!"
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! check if the brackets bound the root region, if not proceed until
|
||
! opposite function signs are found for the brackets
|
||
!
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wl, fl)
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
|
||
keep = (fl * fu > 0.0d+00)
|
||
it = nrmax
|
||
|
||
do while (keep)
|
||
|
||
wl = wu
|
||
fl = fu
|
||
wu = 2.0d+00 * wu
|
||
|
||
call nr_function_srhd_adi_1d(mm, en, dn, wu, fu)
|
||
|
||
it = it - 1
|
||
|
||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||
|
||
end do
|
||
if (it <= 0) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)") "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
|
||
write(*,"(a)" ) "No initial brackets found!"
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! estimate the value of enthalpy close to the root and corresponding u²
|
||
!
|
||
w = wl - fl * (wu - wl) / (fu - fl)
|
||
uu = mm / (w * w - mm)
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! iterate using the Newton-Raphson method in order to find the roots W and |u|²
|
||
! of functions
|
||
!
|
||
do while(keep)
|
||
|
||
! calculate W², (1 + |u|²), and the Lorentz factor, and some repeated
|
||
! expressions
|
||
!
|
||
ww = w * w
|
||
up = 1.0d+00 + uu
|
||
gm = sqrt(up)
|
||
gd = gammaxi * dn
|
||
|
||
! calculate F(W,|u|²) and G(W,|u|²)
|
||
!
|
||
f = (up - gammaxi) * w - up * en + gm * gd
|
||
g = uu * ww - up * mm
|
||
|
||
! calculate dF(W,|u|²)/dW and dF(W,|u|²)/d|u|²
|
||
!
|
||
dfw = up - gammaxi
|
||
dfu = w - en + 0.5d+00 * gd / gm
|
||
|
||
! calculate dG(W,|u|²)/dW and dG(W,|u|²)/d|u|²
|
||
!
|
||
dgw = 2.0d+00 * uu * w
|
||
dgu = ww - mm
|
||
|
||
! invert the Jacobian J = | dF/dW, dF/d|u|² |
|
||
! | dG/dW, dG/d|u|² |
|
||
!
|
||
det = dfw * dgu - dfu * dgw
|
||
|
||
jfw = dgu / det
|
||
jgw = - dfu / det
|
||
jfu = - dgw / det
|
||
jgu = dfw / det
|
||
|
||
! calculate increments dW and d|u|²
|
||
!
|
||
dw = f * jfw + g * jgw
|
||
du = f * jfu + g * jgu
|
||
|
||
! correct W and |u|²
|
||
!
|
||
w = w - dw
|
||
uu = uu - du
|
||
|
||
! check if the new enthalpy gives physical pressure and velocity
|
||
!
|
||
if (w < wl) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
|
||
write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
|
||
info = .false.
|
||
return
|
||
end if
|
||
if (uu < 0.0d+00) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
|
||
write(*,"(a,1x,1es24.16e3)") "Unphysical speed |u|²: ", uu
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! calculate the error
|
||
!
|
||
err = max(abs(dw / w), abs(du))
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! decrease the number of remaining iterations
|
||
!
|
||
it = it - 1
|
||
|
||
end do ! NR iterations
|
||
|
||
! calculate |v|² from |u|²
|
||
!
|
||
vv = uu / (1.0d+00 + uu)
|
||
|
||
! print information about failed convergence or unphysical variables
|
||
!
|
||
if (err >= tol) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||
, "EQUATIONS::nr_iterate_srhd_adi_2dwu()"
|
||
write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srhd_adi_2dwu
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
! ADIABATIC SPECIAL RELATIVITY MAGNETOHYDRODYNAMIC EQUATIONS
|
||
!
|
||
!*******************************************************************************
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine PRIM2CONS_SRMHD_ADI:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine converts primitive variables to their corresponding
|
||
! conservative representation.
|
||
!
|
||
! 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;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine prim2cons_srmhd_adi(q, u, s)
|
||
|
||
! 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) :: vv, bb, vb, vm, vs, ww, wt
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! iterate over all positions
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
! 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
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine prim2cons_srmhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine CONS2PRIM_SRMHD_ADI:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine converts conservative variables to their corresponding
|
||
! primitive representation using an interative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! u - the input array of conservative variables;
|
||
! q - the output array of primitive variables;
|
||
! s - the status flag;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine cons2prim_srmhd_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
|
||
|
||
logical :: info
|
||
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
|
||
|
||
do i = 1, size(u,2)
|
||
|
||
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)
|
||
|
||
call nr_iterate(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
if (info) then
|
||
|
||
vm = 1.0d+00 - vv
|
||
vs = sqrt(vm)
|
||
wt = w + bb
|
||
fc = mb / w
|
||
|
||
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
|
||
q(ivz,i) = (u(imz,i) + fc * u(ibz,i)) / wt
|
||
q(ibx,i) = u(ibx,i)
|
||
q(iby,i) = u(iby,i)
|
||
q(ibz,i) = u(ibz,i)
|
||
q(ibp,i) = u(ibp,i)
|
||
q(ipr,i) = w - en + 0.5d+00 * (bb + (bb * mm - mb * mb) / wt**2)
|
||
|
||
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(1:nf,i)
|
||
call print_message(loc, msg)
|
||
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
|
||
|
||
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 = ", &
|
||
dn, mm, mb, bb, en
|
||
call print_message(loc, msg)
|
||
#endif /* DEBUG */
|
||
end if
|
||
end do
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine cons2prim_srmhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FLUXSPEED_SRMHD_ADI:
|
||
! ------------------------------
|
||
!
|
||
! Subroutine calculates physical fluxes and characteristic speeds from a
|
||
! given equation system.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! q - the input array of primitive variables;
|
||
! u - the input array of conservative variables;
|
||
! f - the output vector of fluxes;
|
||
! c - the output vector of left- and right-going characteristic speeds;
|
||
!
|
||
! References:
|
||
!
|
||
! [1] Mignone, A., Bodo, G.,
|
||
! "An HLLC Riemann solver for relativistic flows -
|
||
! II. Magnetohydrodynamics",
|
||
! Monthly Notices of the Royal Astronomical Society,
|
||
! 2006, Volume 368, Pages 1040-1054
|
||
! [2] van der Holst, B., Keppens, R., Meliani, Z.
|
||
! "A multidimentional grid-adaptive relativistic magnetofluid code",
|
||
! Computer Physics Communications, 2008, Volume 179, Pages 617-627
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine fluxspeed_srmhd_adi(q, u, f, c)
|
||
|
||
use algebra, only : quadratic, quartic
|
||
use helpers, only : print_message
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:) , intent(in) :: q, u
|
||
real(kind=8), dimension(:,:) , intent(out) :: f
|
||
real(kind=8), dimension(:,:), optional, intent(out) :: c
|
||
|
||
integer :: i, nr
|
||
real(kind=8) :: bb, vs
|
||
real(kind=8) :: bx, by, bz, pm, pt
|
||
real(kind=8) :: rh, v1, v2
|
||
real(kind=8) :: ca, cc, c2, gn, rt, zm, zp
|
||
real(kind=8) :: fa, fb, fc, fd, fe, ff, fg
|
||
|
||
real(kind=8), dimension(size(q,2)) :: vv, vm, vb, b2
|
||
|
||
real(kind=8), dimension(5) :: a
|
||
real(kind=8), dimension(4) :: x
|
||
|
||
character(len=80) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::fluxspeed_srmhd_adi()'
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! iterate over all positions
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
! calculate the square of velocity, magnetic field and their scalar product
|
||
!
|
||
vv(i) = sum(q(ivx:ivz,i) * q(ivx:ivz,i))
|
||
bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i))
|
||
vb(i) = sum(q(ivx:ivz,i) * q(ibx:ibz,i))
|
||
|
||
! calculate (1 - |V|²)
|
||
!
|
||
vm(i) = 1.0d+00 - vv(i)
|
||
|
||
! calculate magnetic field components of the magnetic four-vector divided by
|
||
! the Lorentz factor (eq. 3 in [1])
|
||
!
|
||
bx = q(ibx,i) * vm(i) + vb(i) * q(ivx,i)
|
||
by = q(iby,i) * vm(i) + vb(i) * q(ivy,i)
|
||
bz = q(ibz,i) * vm(i) + vb(i) * q(ivz,i)
|
||
|
||
! calculate magnetic and total pressures (eq. 6 in [1])
|
||
!
|
||
b2(i) = bb * vm(i) + vb(i) * vb(i)
|
||
pm = 0.5d+00 * b2(i)
|
||
pt = q(ipr,i) + pm
|
||
|
||
! calculate the relativistic hydrodynamic fluxes (eq. 13 in [1])
|
||
!
|
||
f(idn,i) = u(idn,i) * q(ivx,i)
|
||
f(imx,i) = u(imx,i) * q(ivx,i) - q(ibx,i) * bx + pt
|
||
f(imy,i) = u(imy,i) * q(ivx,i) - q(ibx,i) * by
|
||
f(imz,i) = u(imz,i) * q(ivx,i) - q(ibx,i) * bz
|
||
f(ibx,i) = q(ibp,i)
|
||
f(ibp,i) = cmax2 * q(ibx,i)
|
||
f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i)
|
||
f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i)
|
||
f(ien,i) = u(imx,i) - f(idn,i)
|
||
|
||
end do
|
||
|
||
if (present(c)) then
|
||
|
||
! calculate the characteristic speeds
|
||
!
|
||
do i = 1, size(q,2)
|
||
|
||
! check if the total velocity |V|² is larger than zero
|
||
!
|
||
if (vv(i) > 0.0d+00) then
|
||
|
||
! calculate additional coefficients
|
||
!
|
||
rh = q(idn,i) + q(ipr,i) / gammaxi
|
||
vs = sqrt(vm(i))
|
||
|
||
! check if the normal component of magnetic field Bₓ is larger than zero
|
||
!
|
||
if (abs(q(ibx,i)) > 0.0d+00) then ! Bₓ ≠ 0
|
||
|
||
! prepare parameters for this case
|
||
!
|
||
c2 = adiabatic_index * q(ipr,i) / rh
|
||
v1 = abs(q(ivx,i))
|
||
v2 = v1 * v1
|
||
|
||
fa = rh * (1.0d+00 - c2)
|
||
fb = c2 * (rh - vb(i) * vb(i)) + b2(i)
|
||
fc = sign(1.0d+00, q(ivx,i)) * q(ibx,i) * vs
|
||
fd = c2 * fc * fc
|
||
fe = 1.0d+00 - v2
|
||
ff = c2 * vb(i) * fc
|
||
fg = v1 * vs
|
||
|
||
! prepare polynomial coefficients
|
||
!
|
||
a(5) = fa + fb * vm(i)
|
||
a(4) = 2.0d+00 * (ff * vm(i) + fb * fg)
|
||
a(3) = - fd * vm(i) + 4.0d+00 * ff * fg - fb * fe
|
||
a(2) = - 2.0d+00 * (fd * fg + fe * ff)
|
||
a(1) = fd * fe
|
||
|
||
! call the quartic solver
|
||
!
|
||
nr = quartic(a(1:5), x(1:4))
|
||
|
||
! convert eigenvalues to charasteristic speeds
|
||
!
|
||
x(1:nr) = sign(1.0d+00, q(ivx,i)) * (abs(v1) + x(1:nr) * vs)
|
||
|
||
else ! Bₓ ≠ 0
|
||
|
||
! special case when Bₓ = 0, then the quartic equation reduces to quadratic one
|
||
!
|
||
! prepare parameters for this case
|
||
!
|
||
c2 = adiabatic_index * q(ipr,i) / rh
|
||
cc = (1.0d+00 - c2) / vm(i)
|
||
gn = b2(i) - c2 * vb(i) * vb(i)
|
||
|
||
! prepare polynomial coefficients
|
||
!
|
||
a(3) = rh * (c2 + cc) + gn
|
||
a(2) = - 2.0d+00 * rh * cc * q(ivx,i)
|
||
a(1) = rh * (cc * q(ivx,i)**2 - c2) - gn
|
||
|
||
! solve the quadratic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
end if ! Bx ≠ 0
|
||
|
||
else ! |V|² > 0
|
||
|
||
! special case when |V|² = 0 (Γ = 1), then the quartic equation reduces to
|
||
! bi-quartic one
|
||
!
|
||
! prepare parameters for this case
|
||
!
|
||
rh = q(idn,i) + q(ipr,i) / gammaxi
|
||
vs = sqrt(vm(i))
|
||
rt = rh + b2(i)
|
||
c2 = adiabatic_index * q(ipr,i) / rh
|
||
ca = (q(ibx,i) * vs + vb(i) * q(ivx,i) / vs)**2
|
||
|
||
! prepare polynomial coefficients
|
||
!
|
||
a(3) = 1.0d+00
|
||
a(2) = - ((rh + ca) * c2 + b2(i)) / rt
|
||
a(1) = c2 * ca / rt
|
||
|
||
! solve the bi-quartic equation
|
||
!
|
||
nr = quadratic(a(1:3), x(1:2))
|
||
|
||
! compute the roots
|
||
!
|
||
if (nr > 0) then
|
||
|
||
zm = min(x(1), x(2))
|
||
zp = max(x(1), x(2))
|
||
|
||
if (zm >= 0.0d+00) then
|
||
|
||
zm = sqrt(zm)
|
||
zp = sqrt(zp)
|
||
|
||
x(1) = zp
|
||
x(2) = zm
|
||
x(3) = - zm
|
||
x(4) = - zp
|
||
|
||
nr = 4
|
||
|
||
else
|
||
|
||
if (zp >= 0.0d+00) then
|
||
|
||
zp = sqrt(zp)
|
||
|
||
x(1) = zp
|
||
x(2) = - zp
|
||
|
||
nr = 2
|
||
|
||
else
|
||
|
||
x(:) = 0.0d+00
|
||
|
||
nr = 0
|
||
|
||
end if
|
||
end if
|
||
end if
|
||
|
||
end if ! |V|² > 0
|
||
|
||
! find the minimum and maximum characteristic speeds
|
||
!
|
||
if (nr > 1) then
|
||
|
||
c(1,i) = minval(x(1:nr))
|
||
c(2,i) = maxval(x(1:nr))
|
||
|
||
if (max(abs(c(1,i)), abs(c(2,i))) >= 1.0d+00) then
|
||
call print_message(loc, "Estimation returned unphysical speeds!")
|
||
write(msg,"(a,5(1es24.16))") "A = ", a(1:5)
|
||
call print_message(loc, msg)
|
||
write(msg,"(a,1i2)") "N = ", nr
|
||
call print_message(loc, msg)
|
||
write(msg,"(a,4(1es24.16))") "X = ", x(1:4)
|
||
call print_message(loc, msg)
|
||
end if
|
||
|
||
else
|
||
|
||
! speed estimation failed, so we substitute the minimum and maximum physical
|
||
! speeds equal the speed of light
|
||
!
|
||
c(1,i) = - 1.0d+00
|
||
c(2,i) = 1.0d+00
|
||
|
||
end if
|
||
|
||
end do
|
||
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine fluxspeed_srmhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine GET_MAXIMUM_SPEEDS_SRMHD_ADI:
|
||
! ---------------------------------------
|
||
!
|
||
! Subroutine determines the maximum characteristic speed and eigenvalue
|
||
! in the input array of primitive variables.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! qq - the input array of primitive variables;
|
||
! vm - the maximum physical speed;
|
||
! cm - the maximum eigenvalue;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine get_maximum_speeds_srmhd_adi(qq, vm, cm)
|
||
|
||
use coordinates, only : nb, ne
|
||
|
||
implicit none
|
||
|
||
real(kind=8), dimension(:,:,:,:), intent(in) :: qq
|
||
real(kind=8) , intent(out) :: vm, cm
|
||
|
||
integer :: i, j, k
|
||
real(kind=8) :: vl, vu
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
vm = 0.0d+00
|
||
cm = 1.0d+00
|
||
|
||
#if NDIMS == 2
|
||
k = 1
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
do k = nb, ne
|
||
#endif /* NDIMS == 3 */
|
||
do j = nb, ne
|
||
do i = nb, ne
|
||
vl = minval(qq(ivx:ivz,i,j,k))
|
||
vu = maxval(qq(ivx:ivz,i,j,k))
|
||
vm = max(vm, abs(vl), abs(vu))
|
||
end do
|
||
end do
|
||
#if NDIMS == 3
|
||
end do
|
||
#endif /* NDIMS == 3 */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine get_maximum_speeds_srmhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_VELOCITY_SRMHD_ADI_1D:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine calculates the squared velocity
|
||
!
|
||
! |V|²(W) = (|m|² W² + S² (2 W + |B|²)) / (W² (W² + |B|²)²)
|
||
!
|
||
! and its derivative
|
||
!
|
||
! |V|²'(W) = - 2 (|m|² W³ + 3 S² W² + 3 S² |B|² W + S² |B|⁴)
|
||
! / (W³ (W² + |B|²)³)
|
||
!
|
||
! for a given enthalpy W.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, bb, mb, w - input coefficients for |M|², |B|², M.B, and W,
|
||
! respectively;
|
||
! vv, dv - the values of squared velocity |V|² and its derivative;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv, dv)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), intent(in) :: mm, bb, mb, w
|
||
real(kind=8), intent(out) :: vv
|
||
real(kind=8), optional, intent(out) :: dv
|
||
|
||
! local variables
|
||
!
|
||
real(kind=8) :: ss, ww, www, wt, wt2
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! temporary variables
|
||
!
|
||
ss = mb * mb
|
||
ww = w * w
|
||
www = ww * w
|
||
wt = w + bb
|
||
wt2 = wt * wt
|
||
|
||
! the function and its derivative
|
||
!
|
||
vv = (mm * ww + (w + wt) * ss) / (ww * wt2)
|
||
if (present(dv)) then
|
||
dv = - 2.0d+00 * (mm * www + (3.0d+00 * ww &
|
||
+ (2.0d+00 * w + wt) * bb) * ss) / (www * wt2 * wt)
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_velocity_srmhd_adi_1d
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_PRESSURE_SRMHD_ADI_1D:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine calculates the pressure function
|
||
!
|
||
! P(W) = W - E + ½ |B|² + ½ (|m|² |B|² - S²) / (W + |B|²)²
|
||
!
|
||
! and its derivative
|
||
!
|
||
! P'(W) = 1 - (|m|² |B|² - S²) / (W + |B|²)³
|
||
!
|
||
! for a given enthalpy W.
|
||
!
|
||
! This subroutine is used to find the minimum enthalpy for which the velocity
|
||
! is physical and the pressure is positive.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, bb, mb, en, w - input coefficients for |M|², |B|², M.B, E,
|
||
! and W, respectively;
|
||
! p, dp - the values for the function P(W) and its
|
||
! derivative P'(W);
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_pressure_srmhd_adi_1d(mm, bb, mb, en, w, p, dp)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8) , intent(in) :: mm, bb, mb, en, w
|
||
real(kind=8) , intent(out) :: p
|
||
real(kind=8), optional, intent(out) :: dp
|
||
|
||
! local variables
|
||
!
|
||
real(kind=8) :: wt, wd, ss, fn
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! temporary variables
|
||
!
|
||
wt = w + bb
|
||
wd = wt * wt
|
||
ss = mb * mb
|
||
fn = (mm * bb - ss) / wd
|
||
|
||
! the pressure function and its derivative
|
||
!
|
||
p = w - en + 0.5d+00 * (bb + fn)
|
||
if (present(dp)) dp = 1.0d+00 - fn / wt
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_pressure_srmhd_adi_1d
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_FUNCTION_SRMHD_ADI_1D:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine calculates the energy function
|
||
!
|
||
! F(W) = W - P(W) + ½ |B|² + ½ (|m|² |B|² - S²) / (W + |B|²)² - E
|
||
!
|
||
! and its derivative
|
||
!
|
||
! F'(W) = 1 - dP(W)/dW - (|m|² |B|² - S²) / (W + |B|²)³
|
||
!
|
||
! for a given enthalpy W. It is used to estimate the initial guess.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, bb, mb, en, dn, w - input coefficients for |M|², |B|², M.B, E, D,
|
||
! and W, respectively;
|
||
! f, df - the values of F(W) and its derivative;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn, w
|
||
real(kind=8), intent(out) :: f
|
||
real(kind=8), optional, intent(out) :: df
|
||
|
||
! local variables
|
||
!
|
||
real(kind=8) :: pr, dp, gm2, gm, dg
|
||
real(kind=8) :: ww, wt, wd, ss, sw, ws, fn, dv, ds
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! temporary variables
|
||
!
|
||
ww = w * w
|
||
wt = w + bb
|
||
wd = wt * wt
|
||
ss = mb * mb
|
||
sw = ss / ww
|
||
ws = (w + wt) * sw
|
||
fn = (mm * bb - ss) / wd
|
||
dv = wd - mm - ws
|
||
ds = sqrt(dv)
|
||
|
||
! calculate the Lorentz factor
|
||
!
|
||
gm2 = wd / dv
|
||
gm = wt / ds
|
||
|
||
! calculate the pressure P(W) and energy function F(W)
|
||
!
|
||
pr = gammaxi * (w - gm * dn) / gm2
|
||
f = w - pr - en + 0.5d+00 * (bb + fn)
|
||
|
||
! if desired, calculate the derivatives dP(W)/dW and dF(W)/dW
|
||
!
|
||
if (present(df)) then
|
||
|
||
dg = (1.0d+00 - wt * (wt - sw + ws / w) / dv) / ds
|
||
dp = gammaxi * (1.0d+00 - (2.0d+00 * w / gm - dn) * dg) / gm2
|
||
df = 1.0d+00 - dp - fn / wt
|
||
|
||
end if ! df present
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_function_srmhd_adi_1d
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_INITIAL_BRACKETS_SRMHD_ADI:
|
||
! ----------------------------------------
|
||
!
|
||
! Subroutine finds the initial brackets and initial guess from
|
||
! the positivity condition
|
||
!
|
||
! W³ + (5/2 |B|² - E) W² + 2 (|B|² - E) |B|² W
|
||
! + 1/2 [(|B|⁴ - 2 |B|² E + |m|²) |B|² - S²] > 0
|
||
!
|
||
! coming from the energy equation and
|
||
!
|
||
! W⁴ + 2 |B|² W³ - (|m|² + D² - |B|⁴) W²
|
||
! - (2 S² + D² |B|²) W - (S² + D² |B|²) |B|² > 0
|
||
!
|
||
! coming from the equation of state
|
||
!
|
||
! using analytical. It takes the maximum estimated root as the lower bracket.
|
||
! If the analytical estimation fails, the Newton-Raphson iterative method
|
||
! is used.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, mb - input coefficients for |B|² and m.B, respectively;
|
||
! wl, wu - the lower and upper limits for the enthalpy;
|
||
! wc - the initial root guess;
|
||
! info - the flag is .true. if the initial brackets and guess were found,
|
||
! otherwise it is .false.;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn &
|
||
, wl, wu, wc, fl, fu, info)
|
||
|
||
use algebra, only : cubic_normalized, quartic
|
||
use helpers, only : print_message
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(out) :: wl, wu, wc, fl, fu
|
||
logical , intent(out) :: info
|
||
|
||
logical :: keep
|
||
integer :: it, nr
|
||
real(kind=8) :: dd, ss, ec
|
||
real(kind=8) :: f , df
|
||
real(kind=8) :: dw, err
|
||
|
||
real(kind=8), dimension(5) :: a
|
||
real(kind=8), dimension(4) :: x
|
||
|
||
character(len=80) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::nr_initial_brackets_srmhd_adi()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! calculate temporary variables
|
||
!
|
||
dd = dn * dn
|
||
ss = mb * mb
|
||
ec = en + pmin
|
||
|
||
! set the initial upper bracket
|
||
!
|
||
wu = en + pmin
|
||
|
||
! calculate the cubic equation coefficients for the positivity condition
|
||
! coming from the energy equation; the condition, in fact, finds the minimum
|
||
! enthalphy for which the pressure is equal to pmin
|
||
!
|
||
a(3) = 2.5d+00 * bb - ec
|
||
a(2) = 2.0d+00 * (bb - ec) * bb
|
||
a(1) = 0.5d+00 * (((bb - 2.0d+00 * ec) * bb + mm) * bb - ss)
|
||
|
||
! solve the cubic equation
|
||
!
|
||
nr = cubic_normalized(a(1:3), x(1:3))
|
||
|
||
! if solution was found, use the maximum root as the lower bracket
|
||
!
|
||
if (nr > 0) then
|
||
|
||
wl = x(nr)
|
||
|
||
! calculate the quartic equation coefficients for the positivity condition
|
||
! coming from the pressure equation
|
||
!
|
||
a(5) = 1.0d+00
|
||
a(4) = 2.0d+00 * bb
|
||
a(3) = bb * bb - dd - mm
|
||
a(2) = - 2.0d+00 * (ss + dd * bb)
|
||
a(1) = - (ss + dd * bb) * bb
|
||
|
||
! solve the quartic equation
|
||
!
|
||
nr = quartic(a(1:5), x(1:4))
|
||
|
||
! take the maximum ethalpy from both conditions to guarantee that the pressure
|
||
! obtains from any of those equations is positive
|
||
!
|
||
if (nr > 0) wl = max(wl, x(nr))
|
||
|
||
else ! nr = 0
|
||
|
||
! the root could not be found analytically, so use the iterative solver
|
||
! to find the lower bracket; as the initial guess use the initial upper bracket
|
||
!
|
||
keep = .true.
|
||
it = nrmax
|
||
wl = wu
|
||
do while(keep)
|
||
call nr_pressure_srmhd_adi_1d(mm, bb, mb, ec, wl, f, df)
|
||
dw = f / df
|
||
wl = wl - dw
|
||
err = abs(dw / wl)
|
||
it = it - 1
|
||
keep = (err > tol) .and. it > 0
|
||
end do
|
||
if (it <= 0) then
|
||
call print_message(loc, &
|
||
"Iterative solver failed to find the lower bracket!")
|
||
write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ", &
|
||
dn, mm, mb, bb, en
|
||
call print_message(loc, msg)
|
||
end if
|
||
|
||
end if ! nr > 0
|
||
|
||
! check if the energy function is negative for the lower limit
|
||
!
|
||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wl, fl)
|
||
|
||
if (fl <= 0.0d+00) then
|
||
|
||
! make sure that the upper limit is larger than the lower one
|
||
!
|
||
if (wu <= wl) wu = 2.0d+00 * wl
|
||
|
||
! check if the brackets bound the root region, if not proceed until
|
||
! opposite function signs are found for the brackets
|
||
!
|
||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
|
||
it = nrmax
|
||
keep = fl * fu > 0.0d+00
|
||
do while (keep)
|
||
it = it - 1
|
||
wl = wu
|
||
fl = fu
|
||
wu = 2.0d+00 * wu
|
||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, wu, fu)
|
||
keep = (fl * fu > 0.0d+00) .and. it > 0
|
||
end do
|
||
|
||
! the upper bracket was found, so proceed with determining the root
|
||
!
|
||
if (it > 0 .and. fu >= 0.0d+00) then
|
||
|
||
! estimate the enthalpy value close to the root
|
||
!
|
||
wc = wl - fl * (wu - wl) / (fu - fl)
|
||
|
||
! we have good brackets and guess, so good to go
|
||
!
|
||
info = .true.
|
||
|
||
else ! the upper brack not found
|
||
call print_message(loc, "Could not find the upper bracket!")
|
||
write(msg,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = ", &
|
||
dn, mm, mb, bb, en
|
||
call print_message(loc, msg)
|
||
info = .false.
|
||
|
||
end if
|
||
|
||
else ! the root cannot be found, since it is below the lower bracket
|
||
call print_message(loc, "Positive function for lower bracket!")
|
||
write(msg,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = ", &
|
||
dn, mm, mb, bb, en, wl
|
||
call print_message(loc, msg)
|
||
info = .false.
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_initial_brackets_srmhd_adi
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRMHD_ADI_1DW:
|
||
! -----------------------------------
|
||
!
|
||
! Subroutine finds a root W of equation
|
||
!
|
||
! F(W) = W - P(W) + ½ [(1 + |V|²) |B|² - S² / W²] - E = 0
|
||
!
|
||
! using the Newton-Raphson 1Dw iterative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w , vv - input/output coefficients W and |V|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
! References:
|
||
!
|
||
! Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
|
||
! "Primitive Variable Solvers for Conservative General Relativistic
|
||
! Magnetohydrodynamics",
|
||
! The Astrophysical Journal, 2006, vol. 641, pp. 626-637
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srmhd_adi_1dw(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
use helpers, only : print_message
|
||
|
||
implicit none
|
||
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: f , df, dw
|
||
real(kind=8) :: err
|
||
|
||
character(len=80) :: msg
|
||
|
||
character(len=*), parameter :: loc = 'EQUATIONS::nr_iterate_srmhd_adi_1dw()'
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! find the initial brackets and estimate the initial enthalpy
|
||
!
|
||
call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn &
|
||
, wl, wu, w, fl, fu, info)
|
||
|
||
! continue if brackets found
|
||
!
|
||
if (info) then
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! iterate using the Newton-Raphson method in order to find a root w of the
|
||
! function
|
||
!
|
||
do while(keep)
|
||
|
||
! calculate F(W) and dF(W)/dW
|
||
!
|
||
call nr_function_srmhd_adi_1d(mm, bb, mb, en, dn, w, f, df)
|
||
|
||
! update brackets
|
||
!
|
||
if (f > fl .and. f < 0.0d+00) then
|
||
wl = w
|
||
fl = f
|
||
end if
|
||
if (f < fu .and. f > 0.0d+00) then
|
||
wu = w
|
||
fu = f
|
||
end if
|
||
|
||
! calculate the increment dW, update the solution, and estimate the error
|
||
!
|
||
dw = f / df
|
||
w = w - dw
|
||
err = abs(dw / w)
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! if new W lays out of the brackets, use the bisection method to estimate
|
||
! the new guess
|
||
!
|
||
if (w < wl .or. w > wu) then
|
||
w = 0.5d+00 * (wl + wu)
|
||
end if
|
||
|
||
! decrease the number of remaining iterations
|
||
!
|
||
it = it - 1
|
||
|
||
end do ! NR iterations
|
||
|
||
! let know the user if the convergence failed
|
||
!
|
||
if (err >= tol) then
|
||
call print_message(loc, "Convergence not reached!")
|
||
write(msg,"(a,1x,1es24.16e3)") "Error: ", err
|
||
call print_message(loc, msg)
|
||
end if
|
||
|
||
! calculate |V|² from W
|
||
!
|
||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||
|
||
end if ! correct brackets
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srmhd_adi_1dw
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRMHD_ADI_2DWV:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine finds a root (W, |V|²) of equations
|
||
!
|
||
! F(W,|V|²) = W - P + ½ [(1 + |V|²) |B|² - S² / W²] - E = 0
|
||
! G(W,|V|²) = |V|² (|B|² + W)² - S² (|B|² + 2W) / W² - |M|² = 0
|
||
!
|
||
! using the Newton-Raphson 2D iterative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w , vv - input/output coefficients W and |V|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
! References:
|
||
!
|
||
! Noble, S. C., Gammie, C. F., McKinney, J. C, Del Zanna, L.,
|
||
! "Primitive Variable Solvers for Conservative General Relativistic
|
||
! Magnetohydrodynamics",
|
||
! The Astrophysical Journal, 2006, vol. 641, pp. 626-637
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srmhd_adi_2dwv(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
! local variables
|
||
!
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: vm, vs, wt, mw, wt2
|
||
real(kind=8) :: f, dfw, dfv
|
||
real(kind=8) :: g, dgw, dgv
|
||
real(kind=8) :: det, jfw, jfv, jgw, jgv
|
||
real(kind=8) :: dv, dw
|
||
real(kind=8) :: err
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! find the initial brackets and estimate the initial enthalpy
|
||
!
|
||
call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn &
|
||
, wl, wu, w, fl, fu, info)
|
||
|
||
! if the brackets could not be found, return the lower bracket as the solution
|
||
!
|
||
if (.not. info) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||
write(*,"(a,1x)" ) "The solution lays in unphysical regime."
|
||
write(*,"(a,1x,1es24.16e3)") "Using the lower bracket as solution: ", wl
|
||
|
||
! use the lower bracket, since it guarantees the positive pressure
|
||
!
|
||
w = wl
|
||
|
||
! calculate |V|² from W
|
||
!
|
||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||
|
||
info = .true.
|
||
return
|
||
|
||
end if
|
||
|
||
! and the corresponding |V|²
|
||
!
|
||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! find root with the help of the Newton-Raphson method
|
||
!
|
||
do while(keep)
|
||
|
||
! calculate (S/W)², Wt, Wt²
|
||
!
|
||
mw = (mb / w)**2
|
||
wt = w + bb
|
||
wt2 = wt * wt
|
||
|
||
! prepare (1 - |V|²) and sqrt(1 - |V|²)
|
||
!
|
||
vm = 1.0d+00 - vv
|
||
vs = sqrt(vm)
|
||
|
||
! calculate functions F(W,|V|²) and G(W,|V|²)
|
||
!
|
||
f = w - en - gammaxi * (w * vm - dn * vs) &
|
||
+ 0.5d+00 * ((1.0d+00 + vv) * bb - mw)
|
||
g = vv * wt2 - (wt + w) * mw - mm
|
||
|
||
! calculate derivatives dF(W,|V|²)/dW and dF(W,|V|²)/d|V|²
|
||
!
|
||
dfw = 1.0d+00 - gammaxi * vm + mw / w
|
||
dfv = - gammaxi * (0.5d+00 * dn / vs - w) + 0.5d+00 * bb
|
||
|
||
! calculate derivatives dG(W,|V|²)/dW and dG(W,|V|²)/d|V|²
|
||
!
|
||
dgw = 2.0d+00 * wt * (vv + mw / w)
|
||
dgv = wt2
|
||
|
||
! invert the Jacobian J = | dF/dW, dF/d|V|² |
|
||
! | dG/dW, dG/d|V|² |
|
||
!
|
||
det = dfw * dgv - dfv * dgw
|
||
|
||
jfw = dgv / det
|
||
jgw = - dfv / det
|
||
jfv = - dgw / det
|
||
jgv = dfw / det
|
||
|
||
! calculate increments dW and d|V|²
|
||
!
|
||
dw = f * jfw + g * jgw
|
||
dv = f * jfv + g * jgv
|
||
|
||
! correct W and |V|²
|
||
!
|
||
w = w - dw
|
||
vv = vv - dv
|
||
|
||
! check if the new enthalpy and velocity are physical
|
||
!
|
||
if (w < wl) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
|
||
write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
|
||
info = .false.
|
||
return
|
||
end if
|
||
if (vv < 0.0d+00 .or. vv >= 1.0d+00) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
|
||
write(*,"(a,1x,1es24.16e3)") "Unphysical speed |v|²: ", vv
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! calculate the error
|
||
!
|
||
err = max(abs(dw / w), abs(dv))
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! decrease the number of remaining iterations
|
||
!
|
||
it = it - 1
|
||
|
||
end do ! NR iterations
|
||
|
||
! let know the user if the convergence failed
|
||
!
|
||
if (err >= tol) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwv()"
|
||
write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srmhd_adi_2dwv
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine NR_ITERATE_SRMHD_ADI_2DWU:
|
||
! ------------------------------------
|
||
!
|
||
! Subroutine finds a root (W, |u|²) of equations
|
||
!
|
||
! F(W,|u|²) = W - E - P + ½ [(1 + |u|² / (1 + |u|²)) |B|² - S² / W²] = 0
|
||
! G(W,|u|²) = (|B|² + W)² |u|² / (1 + |u|²) - (2W + |B|²) S² / W² - |M|² = 0
|
||
!
|
||
! using the Newton-Raphson 2D iterative method.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! mm, en - input coefficients for |M|² and E, respectively;
|
||
! bb, bm - input coefficients for |B|² and B.M, respectively;
|
||
! w , vv - input/output coefficients W and |v|²;
|
||
! info - the flag is .true. if the solution was found, otherwise
|
||
! it is .false.;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine nr_iterate_srmhd_adi_2dwu(mm, bb, mb, en, dn, w, vv, info)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! input/output arguments
|
||
!
|
||
real(kind=8), intent(in) :: mm, bb, mb, en, dn
|
||
real(kind=8), intent(inout) :: w, vv
|
||
logical , intent(out) :: info
|
||
|
||
! local variables
|
||
!
|
||
logical :: keep
|
||
integer :: it, cn
|
||
real(kind=8) :: wl, wu, fl, fu
|
||
real(kind=8) :: uu, up, gm
|
||
real(kind=8) :: ss, ww, wt, wt2, wd, wp
|
||
real(kind=8) :: f, dfw, dfu
|
||
real(kind=8) :: g, dgw, dgu
|
||
real(kind=8) :: det, jfw, jfu, jgw, jgu
|
||
real(kind=8) :: dw, du
|
||
real(kind=8) :: err
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
! find the initial brackets and estimate the initial enthalpy
|
||
!
|
||
call nr_initial_brackets_srmhd_adi(mm, bb, mb, en, dn &
|
||
, wl, wu, w, fl, fu, info)
|
||
|
||
! if the brackets could not be found, return the lower bracket as the solution
|
||
!
|
||
if (.not. info) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_1dw()"
|
||
write(*,"(a,1x)" ) "The solution lays in unphysical regime."
|
||
write(*,"(a,1x,1es24.16e3)") "Using the lower bracket as solution: ", wl
|
||
|
||
! use the lower bracket, since it guarantees the positive pressure
|
||
!
|
||
w = wl
|
||
|
||
! calculate |V|² from W
|
||
!
|
||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||
|
||
info = .true.
|
||
return
|
||
|
||
end if
|
||
|
||
! and the corresponding |u|²
|
||
!
|
||
call nr_velocity_srmhd_adi_1d(mm, bb, mb, w, vv)
|
||
uu = vv / (1.0d+00 - vv)
|
||
|
||
! initialize iteration parameters
|
||
!
|
||
info = .true.
|
||
keep = .true.
|
||
it = nrmax
|
||
cn = nrext
|
||
|
||
! find root with the help of the Newton-Raphson method
|
||
!
|
||
do while(keep)
|
||
|
||
! prepare (1 + |u|²) and the Lorentz factor
|
||
!
|
||
up = 1.0d+00 + uu
|
||
gm = sqrt(up)
|
||
|
||
! calculate temporary variables
|
||
!
|
||
ss = mb * mb
|
||
ww = w * w
|
||
wt = w + bb
|
||
wt2 = wt * wt
|
||
wd = w - gm * dn
|
||
wp = wt / up
|
||
|
||
! calculate functions F(W,|u|²) and G(W,|u|²)
|
||
!
|
||
f = w - en - gammaxi * wd / up &
|
||
+ 0.5d+00 * (bb * (1.0d+00 + uu / up) - ss / ww)
|
||
g = wp * wt * uu - (w + wt) * ss / ww - mm
|
||
|
||
! calculate derivatives dF(W,|u|²)/dW and dF(W,|u|²)/d|u|²
|
||
!
|
||
dfw = 1.0d+00 - gammaxi / up + ss / ww / w
|
||
dfu = 0.5d+00 * (gammaxi * (w + wd) + bb) / up**2
|
||
|
||
! calculate derivatives dG(W,|u|²)/dW and dG(W,|u|²)/d|u|²
|
||
!
|
||
dgw = 2.0d+00 * wt * (uu / up + ss / ww / w)
|
||
dgu = wp * wp
|
||
|
||
! invert the Jacobian J = | dF/dW, dF/d|u|² |
|
||
! | dG/dW, dG/d|u|² |
|
||
!
|
||
det = dfw * dgu - dfu * dgw
|
||
|
||
jfw = dgu / det
|
||
jgw = - dfu / det
|
||
jfu = - dgw / det
|
||
jgu = dfw / det
|
||
|
||
! calculate increments dW and d|u|²
|
||
!
|
||
dw = f * jfw + g * jgw
|
||
du = f * jfu + g * jgu
|
||
|
||
! correct W and |u|²
|
||
!
|
||
w = w - dw
|
||
uu = uu - du
|
||
|
||
! check if the new enthalpy and velocity are physical
|
||
!
|
||
if (w < wl) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
|
||
write(*,"(a,1x,2es24.16e3)") "Enthalpy smaller than the limit: ", w, wl
|
||
info = .false.
|
||
return
|
||
end if
|
||
if (uu < 0.0d+00) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "ERROR in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
|
||
write(*,"(a,1x,1es24.16e3)") "Unphysical speed |u|²: ", uu
|
||
info = .false.
|
||
return
|
||
end if
|
||
|
||
! calculate the error
|
||
!
|
||
err = max(abs(dw / w), abs(du))
|
||
|
||
! check the convergence, if the convergence is not reached, iterate until
|
||
! the maximum number of iteration is reached
|
||
!
|
||
if (err < tol) then
|
||
keep = cn > 0
|
||
cn = cn - 1
|
||
else
|
||
keep = it > 0
|
||
end if
|
||
|
||
! decrease the number of remaining iterations
|
||
!
|
||
it = it - 1
|
||
|
||
end do ! NR iterations
|
||
|
||
! calculate |v|² from |u|²
|
||
!
|
||
vv = uu / (1.0d+00 + uu)
|
||
|
||
! let know the user if the convergence failed
|
||
!
|
||
if (err >= tol) then
|
||
write(*,*)
|
||
write(*,"(a,1x,a)" ) "WARNING in" &
|
||
, "EQUATIONS::nr_iterate_srmhd_adi_2dwu()"
|
||
write(*,"(a,1x,1es24.16e3)") "Convergence not reached: ", err
|
||
end if
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine nr_iterate_srmhd_adi_2dwu
|
||
|
||
!===============================================================================
|
||
!
|
||
end module equations
|