Merge branch 'master' into reconnection
This commit is contained in:
commit
c7fb30839f
@ -1,6 +1,6 @@
|
||||
--------------------------------------------------------------------------------
|
||||
# **The AMUN Code**
|
||||
## Copyright (C) 2008-2017 Grzegorz Kowal ##
|
||||
## Copyright (C) 2008-2018 Grzegorz Kowal ##
|
||||
--------------------------------------------------------------------------------
|
||||
|
||||
AMUN is a parallel code to perform numerical simulations in fluid approximation
|
||||
@ -65,8 +65,9 @@ the HDF5 libraries have been installed.
|
||||
|
||||
Compilation
|
||||
===========
|
||||
1. Clone or unpack the code sources downloaded from
|
||||
[Bitbucket](git@bitbucket.org:amunteam/amun-code.git).
|
||||
1. Clone the AMUN source code: `git clone https://bitbucket.org/amunteam/amun-code.git`,
|
||||
or unpack the archive downloaded from page
|
||||
[Downloads](https://bitbucket.org/amunteam/amun-code/downloads/).
|
||||
2. Go to directory **hosts/** and copy file **default** to a new file named
|
||||
exactly as your host name (name returned by command `hostname`).
|
||||
3. Customize your compiler and compilation options in your new host file.
|
||||
|
@ -1,7 +1,7 @@
|
||||
# problem name and parameters
|
||||
#
|
||||
problem = "rt"
|
||||
gamma = 1.4
|
||||
problem = "rayleigh-taylor"
|
||||
gamma = 1.4d+00
|
||||
|
||||
# random number generator parameters
|
||||
#
|
||||
@ -24,33 +24,33 @@ fix_positivity = "off"
|
||||
# mesh parameters
|
||||
#
|
||||
xblocks = 1
|
||||
yblocks = 3
|
||||
yblocks = 2
|
||||
xmin = -2.5d-01
|
||||
xmax = 2.5d-01
|
||||
ymin = -7.5d-01
|
||||
ymax = 7.5d-01
|
||||
zmin = -5.0d-01
|
||||
zmax = 5.0d-01
|
||||
ymin = -5.0d-01
|
||||
ymax = 5.0d-01
|
||||
|
||||
# refinement control
|
||||
#
|
||||
ncells = 16
|
||||
nghosts = 4
|
||||
nghosts = 2
|
||||
minlev = 1
|
||||
maxlev = 6
|
||||
maxlev = 5
|
||||
crefmin = 1.0d-02
|
||||
crefmax = 1.0d-01
|
||||
|
||||
# boundary conditions
|
||||
#
|
||||
xlbndry = "periodic"
|
||||
xubndry = "periodic"
|
||||
ylbndry = "reflecting"
|
||||
yubndry = "reflecting"
|
||||
ylbndry = "hydrostatic"
|
||||
yubndry = "hydrostatic"
|
||||
zlbndry = "periodic"
|
||||
zubndry = "periodic"
|
||||
|
||||
# runtime control parameters
|
||||
#
|
||||
tmax = 1.0d+01
|
||||
tmax = 3.0d+00
|
||||
cfl = 3.0d-01
|
||||
|
||||
# data output control
|
||||
@ -59,5 +59,4 @@ precise_snapshots = "on"
|
||||
snapshot_type = "p"
|
||||
snapshot_interval = 1.0d-01
|
||||
restart_number = -1
|
||||
integrals_interval = 10
|
||||
generate_xdmf = "on"
|
||||
integrals_interval = 10
|
62
problems/sedov-taylor.in
Normal file
62
problems/sedov-taylor.in
Normal file
@ -0,0 +1,62 @@
|
||||
# problem name and parameters
|
||||
#
|
||||
problem = "sedov-taylor"
|
||||
gamma = 1.4d+00
|
||||
|
||||
# physics
|
||||
#
|
||||
equation_system = "hd"
|
||||
equation_of_state = "adi"
|
||||
|
||||
# methods
|
||||
#
|
||||
time_advance = "rk2"
|
||||
riemann_solver = "hll"
|
||||
reconstruction = "tvd"
|
||||
limiter = "mc"
|
||||
fix_positivity = "off"
|
||||
clip_extrema = "off"
|
||||
|
||||
# mesh parameters
|
||||
#
|
||||
xmin = -5.00d-01
|
||||
xmax = 5.00d-01
|
||||
ymin = -5.00d-01
|
||||
ymax = 5.00d-01
|
||||
zmin = -5.00d-01
|
||||
zmax = 5.00d-01
|
||||
|
||||
# refinement control
|
||||
#
|
||||
xblocks = 2
|
||||
yblocks = 2
|
||||
ncells = 16
|
||||
nghosts = 4
|
||||
minlev = 1
|
||||
maxlev = 4
|
||||
crefmin = 2.00d-01
|
||||
crefmax = 8.00d-01
|
||||
epsref = 1.00d-02
|
||||
refinement_variables = "dens pres"
|
||||
|
||||
# boundary conditions
|
||||
#
|
||||
xlbndry = "periodic"
|
||||
xubndry = "periodic"
|
||||
ylbndry = "periodic"
|
||||
yubndry = "periodic"
|
||||
zlbndry = "periodic"
|
||||
zubndry = "periodic"
|
||||
|
||||
# runtime control parameters
|
||||
#
|
||||
tmax = 2.00d-01
|
||||
cfl = 4.00d-01
|
||||
|
||||
# data output control
|
||||
#
|
||||
precise_snapshots = "on"
|
||||
snapshot_type = "p"
|
||||
snapshot_interval = 1.0d-02
|
||||
restart_number = -1
|
||||
integrals_interval = 10
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -130,15 +130,15 @@ module algebra
|
||||
!
|
||||
if (b(2) > 0.0d+00) then
|
||||
tm = - bh - dr
|
||||
x(1) = b(1) / tm
|
||||
x(2) = tm / b(3)
|
||||
else if (b(2) < 0.0d+00) then
|
||||
tm = - bh + dr
|
||||
x(1) = tm / b(3)
|
||||
x(2) = b(1) / tm
|
||||
else if (b(2) < 0.0d+00) then
|
||||
tm = - bh + dr
|
||||
x(1) = b(1) / tm
|
||||
x(2) = tm / b(3)
|
||||
else
|
||||
x(1) = dr / b(3)
|
||||
x(2) = - x(1)
|
||||
x(2) = dr / b(3)
|
||||
x(1) = - x(2)
|
||||
end if
|
||||
|
||||
! update the number of roots
|
||||
@ -303,15 +303,15 @@ module algebra
|
||||
!
|
||||
if (a(2) > 0.0d+00) then
|
||||
tm = bh + dr
|
||||
x(1) = - tm
|
||||
x(2) = - a(1) / tm
|
||||
else if (a(2) < 0.0d+00) then
|
||||
tm = bh - dr
|
||||
x(1) = - a(1) / tm
|
||||
x(2) = - tm
|
||||
else if (a(2) < 0.0d+00) then
|
||||
tm = bh - dr
|
||||
x(1) = - tm
|
||||
x(2) = - a(1) / tm
|
||||
else
|
||||
x(1) = - dr
|
||||
x(2) = dr
|
||||
x(1) = dr
|
||||
x(2) = - dr
|
||||
end if
|
||||
|
||||
! update the number of roots
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -122,7 +122,7 @@ program amun
|
||||
character(len=80) :: fmt, tmp
|
||||
|
||||
real(kind=8) :: tbeg, thrs
|
||||
real(kind=8) :: tm_curr, tm_exec, tm_conv
|
||||
real(kind=8) :: tm_curr, tm_exec, tm_conv, tm_last = 0.0d+00
|
||||
|
||||
#ifdef INTEL
|
||||
! the type of the function SIGNAL should be defined for Intel compiler
|
||||
@ -203,7 +203,7 @@ program amun
|
||||
write (*,"(1x,78('-'))")
|
||||
write (*,"(1x,18('='),17x,a,17x,19('='))") 'A M U N'
|
||||
write (*,"(1x,16('='),4x,a,4x,16('='))") &
|
||||
'Copyright (C) 2008-2017 Grzegorz Kowal'
|
||||
'Copyright (C) 2008-2018 Grzegorz Kowal'
|
||||
write (*,"(1x,18('='),9x,a,9x,19('='))") &
|
||||
'under GNU GPLv3 license'
|
||||
write (*,"(1x,78('-'))")
|
||||
@ -584,15 +584,15 @@ program amun
|
||||
write(*,"(1x,a)" ) "Evolving the system:"
|
||||
write(*,'(4x,a4,5x,a4,11x,a2,12x,a6,7x,a3)') 'step', 'time', 'dt' &
|
||||
, 'blocks', 'ETA'
|
||||
#if defined INTEL || defined PATHSCALE
|
||||
#ifdef INTEL
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1,$)') &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#else /* INTEL | PATHSCALE */
|
||||
#else /* INTEL */
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1)',advance="no") &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#endif /* INTEL | PATHSCALE */
|
||||
#endif /* INTEL */
|
||||
|
||||
end if
|
||||
|
||||
@ -642,30 +642,36 @@ program amun
|
||||
!
|
||||
if (thrs > trun) iterm = 100
|
||||
|
||||
! print progress info to console
|
||||
! print progress info to console, but not too often
|
||||
!
|
||||
if (master) then
|
||||
if (time >= tmax .or. (tm_curr - tm_last) >= 1.0d+00) then
|
||||
|
||||
! calculate days, hours, seconds
|
||||
!
|
||||
ec = int(tm_curr * (tmax - time) / max(1.0d-08, time - tbeg), kind = 4)
|
||||
es = max(0, int(mod(ec, 60)))
|
||||
em = int(mod(ec / 60, 60))
|
||||
eh = int(ec / 3600)
|
||||
ed = int(eh / 24)
|
||||
eh = int(mod(eh, 24))
|
||||
ed = min(9999,ed)
|
||||
ec = int(tm_curr * (tmax - time) / max(1.0d-08, time - tbeg), kind = 4)
|
||||
es = max(0, int(mod(ec, 60)))
|
||||
em = int(mod(ec / 60, 60))
|
||||
eh = int(ec / 3600)
|
||||
ed = int(eh / 24)
|
||||
eh = int(mod(eh, 24))
|
||||
ed = min(9999,ed)
|
||||
|
||||
#if defined INTEL || defined PATHSCALE
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1,$)') &
|
||||
#ifdef INTEL
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1,$)') &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#else /* INTEL | PATHSCALE */
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1)',advance="no") &
|
||||
#else /* INTEL */
|
||||
write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1)',advance="no") &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#endif /* INTEL | PATHSCALE */
|
||||
#endif /* INTEL */
|
||||
|
||||
! update the timestamp
|
||||
!
|
||||
tm_last = tm_curr
|
||||
|
||||
end if
|
||||
end if
|
||||
|
||||
! prepare iterm
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -289,7 +289,7 @@ module gravity
|
||||
! gravitational acceleration constant
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: gacc_const = -1.0d-01
|
||||
real(kind=8), save :: gacc_const = -1.0d+00
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -79,6 +79,10 @@ module interpolations
|
||||
!
|
||||
real(kind=8), save :: sgp = 9.0d+00
|
||||
|
||||
! PPM constant
|
||||
!
|
||||
real(kind=8), save :: ppm_const = 1.25d+00
|
||||
|
||||
! Gaussian process reconstruction coefficients vector
|
||||
!
|
||||
real(kind=8), dimension(:) , allocatable, save :: cgp
|
||||
@ -179,6 +183,7 @@ module interpolations
|
||||
call get_parameter_real ("limo3_rad" , rad )
|
||||
call get_parameter_real ("kappa" , kappa )
|
||||
call get_parameter_real ("kbeta" , kbeta )
|
||||
call get_parameter_real ("ppm_const" , ppm_const )
|
||||
call get_parameter_real ("cfl" , cfl )
|
||||
|
||||
! calculate κ = (1 - ν) / ν
|
||||
@ -215,6 +220,13 @@ module interpolations
|
||||
call print_warning("interpolations:initialize_interpolation" &
|
||||
, "Increase the number of ghost cells (at least 2).")
|
||||
eps = max(1.0d-12, eps)
|
||||
case ("ppm", "PPM")
|
||||
name_rec = "3rd order PPM"
|
||||
interfaces => interfaces_dir
|
||||
reconstruct_states => reconstruct_ppm
|
||||
if (verbose .and. ng < 4) &
|
||||
call print_warning("interpolations:initialize_interpolation" &
|
||||
, "Increase the number of ghost cells (at least 4).")
|
||||
case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z")
|
||||
name_rec = "5th order WENO-Z (Borges et al. 2008)"
|
||||
interfaces => interfaces_dir
|
||||
@ -1648,6 +1660,200 @@ module interpolations
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_PPM:
|
||||
! --------------------------
|
||||
!
|
||||
! Subroutine reconstructs the interface states using the third order
|
||||
! Piecewise Parabolic Method (PPM) with new limiters. This version lacks
|
||||
! the support for flattening important for keeping the spurious oscillations
|
||||
! near strong shocks/discontinuties under control.
|
||||
!
|
||||
! Arguments are described in subroutine reconstruct().
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Colella, P., Woodward, P. R.,
|
||||
! "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
|
||||
! Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
|
||||
! https://dx.doi.org/10.1016/0021-9991(84)90143-8
|
||||
! [2] Colella, P., Sekora, M. D.,
|
||||
! "A limiter for PPM that preserves accuracy at smooth extrema",
|
||||
! Journal of Computational Physics, 2008, vol. 227, pp. 7069-7076,
|
||||
! https://doi.org/10.1016/j.jcp.2008.03.034
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine reconstruct_ppm(n, h, f, fl, fr)
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
integer , intent(in) :: n
|
||||
real(kind=8) , intent(in) :: h
|
||||
real(kind=8), dimension(n), intent(in) :: f
|
||||
real(kind=8), dimension(n), intent(out) :: fl, fr
|
||||
|
||||
! local variables
|
||||
!
|
||||
logical :: cm, cp, ext
|
||||
integer :: i, im1, ip1, im2
|
||||
real(kind=8) :: df2c, df2m, df2p, df2s, df2l
|
||||
real(kind=8) :: dfm, dfp, dcm, dcp
|
||||
real(kind=8) :: alm, alp, amx, sgn, del
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(n) :: fi, df2
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! calculate the high-order interface interpolation; eq. (16) in [2]
|
||||
!
|
||||
fi(1 ) = 5.0d-01 * (f(1 ) + f(2 ))
|
||||
do i = 2, n - 2
|
||||
fi(i) = (7.0d+00 * (f(i) + f(i+1)) - (f(i-1) + f(i+2))) / 1.2d+01
|
||||
end do
|
||||
fi(n-1) = 5.0d-01 * (f(n-1) + f(n ))
|
||||
fi(n ) = f(n)
|
||||
|
||||
! calculate second-order central derivative
|
||||
!
|
||||
df2(1) = 0.0d+00
|
||||
do i = 2, n - 1
|
||||
df2(i) = (f(i+1) + f(i-1)) - 2.0d+00 * f(i)
|
||||
end do
|
||||
df2(n) = 0.0d+00
|
||||
|
||||
! limit the interpolation; eqs. (18) and (19) in [2]
|
||||
!
|
||||
do i = 2, n - 2
|
||||
im1 = i - 1
|
||||
ip1 = i + 1
|
||||
|
||||
if ((f(ip1) - fi(i)) * (fi(i) - f(i)) < 0.0d+00) then
|
||||
df2c = 3.0d+00 * ((f(ip1) + f(i )) - 2.0d+00 * fi(i))
|
||||
df2m = df2(i )
|
||||
df2p = df2(ip1)
|
||||
if (min(df2c, df2m, df2p) * max(df2c, df2m, df2p) > 0.0d+00) then
|
||||
df2l = sign(min(ppm_const * min(abs(df2m), abs(df2p)), abs(df2c)), df2c)
|
||||
else
|
||||
df2l = 0.0d+00
|
||||
end if
|
||||
fi(i) = 5.0d-01 * (f(i) + f(ip1)) - df2l / 6.0d+00
|
||||
end if
|
||||
end do
|
||||
|
||||
! iterate along the vector
|
||||
!
|
||||
do i = 2, n - 1
|
||||
im1 = i - 1
|
||||
ip1 = i + 1
|
||||
|
||||
! limit states if local extremum is detected or the interpolation is not
|
||||
! monotone
|
||||
!
|
||||
alm = fi(im1) - f(i)
|
||||
alp = fi(i ) - f(i)
|
||||
cm = abs(alm) >= 2.0d+00 * abs(alp)
|
||||
cp = abs(alp) >= 2.0d+00 * abs(alm)
|
||||
ext = .false.
|
||||
|
||||
! check if we have local extremum here
|
||||
!
|
||||
if (alm * alp > 0.0d+00) then
|
||||
ext = .true.
|
||||
else if (cm .or. cp) then
|
||||
im2 = max(1, i - 1)
|
||||
|
||||
dfm = fi(im1) - fi(im2)
|
||||
dfp = fi(ip1) - fi(i )
|
||||
dcm = f (i ) - f (im1)
|
||||
dcp = f (ip1) - f (i )
|
||||
if (min(abs(dfm),abs(dfp)) >= min(abs(dcm),abs(dcp))) then
|
||||
ext = dfm * dfp < 0.0d+00
|
||||
else
|
||||
ext = dcm * dcp < 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
! limit states if the local extremum is detected
|
||||
!
|
||||
if (ext) then
|
||||
df2s = 6.0d+00 * (alm + alp)
|
||||
df2m = df2(im1)
|
||||
df2c = df2(i )
|
||||
df2p = df2(ip1)
|
||||
if (min(df2s, df2c, df2m, df2p) &
|
||||
* max(df2s, df2c, df2m, df2p) > 0.0d+00) then
|
||||
df2l = sign(min(ppm_const * min(abs(df2m), abs(df2c), abs(df2p)) &
|
||||
, abs(df2s)), df2s)
|
||||
if (abs(df2l) > 0.0d+00) then
|
||||
alm = alm * df2l / df2s
|
||||
alp = alp * df2l / df2s
|
||||
else
|
||||
alm = 0.0d+00
|
||||
alp = 0.0d+00
|
||||
end if
|
||||
else
|
||||
alm = 0.0d+00
|
||||
alp = 0.0d+00
|
||||
end if
|
||||
else
|
||||
|
||||
! the interpolation is not monotonic so apply additional limiter
|
||||
!
|
||||
if (cp) then
|
||||
sgn = sign(1.0d+00, alp)
|
||||
amx = - 2.5d-01 * alp**2 / (alp + alm)
|
||||
dfp = f(ip1) - f(i)
|
||||
if (sgn * amx >= sgn * dfp) then
|
||||
del = dfp * (dfp - alm)
|
||||
if (del >= 0.0d+00) then
|
||||
alp = - 2.0d+00 * (dfp + sgn * sqrt(del))
|
||||
else
|
||||
alp = - 2.0d+00 * alm
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
|
||||
if (cm) then
|
||||
sgn = sign(1.0d+00, alm)
|
||||
amx = - 2.5d-01 * alm**2 / (alp + alm)
|
||||
dfm = f(im1) - f(i)
|
||||
if (sgn * amx >= sgn * dfm) then
|
||||
del = dfm * (dfm - alp)
|
||||
if (del >= 0.0d+00) then
|
||||
alm = - 2.0d+00 * (dfm + sgn * sqrt(del))
|
||||
else
|
||||
alm = - 2.0d+00 * alp
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
|
||||
! update the states
|
||||
!
|
||||
fr(im1) = f(i) + alm
|
||||
fl(i ) = f(i) + alp
|
||||
|
||||
end do ! i = 2, n
|
||||
|
||||
! update the interpolation of the first and last points
|
||||
!
|
||||
fl(1 ) = fi(1)
|
||||
fr(n-1) = fi(n)
|
||||
fl(n ) = fi(n)
|
||||
fr(n ) = f (n)
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine reconstruct_ppm
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine RECONSTRUCT_WENO5Z:
|
||||
! -----------------------------
|
||||
!
|
||||
|
56
src/io.F90
56
src/io.F90
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -1269,12 +1269,12 @@ module io
|
||||
use coordinates , only : minlev, maxlev, toplev
|
||||
use coordinates , only : nc, ng, in, jn, kn, ir, jr, kr
|
||||
use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax
|
||||
use equations , only : eqsys, eos
|
||||
use equations , only : eqsys, eos, gamma, csnd
|
||||
use error , only : print_error
|
||||
use evolution , only : step, time, dt, dtn
|
||||
use hdf5 , only : hid_t
|
||||
use hdf5 , only : h5gcreate_f, h5gclose_f
|
||||
use mpitools , only : nprocs, nproc
|
||||
use mpitools , only : nprocs, nproc, periodic
|
||||
use random , only : nseeds, get_seeds
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -1290,12 +1290,21 @@ module io
|
||||
integer(hid_t) :: gid
|
||||
integer :: err
|
||||
|
||||
! local vectors
|
||||
!
|
||||
integer, dimension(3) :: per
|
||||
|
||||
! local allocatable arrays
|
||||
!
|
||||
integer(kind=4), dimension(:), allocatable :: seeds
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! store the code name in order to determine the format of data
|
||||
!
|
||||
call write_attribute(fid, 'code' , 'AMUN')
|
||||
call write_attribute(fid, 'version', 'v1.0')
|
||||
|
||||
! create a group to store the global attributes
|
||||
!
|
||||
call h5gcreate_f(fid, 'attributes', gid, err)
|
||||
@ -1314,6 +1323,10 @@ module io
|
||||
|
||||
end if
|
||||
|
||||
! convert periodic(:) to an integer vector
|
||||
!
|
||||
per(:) = merge(1, 0, periodic(:))
|
||||
|
||||
! store string attributes
|
||||
!
|
||||
call write_attribute(gid, 'eqsys' , eqsys )
|
||||
@ -1321,21 +1334,22 @@ module io
|
||||
|
||||
! store the integer attributes
|
||||
!
|
||||
call write_attribute(gid, 'ndims' , NDIMS )
|
||||
call write_attribute(gid, 'last_id', get_last_id())
|
||||
call write_attribute(gid, 'mblocks', get_mblocks())
|
||||
call write_attribute(gid, 'dblocks', get_dblocks())
|
||||
call write_attribute(gid, 'nleafs' , get_nleafs() )
|
||||
call write_attribute(gid, 'ncells' , nc )
|
||||
call write_attribute(gid, 'nghosts', ng )
|
||||
call write_attribute(gid, 'minlev' , minlev )
|
||||
call write_attribute(gid, 'maxlev' , maxlev )
|
||||
call write_attribute(gid, 'toplev' , toplev )
|
||||
call write_attribute(gid, 'nprocs' , nprocs )
|
||||
call write_attribute(gid, 'nproc' , nproc )
|
||||
call write_attribute(gid, 'nseeds' , nseeds )
|
||||
call write_attribute(gid, 'step' , step )
|
||||
call write_attribute(gid, 'isnap' , isnap )
|
||||
call write_attribute(gid, 'ndims' , NDIMS )
|
||||
call write_attribute(gid, 'last_id' , get_last_id())
|
||||
call write_attribute(gid, 'mblocks' , get_mblocks())
|
||||
call write_attribute(gid, 'dblocks' , get_dblocks())
|
||||
call write_attribute(gid, 'nleafs' , get_nleafs() )
|
||||
call write_attribute(gid, 'ncells' , nc )
|
||||
call write_attribute(gid, 'nghosts' , ng )
|
||||
call write_attribute(gid, 'minlev' , minlev )
|
||||
call write_attribute(gid, 'maxlev' , maxlev )
|
||||
call write_attribute(gid, 'toplev' , toplev )
|
||||
call write_attribute(gid, 'nprocs' , nprocs )
|
||||
call write_attribute(gid, 'nproc' , nproc )
|
||||
call write_attribute(gid, 'nseeds' , nseeds )
|
||||
call write_attribute(gid, 'step' , step )
|
||||
call write_attribute(gid, 'isnap' , isnap )
|
||||
call write_attribute(gid, 'periodic', per(:) )
|
||||
|
||||
! store the real attributes
|
||||
!
|
||||
@ -1348,6 +1362,12 @@ module io
|
||||
call write_attribute(gid, 'time', time)
|
||||
call write_attribute(gid, 'dt' , dt )
|
||||
call write_attribute(gid, 'dtn' , dtn )
|
||||
if (eos == 'adi') then
|
||||
call write_attribute(gid, 'gamma', gamma)
|
||||
end if
|
||||
if (eos == 'iso') then
|
||||
call write_attribute(gid, 'csnd' , csnd )
|
||||
end if
|
||||
|
||||
! store the vector attributes
|
||||
!
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
519
src/problems.F90
519
src/problems.F90
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
@ -135,6 +135,9 @@ module problems
|
||||
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
|
||||
|
||||
@ -784,6 +787,470 @@ module problems
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! 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 : pi, pi4, d2r
|
||||
use coordinates, only : im, jm, km
|
||||
use coordinates, only : ax, ay, az, adx, ady, adz, advol
|
||||
use equations , only : prim2cons
|
||||
use equations , only : gamma
|
||||
use equations , only : nv
|
||||
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
use error , only : print_error
|
||||
use parameters , only : get_parameter_real, get_parameter_integer
|
||||
|
||||
! 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, ic, jc, kc
|
||||
real(kind=8) :: xl, yl, zl, xu, yu, zu, rl, ru
|
||||
real(kind=8) :: sn
|
||||
#if NDIMS == 3
|
||||
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, dz, dxh, dyh, dzh, dvol
|
||||
real(kind=8) :: fc_amb, fc_ovr
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(nv,im) :: q, u
|
||||
real(kind=8), dimension(im) :: x
|
||||
real(kind=8), dimension(jm) :: y
|
||||
real(kind=8), dimension(km) :: z
|
||||
|
||||
#if NDIMS == 3
|
||||
! allocatable arrays
|
||||
!
|
||||
real(kind=8), dimension(:), allocatable :: xm, ym, zm
|
||||
real(kind=8), dimension(:), allocatable :: xp, yp, zp
|
||||
#endif /* NDIMS == 3 */
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! quit if no adiabatic equation of state
|
||||
!
|
||||
if (ipr <= 0) then
|
||||
call print_error("problems::setup_problem_sedov_taylor" &
|
||||
, "Only adiabatic equation of state is supported!")
|
||||
stop
|
||||
end if
|
||||
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the problem setup
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare problem constants during the first subroutine call
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! get problem parameters
|
||||
!
|
||||
call get_parameter_real("radius", radius)
|
||||
call get_parameter_real("dens" , dens )
|
||||
call get_parameter_real("pres" , pres )
|
||||
call get_parameter_real("eexp" , eexp )
|
||||
call get_parameter_real("buni" , buni )
|
||||
call get_parameter_real("angle" , angle )
|
||||
|
||||
#if NDIMS == 3
|
||||
! get the fine grid resolution
|
||||
!
|
||||
call get_parameter_integer("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 = (gamma - 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(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
|
||||
y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
|
||||
#if NDIMS == 3
|
||||
z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km)
|
||||
#else /* NDIMS == 3 */
|
||||
z(1:km) = 0.0d+00
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! calculate mesh intervals and areas
|
||||
!
|
||||
dx = adx(pdata%meta%level)
|
||||
dy = ady(pdata%meta%level)
|
||||
dz = adz(pdata%meta%level)
|
||||
dxh = 0.5d+00 * dx
|
||||
dyh = 0.5d+00 * dy
|
||||
#if NDIMS == 3
|
||||
dzh = 0.5d+00 * dz
|
||||
#else /* NDIMS == 3 */
|
||||
dzh = 1.0d+00
|
||||
#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
|
||||
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
|
||||
!
|
||||
do k = 1, km
|
||||
|
||||
#if NDIMS == 3
|
||||
! calculate the corner Z coordinates
|
||||
!
|
||||
zl = abs(z(k)) - dzh
|
||||
zu = abs(z(k)) + dzh
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
do j = 1, jm
|
||||
|
||||
! calculate the corner Y coordinates
|
||||
!
|
||||
yl = abs(y(j)) - dyh
|
||||
yu = abs(y(j)) + dyh
|
||||
|
||||
! sweep along the X coordinate
|
||||
!
|
||||
do i = 1, im
|
||||
|
||||
! 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
|
||||
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
|
||||
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
|
||||
q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb
|
||||
|
||||
end if
|
||||
|
||||
end do ! i = 1, im
|
||||
|
||||
! convert the primitive variables to conservative ones
|
||||
!
|
||||
call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im))
|
||||
|
||||
! copy the conserved variables to the current block
|
||||
!
|
||||
pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im)
|
||||
|
||||
! copy the primitive variables to the current block
|
||||
!
|
||||
pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im)
|
||||
|
||||
end do ! j = 1, jm
|
||||
end do ! k = 1, km
|
||||
|
||||
#if NDIMS == 3
|
||||
! deallocate subgrid coordinates
|
||||
!
|
||||
deallocate(xm, ym, zm)
|
||||
deallocate(xp, yp, zp)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the problems setup
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine setup_problem_sedov_taylor
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine SETUP_PROBLEM_IMPLOSION:
|
||||
! ----------------------------------
|
||||
!
|
||||
@ -1213,6 +1680,14 @@ module problems
|
||||
! 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)
|
||||
@ -1220,9 +1695,10 @@ module problems
|
||||
! include external procedures and variables
|
||||
!
|
||||
use blocks , only : block_data
|
||||
use constants , only : d2r
|
||||
use constants , only : pi2, d2r
|
||||
use coordinates, only : xmin, xmax, xlen
|
||||
use coordinates, only : im, jm, km
|
||||
use coordinates, only : ay, ady
|
||||
use coordinates, only : ax, ay, ady
|
||||
use equations , only : prim2cons
|
||||
use equations , only : nv
|
||||
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
@ -1240,12 +1716,16 @@ module problems
|
||||
|
||||
! default parameter values
|
||||
!
|
||||
real(kind=8), save :: ycut = 0.00d-00
|
||||
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 :: vper = 1.00d-02
|
||||
real(kind=8), save :: gacc = -1.00d-01
|
||||
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
|
||||
@ -1262,9 +1742,8 @@ module problems
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(nv,im) :: q, u
|
||||
real(kind=8), dimension(im) :: x
|
||||
real(kind=8), dimension(im) :: x, yp
|
||||
real(kind=8), dimension(jm) :: y
|
||||
real(kind=8), dimension(km) :: z
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1284,12 +1763,19 @@ module problems
|
||||
call get_parameter_real("dens" , dens )
|
||||
call get_parameter_real("drat" , drat )
|
||||
call get_parameter_real("pres" , pres )
|
||||
call get_parameter_real("lper" , lper )
|
||||
call get_parameter_real("kper" , kper )
|
||||
call get_parameter_real("vper" , vper )
|
||||
call get_parameter_real("hdel" , hdel )
|
||||
call get_parameter_real("gacc" , gacc )
|
||||
call get_parameter_real("buni" , buni )
|
||||
call get_parameter_real("bgui" , bgui )
|
||||
call get_parameter_real("angle" , angle )
|
||||
|
||||
! calculate the density change across the interface
|
||||
!
|
||||
damp = 5.0d-01 * (drat * dens - dens)
|
||||
|
||||
! reset the first execution flag
|
||||
!
|
||||
first = .false.
|
||||
@ -1298,6 +1784,7 @@ module problems
|
||||
|
||||
! prepare block coordinates
|
||||
!
|
||||
x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
|
||||
y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
|
||||
|
||||
! set the ambient density and pressure
|
||||
@ -1324,6 +1811,11 @@ module problems
|
||||
|
||||
end if
|
||||
|
||||
! prepare density perturbation
|
||||
!
|
||||
yp(1:im) = 0.5d+00 * lper * (cos(pi2 * kper * (x(1:im) - xmin) / xlen) &
|
||||
+ cos(pi2 * kper * (xmax - x(1:im)) / xlen)) + ycut
|
||||
|
||||
! iterate over all positions in the YZ plane
|
||||
!
|
||||
do k = 1, km
|
||||
@ -1337,6 +1829,7 @@ module problems
|
||||
else
|
||||
q(idn,1:im) = dens * drat
|
||||
end if
|
||||
q(idn,1:im) = dens + damp * (1.0d+00 + tanh((y(j) - yp(1:im)) / hdel))
|
||||
q(ipr,1:im) = pres + q(idn,1:im) * gacc * y(j)
|
||||
else
|
||||
if (y(j) <= ycut) then
|
||||
@ -1354,9 +1847,11 @@ module problems
|
||||
|
||||
! add a random seed velocity component
|
||||
!
|
||||
do i = 1, im
|
||||
q(ivy,i) = q(ivy,i) + vper * randomn()
|
||||
end do
|
||||
if (vper /= 0.0d+00) then
|
||||
do i = 1, im
|
||||
q(ivy,i) = q(ivy,i) + vper * randomn()
|
||||
end do
|
||||
end if
|
||||
|
||||
! convert the primitive variables to conservative ones
|
||||
!
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2008-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2017-2018 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
|
||||
|
@ -4,7 +4,7 @@
|
||||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||||
!! adaptive mesh.
|
||||
!!
|
||||
!! Copyright (C) 2015-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||||
!! Copyright (C) 2015-2018 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
|
||||
|
Loading…
x
Reference in New Issue
Block a user