gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2014-07-24 12:43:59 -03:00
commit 6011ec0f44
13 changed files with 15172 additions and 6792 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
problems/blast.in
View File

@ -1,13 +1,13 @@
# problem name and parameters
#
problem = "blast"
dens = 1.00d+00
ratio = 1.00d+02
radius = 1.00d-01
csnd = 4.08d-01
gamma = 1.67d-01
buni = 0.00d+00
angle = 0.00d+00
dens = 1.0000d+00
ratio = 1.0000d+02
radius = 1.0000d-01
buni = 0.0000d+00
angle = 4.5000d+01
csnd = 4.0825d-01
gamma = 1.6667d+00
# physics
#
@ -28,19 +28,19 @@ clip_extrema = "off"
rdimx = 2
rdimy = 3
rdimz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.75
ymax = 0.75
zmin = -0.5
zmax = 0.5
xmin = -5.00d-01
xmax = 5.00d-01
ymin = -7.50d-01
ymax = 7.50d-01
zmin = -5.00d-01
zmax = 5.00d-01
# refinement control
#
ncells = 8
nghost = 2
minlev = 1
maxlev = 7
maxlev = 6
crefmin = 2.00d-01
crefmax = 8.00d-01
epsref = 1.00d-02

20
problems/implosion.in
View File

@ -23,31 +23,31 @@ ncells = 8
nghosts = 2
rdimx = 1
rdimy = 1
xmin = -1.00d+00
xmin = 0.00d+00
xmax = 1.00d+00
ymin = -1.00d+00
ymin = 0.00d+00
ymax = 1.00d+00
# refinement control
#
maxlev = 8
maxlev = 7
# boundary conditions
#
#xlbndry = "reflecting"
#xubndry = "reflecting"
#ylbndry = "reflecting"
#yubndry = "reflecting"
xlbndry = "reflecting"
xubndry = "reflecting"
ylbndry = "reflecting"
yubndry = "reflecting"
# runtime control parameters
#
tmax = 1.00d+00
tmax = 4.00d+00
cfl = 4.00d-01
# data output control
#
precise_snapshots = "on"
snapshot_type = "p"
snapshot_interval = 1.0d-02
snapshot_interval = 1.0d-01
restart_number = -1
integrals_interval = 1
integrals_interval = 10

6
problems/kh.in
View File

@ -14,16 +14,13 @@ equation_of_state = "adi"
# methods
#
time_advance = "rk2"
riemann_solver = "hllc"
riemann_solver = "hlld"
reconstruction = "limo3"
limiter = "mc"
fix_positivity = "off"
# mesh parameters
#
rdimx = 1
rdimy = 1
rdimz = 1
xmin = -5.0d-01
xmax = 5.0d-01
ymin = -5.0d-01
@ -49,7 +46,6 @@ zubndry = "periodic"
# runtime control parameters
#
nmax = 1000000
tmax = 1.0d+00
cfl = 4.0d-01

3156
src/blocks.F90
View File

File diff suppressed because it is too large Load Diff

10296
src/boundaries.F90
View File

File diff suppressed because it is too large Load Diff

24
src/coordinates.F90
View File

@ -83,6 +83,11 @@ module coordinates
real(kind=8), save :: zmax = 1.0d+00
real(kind=8), save :: zlen = 1.0d+00
! the domain volume and its inversion
!
real(kind=8), save :: vol = 1.0d+00
real(kind=8), save :: voli = 1.0d+00
! the block coordinates for all levels of refinement
!
real(kind=8), dimension(:,:), allocatable, save :: ax , ay , az
@ -239,6 +244,19 @@ module coordinates
call get_parameter_real ("zmax" , zmax )
#endif /* NDIMS == 3 */
! calculate the domain sizes
!
xlen = xmax - xmin
ylen = ymax - ymin
#if NDIMS == 3
zlen = zmax - zmin
#endif /* NDIMS == 3 */
! calculate the domain volume
!
vol = xlen * ylen * zlen
voli = 1.0d+00 / vol
! allocate space for coordinate variables
!
allocate(ax (toplev, im))
@ -280,10 +298,10 @@ module coordinates
! calculate the cell sizes for each level
!
adx (l) = (xmax - xmin) / (ir * ni)
ady (l) = (ymax - ymin) / (jr * nj)
adx (l) = xlen / (ir * ni)
ady (l) = ylen / (jr * nj)
#if NDIMS == 3
adz (l) = (zmax - zmin) / (kr * nk)
adz (l) = zlen / (kr * nk)
#endif /* NDIMS == 3 */
#if NDIMS == 2
adr (l) = sqrt(adx(l)**2 + ady(l)**2)

322
src/domains.F90
View File

@ -142,9 +142,9 @@ module domains
use blocks , only : metablock_set_leaf, metablock_set_level
use blocks , only : metablock_set_configuration
use blocks , only : metablock_set_coordinates, metablock_set_bounds
use blocks , only : nsides, nfaces
use blocks , only : nsides
use boundaries , only : bnd_type, bnd_periodic
use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax
use coordinates , only : xmin, ymin, zmin, xlen, ylen, zlen
use coordinates , only : ir, jr, kr
! local variables are not implicit by default
@ -167,6 +167,8 @@ module domains
! allocatable arrays
!
integer, dimension(:,:,:), allocatable :: cfg
integer, dimension(:) , allocatable :: im, jm, km
integer, dimension(:) , allocatable :: ip, jp, kp
! local pointer array
!
@ -257,9 +259,9 @@ module domains
!!
! calculate block sizes
!
xl = (xmax - xmin) / ir
yl = (ymax - ymin) / jr
zl = (zmax - zmin) / kr
xl = xlen / ir
yl = ylen / jr
zl = zlen / kr
! fill out block structure fields
!
@ -321,144 +323,208 @@ module domains
!! ASSIGN THE BLOCK NEIGHBORS
!!
! assign boundaries along the X direction
! allocate indices
!
allocate(im(ir), ip(ir))
allocate(jm(jr), jp(jr))
#if NDIMS == 3
allocate(km(kr), kp(kr))
#endif /* NDIMS == 3 */
! generate indices
!
im(:) = cshift((/(i, i = 1, ir)/),-1)
ip(:) = cshift((/(i, i = 1, ir)/), 1)
jm(:) = cshift((/(j, j = 1, jr)/),-1)
jp(:) = cshift((/(j, j = 1, jr)/), 1)
#if NDIMS == 3
km(:) = cshift((/(k, k = 1, kr)/),-1)
kp(:) = cshift((/(k, k = 1, kr)/), 1)
#endif /* NDIMS == 3 */
! check periodicity and reset the edge indices if box is not periodic
!
if (bnd_type(1,1) /= bnd_periodic .or. bnd_type(1,2) /= bnd_periodic) then
im( 1) = 0
ip(ir) = 0
end if
if (bnd_type(2,1) /= bnd_periodic .or. bnd_type(2,2) /= bnd_periodic) then
jm( 1) = 0
jp(jr) = 0
end if
#if NDIMS == 3
if (bnd_type(3,1) /= bnd_periodic .or. bnd_type(3,2) /= bnd_periodic) then
km( 1) = 0
kp(kr) = 0
end if
#endif /* NDIMS == 3 */
! iterate over all initial blocks
!
do k = 1, kr
do j = 1, jr
do i = 1, ir - 1
! assign a pointer
!
pmeta => block_array(i ,j,k)%ptr
! assign neighbor
!
pnext => block_array(i+1,j,k)%ptr
! assign their neighbor pointers
!
do p = 1, nfaces
pmeta%neigh(1,2,p)%ptr => pnext
pnext%neigh(1,1,p)%ptr => pmeta
end do
end do
end do
end do
! if periodic boundary conditions set edge block neighbors
!
if (bnd_type(1,1) == bnd_periodic .and. bnd_type(1,2) == bnd_periodic) then
do k = 1, kr
do j = 1, jr
! assign pointers
!
pmeta => block_array( 1,j,k)%ptr
pnext => block_array(ir,j,k)%ptr
! assign their neighbor pointers
!
do p = 1, nfaces
pmeta%neigh(1,1,p)%ptr => pnext
pnext%neigh(1,2,p)%ptr => pmeta
end do
end do
end do
end if
! assign boundaries along the Y direction
!
do k = 1, kr
do j = 1, jr - 1
do i = 1, ir
! assign a pointer
! assign pmeta with the current block
!
pmeta => block_array(i,j ,k)%ptr
! assign neighbor
!
pnext => block_array(i,j+1,k)%ptr
! assign their neighbor pointers
!
do p = 1, nfaces
pmeta%neigh(2,2,p)%ptr => pnext
pnext%neigh(2,1,p)%ptr => pmeta
end do
end do
end do
end do
! if periodic boundary conditions set edge block neighbors
!
if (bnd_type(2,1) == bnd_periodic .and. bnd_type(2,2) == bnd_periodic) then
do k = 1, kr
do i = 1, ir
! assign pointers
!
pmeta => block_array(i, 1,k)%ptr
pnext => block_array(i,jr,k)%ptr
! assign their neighbor pointers
!
do p = 1, nfaces
pmeta%neigh(2,1,p)%ptr => pnext
pnext%neigh(2,2,p)%ptr => pmeta
end do
end do
end do
end if
pmeta => block_array(i,j,k)%ptr
#if NDIMS == 3
! assign boundaries along the Z direction
! assign face neighbor pointers
!
do k = 1, kr - 1
do j = 1, jr
do i = 1, ir
if (im(i) > 0) then
pmeta%faces(1,1,1,1)%ptr => block_array(im(i),j,k)%ptr
pmeta%faces(1,2,1,1)%ptr => block_array(im(i),j,k)%ptr
pmeta%faces(1,1,2,1)%ptr => block_array(im(i),j,k)%ptr
pmeta%faces(1,2,2,1)%ptr => block_array(im(i),j,k)%ptr
end if
if (ip(i) > 0) then
pmeta%faces(2,1,1,1)%ptr => block_array(ip(i),j,k)%ptr
pmeta%faces(2,2,1,1)%ptr => block_array(ip(i),j,k)%ptr
pmeta%faces(2,1,2,1)%ptr => block_array(ip(i),j,k)%ptr
pmeta%faces(2,2,2,1)%ptr => block_array(ip(i),j,k)%ptr
end if
! assign a pointer
if (jm(j) > 0) then
pmeta%faces(1,1,1,2)%ptr => block_array(i,jm(j),k)%ptr
pmeta%faces(2,1,1,2)%ptr => block_array(i,jm(j),k)%ptr
pmeta%faces(1,1,2,2)%ptr => block_array(i,jm(j),k)%ptr
pmeta%faces(2,1,2,2)%ptr => block_array(i,jm(j),k)%ptr
end if
if (jp(j) > 0) then
pmeta%faces(1,2,1,2)%ptr => block_array(i,jp(j),k)%ptr
pmeta%faces(2,2,1,2)%ptr => block_array(i,jp(j),k)%ptr
pmeta%faces(1,2,2,2)%ptr => block_array(i,jp(j),k)%ptr
pmeta%faces(2,2,2,2)%ptr => block_array(i,jp(j),k)%ptr
end if
if (km(k) > 0) then
pmeta%faces(1,1,1,3)%ptr => block_array(i,j,km(k))%ptr
pmeta%faces(2,1,1,3)%ptr => block_array(i,j,km(k))%ptr
pmeta%faces(1,2,1,3)%ptr => block_array(i,j,km(k))%ptr
pmeta%faces(2,2,1,3)%ptr => block_array(i,j,km(k))%ptr
end if
if (kp(k) > 0) then
pmeta%faces(1,1,2,3)%ptr => block_array(i,j,kp(k))%ptr
pmeta%faces(2,1,2,3)%ptr => block_array(i,j,kp(k))%ptr
pmeta%faces(1,2,2,3)%ptr => block_array(i,j,kp(k))%ptr
pmeta%faces(2,2,2,3)%ptr => block_array(i,j,kp(k))%ptr
end if
#endif /* NDIMS == 3 */
! assign edge neighbor pointers
!
pmeta => block_array(i,j,k )%ptr
#if NDIMS == 2
if (im(i) > 0) then
pmeta%edges(1,1,2)%ptr => block_array(im(i),j,k)%ptr
pmeta%edges(1,2,2)%ptr => block_array(im(i),j,k)%ptr
end if
if (ip(i) > 0) then
pmeta%edges(2,1,2)%ptr => block_array(ip(i),j,k)%ptr
pmeta%edges(2,2,2)%ptr => block_array(ip(i),j,k)%ptr
end if
if (jm(j) > 0) then
pmeta%edges(1,1,1)%ptr => block_array(i,jm(j),k)%ptr
pmeta%edges(2,1,1)%ptr => block_array(i,jm(j),k)%ptr
end if
if (jp(j) > 0) then
pmeta%edges(1,2,1)%ptr => block_array(i,jp(j),k)%ptr
pmeta%edges(2,2,1)%ptr => block_array(i,jp(j),k)%ptr
end if
#endif /* NDIMS == 2 */
#if NDIMS == 3
if (jm(j) > 0 .and. km(k) > 0) then
pmeta%edges(1,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr
pmeta%edges(2,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr
end if
if (jp(j) > 0 .and. km(k) > 0) then
pmeta%edges(1,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr
pmeta%edges(2,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr
end if
if (jm(j) > 0 .and. kp(k) > 0) then
pmeta%edges(1,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr
pmeta%edges(2,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr
end if
if (jp(j) > 0 .and. kp(k) > 0) then
pmeta%edges(1,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr
pmeta%edges(2,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr
end if
! assign neighbor
if (im(i) > 0 .and. km(k) > 0) then
pmeta%edges(1,1,1,2)%ptr => block_array(im(i),j,km(k))%ptr
pmeta%edges(1,2,1,2)%ptr => block_array(im(i),j,km(k))%ptr
end if
if (ip(i) > 0 .and. km(k) > 0) then
pmeta%edges(2,1,1,2)%ptr => block_array(ip(i),j,km(k))%ptr
pmeta%edges(2,2,1,2)%ptr => block_array(ip(i),j,km(k))%ptr
end if
if (im(i) > 0 .and. kp(k) > 0) then
pmeta%edges(1,1,2,2)%ptr => block_array(im(i),j,kp(k))%ptr
pmeta%edges(1,2,2,2)%ptr => block_array(im(i),j,kp(k))%ptr
end if
if (ip(i) > 0 .and. kp(k) > 0) then
pmeta%edges(2,1,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr
pmeta%edges(2,2,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr
end if
if (im(i) > 0 .and. jm(j) > 0) then
pmeta%edges(1,1,1,3)%ptr => block_array(im(i),jm(j),k)%ptr
pmeta%edges(1,1,2,3)%ptr => block_array(im(i),jm(j),k)%ptr
end if
if (ip(i) > 0 .and. jm(j) > 0) then
pmeta%edges(2,1,1,3)%ptr => block_array(ip(i),jm(j),k)%ptr
pmeta%edges(2,1,2,3)%ptr => block_array(ip(i),jm(j),k)%ptr
end if
if (im(i) > 0 .and. jp(j) > 0) then
pmeta%edges(1,2,1,3)%ptr => block_array(im(i),jp(j),k)%ptr
pmeta%edges(1,2,2,3)%ptr => block_array(im(i),jp(j),k)%ptr
end if
if (ip(i) > 0 .and. jp(j) > 0) then
pmeta%edges(2,2,1,3)%ptr => block_array(ip(i),jp(j),k)%ptr
pmeta%edges(2,2,2,3)%ptr => block_array(ip(i),jp(j),k)%ptr
end if
#endif /* NDIMS == 3 */
! assign corner neighbor pointers
!
pnext => block_array(i,j,k+1)%ptr
#if NDIMS == 2
if (im(i) > 0 .and. jm(j) > 0) &
pmeta%corners(1,1)%ptr => block_array(im(i),jm(j),k)%ptr
if (ip(i) > 0 .and. jm(j) > 0) &
pmeta%corners(2,1)%ptr => block_array(ip(i),jm(j),k)%ptr
if (im(i) > 0 .and. jp(j) > 0) &
pmeta%corners(1,2)%ptr => block_array(im(i),jp(j),k)%ptr
if (ip(i) > 0 .and. jp(j) > 0) &
pmeta%corners(2,2)%ptr => block_array(ip(i),jp(j),k)%ptr
#endif /* NDIMS == 2 */
#if NDIMS == 3
if (im(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) &
pmeta%corners(1,1,1)%ptr => block_array(im(i),jm(j),km(k))%ptr
if (ip(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) &
pmeta%corners(2,1,1)%ptr => block_array(ip(i),jm(j),km(k))%ptr
if (im(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) &
pmeta%corners(1,2,1)%ptr => block_array(im(i),jp(j),km(k))%ptr
if (ip(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) &
pmeta%corners(2,2,1)%ptr => block_array(ip(i),jp(j),km(k))%ptr
if (im(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) &
pmeta%corners(1,1,2)%ptr => block_array(im(i),jm(j),kp(k))%ptr
if (ip(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) &
pmeta%corners(2,1,2)%ptr => block_array(ip(i),jm(j),kp(k))%ptr
if (im(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) &
pmeta%corners(1,2,2)%ptr => block_array(im(i),jp(j),kp(k))%ptr
if (ip(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) &
pmeta%corners(2,2,2)%ptr => block_array(ip(i),jp(j),kp(k))%ptr
#endif /* NDIMS == 3 */
end do ! over i
end do ! over j
end do ! over k
! assign their neighbor pointers
! deallocate indices
!
do p = 1, nfaces
pmeta%neigh(3,2,p)%ptr => pnext
pnext%neigh(3,1,p)%ptr => pmeta
end do
end do
end do
end do
! if periodic boundary conditions set edge block neighbors
!
if (bnd_type(3,1) == bnd_periodic .and. bnd_type(3,2) == bnd_periodic) then
do j = 1, jr
do i = 1, ir
! assign pointers
!
pmeta => block_array(i,j, 1)%ptr
pnext => block_array(i,j,kr)%ptr
! assign their neighbor pointers
!
do p = 1, nfaces
pmeta%neigh(3,1,p)%ptr => pnext
pnext%neigh(3,2,p)%ptr => pmeta
end do
end do
end do
end if
deallocate(im, ip)
deallocate(jm, jp)
#if NDIMS == 3
deallocate(km, kp)
#endif /* NDIMS == 3 */
! deallocate the block pointer array

5
src/driver.F90
View File

@ -36,6 +36,7 @@ program amun
!
use blocks , only : initialize_blocks, finalize_blocks, get_nleafs
use boundaries , only : initialize_boundaries, finalize_boundaries
use boundaries , only : boundary_variables
use coordinates , only : initialize_coordinates, finalize_coordinates
use equations , only : initialize_equations, finalize_equations
use evolution , only : initialize_evolution, finalize_evolution
@ -442,6 +443,10 @@ program amun
!
call generate_mesh()
! update boundaries
!
call boundary_variables()
! calculate new timestep
!
call new_time_step(dtnext)

85
src/evolution.F90
View File

@ -257,6 +257,12 @@ module evolution
!
call boundary_variables()
#ifdef DEBUG
! check variables for NaNs
!
call check_variables()
#endif /* DEBUG */
! set all meta blocks to be updated
!
call set_blocks_update(.true.)
@ -818,6 +824,85 @@ module evolution
!-------------------------------------------------------------------------------
!
end subroutine update_variables
#ifdef DEBUG
!
!===============================================================================
!
! subroutine CHECK_VARIABLES:
! --------------------------
!
! Subroutine iterates over all data blocks and converts the conservative
! variables to their primitive representation.
!
!
!===============================================================================
!
subroutine check_variables()
! include external procedures
!
use coordinates , only : im, jm, km
use equations , only : nv, pvars, cvars
! include external variables
!
use blocks , only : block_meta, list_meta
use blocks , only : block_data, list_data
! local variables are not implicit by default
!
implicit none
! local variables
!
integer :: i, j, k, p
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! associate the pointer with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! associate pmeta with the corresponding meta block
!
pmeta => pdata%meta
! check if there are NaNs in primitive variables
!
do k = 1, km
do j = 1, jm
do i = 1, im
do p = 1, nv
if (isnan(pdata%u(p,i,j,k))) then
print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k
end if
if (isnan(pdata%q(p,i,j,k))) then
print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k
end if
end do ! p = 1, nv
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
! assign pointer to the next block
!
pdata => pdata%next
end do
!-------------------------------------------------------------------------------
!
end subroutine check_variables
#endif /* DEBUG */
!===============================================================================
!

217
src/integrals.F90
View File

@ -50,9 +50,11 @@ module integrals
! =================
!
! funit - a file handler to the integrals file;
! sunit - a file handler to the statistics file;
! iintd - the number of steps between subsequent intervals storing;
!
integer(kind=4), save :: funit = 7
integer(kind=4), save :: sunit = 8
integer(kind=4), save :: iintd = 1
! by default everything is private
@ -93,7 +95,6 @@ module integrals
! import external variables and subroutines
!
use error , only : print_error
use mpitools , only : master
use parameters , only : get_parameter_integer
@ -109,7 +110,6 @@ module integrals
! local variables
!
character(len=32) :: fname
logical :: lex, lop
!
!-------------------------------------------------------------------------------
!
@ -136,35 +136,18 @@ module integrals
!
if (master) then
! find the first available file handler to use
!
funit = 6
lex = .false.
lop = .true.
do while(.not. lex .or. lop .and. funit < 100)
funit = funit + 1
inquire(unit = funit, exist = lex, opened = lop)
end do
! check if the file handler could be found
!
if (funit >= 100) then
call print_error('integrals::initialize_integrals' &
, 'Could not find any available file handlers!')
iret = 300
end if
! generate the integrals file name
!
write(fname, "('integrals_',i2.2,'.dat')") irun
! create a new file
! create a new integrals file
!
#ifdef INTEL
open (unit = funit, file = fname, form = 'formatted', status = 'replace' &
, buffered = 'yes')
open (newunit = funit, file = fname, form = 'formatted' &
, status = 'replace', buffered = 'yes')
#else /* INTEL */
open (unit = funit, file = fname, form = 'formatted', status = 'replace')
open (newunit = funit, file = fname, form = 'formatted' &
, status = 'replace')
#endif /* INTEL */
! write the integral file header
@ -174,6 +157,32 @@ module integrals
, 'ener', 'ekin', 'emag', 'eint'
write(funit,"('#')")
! generate the statistics file name
!
write(fname, "('statistics_',i2.2,'.dat')") irun
! create a new statistics file
!
#ifdef INTEL
open (newunit = sunit, file = fname, form = 'formatted' &
, status = 'replace', buffered = 'yes')
#else /* INTEL */
open (newunit = sunit, file = fname, form = 'formatted' &
, status = 'replace')
#endif /* INTEL */
! write the integral file header
!
write(sunit,"('#',a8,23(1x,a18))") 'step', 'time' &
, 'mean(dens)', 'min(dens)', 'max(dens)' &
, 'mean(pres)', 'min(pres)', 'max(pres)' &
, 'mean(velo)', 'min(velo)', 'max(velo)' &
, 'mean(magn)', 'min(magn)', 'max(magn)' &
, 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)' &
, 'mean(Mson)', 'min(Mson)', 'max(Mson)' &
, 'mean(Malf)', 'min(Malf)', 'max(Malf)'
write(sunit,"('#')")
end if ! master
#ifdef PROFILE
@ -214,9 +223,12 @@ module integrals
call start_timer(imi)
#endif /* PROFILE */
! close the integrals file
! close the integrals and statistics files
!
if (master) close(funit)
if (master) then
close(funit)
close(sunit)
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
@ -244,13 +256,17 @@ module integrals
! import external variables and subroutines
!
use blocks , only : block_meta, block_data, list_data
use coordinates , only : ib, ie, jb, je, kb, ke
use coordinates , only : advol
use equations , only : idn, imx, imy, imz, ien, ibx, iby, ibz
use coordinates , only : in, jn, kn, ib, jb, kb, ie, je, ke
use coordinates , only : advol, voli
use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
use equations , only : ien, imx, imy, imz
use equations , only : gamma, csnd
use evolution , only : step, time, dt
use mpitools , only : master
#ifdef MPI
use mpitools , only : reduce_sum_real_array
use mpitools , only : reduce_minimum_real_array
use mpitools , only : reduce_maximum_real_array
#endif /* MPI */
! local variables are not implicit by default
@ -272,7 +288,13 @@ module integrals
! local arrays
!
real(kind=8), dimension(narr) :: arr
real(kind=8), dimension(narr) :: inarr, avarr, mnarr, mxarr
real(kind=8), dimension(in,jn,kn) :: vel, mag, sqd, tmp
! parameters
!
real(kind=8), parameter :: eps = epsilon(1.0d+00)
real(kind=8), parameter :: big = huge(1.0d+00)
!
!-------------------------------------------------------------------------------
!
@ -288,7 +310,25 @@ module integrals
! reset the integrals array
!
arr(:) = 0.0d+00
inarr(:) = 0.0d+00
avarr(:) = 0.0d+00
mnarr(:) = big
mxarr(:) = - big
! reset some statistics if they are not used
!
if (ipr < 1) then
mnarr(2) = 0.0d+00
mxarr(2) = 0.0d+00
end if
if (ibx < 1) then
mnarr(4) = 0.0d+00
mxarr(4) = 0.0d+00
mnarr(5) = 0.0d+00
mxarr(5) = 0.0d+00
mnarr(7) = 0.0d+00
mxarr(7) = 0.0d+00
end if
! associate the pointer with the first block on the data block list
!
@ -305,30 +345,93 @@ module integrals
! sum up density and momenta components
!
arr(1) = arr(1) + sum(pdata%u(idn,ib:ie,jb:je,kb:ke)) * dvol
arr(2) = arr(2) + sum(pdata%u(imx,ib:ie,jb:je,kb:ke)) * dvol
arr(3) = arr(3) + sum(pdata%u(imy,ib:ie,jb:je,kb:ke)) * dvol
arr(4) = arr(4) + sum(pdata%u(imz,ib:ie,jb:je,kb:ke)) * dvol
inarr(1) = inarr(1) + sum(pdata%u(idn,ib:ie,jb:je,kb:ke)) * dvol
inarr(2) = inarr(2) + sum(pdata%u(imx,ib:ie,jb:je,kb:ke)) * dvol
inarr(3) = inarr(3) + sum(pdata%u(imy,ib:ie,jb:je,kb:ke)) * dvol
inarr(4) = inarr(4) + sum(pdata%u(imz,ib:ie,jb:je,kb:ke)) * dvol
! sum up total energy
!
if (ien > 0) then
arr(5) = arr(5) + sum(pdata%u(ien,ib:ie,jb:je,kb:ke)) * dvol
inarr(5) = inarr(5) + sum(pdata%u(ien,ib:ie,jb:je,kb:ke)) * dvol
end if
! sum up kinetic energy
!
arr(6) = arr(6) + sum((pdata%u(imx,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(imy,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(imz,ib:ie,jb:je,kb:ke)**2) &
/ pdata%u(idn,ib:ie,jb:je,kb:ke)) * dvolh
inarr(6) = inarr(6) + sum((pdata%u(imx,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(imy,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(imz,ib:ie,jb:je,kb:ke)**2) &
/ pdata%u(idn,ib:ie,jb:je,kb:ke)) * dvolh
! sum up magnetic energy
!
if (ibx > 0) then
arr(7) = arr(7) + sum(pdata%u(ibx,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(iby,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(ibz,ib:ie,jb:je,kb:ke)**2) * dvolh
inarr(7) = inarr(7) + sum(pdata%u(ibx,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(iby,ib:ie,jb:je,kb:ke)**2 &
+ pdata%u(ibz,ib:ie,jb:je,kb:ke)**2) * dvolh
end if
! get average, minimum and maximum values of density
!
tmp(:,:,:) = pdata%q(idn,ib:ie,jb:je,kb:ke)
avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol
mnarr(1) = min(mnarr(1), minval(tmp(:,:,:)))
mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:)))
! get average, minimum and maximum values of pressure
!
if (ipr > 0) then
tmp(:,:,:) = pdata%q(ipr,ib:ie,jb:je,kb:ke)
avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol
mnarr(2) = min(mnarr(2), minval(tmp(:,:,:)))
mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:)))
end if
! get average, minimum and maximum values of velocity amplitude
!
vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,ib:ie,jb:je,kb:ke)**2, 1))
avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol
mnarr(3) = min(mnarr(3), minval(vel(:,:,:)))
mxarr(3) = max(mxarr(3), maxval(vel(:,:,:)))
! get average, minimum and maximum values of magnetic field amplitude, and
! divergence potential
!
if (ibx > 0) then
mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,ib:ie,jb:je,kb:ke)**2, 1))
avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol
mnarr(4) = min(mnarr(4), minval(mag(:,:,:)))
mxarr(4) = max(mxarr(4), maxval(mag(:,:,:)))
tmp(:,:,:) = pdata%q(ibp,ib:ie,jb:je,kb:ke)
avarr(5) = avarr(5) + sum(tmp(:,:,:)) * dvol
mnarr(5) = min(mnarr(5), minval(tmp(:,:,:)))
mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:)))
end if
! calculate square root of density
!
sqd(:,:,:) = sqrt(pdata%q(idn,ib:ie,jb:je,kb:ke))
! get average, minimum and maximum values of sonic Mach number
!
if (ipr > 0) then
tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) &
/ sqrt(gamma * pdata%q(ipr,ib:ie,jb:je,kb:ke))
else
tmp(:,:,:) = vel(:,:,:) / csnd
end if
avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol
mnarr(6) = min(mnarr(6), minval(tmp(:,:,:)))
mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:)))
! get average, minimum and maximum values of Alfvénic Mach number
!
if (ibx > 0) then
tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:))
avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol
mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
end if
! associate the pointer with the next block on the list
@ -340,16 +443,36 @@ module integrals
#ifdef MPI
! sum the integral array from all processes
!
call reduce_sum_real_array(narr, arr(:), iret)
call reduce_sum_real_array(narr, inarr(:), iret)
! reduce average, minimum and maximum values
!
call reduce_sum_real_array(narr, avarr(:), iret)
call reduce_minimum_real_array(narr, mnarr(:), iret)
call reduce_maximum_real_array(narr, mxarr(:), iret)
#endif /* MPI */
! calculate the internal energy
!
if (ien > 0) arr(8) = arr(5) - arr(6) - arr(7)
if (ien > 0) inarr(8) = inarr(5) - inarr(6) - inarr(7)
! write down the integrals to the integrals file
! normalize the averages by the volume of domain
!
if (master) write(funit,"(i9,10(1x,1e18.8))") step, time, dt, arr(1:8)
avarr(:) = avarr(:) * voli
! write down the integrals and statistics to appropriate files
!
if (master) then
write(funit,"(i9,10(1x,1e18.8))") step, time, dt, inarr(1:8)
write(sunit,"(i9,23(1x,1e18.8))") step, time &
, avarr(1), mnarr(1), mxarr(1) &
, avarr(2), mnarr(2), mxarr(2) &
, avarr(3), mnarr(3), mxarr(3) &
, avarr(4), mnarr(4), mxarr(4) &
, avarr(5), mnarr(5), mxarr(5) &
, avarr(6), mnarr(6), mxarr(6) &
, avarr(7), mnarr(7), mxarr(7)
end if
#ifdef PROFILE
! stop accounting time for the integrals storing

6291
src/io.F90
View File

File diff suppressed because it is too large Load Diff

1302
src/mesh.F90
View File

File diff suppressed because it is too large Load Diff

212
src/problems.F90
View File

@ -232,12 +232,12 @@ module problems
use blocks , only : block_data
use constants , only : d2r
use coordinates, only : im, jm, km
use coordinates, only : ax, ay, az, adx, ady, adz
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 parameters , only : get_parameter_real
use parameters , only : get_parameter_real, get_parameter_integer
! local variables are not implicit by default
!
@ -249,12 +249,15 @@ module problems
! 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
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
!
@ -265,10 +268,19 @@ module problems
! local variables
!
integer :: i, j, k
integer :: i, j, k, ic, jc, kc
real(kind=8) :: xl, yl, zl, xu, yu, zu, rl, ru
real(kind=8) :: xb, yb, xt, yt
real(kind=8) :: dx, dy, dz, dxh, dyh, dzh, daxy
#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) :: sn, ph
#endif /* NDIMS == 3 */
real(kind=8) :: dx, dy, dz, dxh, dyh, dzh, dvol
real(kind=8) :: fc_amb, fc_ovr
! local arrays
@ -277,6 +289,13 @@ module problems
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 */
!
!-------------------------------------------------------------------------------
!
@ -299,6 +318,16 @@ module problems
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 overdense and ambient region densities
!
dn_amb = dens
@ -346,7 +375,33 @@ module problems
#else /* NDIMS == 3 */
dzh = 1.0d+00
#endif /* NDIMS == 3 */
daxy = dx * dy
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
!
@ -439,41 +494,122 @@ module problems
else
#if NDIMS == 3
! in 3D simply set the ambient values since the integration is more complex
! 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
! set the ambient region density
! update the integration factor depending on the subcell position
!
q(idn,i) = dn_amb
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
! set the ambient medium pressure
!
if (ipr > 0) q(ipr,i) = pr_amb
end do ! ic = 1, nsubgrid
end do ! jc = 1, nsubgrid
end do ! kc = 1, nsubgrid
#else /* NDIMS == 3 */
! calculate the bounds of area integration
! calculate the distance of remaining corners
!
xb = max(xl, sqrt(max(0.0d+00, r2 - yu * yu)))
xt = min(xu, sqrt(max(0.0d+00, r2 - yl * yl)))
yb = max(yl, sqrt(max(0.0d+00, r2 - xu * xu)))
yt = min(yu, sqrt(max(0.0d+00, r2 - xl * xl)))
rlu = xl * xl + yu * yu
rul = xu * xu + yl * yl
! integrate the area below the circle within the current cell for both
! functions, y = f(x) and x = g(y), and then average them to be sure that we
! are getting the ideal symmetry
! separate in the cases of which corners lay inside, and which outside
! the radius
!
fc_ovr = 0.5d+00 * (r2 * (asin(xt / radius) - asin(xb / radius)) &
+ (xt * yb - xb * yt)) - yl * (xt - xb)
fc_ovr = fc_ovr + (xb - xl) * dy
if (min(rlu, rul) >= r2) then
fc_amb = 0.5d+00 * (r2 * (asin(yt / radius) - asin(yb / radius)) &
+ (yt * xb - yb * xt)) - xl * (yt - yb)
fc_amb = fc_amb + (yb - yl) * dx
! 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
fc_ovr = 0.5d+00 * (fc_ovr + fc_amb)
! 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 / daxy
fc_ovr = fc_ovr / dvol
fc_amb = 1.0d+00 - fc_ovr
! integrate the density over the edge cells
@ -483,7 +619,6 @@ module problems
! integrate the pressure over the edge cells
!
if (ipr > 0) q(ipr,i) = fc_ovr * pr_ovr + fc_amb * pr_amb
#endif /* NDIMS == 3 */
end if
@ -504,6 +639,13 @@ module problems
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
!