!!****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