[FLASH-BUGS] FPE in 3D mhd setup
Sean Matt
matt at physics.mcmaster.ca
Fri Feb 21 12:58:11 CST 2003
Hello,
We have been having a few problems getting a new 3D mhd test
problem setup to run. We am getting a floating point exception (FPE)
crash, as well as bad results (before the crash). We are running on Tru64
(OSF1 V5.1) using a Compaq fortran compiler (V5.5-1877-48BBF).
The test is the spin down of a magnetic rotating disk via the
transport of torsional alfven waves along field lines threading the disk.
Our setup files are attached to this email; please let me know if you have
a problem receiving them. There is an analytical solution for this
problem (see references listed in Config file), and we think it will be an
excellent test of the FLASH mhd.
Earlier, we have run the problem using the /hydro module, and with
no B fields. The initial setup is a dense, rotating cylinder in a cold
ambient medium. The hydro setup seems to work fine, and as the system
evolves the disk spins apart, as it should (there is nothing keeping it
together).
Currently, we have switched to the /hydro/mhd module, but
initialize magx, magy, magz to zero. When we run, we get a FPE. The
error occurs after 13 timesteps, in the current setup. If we initialize
magx, magy, magz to 1.0e-10 or to 1.0, the FPE still occurs, though
usually at a slightly different timestep.
Using dbx, we've been able to track down the FPE. It occurs in
the routine /source/hydro/mhd/mhd_fluxes.F90 (around line 76) at the
statement:
s = sqrt(Um(idn,i)/Up(idn,i-1))
We found, that at the time of the crash, one of the values of Um and Up
was positive, while the other was negative. This results in the attempted
calculation of the square root of a negative number, hence the FPE. Um
and Up are passed in from mhd_sweep.F90 just after being initialized by
the routine mhd_interpolate. We gather that Um and Up are
"time-interpolated data" at the left and right interfaces of the cell.
Since the index "idn" refers to the gas density, we're not sure that these
should ever be negative (though we are only just learning about the 8-wave
solver). We were unable to track the problem further. This is our first
major concern.
Second (this may or may not be related to the above), we have
looked at output data prior to the FPE crash at 13 timesteps. The data
doesn't look very good (esp. when compared to the hydro results of the
same setup at the same time). In particular, the density in the cylinder,
which is initially constant, quickly takes on a "mottled" appearance
(meaning that the values seem to fluxuate in space about the correct
value). The pressure in the cylinder also has the wrong value and is also
mottled. It is possible that whatever causes the FPE, also causes this
effect. We really have no idea. If this is an actual bug, it would be
crucial to fix.
We would appreciate any insight you may have about these problems.
Please let us know if you require more information.
-Sean
-------------- next part --------------
# Configuration file for the magnetic rotator problem
# Mouschovias, T. Ch., & Paleologou, E. V. 1980, ApJ, 237, 877
# Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791
# Stone, J. M., Hawley, J. F., Evans, C. R., & Norman, M. L. 1992,
# ApJ, 388, 415
REQUIRES driver
REQUIRES materials/eos/gamma
REQUIRES hydro/mhd
#REQUIRES hydro/mhd/divb_project
#REQUIRES hydro/
# Add custom variables below
D rho_disk Disk density
D rho_amb Ambient density
D B_z Vertical magnetic field strength
D z_disk Half thickness of disk
D r_disk Truncation radius of disk
D omega_disk Angular rotation rate of disk
D const_press Pressure everywhere
D xctr Center for the angular momentum calculation
D yctr Center for the angular momentum calculation
D zctr Center for the angular momentum calculation
PARAMETER rho_disk REAL 10.0
PARAMETER rho_amb REAL 1.0
PARAMETER B_z REAL 0.282095
PARAMETER z_disk REAL 1.0
PARAMETER r_disk REAL 5.0
PARAMETER omega_disk REAL 1.0
PARAMETER const_press REAL 0.00006
PARAMETER xctr REAL 0.0
PARAMETER yctr REAL 0.0
PARAMETER zctr REAL 0.0
-------------- next part --------------
# Runtime parameters for the magnetic rotator problem.
# Size of computational volume
Nblockx = 1
Nblocky = 1
Nblockz = 1
xmin = -15.
xmax = 15.
ymin = -15.
ymax = 15.
zmin = 0.
zmax = 30.
# Boundary conditions
xl_boundary_type = "outflow"
xr_boundary_type = "outflow"
yl_boundary_type = "outflow"
yr_boundary_type = "outflow"
zl_boundary_type = "reflecting"
zr_boundary_type = "outflow"
# Parameters for initial model
rho_disk = 10.0
rho_amb = 1.0
#B_z = 0.282095
B_z = 0.0
z_disk = 1.0
r_disk = 5.0
omega_disk = 1.0
const_press = 0.00006
# Specific heats ratio
gamma = 1.6666666667
# Turn off dual energy formulation.
eint_switch = 0.0
# Simulation (grid, time, I/O) parameters
dtini = 1.e-3
dtmin = 1.e-3
lrefine_min = 2
lrefine_max = 4
restart = .false.
CPNumber = 0
basenm = "MagRot_l5_"
log_file = "MagRot_l5.log"
run_comment = "Magnetic Rotator Problem"
plot_var_1 = "dens"
plot_var_2 = "pres"
plot_var_3 = "velx"
plot_var_4 = "vely"
plot_var_5 = "velz"
plot_var_6 = "magx"
plot_var_7 = "magy"
plot_var_8 = "magz"
tplot = 1.0
trstrt = 5.0
nrstrt = 9999999
nend = 9999999
wall_clock_checkpoint = 9999999
tmax = 5.0
# MHD parameters
cfl = 1.0
killdivb = .FALSE.
-------------- next part --------------
!******************************************************************************
!
! Subroutine: init_block()
!
! Description: Initializes fluid data for a specified block.
! This version sets up the magnetic rotator MHD problem.
!
! References: Mouschovias, T. Ch., & Paleologou, E. V. 1980, ApJ, 237, 877
! Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791
! Stone, J. M., Hawley, J. F., Evans, C. R., & Norman, M. L.
! 1992 , ApJ, 388, 415
!
! Parameters: block_no The number of the block to initialize
!
!
subroutine init_block (block_no)
!
!==============================================================================
use logfile, ONLY: write_logfile
use physical_constants
use runtime_parameters, ONLY : get_parm_from_context
use dBase, ONLY: iu_bnd, ju_bnd, ku_bnd, &
& ndim, &
& nguard, ionmax, &
& dBaseKeyNumber, dBaseSpecies, &
& dBaseGetCoords, dBasePutData, &
& dBasePropertyInteger
implicit none
integer, INTENT(in) :: block_no
integer, save :: MyPE, MasterPE
integer, save :: idens, ipres, iener, igamc, igame
integer, save :: ivelx, ively, ivelz
integer, save :: imagx, imagy, imagz
integer, save :: iPoint
integer, save :: iXcoord, iYcoord, iZcoord
integer, save :: izn, iznl, iznr
real, save :: const_PI
real, save :: gamma
real, save :: rho_disk, rho_amb, B_z, z_disk, omega_disk, const_press
real, save :: r_disk
real :: rho, pres, eng, gamc, game
real :: vx, vy, vz, bx, by, bz
integer, parameter :: q = max(iu_bnd,ju_bnd,ku_bnd)
real, DIMENSION(q) :: x, y, z, zL, zR
real :: xx, yy, zz, zzL, zzR, cyl_r, theta
integer :: i, j, k, n
logical, save:: first_call = .true.
!==============================================================================
!------------------------------------------------------------------------------
! Compute direction cosines for the shock normal.
if (first_call) then
MyPE = dBasePropertyInteger('MyProcessor')
MasterPE = dBasePropertyInteger('MasterProcessor')
call get_constant_from_db("pi", const_PI)
idens = dBaseKeyNumber('dens')
ivelx = dBaseKeyNumber('velx')
ively = dBaseKeyNumber('vely')
ivelz = dBaseKeyNumber('velz')
imagx = dBaseKeyNumber('magx')
imagy = dBaseKeyNumber('magy')
imagz = dBaseKeyNumber('magz')
iener = dBaseKeyNumber('ener')
ipres = dBaseKeyNumber('pres')
igamc = dBaseKeyNumber('gamc')
igame = dBaseKeyNumber('game')
iPoint = dBaseKeyNumber('Point')
iXcoord = dBaseKeyNumber('xCoord')
iYcoord = dBaseKeyNumber('yCoord')
iZcoord = dBaseKeyNumber('zCoord')
izn = dBaseKeyNumber('zn' )
! iznl = dBaseKeyNumber('znl' )
! iznr = dBaseKeyNumber('znr' )
call get_parm_from_context('rho_disk', rho_disk)
call get_parm_from_context('rho_amb', rho_amb)
call get_parm_from_context('B_z', B_z)
call get_parm_from_context('z_disk', z_disk)
call get_parm_from_context('r_disk', r_disk)
call get_parm_from_context('omega_disk', omega_disk)
call get_parm_from_context('const_press', const_press)
call get_parm_from_context('gamma', gamma)
!------------------------------------------------------------------------------
! Write a message to stdout describing the problem setup.
if (MyPE .eq. MasterPE) then
call write_logfile("flash: initializing the magnetic rotator problem.")
write (*,*)
write (*,*) 'Flash: initializing the magnetic rotator problem.'
write (*,*)
write (*,1) 'rho_disk = ', rho_disk, 'rho_amb = ', rho_amb
write (*,1) 'B_z = ', B_z, 'z_disk = ', z_disk, &
& 'omega_disk = ', omega_disk
write (*,2) 'gamma = ', gamma, 'ndim = ', ndim
write (*,*)
1 format (1X, 1P, 5(A7, E13.6, :, 1X))
2 format (1X, 1P, 2(A7, E13.6, 1X), A7, I13)
endif
first_call = .false.
endif
!------------------------------------------------------------------------------
! Loop through the block and initialize all values.
x(:) = 0.0
y(:) = 0.0
z(:) = 0.0
! zL(:) = 0.0
! zR(:) = 0.0
! If simulation is 3D, initialze, with rotation axis along z.
if (ndim == 3) then
!***** Get coordinates. I think this gets the center cell values.
call dBaseGetCoords(izn, iXcoord, block_no, x)
call dBaseGetCoords(izn, iYcoord, block_no, y)
call dBaseGetCoords(izn, iZcoord, block_no, z)
! call dBaseGetCoords(iznl, iZcoord, block_no, zL)
! call dBaseGetCoords(iznr, iZcoord, block_no, zR)
!***** Loop through the block.
do k = 1, ku_bnd
zz = z(k)
! zzL= zL(k)
! zzR= zR(k)
do j = 1, ju_bnd
yy = y(j)
do i = 1, iu_bnd
xx = x (i)
!***** Calculate the cylindrical r and azimuthal angle.
cyl_r = sqrt(xx*xx + yy*yy)
if (cyl_r > 0.0) then
theta = ACOS(xx/cyl_r)
if (yy < 0.0) theta = const_PI*2.0 - theta
else
theta = 0.0
endif
!***** Set initial values.
game = gamma
gamc = gamma
rho = rho_amb
pres = const_press
vx = 0.0
vy = 0.0
vz = 0.0
bx = 0.0
by = 0.0
bz = B_z
!***** Create the rotating disk.
if ((ABS(zz) < z_disk) .and. (cyl_r < r_disk)) then
rho = rho_disk
vx =-omega_disk * cyl_r * SIN(theta)
vy = omega_disk * cyl_r * COS(theta)
! elseif ((((ABS(zzL) < z_disk) .and. (ABS(zzR) > z_disk)) .or. &
! & ((ABS(zzL) > z_disk) .and. (ABS(zzR) < z_disk))) &
! & .and. (cyl_r < r_disk)) then
!
! rho = 0.25 * rho_disk
endif
eng = 0.5 * (vx*vx + vy*vy + vz*vz) + &
& pres / (gamma - 1.0) / rho
call dBasePutData(idens, iPoint, i, j, k, block_no, rho)
call dBasePutData(ipres, iPoint, i, j, k, block_no, pres)
call dBasePutData(iener, iPoint, i, j, k, block_no, eng)
call dBasePutData(ivelx, iPoint, i, j, k, block_no, vx )
call dBasePutData(ively, iPoint, i, j, k, block_no, vy )
call dBasePutData(ivelz, iPoint, i, j, k, block_no, vz )
call dBasePutData(imagx, iPoint, i, j, k, block_no, bx )
call dBasePutData(imagy, iPoint, i, j, k, block_no, by )
call dBasePutData(imagz, iPoint, i, j, k, block_no, bz )
call dBasePutData(igamc, iPoint, i, j, k, block_no, gamc)
call dBasePutData(igame, iPoint, i, j, k, block_no, game)
do n=1,ionmax
call dBasePutData(dBaseSpecies(n),iPoint,i,j,k,block_no, 1.0)
enddo
enddo
enddo
enddo
! If simulation is 1D, print error message.
elseif (ndim < 2) then
if (MyPE .eq. MasterPE) then
write(*,*) 'Simulation must have at least 2 dimensions!'
endif
! If simulation is 2D, initialize axisymmetric setup, using x,y as r,z.
!***** AXISYMMETRY DOES NOT WORK !!!! *****
elseif (ndim == 2) then
if (MyPE .eq. MasterPE) then
write(*,*) 'Simulation does not work in 2 dimensions!'
endif
!***** Get coordinates. I think (hope) this gets the center cell
!***** values.
call dBaseGetCoords(izn, iXcoord, block_no, x)
call dBaseGetCoords(izn, iYcoord, block_no, y)
!***** Loop through the block.
do k = 1, ku_bnd
do j = 1, ju_bnd
yy = y(j)
do i = 1, iu_bnd
xx = x (i)
!***** Set initial values.
game = gamma
gamc = gamma
rho = rho_amb
pres = const_press
vx = 0.0
vy = 0.0
vz = 0.0
bx = 0.0
by = B_z
bz = 0.0
!***** Create the disk.
if (yy .lt. z_disk) then
rho = rho_disk
vz = omega_disk * xx
endif
eng = 0.5 * (vx*vx + vy*vy + vz*vz) + &
& pres / (gamma - 1.0) / rho
call dBasePutData(idens, iPoint, i, j, k, block_no, rho)
call dBasePutData(ipres, iPoint, i, j, k, block_no, pres)
call dBasePutData(iener, iPoint, i, j, k, block_no, eng)
call dBasePutData(ivelx, iPoint, i, j, k, block_no, vx )
call dBasePutData(ively, iPoint, i, j, k, block_no, vy )
call dBasePutData(ivelz, iPoint, i, j, k, block_no, vz )
call dBasePutData(imagx, iPoint, i, j, k, block_no, bx )
call dBasePutData(imagy, iPoint, i, j, k, block_no, by )
call dBasePutData(imagz, iPoint, i, j, k, block_no, bz )
call dBasePutData(igamc, iPoint, i, j, k, block_no, gamc)
call dBasePutData(igame, iPoint, i, j, k, block_no, game)
do n=1,ionmax
call dBasePutData(dBaseSpecies(n),iPoint,i,j,k,block_no, 1.0)
enddo
enddo
enddo
enddo
else
if (MyPE .eq. MasterPE) then
write(*,*) 'Simulations must have 1, 2, or 3 dimensions!'
endif
endif
!==============================================================================
end subroutine init_block
-------------- next part --------------
!!****f* source/io/wr_integrals
!!
!! NAME
!! wr_integrals
!!
!! SYNOPSIS
!! call wr_integrals(isfirst, simtime)
!! call wr_integrals(integer, real)
!!
!! DESCRIPTION
!! Compute the values of integral quantities (eg. total energy)
!! and write them to an ASCII file. If this is the initial step,
!! create the file and write a header to it before writing the
!! data.
!!
!! Presently, this supports 1, 2, and 3-d Cartesian geometry and 2-d
!! cylindrical geometry (r,z). More geometries can be added by modifying
!! the volume of each zone (dvol) computed in SumLocalIntegrals.
!!
!!***
subroutine wr_integrals (isfirst, simtime)
!=============================================================================
use dBase, ONLY : &
& nxb, nyb, nzb, ndim, &
& nguard, k2d, k3d, &
& dBasePropertyInteger, dBaseLocalBlockCount, dBaseNodeType
use runtime_parameters, ONLY: get_parm_from_context, &
& global_parm_context
implicit none
#include "mpif.h"
integer :: isfirst
integer, parameter :: nsum = 11 ! Number of summed quantities;
! coordinate w/SumLocalIntegrals
real :: simtime, sum(nsum), lsum(nsum)
integer :: i, j
integer :: imin, imax, jmin, jmax, kmin, kmax
integer :: funit = 99
integer :: error
character (len=80) :: fname = "flash.dat"
integer :: lnblocks
integer :: nodetype
logical, save :: restart
! store the master and current processor numbers
integer, save :: MyPE, MasterPE
logical, save :: firstCall = .TRUE.
!=============================================================================
if (firstCall) then
! get the processor information
MyPE = dBasePropertyInteger("MyProcessor")
MasterPE = dBasePropertyInteger("MasterProcessor")
! get the runtime parameters
call get_parm_from_context(global_parm_context, &
& 'restart', restart)
firstCall = .FALSE.
endif
! get the number of blocks on this processor
lnblocks = dbaseLocalBlockCount()
! Sum quantities over all locally held leaf-node blocks.
sum = 0.
lsum = 0.
imin = nguard + 1 ! Determine index ranges of
imax = imin + nxb - 1 ! cells in each block.
jmin = nguard*k2d + 1
jmax = jmin + nyb - 1
kmin = nguard*k3d + 1
kmax = kmin + nzb - 1
do i = 1, lnblocks
nodetype = dBaseNodeType(i)
if (nodetype == 1) then
call SumLocalIntegrals (lsum, i, imin, imax, jmin, jmax, &
& kmin, kmax)
endif
enddo
!=============================================================================
! Now the master PE sums the local contributions from all of
! the processors and writes the total to a file.
call MPI_Reduce (lsum, sum, nsum, MPI_Double_Precision, MPI_Sum, &
& MasterPE, MPI_Comm_World, error)
if (MyPE == MasterPE) then
! create the file from scratch if it is a not a restart simulation,
! otherwise append to the end of the file
if (isfirst == 0) then
open (funit, file=fname, position='APPEND')
else
if (.NOT. restart) then
open (funit, file=fname)
call WriteIntHeader (funit)
else
open (funit, file=fname, position='APPEND')
write (funit, 11)
11 format('# simulation restarted')
endif
endif
write (funit, 10) simtime, sum ! Write the sums to the file.
10 format (1X, 1P, 50(E25.18, :, 1X))
close (funit) ! Close the file.
endif
call MPI_Barrier (MPI_Comm_World, error)
!=============================================================================
return
end
!*****************************************************************************
! Routine: SumLocalIntegrals()
! Description: Accumulate values of the energy, angular momentum, etc. from
! a given range of cells in a given block. The results are
! stored in a given vector (lsum), which is assumed to have been
! initialized by the calling routine.
!
subroutine SumLocalIntegrals (lsum, block_no, imin, imax, jmin, &
& jmax, kmin, kmax)
!=============================================================================
use physical_constants, ONLY: get_constant_from_db
use runtime_parameters, ONLY: get_parm_from_context, &
& GLOBAL_PARM_CONTEXT
use dBase, ONLY: nxb, nyb, nzb, mdim, ndim, nguard, k2d, k3d, &
& geom_cartesian, geom_cylrad, &
& dBaseBlockSize, &
& dBaseBlockCoord, &
& dBaseGetCoords, &
& dBaseKeyNumber, &
& dBaseGetDataPtrSingleBlock
implicit none
integer :: block_no
integer :: imin, imax, jmin, jmax, kmin, kmax
real :: lsum(*)
integer, save :: idens, ivelx, ively, ivelz, ipres, iener, &
& igame, igpot
integer, save :: igeomx, igeomy, igeomz
integer, save :: iXcoord
integer, save :: iznl, iznr
integer, parameter :: q = max(nxb+2*nguard, &
& nyb+2*nguard*k2d, &
& nzb+2*nguard*k3d)
real, dimension(q) :: xl, xr
real :: size(mdim), coord(mdim)
integer :: i, j, k
real :: dvol, delx, dely, delz
real :: xdist, ydist, zdist
real, save :: pi
real, save :: xctr, yctr, zctr
real, DIMENSION(:,:,:,:), POINTER :: solnData
logical, save :: firstCall = .TRUE.
!=============================================================================
if (firstCall) then
idens = dBaseKeyNumber('dens')
ivelx = dBaseKeyNumber('velx')
ively = dBaseKeyNumber('vely')
ivelz = dBaseKeyNumber('velz')
ipres = dBaseKeyNumber('pres')
iener = dBaseKeyNumber('ener')
igame = dBaseKeyNumber('game')
igpot = dBaseKeyNumber('gpot')
iXcoord = dBaseKeyNumber('xCoord')
iznl = dBaseKeyNumber('znl')
iznr = dBaseKeyNumber('znr')
! get the center coordinates
! should be improved by computing the center of mass system (cms)
call get_parm_from_context("xctr", xctr)
call get_parm_from_context("yctr", yctr)
call get_parm_from_context("zctr", zctr)
! get the geometry information for the computation domain
call get_parm_from_context(GLOBAL_PARM_CONTEXT, &
& 'igeomx', igeomx)
call get_parm_from_context(GLOBAL_PARM_CONTEXT, &
& 'igeomy', igeomy)
call get_parm_from_context(GLOBAL_PARM_CONTEXT, &
& 'igeomz', igeomz)
! grab the necessary physical constants
call get_constant_from_db("pi", pi)
firstCall = .FALSE.
endif
! Hydro variables must be defined, else we exit without adding to local sums.
if ((idens < 0) .or. &
& (ipres < 0) .or. &
& (ivelx < 0) .or. &
& (ively < 0) .or. &
& (ivelz < 0) .or. &
& (iener < 0) .or. &
& (igame < 0)) &
& return
! get the size of the current block, and compute the size and volume of a
! cell
!
! deal with the different geometries here.
size(:) = dBaseBlockSize(block_no)
coord(:) = dBaseBlockCoord(block_no)
delx = size(1) / real(nxb)
dely = size(2) / real(nyb)
delz = size(3) / real(nzb)
call dBaseGetCoords(iznl, iXcoord, block_no, xl)
! get a pointer to the current block of data
solnData => dBaseGetDataPtrSingleBlock(block_no)
!-----------------------------------------------------------------------------
! Cartesian geometries
!-----------------------------------------------------------------------------
! for Cartesian geometries, the volume of a cell is constant within a block,
! so we can compute it outside of the loops over all the zones in a block
if (igeomx == 0 .AND. igeomy == 0 .AND. igeomz == 0) then
dvol = delx
if (ndim >= 2) dvol = dvol * dely
if (ndim == 3) dvol = dvol * delz
! Sum contributions from the indicated range of cells.
do k = kmin, kmax
zdist = (coord(3) - zctr) + real(k - 1 - kmax/2)*delz
do j = jmin, jmax
ydist = (coord(2) - yctr) + real(j - 1 - jmax/2)*dely
do i = imin, imax
xdist = (coord(1) - xctr) + real(i - 1 - imax/2)*delx
! mass
lsum(1) = lsum(1) + solnData(idens,i,j,k)*dvol
! momentum
lsum(2) = lsum(2) + solnData(idens,i,j,k) * &
& solnData(ivelx,i,j,k)*dvol
lsum(3) = lsum(3) + solnData(idens,i,j,k) * &
& solnData(ively,i,j,k)*dvol
lsum(4) = lsum(4) + solnData(idens,i,j,k) * &
& solnData(ivelz,i,j,k)*dvol
! total energy
lsum(5) = lsum(5) + solnData(iener,i,j,k) * &
& solnData(idens,i,j,k)*dvol
! kinetic energy
lsum(6) = lsum(6) + 0.5*solnData(idens,i,j,k) * &
& (solnData(ivelx,i,j,k)**2+ &
& solnData(ively,i,j,k)**2+ &
& solnData(ivelz,i,j,k)**2)*dvol
! internal energy
lsum(7) = lsum(7) + solnData(ipres,i,j,k) / &
& (solnData(igame,i,j,k)-1.)*dvol
! potential energy
if (igpot >= 0) &
& lsum(8) = lsum(8) + 0.5*solnData(idens,i,j,k) * &
& solnData(igpot,i,j,k)*dvol
! x-angular momentum
lsum(9) = lsum(9) + solnData(idens,i,j,k) &
& * dvol * ( ydist * solnData(ivelz,i,j,k) &
& - zdist * solnData(ively,i,j,k) )
! y-angular momentum
lsum(10) = lsum(10) + solnData(idens,i,j,k) &
& * dvol * ( zdist * solnData(ivelx,i,j,k) &
& - xdist * solnData(ivelz,i,j,k) )
! z-angular momentum
lsum(11) = lsum(11) + solnData(idens,i,j,k) &
& * dvol * ( xdist * solnData(ively,i,j,k) &
- ydist * solnData(ivelx,i,j,k) )
enddo
enddo
enddo
!-----------------------------------------------------------------------------
! 2-d cylindrical geometry
!-----------------------------------------------------------------------------
elseif (igeomx == geom_cylrad .AND. &
& igeomy == geom_cartesian .AND. &
& ndim == 2) then
xl(:) = 0.0
xr(:) = 0.0
call dBaseGetCoords(iznl, iXcoord, block_no, xl)
call dBaseGetCoords(iznr, iXcoord, block_no, xr)
! Sum contributions from the indicated range of cells.
do k = kmin, kmax
do j = jmin, jmax
do i = imin, imax
! now the volume factor is changing from zone to zone within the block,
! so we need to compute it for each zone separately
!
! the volume here is annular
dvol = 2.0*pi*(xr(i)**2 - xl(i)**2)*dely
! mass
lsum(1) = lsum(1) + solnData(idens,i,j,k)*dvol
! momentum
lsum(2) = lsum(2) + solnData(idens,i,j,k) * &
& solnData(ivelx,i,j,k)*dvol
lsum(3) = lsum(3) + solnData(idens,i,j,k) * &
& solnData(ively,i,j,k)*dvol
lsum(4) = lsum(4) + solnData(idens,i,j,k) * &
& solnData(ivelz,i,j,k)*dvol
! total energy
lsum(5) = lsum(5) + solnData(iener,i,j,k) * &
& solnData(idens,i,j,k)*dvol
! kinetic energy
lsum(6) = lsum(6) + 0.5*solnData(idens,i,j,k) * &
& (solnData(ivelx,i,j,k)**2+ &
& solnData(ively,i,j,k)**2+ &
& solnData(ivelz,i,j,k)**2)*dvol
! internal energy
lsum(7) = lsum(7) + solnData(ipres,i,j,k) / &
& (solnData(igame,i,j,k)-1.)*dvol
! potential energy
if (igpot >= 0) &
& lsum(8) = lsum(8) + 0.5*solnData(idens,i,j,k) * &
& solnData(igpot,i,j,k)*dvol
! lsum(9) = lsum(9) + 0. ! <---- NEED ANGULAR MOMENTUM
! lsum(10) = lsum(10) + 0. ! <---- NEED ANGULAR MOMENTUM
! lsum(11) = lsum(11) + 0. ! <---- NEED ANGULAR MOMENTUM
enddo
enddo
enddo
else
! the geometry is not supported. We should report this to the .dat file,
! but for now, just return
return
endif
return
end
!****************************************************************************
!
! Routine: WriteIntHeader()
!
! Description: Writes header information to the integral output file. Should
! coordinate the labels used here with the values computed in
! SumLocalIntegrals().
subroutine WriteIntHeader (funit)
!=============================================================================
integer :: funit
!=============================================================================
write (funit, 10) '#time ', 'mass', 'x-momentum', &
& 'y-momentum', &
& 'z-momentum', 'energy', 'kinetic', 'internal', &
& 'potential', 'x-angular momentum', &
& 'y-angular momentum', 'z-angular momentum'
10 format (50(A25, :, 1X))
!=============================================================================
return
end
More information about the flash-bugs
mailing list