2484 lines
64 KiB
Fortran
2484 lines
64 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-2024 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: PROBLEMS
|
|
!!
|
|
!! This module handles the initialization of various test and research
|
|
!! problems.
|
|
!!
|
|
!!
|
|
!!******************************************************************************
|
|
!
|
|
module problems
|
|
|
|
implicit none
|
|
|
|
character(len=64), save :: problem_name = "none"
|
|
|
|
private
|
|
|
|
public :: initialize_problems, finalize_problems
|
|
public :: problem_name
|
|
|
|
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
!
|
|
contains
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PUBLIC SUBROUTINES *****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine INITIALIZE_PROBLEMS:
|
|
! ------------------------------
|
|
!
|
|
! Subroutine prepares module PROBLEMS.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! problem - the problem name
|
|
! rcount - the run count for restarted jobs
|
|
! verbose - a logical flag turning the information printing;
|
|
! status - an integer flag for error return value;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine initialize_problems(problem, rcount, verbose, status)
|
|
|
|
use mesh , only : setup_problem
|
|
use parameters , only : get_parameter
|
|
|
|
implicit none
|
|
|
|
character(len=64), intent(in) :: problem
|
|
integer , intent(in) :: rcount
|
|
logical , intent(in) :: verbose
|
|
integer , intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
problem_name = problem
|
|
|
|
select case(trim(problem))
|
|
|
|
case("riemann")
|
|
setup_problem => setup_problem_riemann
|
|
|
|
case("blast")
|
|
setup_problem => setup_problem_blast
|
|
|
|
case("st", "sedov-taylor", "ST", "Sedov-Taylor")
|
|
setup_problem => setup_problem_sedov_taylor
|
|
|
|
case("implosion")
|
|
setup_problem => setup_problem_implosion
|
|
|
|
case("kh", "kelvinhelmholtz", "kelvin-helmholtz", "Kelvin-Helmholtz")
|
|
setup_problem => setup_problem_kelvin_helmholtz
|
|
|
|
case("rt", "rayleightaylor", "rayleigh-taylor")
|
|
setup_problem => setup_problem_rayleigh_taylor
|
|
|
|
case("ot", "orszag-tang", "Orszag-Tang")
|
|
setup_problem => setup_problem_orszag_tang
|
|
|
|
case("current-sheet", "current_sheet")
|
|
setup_problem => setup_problem_current_sheet
|
|
|
|
case("tearing-mode", "tearing_mode", "tearing")
|
|
setup_problem => setup_problem_tearing
|
|
|
|
case("turbulence")
|
|
setup_problem => setup_problem_turbulence
|
|
|
|
case("stellar_wind", "stellar-wind", "wind")
|
|
setup_problem => setup_problem_wind
|
|
|
|
case default
|
|
|
|
end select
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine initialize_problems
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine FINALIZE_PROBLEMS:
|
|
! ----------------------------
|
|
!
|
|
! Subroutine releases memory used by the module.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - an integer flag for error return value;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine finalize_problems(status)
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine finalize_problems
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PRIVATE SUBROUTINES ****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_RIEMANN:
|
|
! --------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the general Riemann problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_riemann(pdata)
|
|
|
|
use blocks , only : block_data
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, adx
|
|
use equations , only : prim2cons
|
|
use equations , only : nv
|
|
use equations , only : qpbnd
|
|
|
|
implicit none
|
|
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
integer :: p, i, j, k
|
|
real(kind=8) :: xl, xr
|
|
real(kind=8) :: dx, dxh
|
|
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare block coordinates
|
|
!
|
|
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
|
|
|
! calculate mesh intervals and areas
|
|
!
|
|
dx = adx(pdata%meta%level)
|
|
dxh = 0.5d+00 * dx
|
|
|
|
! set the left and right states of the primitive variables
|
|
!
|
|
do i = 1, nn
|
|
xl = x(i) - dxh
|
|
xr = x(i) + dxh
|
|
|
|
if (xr <= 0.0d+00) then
|
|
do p = 1, nv
|
|
q(p,i) = qpbnd(p,1,1)
|
|
end do
|
|
else if (xl >= 0.0d+00) then
|
|
do p = 1, nv
|
|
q(p,i) = qpbnd(p,1,2)
|
|
end do
|
|
else
|
|
do p = 1, nv
|
|
q(p,i) = (xr * qpbnd(p,1,2) - xl * qpbnd(p,1,1)) / dx
|
|
end do
|
|
end if
|
|
end do ! i = 1, im
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! iterate over all positions in the YZ plane
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_riemann
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_BLAST:
|
|
! ------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the spherical blast problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_blast(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : d2r
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay, adx, ady, advol
|
|
#if NDIMS == 3
|
|
use coordinates, only : az, adz
|
|
#endif /* NDIMS == 3 */
|
|
use equations , only : prim2cons
|
|
use equations , only : adiabatic_index
|
|
use equations , only : nv
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use parameters , only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: ratio = 1.00d+02
|
|
real(kind=8), save :: radius = 1.00d-01
|
|
real(kind=8), save :: csnd = 4.0824829046386301635d-01
|
|
real(kind=8), save :: buni = 1.00d+00
|
|
real(kind=8), save :: angle = 4.50d+01
|
|
#if NDIMS == 3
|
|
integer , save :: nsubgrid = 10
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! local saved parameters
|
|
!
|
|
logical , save :: first = .true.
|
|
real(kind=8), save :: dn_amb, dn_ovr
|
|
real(kind=8), save :: pr_amb, pr_ovr
|
|
real(kind=8), save :: bx, by
|
|
real(kind=8), save :: r2
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
#if NDIMS == 3
|
|
integer :: ic, jc, kc
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: xl, yl, xu, yu, rl, ru
|
|
real(kind=8) :: dx, dy, dxh, dyh, dvol
|
|
real(kind=8) :: sn
|
|
#if NDIMS == 3
|
|
real(kind=8) :: zl, zu, dz, dzh
|
|
real(kind=8) :: xb, yb, zb
|
|
real(kind=8) :: xt, yt, zt
|
|
real(kind=8) :: fc_inc
|
|
#else /* NDIMS == 3 */
|
|
real(kind=8) :: rlu, rul
|
|
real(kind=8) :: xb, yb
|
|
real(kind=8) :: xt, yt
|
|
real(kind=8) :: ph
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: fc_amb, fc_ovr
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y
|
|
#if NDIMS == 3
|
|
real(kind=8), dimension(nn) :: z
|
|
|
|
! allocatable arrays
|
|
!
|
|
real(kind=8), dimension(:), allocatable :: xm, ym, zm
|
|
real(kind=8), dimension(:), allocatable :: xp, yp, zp
|
|
#endif /* NDIMS == 3 */
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare problem constants during the first subroutine call
|
|
!
|
|
if (first) then
|
|
|
|
! get problem parameters
|
|
!
|
|
call get_parameter("dens" , dens )
|
|
call get_parameter("ratio" , ratio )
|
|
call get_parameter("radius" , radius)
|
|
call get_parameter("sound_speed", csnd )
|
|
call get_parameter("buni" , buni )
|
|
call get_parameter("angle" , angle )
|
|
|
|
#if NDIMS == 3
|
|
! get the fine grid resolution
|
|
!
|
|
call get_parameter("nsubgrid", nsubgrid)
|
|
|
|
! correct subgrid resolution if necessary
|
|
!
|
|
nsubgrid = max(1, nsubgrid)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate the overdense and ambient region densities
|
|
!
|
|
dn_amb = dens
|
|
if (ipr > 0) then
|
|
dn_ovr = dn_amb
|
|
|
|
! calculate parallel and perpendicular pressures from sound speeds
|
|
!
|
|
pr_amb = dens * csnd * csnd / adiabatic_index
|
|
pr_ovr = pr_amb * ratio
|
|
|
|
else
|
|
dn_ovr = dn_amb * ratio
|
|
end if
|
|
|
|
! calculate initial uniform field components
|
|
!
|
|
if (ibx > 0) then
|
|
sn = sin(d2r * angle)
|
|
bx = buni * sqrt(1.0d+00 - sn * sn)
|
|
by = buni * sn
|
|
end if
|
|
|
|
! calculate the square of radius
|
|
!
|
|
r2 = radius * radius
|
|
|
|
! reset the first execution flag
|
|
!
|
|
first = .false.
|
|
|
|
end if ! first call
|
|
|
|
! prepare block coordinates
|
|
!
|
|
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 */
|
|
|
|
! calculate mesh intervals and areas
|
|
!
|
|
dx = adx(pdata%meta%level)
|
|
dy = ady(pdata%meta%level)
|
|
#if NDIMS == 3
|
|
dz = adz(pdata%meta%level)
|
|
#endif /* NDIMS == 3 */
|
|
dxh = 0.5d+00 * dx
|
|
dyh = 0.5d+00 * dy
|
|
#if NDIMS == 3
|
|
dzh = 0.5d+00 * dz
|
|
#endif /* NDIMS == 3 */
|
|
dvol = advol(pdata%meta%level)
|
|
|
|
#if NDIMS == 3
|
|
! allocate subgrid coordinates
|
|
!
|
|
allocate(xm(nsubgrid), ym(nsubgrid), zm(nsubgrid))
|
|
allocate(xp(nsubgrid), yp(nsubgrid), zp(nsubgrid))
|
|
|
|
! and generate them
|
|
!
|
|
xm(:) = (1.0d+00 * (/(i, i = 0, nsubgrid - 1)/)) / nsubgrid
|
|
ym(:) = xm(:)
|
|
zm(:) = xm(:)
|
|
xm(:) = xm(:) * dx
|
|
ym(:) = ym(:) * dy
|
|
zm(:) = zm(:) * dz
|
|
xp(:) = (1.0d+00 * (/(i, i = 1, nsubgrid )/)) / nsubgrid
|
|
yp(:) = xp(:)
|
|
zp(:) = xp(:)
|
|
xp(:) = xp(:) * dx
|
|
yp(:) = yp(:) * dy
|
|
zp(:) = zp(:) * dz
|
|
|
|
! calculate the factor increment for the given subgrid
|
|
!
|
|
fc_inc = dvol / nsubgrid**3
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! set the ambient density and pressure
|
|
!
|
|
q(idn,:) = dn_amb
|
|
if (ipr > 0) q(ipr,:) = pr_amb
|
|
|
|
! reset velocity components
|
|
!
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
! if magnetic field is present, set it to be uniform with the desired strength
|
|
! and orientation
|
|
!
|
|
if (ibx > 0) then
|
|
|
|
q(ibx,:) = bx
|
|
q(iby,:) = by
|
|
q(ibz,:) = 0.0d+00
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
end if
|
|
|
|
! iterate over all positions in the YZ plane
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
|
|
! calculate the corner Z coordinates
|
|
!
|
|
zl = abs(z(k)) - dzh
|
|
zu = abs(z(k)) + dzh
|
|
#endif /* NDIMS == 3 */
|
|
|
|
do j = 1, nn
|
|
|
|
! calculate the corner Y coordinates
|
|
!
|
|
yl = abs(y(j)) - dyh
|
|
yu = abs(y(j)) + dyh
|
|
|
|
! sweep along the X coordinate
|
|
!
|
|
do i = 1, nn
|
|
|
|
! calculate the corner X coordinates
|
|
!
|
|
xl = abs(x(i)) - dxh
|
|
xu = abs(x(i)) + dxh
|
|
|
|
! calculate the minimum and maximum corner distances from the origin
|
|
!
|
|
#if NDIMS == 3
|
|
rl = xl * xl + yl * yl + zl * zl
|
|
ru = xu * xu + yu * yu + zu * zu
|
|
#else /* NDIMS == 3 */
|
|
rl = xl * xl + yl * yl
|
|
ru = xu * xu + yu * yu
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! set the initial density and pressure in cells laying completely within
|
|
! the blast radius
|
|
!
|
|
if (ru <= r2) then
|
|
|
|
! set the overpressure region density
|
|
!
|
|
q(idn,i) = dn_ovr
|
|
|
|
! set the overpressure region pressure
|
|
!
|
|
if (ipr > 0) q(ipr,i) = pr_ovr
|
|
|
|
! set the initial pressure in the cell completely outside the radius
|
|
!
|
|
else if (rl >= r2) then
|
|
|
|
! set the ambient region density
|
|
!
|
|
q(idn,i) = dn_amb
|
|
|
|
! set the ambient medium pressure
|
|
!
|
|
if (ipr > 0) q(ipr,i) = pr_amb
|
|
|
|
! integrate density or pressure in cells which are crossed by the circule with
|
|
! the given radius
|
|
!
|
|
else
|
|
|
|
#if NDIMS == 3
|
|
! interpolate the factor using subgrid
|
|
!
|
|
fc_ovr = 0.0d+00
|
|
do kc = 1, nsubgrid
|
|
zb = (zl + zm(kc))**2
|
|
zt = (zl + zp(kc))**2
|
|
do jc = 1, nsubgrid
|
|
yb = (yl + ym(jc))**2
|
|
yt = (yl + yp(jc))**2
|
|
do ic = 1, nsubgrid
|
|
xb = (xl + xm(ic))**2
|
|
xt = (xl + xp(ic))**2
|
|
|
|
! update the integration factor depending on the subcell position
|
|
!
|
|
if ((xt + yt + zt) <= r2) then
|
|
fc_ovr = fc_ovr + fc_inc
|
|
else if ((xb + yb + zb) < r2) then
|
|
fc_ovr = fc_ovr + 0.5d+00 * fc_inc
|
|
end if
|
|
|
|
end do ! ic = 1, nsubgrid
|
|
end do ! jc = 1, nsubgrid
|
|
end do ! kc = 1, nsubgrid
|
|
#else /* NDIMS == 3 */
|
|
|
|
! calculate the distance of remaining corners
|
|
!
|
|
rlu = xl * xl + yu * yu
|
|
rul = xu * xu + yl * yl
|
|
|
|
! separate in the cases of which corners lay inside, and which outside
|
|
! the radius
|
|
!
|
|
if (min(rlu, rul) >= r2) then
|
|
|
|
! only one cell corner inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xb = sqrt(r2 - yl**2) - xl
|
|
yb = sqrt(r2 - xl**2) - yl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(xb**2 + yb**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * (xb * yb + (ph - sn) * r2)
|
|
|
|
else if (rlu >= r2) then
|
|
|
|
! two lower corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
yb = sqrt(r2 - xl**2) - yl
|
|
yt = sqrt(r2 - xu**2) - yl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(dx**2 + (yt - yb)**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * ((yt + yb) * dx + (ph - sn) * r2)
|
|
|
|
else if (rul >= r2) then
|
|
|
|
! two left corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xb = sqrt(r2 - yl**2) - xl
|
|
xt = sqrt(r2 - yu**2) - xl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt((xt - xb)**2 + dy**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * ((xt + xb) * dy + (ph - sn) * r2)
|
|
|
|
else
|
|
|
|
! three corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xt = xu - sqrt(r2 - yu**2)
|
|
yt = yu - sqrt(r2 - xu**2)
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(xt**2 + yt**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = dvol - 0.5d+00 * (xt * yt - (ph - sn) * r2)
|
|
|
|
end if
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! normalize coefficients
|
|
!
|
|
fc_ovr = fc_ovr / dvol
|
|
fc_amb = 1.0d+00 - fc_ovr
|
|
|
|
! integrate the density over the edge cells
|
|
!
|
|
q(idn,i) = fc_ovr * dn_ovr + fc_amb * dn_amb
|
|
|
|
! integrate the pressure over the edge cells
|
|
!
|
|
if (ipr > 0) q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb
|
|
|
|
end if
|
|
|
|
end do ! i = 1, nn
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
|
|
! deallocate subgrid coordinates
|
|
!
|
|
deallocate(xm, ym, zm)
|
|
deallocate(xp, yp, zp)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_blast
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_SEDOV_TAYLOR:
|
|
! -------------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the spherical Sedov-Taylor
|
|
! blast problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
! References:
|
|
!
|
|
! [1] Almgren, A. S. et al.,
|
|
! "CASTRO: A New Compressible Astrophysical Solver.
|
|
! I. Hydrodynamics and Self-Gravity",
|
|
! The Astrophysical Journal, 2010, vol. 715, pp. 1221-1238,
|
|
! http://dx.doi.org/10.1088/0004-637X/715/2/1221
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_sedov_taylor(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : d2r
|
|
#if NDIMS == 3
|
|
use constants , only : pi4
|
|
#else /* NDIMS == 3 */
|
|
use constants , only : pi
|
|
#endif /* NDIMS == 3 */
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay, adx, ady, advol
|
|
#if NDIMS == 3
|
|
use coordinates, only : az, adz
|
|
#endif /* NDIMS == 3 */
|
|
use equations , only : prim2cons
|
|
use equations , only : adiabatic_index
|
|
use equations , only : nv
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use parameters , only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: radius = 1.00d-02
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: pres = 1.00d-05
|
|
real(kind=8), save :: eexp = 1.00d+00
|
|
real(kind=8), save :: buni = 0.00d+00
|
|
real(kind=8), save :: angle = 0.00d+00
|
|
#if NDIMS == 3
|
|
integer , save :: nsubgrid = 10
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! local saved parameters
|
|
!
|
|
logical , save :: first = .true.
|
|
real(kind=8), save :: dn_amb, dn_ovr
|
|
real(kind=8), save :: pr_amb, pr_ovr
|
|
real(kind=8), save :: bx, by
|
|
real(kind=8), save :: r2
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
#if NDIMS == 3
|
|
integer :: ic, jc, kc
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: xl, yl, xu, yu, rl, ru
|
|
real(kind=8) :: sn
|
|
#if NDIMS == 3
|
|
real(kind=8) :: zl, zu, dz, dzh
|
|
real(kind=8) :: xb, yb, zb
|
|
real(kind=8) :: xt, yt, zt
|
|
real(kind=8) :: fc_inc
|
|
#else /* NDIMS == 3 */
|
|
real(kind=8) :: rlu, rul
|
|
real(kind=8) :: xb, yb
|
|
real(kind=8) :: xt, yt
|
|
real(kind=8) :: ph
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: dx, dy, dxh, dyh, dvol
|
|
real(kind=8) :: fc_amb, fc_ovr
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y
|
|
#if NDIMS == 3
|
|
real(kind=8), dimension(nn) :: z
|
|
|
|
! allocatable arrays
|
|
!
|
|
real(kind=8), dimension(:), allocatable :: xm, ym, zm
|
|
real(kind=8), dimension(:), allocatable :: xp, yp, zp
|
|
#endif /* NDIMS == 3 */
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare problem constants during the first subroutine call
|
|
!
|
|
if (first) then
|
|
|
|
! get problem parameters
|
|
!
|
|
call get_parameter("radius", radius)
|
|
call get_parameter("dens" , dens )
|
|
call get_parameter("pres" , pres )
|
|
call get_parameter("eexp" , eexp )
|
|
call get_parameter("buni" , buni )
|
|
call get_parameter("angle" , angle )
|
|
|
|
#if NDIMS == 3
|
|
! get the fine grid resolution
|
|
!
|
|
call get_parameter("nsubgrid", nsubgrid)
|
|
|
|
! correct subgrid resolution if necessary
|
|
!
|
|
nsubgrid = max(1, nsubgrid)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate the volume of the injection region
|
|
!
|
|
#if NDIMS == 3
|
|
dvol = pi4 * radius**3 / 3.0d+00
|
|
#else /* NDIMS == 3 */
|
|
dvol = pi * radius**2
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate the overdense and ambient region densities and pressures
|
|
!
|
|
dn_amb = dens
|
|
dn_ovr = dn_amb
|
|
pr_amb = pres
|
|
pr_ovr = (adiabatic_index - 1.0d+00) * eexp / dvol
|
|
|
|
! calculate initial uniform field components
|
|
!
|
|
if (ibx > 0) then
|
|
sn = sin(d2r * angle)
|
|
bx = buni * sqrt(1.0d+00 - sn * sn)
|
|
by = buni * sn
|
|
end if
|
|
|
|
! calculate the square of radius
|
|
!
|
|
r2 = radius * radius
|
|
|
|
! reset the first execution flag
|
|
!
|
|
first = .false.
|
|
|
|
end if ! first call
|
|
|
|
! prepare block coordinates
|
|
!
|
|
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 */
|
|
|
|
! calculate mesh intervals and areas
|
|
!
|
|
dx = adx(pdata%meta%level)
|
|
dy = ady(pdata%meta%level)
|
|
#if NDIMS == 3
|
|
dz = adz(pdata%meta%level)
|
|
#endif /* NDIMS == 3 */
|
|
dxh = 0.5d+00 * dx
|
|
dyh = 0.5d+00 * dy
|
|
#if NDIMS == 3
|
|
dzh = 0.5d+00 * dz
|
|
#endif /* NDIMS == 3 */
|
|
dvol = advol(pdata%meta%level)
|
|
|
|
#if NDIMS == 3
|
|
! allocate subgrid coordinates
|
|
!
|
|
allocate(xm(nsubgrid), ym(nsubgrid), zm(nsubgrid))
|
|
allocate(xp(nsubgrid), yp(nsubgrid), zp(nsubgrid))
|
|
|
|
! and generate them
|
|
!
|
|
xm(:) = (1.0d+00 * (/(i, i = 0, nsubgrid - 1)/)) / nsubgrid
|
|
ym(:) = xm(:)
|
|
zm(:) = xm(:)
|
|
xm(:) = xm(:) * dx
|
|
ym(:) = ym(:) * dy
|
|
zm(:) = zm(:) * dz
|
|
xp(:) = (1.0d+00 * (/(i, i = 1, nsubgrid )/)) / nsubgrid
|
|
yp(:) = xp(:)
|
|
zp(:) = xp(:)
|
|
xp(:) = xp(:) * dx
|
|
yp(:) = yp(:) * dy
|
|
zp(:) = zp(:) * dz
|
|
|
|
! calculate the factor increment for the given subgrid
|
|
!
|
|
fc_inc = dvol / nsubgrid**3
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! set density and pressure of the ambient
|
|
!
|
|
q(idn,:) = dn_amb
|
|
if (ipr > 0) q(ipr,:) = pr_amb
|
|
|
|
! reset velocity components
|
|
!
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
! if magnetic field is present, set it to be uniform with the desired strength
|
|
! and orientation
|
|
!
|
|
if (ibx > 0) then
|
|
|
|
q(ibx,:) = bx
|
|
q(iby,:) = by
|
|
q(ibz,:) = 0.0d+00
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
end if
|
|
|
|
! iterate over all positions in the YZ plane
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
|
|
! calculate the corner Z coordinates
|
|
!
|
|
zl = abs(z(k)) - dzh
|
|
zu = abs(z(k)) + dzh
|
|
#endif /* NDIMS == 3 */
|
|
|
|
do j = 1, nn
|
|
|
|
! calculate the corner Y coordinates
|
|
!
|
|
yl = abs(y(j)) - dyh
|
|
yu = abs(y(j)) + dyh
|
|
|
|
! sweep along the X coordinate
|
|
!
|
|
do i = 1, nn
|
|
|
|
! calculate the corner X coordinates
|
|
!
|
|
xl = abs(x(i)) - dxh
|
|
xu = abs(x(i)) + dxh
|
|
|
|
! calculate the minimum and maximum corner distances from the origin
|
|
!
|
|
#if NDIMS == 3
|
|
rl = xl * xl + yl * yl + zl * zl
|
|
ru = xu * xu + yu * yu + zu * zu
|
|
#else /* NDIMS == 3 */
|
|
rl = xl * xl + yl * yl
|
|
ru = xu * xu + yu * yu
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! set the initial density and pressure in cells laying completely within
|
|
! the blast radius
|
|
!
|
|
if (ru <= r2) then
|
|
|
|
! set density and pressure for the overpressure region
|
|
!
|
|
q(idn,i) = dn_ovr
|
|
if (ipr > 0) q(ipr,i) = pr_ovr
|
|
|
|
! set the initial pressure in the cell completely outside the radius
|
|
!
|
|
else if (rl >= r2) then
|
|
|
|
! set density and pressure of the ambient
|
|
!
|
|
q(idn,i) = dn_amb
|
|
if (ipr > 0) q(ipr,i) = pr_amb
|
|
|
|
! integrate density or pressure in cells which are crossed by the circule with
|
|
! the given radius
|
|
!
|
|
else
|
|
|
|
#if NDIMS == 3
|
|
! interpolate the factor using subgrid
|
|
!
|
|
fc_ovr = 0.0d+00
|
|
do kc = 1, nsubgrid
|
|
zb = (zl + zm(kc))**2
|
|
zt = (zl + zp(kc))**2
|
|
do jc = 1, nsubgrid
|
|
yb = (yl + ym(jc))**2
|
|
yt = (yl + yp(jc))**2
|
|
do ic = 1, nsubgrid
|
|
xb = (xl + xm(ic))**2
|
|
xt = (xl + xp(ic))**2
|
|
|
|
! update the integration factor depending on the subcell position
|
|
!
|
|
if ((xt + yt + zt) <= r2) then
|
|
fc_ovr = fc_ovr + fc_inc
|
|
else if ((xb + yb + zb) < r2) then
|
|
fc_ovr = fc_ovr + 0.5d+00 * fc_inc
|
|
end if
|
|
|
|
end do ! ic = 1, nsubgrid
|
|
end do ! jc = 1, nsubgrid
|
|
end do ! kc = 1, nsubgrid
|
|
#else /* NDIMS == 3 */
|
|
|
|
! calculate the distance of remaining corners
|
|
!
|
|
rlu = xl * xl + yu * yu
|
|
rul = xu * xu + yl * yl
|
|
|
|
! separate in the cases of which corners lay inside, and which outside
|
|
! the radius
|
|
!
|
|
if (min(rlu, rul) >= r2) then
|
|
|
|
! only one cell corner inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xb = sqrt(r2 - yl**2) - xl
|
|
yb = sqrt(r2 - xl**2) - yl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(xb**2 + yb**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * (xb * yb + (ph - sn) * r2)
|
|
|
|
else if (rlu >= r2) then
|
|
|
|
! two lower corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
yb = sqrt(r2 - xl**2) - yl
|
|
yt = sqrt(r2 - xu**2) - yl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(dx**2 + (yt - yb)**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * ((yt + yb) * dx + (ph - sn) * r2)
|
|
|
|
else if (rul >= r2) then
|
|
|
|
! two left corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xb = sqrt(r2 - yl**2) - xl
|
|
xt = sqrt(r2 - yu**2) - xl
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt((xt - xb)**2 + dy**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = 0.5d+00 * ((xt + xb) * dy + (ph - sn) * r2)
|
|
|
|
else
|
|
|
|
! three corners inside the radius
|
|
!
|
|
! calculate middle coordinates of the radius-edge crossing point
|
|
!
|
|
xt = xu - sqrt(r2 - yu**2)
|
|
yt = yu - sqrt(r2 - xu**2)
|
|
|
|
! calculate the sin(½φ), φ, and sin(φ)
|
|
!
|
|
sn = 0.5d+00 * sqrt(xt**2 + yt**2) / radius
|
|
ph = 2.0d+00 * asin(sn)
|
|
sn = sin(ph)
|
|
|
|
! calculate the area of cell intersection with the radius
|
|
!
|
|
fc_ovr = dvol - 0.5d+00 * (xt * yt - (ph - sn) * r2)
|
|
|
|
end if
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! normalize coefficients
|
|
!
|
|
fc_ovr = fc_ovr / dvol
|
|
fc_amb = 1.0d+00 - fc_ovr
|
|
|
|
! integrate density and pressure over the edge cells
|
|
!
|
|
q(idn,i) = fc_ovr * dn_ovr + fc_amb * dn_amb
|
|
if (ipr > 0) q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb
|
|
|
|
end if
|
|
|
|
end do ! i = 1, nn
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
|
|
! deallocate subgrid coordinates
|
|
!
|
|
deallocate(xm, ym, zm)
|
|
deallocate(xp, yp, zp)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_sedov_taylor
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_IMPLOSION:
|
|
! ----------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the implosion problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_implosion(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data, ndims
|
|
use constants , only : d2r
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay, adx, ady
|
|
#if NDIMS == 3
|
|
use coordinates, only : az, adz
|
|
#endif /* NDIMS == 3 */
|
|
use equations , only : prim2cons
|
|
use equations , only : nv
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use parameters , only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: sline = 1.50d-01
|
|
real(kind=8), save :: adens = 1.00d+00
|
|
real(kind=8), save :: apres = 1.00d+00
|
|
real(kind=8), save :: drat = 1.25d-01
|
|
real(kind=8), save :: prat = 1.40d-01
|
|
real(kind=8), save :: buni = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
real(kind=8), save :: angle = 0.00d+00
|
|
|
|
! local saved parameters
|
|
!
|
|
logical , save :: first = .true.
|
|
real(kind=8), save :: odens = 1.25d-01
|
|
real(kind=8), save :: opres = 1.40d-01
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
real(kind=8) :: rl, ru, dx, dy, dxh, dyh, ds, dl, dr
|
|
#if NDIMS == 3
|
|
real(kind=8) :: dz, dzh
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: sn, cs
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, xl, xu
|
|
real(kind=8), dimension(nn) :: y, yl, yu
|
|
#if NDIMS == 3
|
|
real(kind=8), dimension(nn) :: z, zl, zu
|
|
#endif /* NDIMS == 3 */
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare problem constants during the first subroutine call
|
|
!
|
|
if (first) then
|
|
|
|
! get problem parameters
|
|
!
|
|
call get_parameter("shock_line" , sline )
|
|
call get_parameter("ambient_density" , adens )
|
|
call get_parameter("ambient_pressure", apres )
|
|
call get_parameter("density_ratio" , drat )
|
|
call get_parameter("pressure_ratio" , prat )
|
|
call get_parameter("buni" , buni )
|
|
call get_parameter("bgui" , bgui )
|
|
call get_parameter("angle" , angle )
|
|
|
|
! calculate parameters
|
|
!
|
|
odens = drat * adens
|
|
opres = prat * apres
|
|
|
|
! reset the first execution flag
|
|
!
|
|
first = .false.
|
|
|
|
end if ! first call
|
|
|
|
! prepare block coordinates
|
|
!
|
|
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 */
|
|
|
|
! calculate mesh intervals and areas
|
|
!
|
|
dx = adx(pdata%meta%level)
|
|
dy = ady(pdata%meta%level)
|
|
#if NDIMS == 3
|
|
dz = adz(pdata%meta%level)
|
|
#endif /* NDIMS == 3 */
|
|
dxh = 0.5d+00 * dx
|
|
dyh = 0.5d+00 * dy
|
|
#if NDIMS == 3
|
|
dzh = 0.5d+00 * dz
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate edge coordinates
|
|
!
|
|
xl(:) = abs(x(:)) - dxh
|
|
xu(:) = abs(x(:)) + dxh
|
|
yl(:) = abs(y(:)) - dyh
|
|
yu(:) = abs(y(:)) + dyh
|
|
#if NDIMS == 3
|
|
zl(:) = abs(z(:)) - dzh
|
|
zu(:) = abs(z(:)) + dzh
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! reset velocity components
|
|
!
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
! if magnetic field is present, set it to be uniform with the desired strength
|
|
! and orientation
|
|
!
|
|
if (ibx > 0) then
|
|
|
|
! calculate the orientation angles
|
|
!
|
|
sn = sin(d2r * angle)
|
|
cs = sqrt(1.0d+00 - sn * sn)
|
|
|
|
! set magnetic field components
|
|
!
|
|
q(ibx,:) = buni * cs
|
|
q(iby,:) = buni * sn
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
end if
|
|
|
|
! iterate over all positions
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
do i = 1, nn
|
|
|
|
! calculate the distance from the origin
|
|
!
|
|
#if NDIMS == 3
|
|
rl = xl(i) + yl(j) + zl(k)
|
|
ru = xu(i) + yu(j) + zu(k)
|
|
#else /* NDIMS == 3 */
|
|
rl = xl(i) + yl(j)
|
|
ru = xu(i) + yu(j)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! initialize density and pressure
|
|
!
|
|
if (ru <= sline) then
|
|
q(idn,i) = odens
|
|
if (ipr > 0) q(ipr,i) = opres
|
|
else if (rl >= sline) then
|
|
q(idn,i) = adens
|
|
if (ipr > 0) q(ipr,i) = apres
|
|
else
|
|
ds = (sline - rl) / dx
|
|
if (ds <= 1.0d+00) then
|
|
dl = 5.0d-01 * ds**ndims
|
|
dr = 1.0d+00 - dl
|
|
else
|
|
ds = (ru - sline) / dx
|
|
dr = 5.0d-01 * ds**ndims
|
|
dl = 1.0d+00 - dr
|
|
end if
|
|
|
|
q(idn,i) = adens * dl + odens * dr
|
|
if (ipr > 0) q(ipr,i) = apres * dl + opres * dr
|
|
end if
|
|
|
|
end do ! i = 1, im
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_implosion
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_KELVIN_HELMHOLTZ:
|
|
! -----------------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the Kelvin-Helmholtz instability
|
|
! problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_kelvin_helmholtz(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : d2r
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ay, ady
|
|
use equations , only : prim2cons
|
|
use equations , only : nv
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
|
|
use parameters , only : get_parameter
|
|
use random , only : randsym
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: ycut = 2.50d-01
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: drat = 2.00d+00
|
|
real(kind=8), save :: pres = 2.50d+00
|
|
real(kind=8), save :: vamp = 1.00d+00
|
|
real(kind=8), save :: vper = 1.00d-02
|
|
real(kind=8), save :: buni = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
real(kind=8), save :: angle = 0.00d+00
|
|
|
|
! local saved parameters
|
|
!
|
|
logical , save :: first = .true.
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
real(kind=8) :: yl, yu, dy, dyh
|
|
real(kind=8) :: sn, cs
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: y
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare problem constants during the first subroutine call
|
|
!
|
|
if (first) then
|
|
|
|
! get problem parameters
|
|
!
|
|
call get_parameter("ycut" , ycut )
|
|
call get_parameter("dens" , dens )
|
|
call get_parameter("drat" , drat )
|
|
call get_parameter("pres" , pres )
|
|
call get_parameter("vamp" , vamp )
|
|
call get_parameter("vper" , vper )
|
|
call get_parameter("buni" , buni )
|
|
call get_parameter("bgui" , bgui )
|
|
call get_parameter("angle" , angle )
|
|
|
|
! reset the first execution flag
|
|
!
|
|
first = .false.
|
|
|
|
end if ! first call
|
|
|
|
! prepare block coordinates
|
|
!
|
|
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
|
|
|
! calculate mesh intervals and areas
|
|
!
|
|
dy = ady(pdata%meta%level)
|
|
dyh = 0.5d+00 * dy
|
|
|
|
! set the ambient density and pressure
|
|
!
|
|
q(idn,:) = dens
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
|
|
! if magnetic field is present, set it to be uniform with the desired strength
|
|
! and orientation
|
|
!
|
|
if (ibx > 0) then
|
|
|
|
! calculate the orientation angles
|
|
!
|
|
sn = sin(d2r * angle)
|
|
cs = sqrt(1.0d+00 - sn * sn)
|
|
|
|
! set magnetic field components
|
|
!
|
|
q(ibx,:) = buni * cs
|
|
q(iby,:) = buni * sn
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
end if
|
|
|
|
! iterate over all positions in the YZ plane
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
! calculate the corner Y coordinates
|
|
!
|
|
yl = abs(y(j)) - dyh
|
|
yu = abs(y(j)) + dyh
|
|
|
|
! set the primitive variables for two regions
|
|
!
|
|
if (yu <= ycut) then
|
|
q(idn,:) = dens * drat
|
|
q(ivx,:) = vamp
|
|
if (isl > 0) q(isl,:) = 1.0d+00
|
|
else if (yl >= ycut) then
|
|
q(idn,:) = dens
|
|
q(ivx,:) = - vamp
|
|
if (isl > 0) q(isl,:) = 0.0d+00
|
|
else
|
|
q(idn,:) = dens * ((yu - ycut) * drat + (ycut - yl)) / dy
|
|
q(ivx,:) = vamp * ((yu - ycut) - (ycut - yl)) / dy
|
|
if (isl > 0) q(isl,:) = (yu - ycut) / dy
|
|
end if
|
|
|
|
! reset remaining velocity components
|
|
!
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
! set the pressure
|
|
!
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
|
|
! add a random seed velocity component
|
|
!
|
|
do i = 1, nn
|
|
q(ivx,i) = q(ivx,i) + vper * randsym()
|
|
q(ivy,i) = q(ivy,i) + vper * randsym()
|
|
#if NDIMS == 3
|
|
q(ivz,i) = q(ivz,i) + vper * randsym()
|
|
#endif /* NDIMS == 3 */
|
|
end do
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:), .true.)
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_kelvin_helmholtz
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_RAYLEIGH_TAYLOR:
|
|
! ----------------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the Rayleigh-Taylor instability
|
|
! problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
! References:
|
|
!
|
|
! [1] Almgren, A. S. et al.,
|
|
! "CASTRO: A New Compressible Astrophysical Solver.
|
|
! I. Hydrodynamics and Self-Gravity",
|
|
! The Astrophysical Journal, 2010, vol. 715, pp. 1221-1238,
|
|
! http://dx.doi.org/10.1088/0004-637X/715/2/1221
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_rayleigh_taylor(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : pi2, d2r
|
|
use coordinates, only : xmin, xmax, xlen
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay
|
|
use equations , only : prim2cons
|
|
use equations , only : nv
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
|
|
use equations , only : csnd2
|
|
use parameters , only : get_parameter
|
|
use random , only : randsym
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: drat = 2.00d+00
|
|
real(kind=8), save :: damp = 5.00d-01
|
|
real(kind=8), save :: pres = 5.00d+00
|
|
real(kind=8), save :: ycut = 0.00d+00
|
|
real(kind=8), save :: vper = 0.00d+00
|
|
real(kind=8), save :: lper = 1.00d-02
|
|
real(kind=8), save :: kper = 1.00d+00
|
|
real(kind=8), save :: hdel = 5.00d-03
|
|
real(kind=8), save :: gacc = -1.00d+00
|
|
real(kind=8), save :: buni = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
real(kind=8), save :: angle = 0.00d+00
|
|
|
|
! local saved parameters
|
|
!
|
|
logical , save :: first = .true.
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
real(kind=8) :: sn, cs
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y, yp
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! prepare problem constants during the first subroutine call
|
|
!
|
|
if (first) then
|
|
|
|
! get problem parameters
|
|
!
|
|
call get_parameter("ycut" , ycut )
|
|
call get_parameter("dens" , dens )
|
|
call get_parameter("drat" , drat )
|
|
call get_parameter("pres" , pres )
|
|
call get_parameter("lper" , lper )
|
|
call get_parameter("kper" , kper )
|
|
call get_parameter("vper" , vper )
|
|
call get_parameter("hdel" , hdel )
|
|
call get_parameter("gacc" , gacc )
|
|
call get_parameter("buni" , buni )
|
|
call get_parameter("bgui" , bgui )
|
|
call get_parameter("angle" , angle )
|
|
|
|
! calculate the density change across the interface
|
|
!
|
|
damp = 5.0d-01 * (drat * dens - dens)
|
|
|
|
! reset the first execution flag
|
|
!
|
|
first = .false.
|
|
|
|
end if ! first call
|
|
|
|
! prepare block coordinates
|
|
!
|
|
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
|
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
|
|
|
! set the ambient density and pressure
|
|
!
|
|
q(idn,:) = dens
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
|
|
! if magnetic field is present, set it to be uniform with the desired strength
|
|
! and orientation
|
|
!
|
|
if (ibx > 0) then
|
|
|
|
! calculate the orientation angles
|
|
!
|
|
sn = sin(d2r * angle)
|
|
cs = sqrt(1.0d+00 - sn * sn)
|
|
|
|
! set magnetic field components
|
|
!
|
|
q(ibx,:) = buni * cs
|
|
q(iby,:) = buni * sn
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
end if
|
|
|
|
! prepare density perturbation
|
|
!
|
|
yp(:) = 0.5d+00 * lper * (cos(pi2 * kper * (x(:) - xmin) / xlen) &
|
|
+ cos(pi2 * kper * (xmax - x(:)) / xlen)) + ycut
|
|
|
|
! iterate over all positions in the YZ plane
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
! set density and pressure
|
|
!
|
|
if (ipr > 0) then
|
|
if (y(j) <= ycut) then
|
|
q(idn,:) = dens
|
|
else
|
|
q(idn,:) = dens * drat
|
|
end if
|
|
q(idn,:) = dens + damp * (1.0d+00 + tanh((y(j) - yp(:)) / hdel))
|
|
q(ipr,:) = pres + q(idn,:) * gacc * y(j)
|
|
else
|
|
if (y(j) <= ycut) then
|
|
q(idn,:) = dens * exp(gacc * y(j) / csnd2)
|
|
else
|
|
q(idn,:) = dens * drat * exp(gacc * y(j) / csnd2)
|
|
end if
|
|
end if
|
|
if (isl > 0) then
|
|
if (y(j) <= ycut) then
|
|
q(isl,:) = 0.0d+00
|
|
else
|
|
q(isl,:) = 1.0d+00
|
|
end if
|
|
end if
|
|
|
|
! reset the velocity components
|
|
!
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
! add a random seed velocity component
|
|
!
|
|
if (abs(vper) > 0.0d+00) then
|
|
do i = 1, nn
|
|
q(ivy,i) = q(ivy,i) + vper * randsym()
|
|
end do
|
|
end if
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:), .true.)
|
|
|
|
! copy the conserved variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
|
|
! copy the primitive variables to the current block
|
|
!
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_rayleigh_taylor
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_ORSZAG_TANG:
|
|
! ------------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the Orszag-Tang test
|
|
! problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
! References:
|
|
!
|
|
! [1] Orszag, S. A., Tang, C. -M.,
|
|
! "Small-scale structure of two-dimensional magnetohydrodynamic
|
|
! turbulence",
|
|
! Journal of Fluid Mechanics, 1979, 90, pp. 129-143,
|
|
! https://doi.org/10.1017/S002211207900210X
|
|
! [2] Mignone, A.,
|
|
! "A simple and accurate Riemann solver for isothermal MHD",
|
|
! Journal of Computational Physics, 2007, 225, pp. 1427-1441,
|
|
! https://doi.org/10.1016/j.jcp.2007.01.033
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_orszag_tang(pdata)
|
|
|
|
use blocks , only : block_data
|
|
use constants , only : pi2
|
|
use coordinates, only : nn => bcells, ax, ay
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use equations , only : nv, magnetized, csnd, adiabatic_index
|
|
use equations , only : prim2cons
|
|
use helpers , only : print_message
|
|
use parameters , only : get_parameter
|
|
|
|
implicit none
|
|
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: pres = 1.00d+00
|
|
real(kind=8), save :: vamp = 1.00d+00
|
|
real(kind=8), save :: bamp = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
|
|
logical, save :: first = .true.
|
|
|
|
integer :: j, k
|
|
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y
|
|
|
|
character(len=*), parameter :: loc = 'PROBLEMS::setup_problem_orszag_tang()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
if (first) then
|
|
|
|
call get_parameter("dens", dens)
|
|
call get_parameter("vamp", vamp)
|
|
call get_parameter("bamp", bamp)
|
|
call get_parameter("bgui", bgui)
|
|
|
|
if (ipr > 0) then
|
|
pres = dens / adiabatic_index
|
|
bamp = bamp / adiabatic_index
|
|
else
|
|
vamp = csnd * vamp
|
|
bamp = csnd * bamp
|
|
end if
|
|
|
|
first = .false.
|
|
end if
|
|
|
|
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
|
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
|
|
|
q(idn,:) = dens
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
q(ivy,:) = - vamp * sin(pi2 * x(:))
|
|
q(ivz,:) = 0.0d+00
|
|
|
|
if (magnetized) then
|
|
q(iby,:) = bamp * sin(pi2 * x(:) * 2.0d+00)
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
q(ivx,:) = vamp * sin(pi2 * y(j))
|
|
q(ibx,:) = bamp * sin(pi2 * y(j))
|
|
|
|
call prim2cons(q(:,:), u(:,:), .true.)
|
|
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do
|
|
#if NDIMS == 3
|
|
end do
|
|
#endif /* NDIMS == 3 */
|
|
else
|
|
call print_message(loc, "Orszag-Tang problem requires MHD!")
|
|
end if
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_orszag_tang
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_CURRENT_SHEET:
|
|
! --------------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the current sheet test problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_current_sheet(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : pi2
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay
|
|
use equations , only : prim2cons
|
|
use equations , only : nv, magnetized
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use parameters , only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: xpos = 2.50d-01
|
|
real(kind=8), save :: xdel = 1.00d-08
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: beta = 1.00d-01
|
|
real(kind=8), save :: vamp = 1.00d-01
|
|
real(kind=8), save :: bamp = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
|
|
! local saved parameters
|
|
!
|
|
logical, save :: first = .true.
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
real(kind=8) :: t
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y, vx
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! retrieve the problem parameters during the first subroutine call
|
|
!
|
|
if (first) then
|
|
call get_parameter("position" , xpos)
|
|
call get_parameter("thickness" , xdel)
|
|
call get_parameter("density" , dens)
|
|
call get_parameter("beta" , beta)
|
|
call get_parameter("vel_amplitude", vamp)
|
|
call get_parameter("mag_amplitude", bamp)
|
|
call get_parameter("mag_guide" , bgui)
|
|
|
|
first = .false.
|
|
end if
|
|
|
|
! prepare the block coordinates
|
|
!
|
|
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
|
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
|
|
|
! prepare the vector of velocity perturbation (varying along Y)
|
|
!
|
|
vx(:) = vamp * sin(pi2 * y(:))
|
|
|
|
! set the fields, along X, which do not vary along the YZ planes
|
|
!
|
|
q(idn,:) = dens
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
if (ipr > 0) q(ipr,:) = 0.5d+00 * beta
|
|
|
|
! if magnetic field is present, set its antiparallel configuration
|
|
!
|
|
if (magnetized) then
|
|
|
|
! set magnetic field components
|
|
!
|
|
q(ibx,:) = 0.0d+00
|
|
do i = 1, nn
|
|
t = (xpos - abs(x(i))) / xdel
|
|
q(iby,i) = sign(bamp, t) * min(abs(t), 1.0d+00)
|
|
end do
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
|
|
! correct the gas pressure due to the variance of magnetic pressure
|
|
!
|
|
if (ipr > 0) q(ipr,:) = q(ipr,:) + 0.5d+00 * (bamp**2 - q(iby,:)**2)
|
|
|
|
end if
|
|
|
|
! iterate over all positions in the YZ planes
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
! set the velocity perturbation
|
|
!
|
|
q(ivx,:) = vx(j)
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! copy the conserved and primitive variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_current_sheet
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_TEARING:
|
|
! --------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the resistive tearing instability
|
|
! test problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
! References:
|
|
!
|
|
! [1] Lapenta, G.,
|
|
! "Self-Feeding Turbulent Magnetic Reconnection on Macroscopic Scales",
|
|
! Physical Review Letters, 2008, vol. 100, id. 235001,
|
|
! http://dx.doi.org/10.1103/PhysRevLett.100.235001
|
|
!
|
|
! [2] Landi, S., Del Zanna, L., Papini, E., Pucci, F., and Velli, M.
|
|
! "Resistive Magnetohydrodynamics Simulations of the Ideal Tearing Mode",
|
|
! The Astrophysical Journal, 2015, vol. 806, id. 131,
|
|
! http://dx.doi.org/10.1088/0004-637X/806/1/131
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_tearing(pdata)
|
|
|
|
! include external procedures and variables
|
|
!
|
|
use blocks , only : block_data
|
|
use constants , only : pi2
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay, adx, ady
|
|
use equations , only : prim2cons
|
|
use equations , only : nv, magnetized, csnd2
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use parameters , only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! input arguments
|
|
!
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
! default parameter values
|
|
!
|
|
real(kind=8), save :: beta = 1.00d-01
|
|
real(kind=8), save :: delta = 1.00d-08
|
|
real(kind=8), save :: eta = 1.00d-08
|
|
real(kind=8), save :: zeta = 0.00d+00
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: pres = 5.00d-02
|
|
real(kind=8), save :: pmag = 1.00d+00
|
|
real(kind=8), save :: bamp = 1.00d+00
|
|
real(kind=8), save :: bnor = 0.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
real(kind=8), save :: bper = 1.00d-02
|
|
real(kind=8), save :: vper = 1.00d-03
|
|
real(kind=8), save :: kper = 1.00d+00
|
|
real(kind=8), save :: ss = 1.00d+00
|
|
|
|
! local saved parameters
|
|
!
|
|
logical, save :: first = .true.
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k
|
|
real(kind=8) :: dx, dy, t
|
|
|
|
! local arrays
|
|
!
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: xc, xl, xu, yl, yu, xi
|
|
real(kind=8), dimension(nn) :: vx, vy, bx, by, b0
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! retrieve the problem parameters during the first subroutine call
|
|
!
|
|
if (first) then
|
|
call get_parameter("resistivity", eta)
|
|
call get_parameter("beta" , beta)
|
|
call get_parameter("zeta" , zeta)
|
|
call get_parameter("density" , dens)
|
|
call get_parameter("bamp" , bamp)
|
|
call get_parameter("bnor" , bnor)
|
|
call get_parameter("bgui" , bgui)
|
|
call get_parameter("bper" , bper)
|
|
call get_parameter("vper" , vper)
|
|
call get_parameter("kper" , kper)
|
|
|
|
ss = abs(bamp) / (sqrt(dens) * max(1.0d-16, eta))
|
|
delta = ss**(-1.0d+00/3.0d+00)
|
|
if (abs(bper) > 0.0d+00) vper = 0.0d+00
|
|
pres = 0.5d+00 * beta
|
|
pmag = 0.5d+00 * (bamp**2 + bnor**2 + bgui**2)
|
|
|
|
first = .false.
|
|
end if
|
|
|
|
! prepare the block coordinates
|
|
!
|
|
dx = adx(pdata%meta%level)
|
|
dy = ady(pdata%meta%level)
|
|
xc(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
|
xl(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - 0.5d+00 * dx
|
|
xu(:) = xl(:) + dx
|
|
yl(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - 0.5d+00 * dy
|
|
yu(:) = yl(:) + dy
|
|
|
|
! prepare the vector of velocity perturbation (varying along Y)
|
|
!
|
|
if (abs(vper) > 0.0d+00) then
|
|
xi(:) = xc(:) * sqrt(ss)
|
|
vx(:) = vper * tanh(xi(:)) * exp(- xi(:)**2)
|
|
vy(:) = vper * (2.0d+00 * xi(:) * tanh(xi(:)) - 1.0d+00 / cosh(xi(:))**2)&
|
|
* exp(- xi(:)**2) * sqrt(ss) / (pi2 * kper)
|
|
end if
|
|
if (abs(bper) > 0.0d+00) then
|
|
bx(:) = bper * (sin(pi2 * xl(:)) - sin(pi2 * xu(:))) / (pi2 * dx)
|
|
by(:) = bper * (cos(pi2 * xl(:)) - cos(pi2 * xu(:))) / (pi2 * dx * kper)
|
|
end if
|
|
|
|
! set the fields, along X, which do not vary along the YZ planes
|
|
!
|
|
q(idn,:) = dens
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
|
|
! if magnetic field is present, set its antiparallel configuration
|
|
!
|
|
if (magnetized) then
|
|
|
|
! set magnetic field components
|
|
!
|
|
q(ibx,:) = bnor
|
|
do i = 1, nn
|
|
t = delta * (log(cosh(xu(i) / delta)) &
|
|
- log(cosh(xl(i) / delta))) / dx
|
|
q(iby,i) = bamp * max(-1.0d+00, min(1.0d+00, t))
|
|
q(ibz,i) = zeta * sqrt(bamp**2 - q(iby,i)**2) + bgui
|
|
end do
|
|
q(ibp,:) = 0.0d+00
|
|
b0(:) = q(iby,:)
|
|
|
|
! correct the gas pressure due to the variance of magnetic pressure
|
|
!
|
|
if (ipr > 0) then
|
|
q(ipr,:) = pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1)
|
|
else
|
|
q(idn,:) = (pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1)) / csnd2
|
|
end if
|
|
|
|
end if
|
|
|
|
! iterate over all positions in the YZ planes
|
|
!
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
|
|
! set the perturbation
|
|
!
|
|
if (abs(vper) > 0.0d+00) then
|
|
q(ivx,:) = vx(:) * (sin(pi2 * kper * yl(j)) &
|
|
- sin(pi2 * kper * yu(j))) / (pi2 * kper * dy)
|
|
q(ivy,:) = vy(:) * (cos(pi2 * kper * yu(j)) &
|
|
- cos(pi2 * kper * yl(j))) / (pi2 * kper * dy)
|
|
end if
|
|
if (abs(bper) > 0.0d+00) then
|
|
q(ibx,:) = bx(:) * (cos(pi2 * kper * yu(j)) &
|
|
- cos(pi2 * kper * yl(j))) / (pi2 * kper * dy) &
|
|
+ bnor
|
|
q(iby,:) = by(:) * (sin(pi2 * kper * yl(j)) &
|
|
- sin(pi2 * kper * yu(j))) / (pi2 * kper * dy) &
|
|
+ b0(:)
|
|
end if
|
|
|
|
! convert the primitive variables to conservative ones
|
|
!
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
! copy the conserved and primitive variables to the current block
|
|
!
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
|
|
end do ! j = 1, nn
|
|
#if NDIMS == 3
|
|
end do ! k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_tearing
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_TURBULENCE:
|
|
! -----------------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the turbulence problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_turbulence(pdata)
|
|
|
|
use blocks , only : block_data
|
|
use coordinates, only : nn => bcells
|
|
use equations , only : nv, magnetized, adiabatic_index, csnd2
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use equations , only : prim2cons
|
|
use parameters , only : get_parameter
|
|
|
|
implicit none
|
|
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
real(kind=8), save :: beta = 1.00d+00
|
|
real(kind=8), save :: dens = 1.00d+00
|
|
real(kind=8), save :: bamp = 1.00d+00
|
|
real(kind=8), save :: bgui = 0.00d+00
|
|
real(kind=8), save :: pres = 5.00d-02
|
|
real(kind=8), save :: pmag = 1.00d+00
|
|
|
|
logical, save :: first = .true.
|
|
|
|
integer :: j, k
|
|
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
if (first) then
|
|
call get_parameter("density", dens)
|
|
|
|
pres = csnd2 * dens
|
|
if (ipr > 0) pres = pres / adiabatic_index
|
|
|
|
if (magnetized) then
|
|
call get_parameter("bamp", bamp)
|
|
call get_parameter("bgui", bgui)
|
|
|
|
pmag = 0.5d+00 * (bamp**2 + bgui**2)
|
|
beta = pres / pmag
|
|
end if
|
|
|
|
first = .false.
|
|
end if
|
|
|
|
q(idn,:) = dens
|
|
q(ivx,:) = 0.0d+00
|
|
q(ivy,:) = 0.0d+00
|
|
q(ivz,:) = 0.0d+00
|
|
if (ipr > 0) q(ipr,:) = pres
|
|
|
|
if (magnetized) then
|
|
q(ibx,:) = bamp
|
|
q(iby,:) = 0.0d+00
|
|
q(ibz,:) = bgui
|
|
q(ibp,:) = 0.0d+00
|
|
end if
|
|
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
end do
|
|
#if NDIMS == 3
|
|
end do
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_turbulence
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_PROBLEM_WIND:
|
|
! -----------------------------
|
|
!
|
|
! Subroutine sets the initial conditions for the stellar wind problem.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pdata - pointer to the datablock structure of the currently initialized
|
|
! block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_problem_wind(pdata)
|
|
|
|
use blocks , only : block_data
|
|
use coordinates, only : nn => bcells
|
|
use coordinates, only : ax, ay
|
|
#if NDIMS == 3
|
|
use coordinates, only : az
|
|
#endif /* NDIMS == 3 */
|
|
use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2
|
|
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
|
use equations , only : prim2cons
|
|
use parameters , only : get_parameter
|
|
|
|
implicit none
|
|
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
|
|
real(kind=8), save :: radius = 1.00d-01
|
|
real(kind=8), save :: dens = 1.00d+06
|
|
real(kind=8), save :: sonic = 3.00d+00
|
|
real(kind=8), save :: wind = 5.00d-01
|
|
real(kind=8), save :: pres = 5.00d-02
|
|
real(kind=8), save :: bamp = 1.00d-03
|
|
|
|
logical, save :: first = .true.
|
|
|
|
integer :: i, j, k
|
|
real(kind=8) :: r, r2, s, w
|
|
|
|
real(kind=8), dimension(nv,nn) :: q, u
|
|
real(kind=8), dimension(nn) :: x, y
|
|
#if NDIMS == 3
|
|
real(kind=8), dimension(nn) :: z
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
if (first) then
|
|
call get_parameter("radius" , radius)
|
|
call get_parameter("density" , dens )
|
|
call get_parameter("sonic" , sonic )
|
|
call get_parameter("wind_speed", wind )
|
|
|
|
if (ipr > 0) then
|
|
pres = (wind / sonic)**2 * dens / adiabatic_index
|
|
else
|
|
csnd = wind / sonic
|
|
csnd2 = csnd * csnd
|
|
end if
|
|
|
|
if (magnetized) &
|
|
call get_parameter("bamp", bamp)
|
|
|
|
first = .false.
|
|
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 */
|
|
|
|
#if NDIMS == 3
|
|
do k = 1, nn
|
|
#else /* NDIMS == 3 */
|
|
k = 1
|
|
#endif /* NDIMS == 3 */
|
|
do j = 1, nn
|
|
do i = 1, nn
|
|
#if NDIMS == 3
|
|
r2 = x(i)**2 + y(j)**2 + z(k)**2
|
|
#else /* NDIMS == 3 */
|
|
r2 = x(i)**2 + y(j)**2
|
|
#endif /* NDIMS == 3 */
|
|
r = max(sqrt(r2), 1.0d-16)
|
|
s = max(1.0d+00, r / radius)**(NDIMS - 1)
|
|
|
|
q(idn,i) = dens / s
|
|
if (ipr > 0) q(ipr,i) = pres / s**adiabatic_index
|
|
|
|
q(ivx,i) = wind * x(i) / r
|
|
q(ivy,i) = wind * y(j) / r
|
|
#if NDIMS == 3
|
|
q(ivz,i) = wind * z(k) / r
|
|
#else /* NDIMS == 3 */
|
|
q(ivz,i) = 0.0d+00
|
|
#endif /* NDIMS == 3 */
|
|
|
|
if (magnetized) then
|
|
w = max(1.0d+00, r / radius)**(5.0d-1 * NDIMS)
|
|
|
|
q(ibx,i) = 0.0d+00
|
|
q(iby,i) = 0.0d+00
|
|
q(ibz,i) = bamp / w
|
|
q(ibp,i) = 0.0d+00
|
|
end if
|
|
end do
|
|
|
|
call prim2cons(q(:,:), u(:,:))
|
|
|
|
pdata%q(:,:,j,k) = q(:,:)
|
|
pdata%u(:,:,j,k) = u(:,:)
|
|
end do
|
|
#if NDIMS == 3
|
|
end do
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_problem_wind
|
|
|
|
!===============================================================================
|
|
!
|
|
end module problems
|