diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 6555e24d79..52c6bc945b 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -12,7 +12,7 @@ contains subroutine set_up_cmfd() - use global, only: cmfd + use global, only: cmfd,cmfd_coremap use cmfd_header, only: allocate_cmfd use cmfd_output, only: neutron_balance,write_cmfd_hdf5 @@ -26,7 +26,7 @@ contains call neutron_balance() ! check for core map - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then call set_coremap() end if @@ -117,7 +117,7 @@ contains ! detect zero flux if ((flux - 0.0D0) < 1.0D-10) then - if (.not. allocated(cmfd%coremap)) then + if (.not. cmfd_coremap) then write(*,*) 'Fatal: detected zero flux without coremap' stop else @@ -238,7 +238,7 @@ contains subroutine compute_dtilde() - use global, only: cmfd + use global, only: cmfd,cmfd_coremap ! local variables integer :: nx ! maximum number of cells in x direction @@ -323,7 +323,7 @@ contains neig_hxyz = cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3)) ! check for fuel-reflector interface - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == & & 99999 .and. cmfd % coremap(i,j,k) /= 99999) then @@ -374,7 +374,7 @@ contains subroutine compute_dhat() - use global, only:cmfd + use global, only:cmfd,cmfd_coremap ! local variables integer :: nx ! maximum number of cells in x direction @@ -458,7 +458,7 @@ contains product(cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3))) ! check for fuel-reflector interface - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == & & 99999 .and. cmfd % coremap(i,j,k) /= 99999) then diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index f048134a1a..8d198fb5f6 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -51,10 +51,11 @@ contains allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), & & cmfd % indices(3))) cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3))) + cmfd_coremap = .TRUE. end if ! check for core map activation by printing note - if (allocated(cmfd % coremap)) print *,"Core Map Overlay Activated" + if (cmfd_coremap) print *,"Core Map Overlay Activated" ! create tally objects call create_cmfd_tally() @@ -301,7 +302,7 @@ contains subroutine read_cmfd_hdf5() use cmfd_header, only: allocate_cmfd - use global, only: cmfd + use global, only: cmfd,cmfd_coremap #ifdef HDF5 use global, only: hdf5_output_file,hdf5_err @@ -323,6 +324,8 @@ contains integer :: nz ! number of mesh cells in z direction integer :: ng ! number of energy groups + integer :: core_map_int + ! open output file call hdf5_open_output() @@ -407,17 +410,51 @@ contains call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%hxyz,dim4,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) - ! read coremap to cmfd object - call h5dopen_f(hdf5_output_file,"cmfd/coremap",dataset_id,hdf5_err) - dim3 = (/nx,ny,nz/) - call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,hdf5_err) + ! read in core_map logical + call h5dopen_f(hdf5_output_file,"cmfd/coremap_active",dataset_id,hdf5_err) + dim1 = (/1/) + call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,core_map_int,dim1,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) - ! read mat_dim to cmfd object - call h5dopen_f(hdf5_output_file,"cmfd/mat_dim",dataset_id,hdf5_err) - dim1 = (/1/) - call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,hdf5_err) - call h5dclose_f(dataset_id,hdf5_err) + ! now set coremap depending on logical + select case(core_map_int) + + ! core map is not active + case(0) + + ! set logical to false + cmfd_coremap = .FALSE. + + ! core map is active + case(1) + + ! set logical to true + cmfd_coremap = .TRUE. + + ! allocate coremap in cmfd obj + allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), & + & cmfd % indices(3))) + + + ! read coremap to cmfd object + call h5dopen_f(hdf5_output_file,"cmfd/coremap",dataset_id,hdf5_err) + dim3 = (/nx,ny,nz/) + call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + ! read mat_dim to cmfd object + call h5dopen_f(hdf5_output_file,"cmfd/mat_dim",dataset_id,hdf5_err) + dim1 = (/1/) + call h5dread_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + ! something is wrong + case default + + write(*,*) 'FATAL ==> Could not detect core map in hdf5 output file' + stop + + end select ! close output file call hdf5_close_output() diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index 774b19b251..e0518b0b99 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -48,7 +48,7 @@ contains subroutine get_M_indices(this) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap type(loss_operator) :: this @@ -62,7 +62,7 @@ contains this%nnz = 7 + ng - 1 ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then this%n = cmfd % mat_dim * ng else this%n = nx*ny*nz*ng @@ -76,7 +76,7 @@ contains subroutine build_loss_matrix(this) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap type(loss_operator) :: this @@ -123,7 +123,7 @@ contains XLOOP: do i = 1,nx ! check if not including reflector - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! check if at a reflector if (cmfd % coremap(i,j,k) == 99999) then @@ -170,7 +170,7 @@ contains if (bound(l) /= nxyz(xyz_idx,dir_idx)) then ! check for core map - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! check that neighbor is not reflector if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= & @@ -272,7 +272,7 @@ contains function get_matrix_idx(g,i,j,k) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap ! arguments integer :: get_matrix_idx ! the index location in matrix @@ -285,7 +285,7 @@ contains integer :: nidx ! index in matrix ! check if coremap is used - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! get idx from core map nidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) diff --git a/src/cmfd_output.F90 b/src/cmfd_output.F90 index 77055c62ee..f1136ff849 100644 --- a/src/cmfd_output.F90 +++ b/src/cmfd_output.F90 @@ -112,7 +112,7 @@ contains subroutine write_cmfd_hdf5() - use global, only: cmfd + use global, only: cmfd,cmfd_coremap #ifdef HDF5 use hdf5 @@ -144,8 +144,8 @@ contains ! write indices from cmfd object dim1 = (/4/) call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/indices",H5T_NATIVE_INTEGER,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/indices",H5T_NATIVE_INTEGER, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%indices,dim1,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -153,8 +153,8 @@ contains ! write totalxs from cmfd object dim4 = (/ng,nx,ny,nz/) call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/totalxs",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/totalxs",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%totalxs,dim4,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -162,8 +162,8 @@ contains ! write p1scattxs from cmfd object dim4 = (/ng,nx,ny,nz/) call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/p1scattxs",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/p1scattxs",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%p1scattxs,dim4,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -171,8 +171,8 @@ contains ! write scattxs from cmfd object dim5 = (/ng,ng,nx,ny,nz/) call h5screate_simple_f(5,dim5,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/scattxs",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/scattxs",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%scattxs,dim5,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -180,8 +180,8 @@ contains ! write nfissxs from cmfd object dim5 = (/ng,ng,nx,ny,nz/) call h5screate_simple_f(5,dim5,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/nfissxs",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/nfissxs",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%nfissxs,dim5,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -189,8 +189,8 @@ contains ! write diffcof from cmfd object dim4 = (/ng,nx,ny,nz/) call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/diffcof",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/diffcof",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%diffcof,dim4,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -198,8 +198,8 @@ contains ! write current from cmfd object dim5 = (/12,ng,nx,ny,nz/) call h5screate_simple_f(5,dim5,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/current",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/current",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%current,dim5,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -207,8 +207,8 @@ contains ! write flux from cmfd object dim4 = (/ng,nx,ny,nz/) call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/flux",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/flux",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%flux,dim4,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -216,8 +216,8 @@ contains ! write dtilde from cmfd object dim5 = (/6,ng,nx,ny,nz/) call h5screate_simple_f(5,dim5,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/dtilde",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/dtilde",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dtilde,dim5,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -225,8 +225,8 @@ contains ! write dhat from cmfd object dim5 = (/6,ng,nx,ny,nz/) call h5screate_simple_f(5,dim5,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/dhat",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/dhat",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%dhat,dim5,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -234,8 +234,8 @@ contains ! write albedo from cmfd object dim1 = (/6/) call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/albedo",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/albedo",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%albedo,dim1,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) @@ -243,29 +243,54 @@ contains ! write hxyz from cmfd object dim4 = (/3,nx,ny,nz/) call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/hxyz",H5T_NATIVE_DOUBLE,dataspace_id, & - & dataset_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/hxyz",H5T_NATIVE_DOUBLE, & + & dataspace_id,dataset_id,hdf5_err) call h5dwrite_f(dataset_id,H5T_NATIVE_DOUBLE,cmfd%hxyz,dim4,hdf5_err) call h5sclose_f(dataspace_id,hdf5_err) call h5dclose_f(dataset_id,hdf5_err) - ! write coremap from cmfd object - dim3 = (/nx,ny,nz/) - call h5screate_simple_f(3,dim3,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/coremap",H5T_NATIVE_INTEGER,dataspace_id, & - & dataset_id,hdf5_err) - call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,hdf5_err) - call h5sclose_f(dataspace_id,hdf5_err) - call h5dclose_f(dataset_id,hdf5_err) + ! only write the following if core map is active + if (cmfd_coremap) then - ! write mat_dim from cmfd object - dim1 = (/1/) - call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) - call h5dcreate_f(hdf5_output_file,"cmfd/mat_dim",H5T_NATIVE_INTEGER,dataspace_id, & - & dataset_id,hdf5_err) - call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,hdf5_err) - call h5sclose_f(dataspace_id,hdf5_err) - call h5dclose_f(dataset_id,hdf5_err) + ! write coremap from cmfd object + dim3 = (/nx,ny,nz/) + call h5screate_simple_f(3,dim3,dataspace_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/coremap",H5T_NATIVE_INTEGER, & + & dataspace_id,dataset_id,hdf5_err) + call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%coremap,dim3,hdf5_err) + call h5sclose_f(dataspace_id,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + ! write mat_dim from cmfd object + dim1 = (/1/) + call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/mat_dim",H5T_NATIVE_INTEGER, & + & dataspace_id,dataset_id,hdf5_err) + call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,cmfd%mat_dim,dim1,hdf5_err) + call h5sclose_f(dataspace_id,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + ! write core map logical to true (1) + dim1 = (/1/) + call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/coremap_active", & + & H5T_NATIVE_INTEGER,dataspace_id,dataset_id,hdf5_err) + call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,1,dim1,hdf5_err) + call h5sclose_f(dataspace_id,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + else + + ! write core map logical to false (0) + dim1 = (/1/) + call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err) + call h5dcreate_f(hdf5_output_file,"cmfd/coremap_active", & + & H5T_NATIVE_INTEGER,dataspace_id,dataset_id,hdf5_err) + call h5dwrite_f(dataset_id,H5T_NATIVE_INTEGER,0,dim1,hdf5_err) + call h5sclose_f(dataspace_id,hdf5_err) + call h5dclose_f(dataset_id,hdf5_err) + + end if ! close the CMFD group call h5gclose_f(group_id, hdf5_err) @@ -278,7 +303,7 @@ contains subroutine write_cmfd_vtk() - use global, only: cmfd,meshes + use global, only: cmfd,cmfd_coremap,meshes use mesh_header, only: StructuredMesh use vtk_writer @@ -339,7 +364,7 @@ contains XLOOP: do i = 1,nx ! check for non accelerated region - if (allocated(cmfd%coremap)) then + if (cmfd_coremap) then if (cmfd%coremap(i,j,k) == 99999) then cycle end if @@ -405,7 +430,7 @@ contains function get_matrix_idx(g,i,j,k) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap ! arguments integer :: get_matrix_idx ! the index location in matrix @@ -418,7 +443,7 @@ contains integer :: nidx ! index in matrix ! check if coremap is used - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! get idx from core map nidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) diff --git a/src/cmfd_prod_operator.F90 b/src/cmfd_prod_operator.F90 index 77a594be9a..d75caa67d2 100644 --- a/src/cmfd_prod_operator.F90 +++ b/src/cmfd_prod_operator.F90 @@ -48,7 +48,7 @@ contains subroutine get_F_indices(this) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap type(prod_operator) :: this @@ -62,7 +62,7 @@ contains this%nnz = 7 + ng - 1 ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then this%n = cmfd % mat_dim * ng else this%n = nx*ny*nz*ng @@ -76,7 +76,7 @@ contains subroutine build_prod_matrix(this) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap type(prod_operator) :: this @@ -103,7 +103,7 @@ contains XLOOP: do i = 1,nx ! check if not including reflector - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! check if at a reflector if (cmfd % coremap(i,j,k) == 99999) then @@ -152,7 +152,7 @@ contains function get_matrix_idx(g,i,j,k) - use global, only: cmfd + use global, only: cmfd,cmfd_coremap ! arguments integer :: get_matrix_idx ! the index location in matrix @@ -165,7 +165,7 @@ contains integer :: nidx ! index in matrix ! check if coremap is used - if (allocated(cmfd % coremap)) then + if (cmfd_coremap) then ! get idx from core map nidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) diff --git a/src/global.F90 b/src/global.F90 index 7224dfe428..0f0661eecc 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -211,6 +211,9 @@ module global ! Flag for CMFD only logical :: cmfd_only = .FALSE. + ! Flag for coremap accelerator + logical :: cmfd_coremap = .FALSE. + contains !===============================================================================