diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 821f2d3306..c751ab6399 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -15,18 +15,24 @@ ace_header.o: endf_header.o cmfd_data.o: cmfd_output.o cmfd_data.o: datatypes.o +cmfd_data.o: error.o cmfd_data.o: global.o cmfd_data.o: mesh.o +cmfd_data.o: mesh_header.o cmfd_data.o: string.o cmfd_data.o: xml-fortran/templates/cmfd_t.o cmfd_execute.o: cmfd_data.o +cmfd_output.o : cmfd_output.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_output.o: mesh_header.o +cmfd_output.o: vtk_writer.o + cmfd_power_solver.o: cmfd_loss_operator.o cmfd_power_solver.o: cmfd_prod_operator.o diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 247c4adb20..29160f5ec2 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -1,10 +1,11 @@ module cmfd_execute - use cmfd_data + use cmfd_data, only: set_up_cmfd + use cmfd_output, only: write_cmfd_vtk use cmfd_power_solver, only: cmfd_power_execute use cmfd_slepc_solver, only: cmfd_slepc_execute use cmfd_snes_solver, only: cmfd_snes_execute - use global + use global, only: cmfd implicit none @@ -38,6 +39,9 @@ contains ! finalize slepc call SlepcFinalize(ierr) + ! write vtk file + call write_cmfd_vtk() + end subroutine execute_cmfd end module cmfd_execute diff --git a/src/cmfd_output.F90 b/src/cmfd_output.F90 index e88fcbb3d3..daf15a2a59 100644 --- a/src/cmfd_output.F90 +++ b/src/cmfd_output.F90 @@ -2,7 +2,12 @@ module cmfd_output implicit none private - public :: write_cmfd_hdf5 + public :: write_cmfd_hdf5,write_cmfd_vtk + + 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 contains @@ -32,11 +37,6 @@ contains 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) @@ -177,4 +177,167 @@ contains end subroutine write_cmfd_hdf5 +!=============================================================================== +! WRITE_CMFD_VTK outputs mesh data in vtk file for viewing +!=============================================================================== + + subroutine write_cmfd_vtk() + + use global, only: cmfd,meshes + use mesh_header, only: StructuredMesh + use vtk_writer + + integer :: i ! x loop counter + integer :: j ! y loop counter + integer :: k ! z loop counter + integer :: g ! group counter + 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) + 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_cmfd_vtk + +!=============================================================================== +! GET_MATRIX_IDX takes (x,y,z,g) indices and computes location in matrix +!=============================================================================== + + function get_matrix_idx(g,i,j,k) + + use global, only: cmfd + + ! 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 + + ! 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 + end module cmfd_output