diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 68375f03ba..baef724227 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -33,42 +33,49 @@ contains ! initialize mpi communicator call petsc_init_mpi() + ! filter procs if (rank < 6) then - ! initialize slepc/petsc (communicates to world) - if(current_cycle == n_inactive + 1) call SlepcInitialize & - & (PETSC_NULL_CHARACTER,ierr) + ! initialize slepc/petsc (communicates to world) + if(current_cycle == n_inactive + 1) call SlepcInitialize & + & (PETSC_NULL_CHARACTER,ierr) - ! set global variable for number of procs in cmfd calc - call MPI_COMM_SIZE(PETSC_COMM_WORLD,n_procs_cmfd,ierr) + ! set global variable for number of procs in cmfd calc + call MPI_COMM_SIZE(PETSC_COMM_WORLD,n_procs_cmfd,ierr) - ! only run if master process - if (master) then + ! only run if master process + if (master) then - ! begin timer - call timer_start(time_cmfd) + ! begin timer + call timer_start(time_cmfd) - ! set up cmfd - if(.not. cmfd_only) call set_up_cmfd() + ! set up cmfd + if(.not. cmfd_only) call set_up_cmfd() - end if + end if - ! broadcast cmfd object to all procs - call cmfd_bcast() + ! broadcast cmfd object to all procs + call cmfd_bcast() - ! execute snes solver - call cmfd_snes_execute() + ! execute snes solver + call cmfd_snes_execute() + + ! only run if master process + if (master) then + + ! stop timer + call timer_stop(time_cmfd) - ! stop timer - call timer_stop(time_cmfd) + ! write vtk file + !if(.not. cmfd_only) call write_cmfd_vtk() - ! write vtk file - !if(.not. cmfd_only) call write_cmfd_vtk() + ! compute fission source for reweighting + call calc_fission_source() - call calc_fission_source() + end if - ! finalize slepc - if (current_cycle == n_cycles) call SlepcFinalize(ierr) + ! finalize slepc + if (current_cycle == n_cycles) call SlepcFinalize(ierr) end if @@ -170,7 +177,7 @@ contains subroutine calc_fission_source() - use global, only: cmfd,cmfd_coremap + use global, only: cmfd,cmfd_coremap,entropy_on ! local variables integer :: nx ! maximum number of cells in x direction @@ -184,6 +191,7 @@ contains integer :: idx ! index in vector real(8) :: hxyz(3) ! cell dimensions of current ijk cell real(8) :: vol ! volume of cell + real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy ! get maximum of spatial and group indices nx = cmfd%indices(1) @@ -230,6 +238,28 @@ contains ! normalize source such that it sums to 1.0 cmfd%source = cmfd%source/sum(cmfd%source) + ! compute entropy + if (entropy_on) then + + ! allocate tmp array + if (.not.allocated(source)) allocate(source(ng,nx,ny,nz)) + + ! initialize the source + source = 0.0_8 + + ! compute log + where (cmfd%source > 0.0_8) + source = cmfd%source*log(cmfd%source)/log(2.0) + end where + + ! sum that source + cmfd%entropy = -1.0_8*sum(source) + + ! deallocate tmp array + if (allocated(source)) deallocate(source) + + end if + end subroutine calc_fission_source !=============================================================================== diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index d6965af23f..8298bcae9f 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -61,7 +61,7 @@ module cmfd_header real(8), allocatable :: resnb(:,:,:,:) ! entropy calculation - real(8) :: entropy + real(8) :: entropy = 0.0_8 end type cmfd_obj diff --git a/src/intercycle.F90 b/src/intercycle.F90 index f1301920b9..2861f0ab76 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -481,7 +481,7 @@ contains if (entropy_on) then if (cmfd_on) then write(UNIT=OUTPUT_UNIT, FMT=104) current_cycle, k_cycle, & - entropy, keff, keff_std, cmfd%keff + entropy, keff, keff_std, cmfd%keff, cmfd%entropy else write(UNIT=OUTPUT_UNIT, FMT=103) current_cycle, k_cycle, & entropy, keff, keff_std @@ -522,7 +522,7 @@ contains 101 format (2X,I5,2X,F8.5,5X,F8.5," +/-",F8.5) 102 format (2X,I5,2X,F8.5,3X,F8.5) 103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5) -104 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5,2X,F8.5) +104 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5,2X,F8.5,4X,F8.5) 105 format (2X,I5,2X,F8.5,5X,F8.5," +/-",F8.5,2X,F8.5) end subroutine calculate_keff diff --git a/src/main.F90 b/src/main.F90 index 31f33c5404..f2fa99055b 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -56,9 +56,9 @@ contains ! Display column titles if (entropy_on) then if (cmfd_on) then - message = " Cycle k(cycle) Entropy Average k CMFD k" + message = " Cycle k(cycle) Entropy Average k CMFD k CMFD Ent." call write_message(1) - message = " ===== ======== ======= =================== ==========" + message = " ===== ======== ======= =================== ========== =========" call write_message(1) else message = " Cycle k(cycle) Entropy Average k"