deleted cmfd utils module, added new files to objects and dependencies, took debug default off of makefile

This commit is contained in:
Bryan Herman 2012-01-27 20:32:09 -05:00
parent 30cdfc0ca2
commit 29d9d7c45d
4 changed files with 31 additions and 1007 deletions

View file

@ -13,22 +13,30 @@ ace.o: string.o
ace_header.o: constants.o
ace_header.o: endf_header.o
cmfd_execute.o: cmfd_utils.o
cmfd_execute.o: global.o
cmfd_execute.o: mesh.o
cmfd_execute.o: mesh_header.o
cmfd_execute.o: tally_header.o
cmfd_execute.o: timing.o
cmfd_data.o: datatypes.o
cmfd_data.o: global.o
cmfd_data.o: mesh.o
cmfd_data.o: string.o
cmfd_data.o: xml-fortran/templates/cmfd_t.o
cmfd_utils.o: constants.o
cmfd_utils.o: datatypes.o
cmfd_utils.o: global.o
cmfd_utils.o: hdf5_interface.o
cmfd_utils.o: mesh.o
cmfd_utils.o: mesh_header.o
cmfd_utils.o: string.o
cmfd_utils.o: vtk_writer.o
cmfd_utils.o: xml-fortran/templates/cmfd_t.o
cmfd_execute.o: cmfd_data.o
cmfd_execute.o: cmfd_power_solver.o
cmfd_execute.o: cmfd_slepc_solver.o
cmfd_execute.o: cmfd_snes_solver.o
cmfd_loss_operator.o: global.o
cmfd_power_solver.o: cmfd_loss_operator.o
cmfd_power_solver.o: cmfd_prod_operator.o
cmfd_prod_operator.o: global.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: cmfd_loss_operator.o
cmfd_snes_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: cmfd_slepc_solver.o
cross_section.o: ace_header.o
cross_section.o: constants.o
@ -116,7 +124,7 @@ initialize.o: string.o
initialize.o: tally.o
initialize.o: timing.o
input_xml.o: cmfd_utils.o
input_xml.o: cmfd_data.o
input_xml.o: constants.o
input_xml.o: datatypes.o
input_xml.o: error.o

View file

@ -16,7 +16,7 @@ include OBJECTS
#===============================================================================
COMPILER = petsc
DEBUG = yes
DEBUG = no
PROFILE = no
OPTIMIZE = no
USE_MPI = no

View file

@ -2,9 +2,14 @@ objects = \
ace.o \
ace_header.o \
bank_header.o \
cmfd_data.o \
cmfd_execute.o \
cmfd_header.o \
cmfd_utils.o \
cmfd_loss_operator.o \
cmfd_power_solver.o \
cmfd_prod_operator.o \
cmfd_slepc_solver.o \
cmfd_snes_solver.o \
cross_section.o \
datatypes.o \
datatypes_header.o \

View file

@ -1,989 +0,0 @@
module cmfd_utils
use cmfd_header
use constants
use datatypes, only: dict_add_key, dict_get_key
use error, only: fatal_error, warning
use global
use mesh, only: mesh_indices_to_bin
use mesh_header, only: StructuredMesh
use string
use tally_header, only: TallyObject, TallyScore
implicit none
contains
!===============================================================================
! READ_INPUT reads the CMFD input file and organizes it into a data structure
!===============================================================================
subroutine read_cmfd_xml()
use xml_data_cmfd_t
integer :: ng=1 ! number of energy groups (default 1)
integer :: n_words ! number of words read
logical :: file_exists ! does cmfd.xml exist?
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: words(MAX_WORDS)
! read cmfd infput file
filename = "cmfd.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
write(*,*) "Cannot perform CMFD"
STOP
end if
! parse cmfd.xml file
call read_xml_file_cmfd_t(filename)
! set spatial dimensions in cmfd object
cmfd % indices(1:3) = mesh_ % dimension(1:3) ! sets spatial dimensions
! get number of energy groups
if (len_trim(mesh_ % energy) > 0) then
call split_string(mesh_ % energy, words, n_words)
ng = n_words - 1
end if
cmfd % indices(4) = ng ! sets energy group dimension
! set global albedo
cmfd % albedo = mesh_ % albedo
! get acceleration map
if (associated(mesh_ % map)) then
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
& cmfd % indices(3)))
cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3)))
end if
! check for core map activation by printing note
if (allocated(cmfd % coremap)) print *,"Core Map Overlay Activated"
! create tally objects
call create_cmfd_tally()
end subroutine read_cmfd_xml
!===============================================================================
! ALLOCATE_CMFD allocates all of the space for the cmfd object based on tallies
!===============================================================================
subroutine allocate_cmfd()
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! allocate flux, cross sections and diffusion coefficient
if (.not. allocated(cmfd % flux)) allocate(cmfd % flux(ng,nx,ny,nz))
if (.not. allocated(cmfd % totalxs)) allocate(cmfd % totalxs(ng,nx,ny,nz))
if (.not. allocated(cmfd % p1scattxs)) allocate(cmfd % p1scattxs(ng,nx,ny,nz))
if (.not. allocated(cmfd % scattxs)) allocate(cmfd % scattxs(ng,ng,nx,ny,nz))
if (.not. allocated(cmfd % nfissxs)) allocate(cmfd % nfissxs(ng,ng,nx,ny,nz))
if (.not. allocated(cmfd % diffcof)) allocate(cmfd % diffcof(ng,nx,ny,nz))
! allocate dtilde and dhat
if (.not. allocated(cmfd % dtilde)) allocate(cmfd % dtilde(6,ng,nx,ny,nz))
if (.not. allocated(cmfd % dhat)) allocate(cmfd % dhat(6,ng,nx,ny,nz))
! allocate dimensions for each box (here for general case)
if (.not. allocated(cmfd % hxyz)) allocate(cmfd % hxyz(3,nx,ny,nz))
! allocate cmfd fission source pdf
!allocate( cmfd % sourcepdf(ng,nx,ny,nz) )
! allocate surface currents
if (.not. allocated(cmfd % current)) allocate(cmfd % current(12,ng,nx,ny,nz))
! allocate for coremap
if (cmfd_only) allocate(cmfd % coremap(nx,ny,nz))
end subroutine allocate_cmfd
!===============================================================================
! GET_MATRIX_IDX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
function get_matrix_idx(g,i,j,k,ng,nx,ny)
! arguments
integer :: get_matrix_idx ! the index location in matrix
integer :: i ! current x index
integer :: j ! current y index
integer :: k ! current z index
integer :: g ! current group index
integer :: ng ! max energy groups
integer :: nx ! maximum cells in x direction
integer :: ny ! maximum cells in y direction
! local variables
integer :: nidx ! index in matrix
! check if coremap is used
if (allocated(cmfd % coremap)) then
! get idx from core map
nidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! compute index
nidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
! record value to function
get_matrix_idx = nidx
end function get_matrix_idx
!===============================================================================
! PRINT_CMFD is a test routine to check if info from tally is being accessed
!===============================================================================
subroutine print_cmfd()
integer :: bins(TALLY_TYPES) ! bin for tally_types, for filters
integer :: ijk(3) ! indices for mesh cell where tally is
integer :: score_index ! index in tally score to get value
real(8) :: tally_val ! value of tally being extracted
type(TallyObject), pointer :: t ! pointer for a tally object
type(StructuredMesh), pointer :: m ! pointer for mesh object
! associate pointers with objects
t => tallies(3)
m => meshes(t % mesh)
! set all bins to 1
bins = 1
! get mesh indices, first we will first force to 1,1,1
! ijk = (/ 1, 1, 1 /)
! apply filters, here we will just try a mesh filter first
! bins(T_MESH) = mesh_indices_to_bin(m,ijk)
! calculate score index from bins
! score_index = sum((bins - 1) * t%stride) + 1
! get value from tally object
! tally_val = t%scores(score_index,2)%val
! write value to file
! write(7,*) "Tally value is:",tally_val
! Left Surface
ijk = (/ 1-1, 1, 1 /)
score_index = sum(t % stride(1:3) * ijk) + IN_RIGHT
print *, "Outgiong Current from Left", t % scores(score_index,1) % val
score_index = sum(t % stride(1:3) * ijk) + OUT_RIGHT
print *, "Incoming Current from Left", t % scores(score_index,1) % val
end subroutine print_cmfd
!===============================================================================
! CREATE_CMFD_TALLY creates the tally object for OpenMC to process for CMFD
! accleration.
! There are 3 tally types:
! 1: Only an energy in filter-> flux,total,p1 scatter
! 2: Energy in and energy out filter-> nu-scatter,nu-fission
! 3: Surface current
!===============================================================================
subroutine create_cmfd_tally()
use xml_data_cmfd_t
integer :: i ! loop counter
integer :: j ! loop counter
integer :: id ! user-specified identifier
integer :: index ! index in mesh array
integer :: n ! size of arrays in mesh specification
integer :: ng=1 ! number of energy groups (default 1)
integer :: n_words ! number of words read
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: words(MAX_WORDS)
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
! parse cmfd.xml file
filename = trim(path_input) // "cmfd.xml"
call read_xml_file_cmfd_t(filename)
! allocate mesh
n_meshes = 1
allocate(meshes(n_meshes))
m => meshes(1)
! set mesh id
m % id = 1
! set mesh type to rectangular
m % type = LATTICE_RECT
! determine number of dimensions for mesh
n = size(mesh_ % dimension)
if (n /= 2 .and. n /= 3) then
message = "Mesh must be two or three dimensions."
call fatal_error()
end if
m % n_dimension = n
! allocate attribute arrays
allocate(m % dimension(n))
allocate(m % origin(n))
allocate(m % width(n))
allocate(m % upper_right(n))
! read dimensions in each direction
m % dimension = mesh_ % dimension
! read mesh origin location
if (m % n_dimension /= size(mesh_ % origin)) then
message = "Number of entries on <origin> must be the same as " // &
"the number of entries on <dimension>."
call fatal_error()
end if
m % origin = mesh_ % origin
! read mesh widths
if (size(mesh_ % width) /= size(mesh_ % origin)) then
message = "Number of entries on <width> must be the same as " // &
"the number of entries on <origin>."
call fatal_error()
end if
m % width = mesh_ % width
! set upper right coordinate
m % upper_right = m % origin + m % dimension * m % width
! add mesh to dictionary
call dict_add_key(mesh_dict, m % id, 1)
! allocate tallies
n_tallies = 3
allocate(tallies(n_tallies))
! begin loop around tallies
do i = 1,n_tallies
t => tallies(i)
! allocate arrays for number of bins and stride in scores array
allocate(t % n_bins(TALLY_TYPES))
allocate(t % stride(TALLY_TYPES))
! initialize number of bins and stride
t % n_bins = 0
t % stride = 0
! record tally id which is equivalent to loop number
t % id = i
! set mesh filter mesh id = 1
t % mesh = 1
m => meshes(1)
t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
! read and set incoming energy mesh filter
if (len_trim(mesh_ % energy) > 0) then
call split_string(mesh_ % energy,words,n_words)
ng = n_words
allocate(t % energy_in(n_words))
do j = 1,n_words
t % energy_in(j) = str_to_real(words(j))
end do
t % n_bins(T_ENERGYIN) = n_words - 1
end if
if (i == 1) then
! allocate macro reactions
allocate(t % macro_bins(3))
t % n_macro_bins = 3
! set macro_bins
t % macro_bins(1) % scalar = MACRO_FLUX
t % macro_bins(2) % scalar = MACRO_TOTAL
t % macro_bins(3) % scalar = MACRO_SCATTER_1
else if (i == 2) then
! read and set outgoing energy mesh filter
if (len_trim(mesh_ % energy) > 0) then
call split_string(mesh_ % energy, words, n_words)
allocate(t % energy_out(n_words))
do j = 1, n_words
t % energy_out(j) = str_to_real(words(j))
end do
t % n_bins(T_ENERGYOUT) = n_words - 1
end if
! allocate macro reactions
allocate(t % macro_bins(2))
t % n_macro_bins = 2
! set macro_bins
t % macro_bins(1) % scalar = MACRO_NU_SCATTER
t % macro_bins(2) % scalar = MACRO_NU_FISSION
else if (i == 3) then
! allocate macro reactions
allocate(t % macro_bins(1))
t % n_macro_bins = 1
! set macro bins
t % macro_bins(1) % scalar = MACRO_CURRENT
t % surface_current = .true.
! since the number of bins for the mesh filter was already set
! assuming it was a flux tally, we need to adjust the number of
! bins
t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
! get pointer to mesh
id = t % mesh
index = dict_get_key(mesh_dict, id)
m => meshes(index)
! we need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
if (size(m % dimension) == 2) then
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
& product(m % dimension + 1) * 4
elseif (size(m % dimension) == 3) then
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
product(m % dimension + 1) * 6
end if
end if
end do
end subroutine create_cmfd_tally
!===============================================================================
! NEUTRON_BALANCE writes a file that contains n. bal. info for all cmfd mesh
!===============================================================================
subroutine neutron_balance()
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for g
integer :: h ! iteration counter for outgoing groups
integer :: l ! iteration counter for leakage
integer :: io_error ! error for opening file unit
real(8) :: leakage ! leakage term in neutron balance
real(8) :: interactions ! total number of interactions in balance
real(8) :: scattering ! scattering term in neutron balance
real(8) :: fission ! fission term in neutron balance
real(8) :: res ! residual of neutron balance
character(MAX_FILE_LEN) :: filename
character(30) :: label
! open cmfd file for output
filename = "cmfd.out"
open(FILE=filename, UNIT=UNIT_CMFD, STATUS='replace', ACTION='write', &
IOSTAT=io_error)
! extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! begin loop around space and energy groups
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
GROUPG: do g = 1,ng
! get leakage
leakage = 0.0
LEAK: do l = 1,3
leakage = leakage + ((cmfd % current(4*l,g,i,j,k) - &
& cmfd % current(4*l-1,g,i,j,k))) - &
& ((cmfd % current(4*l-2,g,i,j,k) - &
& cmfd % current(4*l-3,g,i,j,k)))
end do LEAK
! interactions
interactions = cmfd % totalxs(g,i,j,k) * cmfd % flux(g,i,j,k)
! get scattering and fission
scattering = 0.0
fission = 0.0
GROUPH: do h = 1,ng
scattering = scattering + cmfd % scattxs(h,g,i,j,k) * &
& cmfd % flux(h,i,j,k)
fission = fission + cmfd % nfissxs(h,g,i,j,k) * &
& cmfd % flux(h,i,j,k)
end do GROUPH
! compute residual
res = leakage + interactions - scattering - (ONE/keff)*fission
! write output
label = "MESH (" // trim(int4_to_str(i)) // ". " // &
& trim(int4_to_str(j)) // ", " // trim(int4_to_str(k)) // &
& ") GROUP " // trim(int4_to_str(g))
write(UNIT=UNIT_CMFD, FMT='(A,T35,A)') label, &
& trim(real_to_str(res))
end do GROUPG
end do XLOOP
end do YLOOP
end do ZLOOP
! close file
close(UNIT=UNIT_CMFD)
end subroutine neutron_balance
!===============================================================================
! SET_COREMAP is a routine that sets the core mapping information
!===============================================================================
subroutine set_coremap()
integer :: kount=1 ! counter for unique fuel assemblies
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
! extract spatial indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
! count how many fuel assemblies exist
cmfd % mat_dim = sum(cmfd % coremap - 1)
! begin loops over spatial indices
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
! check for reflector
if (cmfd % coremap(i,j,k) == 1) then
! reset value to 99999
cmfd % coremap(i,j,k) = 99999
else
! must be a fuel --> give unique id number
cmfd % coremap(i,j,k) = kount
kount = kount + 1
end if
end do XLOOP
end do YLOOP
end do ZLOOP
end subroutine set_coremap
!===============================================================================
! GET_REFLECTOR_ALBEDO is a function that calculates the albedo to the reflector
!===============================================================================
function get_reflector_albedo(l,g,i,j,k)
! function variable
real(8) :: get_reflector_albedo ! reflector albedo
! local variable
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: shift_idx ! parameter to shift index by +1 or -1
real(8) :: current(12) ! partial currents for all faces of mesh cell
real(8) :: albedo ! the albedo
! get partial currents from object
current = cmfd%current(:,g,i,j,k)
! define xyz and +/- indices
shift_idx = -2*mod(l,2) + 1 ! shift neig by -1 or +1
! calculate albedo
albedo = (current(2*l-1)/current(2*l))**(shift_idx)
! assign to function variable
get_reflector_albedo = albedo
end function get_reflector_albedo
!===============================================================================
! WRITE_HDF5 writes an hdf5 output file with the cmfd object for restarts
!===============================================================================
subroutine write_hdf5()
use hdf5
! character(LEN=7), parameter :: filename = "cmfd.h5" ! File name
character(LEN=4), parameter :: grpname = "cmfd" ! Group name
! integer(HID_T) :: file_id ! File identifier
integer(HID_T) :: group_id ! Group identifier
integer(HID_T) :: dataspace_id ! Data space identifier
integer(HID_T) :: dataset_id ! Dataset identifier
integer :: error ! Error flag
integer(HSIZE_T), dimension(1) :: dim1 ! vector for hdf5 dimensions
integer(HSIZE_T), dimension(3) :: dim3 ! vector for hdf5 dimensions
integer(HSIZE_T), dimension(4) :: dim4 ! vector for hdf5 dimensions
integer(HSIZE_T), dimension(5) :: dim5 ! vector for hdf5 dimensions
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! initialize FORTRAN interface.
! call h5open_f(error)
! create a new file using default properties.
! call h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, error)
! create the CMFD group
call h5gcreate_f(hdf5_output_file, grpname, group_id, error)
! write indices from cmfd object
dim1 = (/4/)
call h5screate_simple_f(1,dim1,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/indices",H5T_NATIVE_INTEGER,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%indices,dim1,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write totalxs from cmfd object
dim4 = (/ng,nx,ny,nz/)
call h5screate_simple_f(4,dim4,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/totalxs",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%totalxs,dim4,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write p1scattxs from cmfd object
dim4 = (/ng,nx,ny,nz/)
call h5screate_simple_f(4,dim4,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/p1scattxs",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%p1scattxs,dim4,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write scattxs from cmfd object
dim5 = (/ng,ng,nx,ny,nz/)
call h5screate_simple_f(5,dim5,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/scattxs",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%scattxs,dim5,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write nfissxs from cmfd object
dim5 = (/ng,ng,nx,ny,nz/)
call h5screate_simple_f(5,dim5,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/nfissxs",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%nfissxs,dim5,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write diffcof from cmfd object
dim4 = (/ng,nx,ny,nz/)
call h5screate_simple_f(4,dim4,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/diffcof",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%diffcof,dim4,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write current from cmfd object
dim5 = (/12,ng,nx,ny,nz/)
call h5screate_simple_f(5,dim5,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/current",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%current,dim5,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write flux from cmfd object
dim4 = (/ng,nx,ny,nz/)
call h5screate_simple_f(4,dim4,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/flux",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%flux,dim4,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write dtilde from cmfd object
dim5 = (/6,ng,nx,ny,nz/)
call h5screate_simple_f(5,dim5,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/dtilde",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dtilde,dim5,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write dhat from cmfd object
dim5 = (/6,ng,nx,ny,nz/)
call h5screate_simple_f(5,dim5,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/dhat",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dhat,dim5,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write albedo from cmfd object
dim1 = (/6/)
call h5screate_simple_f(1,dim1,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/albedo",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%albedo,dim1,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write hxyz from cmfd object
dim4 = (/3,nx,ny,nz/)
call h5screate_simple_f(4,dim4,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/hxyz",H5T_NATIVE_DOUBLE,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%hxyz,dim4,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write coremap from cmfd object
dim3 = (/nx,ny,nz/)
call h5screate_simple_f(3,dim3,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/coremap",H5T_NATIVE_INTEGER,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! write mat_dim from cmfd object
dim1 = (/1/)
call h5screate_simple_f(1,dim1,dataspace_id,error)
call h5dcreate_f(hdf5_output_file,"cmfd/mat_dim",H5T_NATIVE_INTEGER,dataspace_id, &
& dataset_id,error)
call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,error)
call h5sclose_f(dataspace_id,error)
call h5dclose_f(dataset_id,error)
! close the CMFD group
call h5gclose_f(group_id, error)
! terminate access to the file.
! call h5fclose_f(hdf5_output_file, error)
! close FORTRAN interface.
! call h5close_f(error)
end subroutine write_hdf5
!===============================================================================
! READ_HDF5 writes an hdf5 output file with the cmfd object for restarts
!===============================================================================
subroutine read_hdf5()
use hdf5
use hdf5_interface, only: hdf5_open_output, hdf5_close_output
! integer(HID_T) :: file_id ! File identifier
integer(HID_T) :: dataset_id ! Dataset identifier
integer :: error ! Error flag
integer(HSIZE_T), dimension(1) :: dim1
integer(HSIZE_T), dimension(3) :: dim3
integer(HSIZE_T), dimension(4) :: dim4
integer(HSIZE_T), dimension(5) :: dim5
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! open output file
call hdf5_open_output()
! read indices to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/indices",dataset_id,error)
dim1 = (/4/)
call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%indices,dim1,error)
call h5dclose_f(dataset_id,error)
! get indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! allocate cmfd object
call allocate_cmfd()
! read totalxs to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/totalxs",dataset_id,error)
dim4 = (/ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%totalxs,dim4,error)
call h5dclose_f(dataset_id,error)
! read p1scattxs to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/p1scattxs",dataset_id,error)
dim4 = (/ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%p1scattxs,dim4,error)
call h5dclose_f(dataset_id,error)
! read scattxs to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/scattxs",dataset_id,error)
dim5 = (/ng,ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%scattxs,dim5,error)
call h5dclose_f(dataset_id,error)
! read scattxs to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/nfissxs",dataset_id,error)
dim5 = (/ng,ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%nfissxs,dim5,error)
call h5dclose_f(dataset_id,error)
! read diffcof to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/diffcof",dataset_id,error)
dim4 = (/ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%diffcof,dim4,error)
call h5dclose_f(dataset_id,error)
! read current to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/current",dataset_id,error)
dim5 = (/12,ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%current,dim5,error)
call h5dclose_f(dataset_id,error)
! read flux to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/flux",dataset_id,error)
dim4 = (/ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%flux,dim4,error)
call h5dclose_f(dataset_id,error)
! read dtilde to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/dtilde",dataset_id,error)
dim5 = (/6,ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dtilde,dim5,error)
call h5dclose_f(dataset_id,error)
! read dhat to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/dhat",dataset_id,error)
dim5 = (/6,ng,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dhat,dim5,error)
call h5dclose_f(dataset_id,error)
! read albedo to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/albedo",dataset_id,error)
dim1 = (/6/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%albedo,dim1,error)
call h5dclose_f(dataset_id,error)
! read hxyz to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/hxyz",dataset_id,error)
dim4 = (/3,nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%hxyz,dim4,error)
call h5dclose_f(dataset_id,error)
! read coremap to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/coremap",dataset_id,error)
dim3 = (/nx,ny,nz/)
call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,error)
call h5dclose_f(dataset_id,error)
! read mat_dim to cmfd object
call h5dopen_f(hdf5_output_file,"cmfd/mat_dim",dataset_id,error)
dim1 = (/1/)
call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,error)
call h5dclose_f(dataset_id,error)
! close output file
call hdf5_close_output()
end subroutine read_hdf5
!===============================================================================
! WRITE_PARAVIEW_VTK outputs mesh data in vtk file for viewing
!===============================================================================
subroutine write_vtk()
use vtk_writer
integer :: i ! x loop counter
integer :: j ! y loop counter
integer :: k ! z loop counter
integer :: g ! group counter
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: n_idx ! index in eigenvector
real(8) :: x_m ! -x coordinate
real(8) :: x_p ! +x coordinate
real(8) :: y_m ! -y coordinate
real(8) :: y_p ! +y coordinate
real(8) :: z_m ! -z coordinate
real(8) :: z_p ! +z coordinate
real(8) :: x_uns(1:8) ! array of x points
real(8) :: y_uns(1:8) ! array of y points
real(8) :: z_uns(1:8) ! array of z points
type(StructuredMesh), pointer :: m => null() ! pointer to mesh
! vtk specific variables
integer :: E_IO ! error code
integer :: nn ! number of nodes
integer :: nc ! number of cells
integer :: con(8) ! connectivity vector
integer :: off(1:1) ! offset, number of nodes in cell
integer(1) :: cell_id(1:1) ! cell type
real(8) :: real_buffer(1:1) ! real data buffer 8-byte
character(len=40) :: varname ! name of output variable
character(len=3) :: str_g ! string for energy group #
! extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! point to mesh object
m => meshes(1)
! set up vtk file
E_IO = VTK_INI_XML(output_format = 'ASCII', &
& filename = 'cmfd_unst.vtu', &
& mesh_topology = 'UnstructuredGrid')
! set vtk parameters
nn = 8
nc = 1
con = (/0,1,2,3,4,5,6,7/)
off = (/8/)
cell_id = (/11/)
! begin loop to construct mesh
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
! check for non accelerated region
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == 99999) then
cycle
end if
end if
! calculate all coordinates
x_m = dble(i - 1)*m%width(1) + m%origin(1)
x_p = dble(i)*m%width(1) + m%origin(1)
y_m = dble(j - 1)*m%width(2) + m%origin(2)
y_p = dble(j)*m%width(2) + m%origin(2)
z_m = dble(k - 1)*m%width(3) + m%origin(3)
z_p = dble(k)*m%width(3) + m%origin(3)
! set up points arrays
x_uns = (/x_m,x_p,x_m,x_p,x_m,x_p,x_m,x_p/)
y_uns = (/y_m,y_m,y_p,y_p,y_m,y_m,y_p,y_p/)
z_uns = (/z_m,z_m,z_m,z_m,z_p,z_p,z_p,z_p/)
! set up geometry piece
E_IO = VTK_GEO_XML(nn,nc,x_uns,y_uns,z_uns)
! open data block in vtk file
E_IO = VTK_DAT_XML('cell','open')
! loop around energy
GROUP: do g = 1,ng
! convert group int to str
write(str_g,'(I3)') g
! write out flux
n_idx = get_matrix_idx(g,i,j,k,ng,nx,ny)
real_buffer = (/cmfd%phi(n_idx)/)
varname = 'flux_'//trim(adjustl(str_g))
E_IO = VTK_VAR_XML(nc,varname,real_buffer)
end do GROUP
! close data block in vtk file
E_IO = VTK_DAT_XML('cell','close')
! write out connectivity
E_IO = VTK_CON_XML(nc,con,off,cell_id)
! close geometry piece
E_IO = VTK_GEO_XML()
end do XLOOP
end do YLOOP
end do ZLOOP
! close vtk file
E_IO = VTK_END_XML()
end subroutine write_vtk
end module cmfd_utils