From e0397ef15ff0226e516ecacbe910ae259c54d124 Mon Sep 17 00:00:00 2001
From: Bryan Herman <bherman@mit.edu>
Date: Sat, 15 Mar 2014 08:18:07 -0400
Subject: [PATCH] added keff balance

---
 src/cmfd_data.F90   | 32 +++++++++++++++++++++++++-------
 src/cmfd_header.F90 |  3 +++
 2 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index 269f8b048..3b46e6de6 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -108,6 +108,8 @@ contains
     cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width
     cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width
 
+    cmfd % keff_bal = ZERO
+
    ! Begin loop around tallies
    TAL: do ital = 1, n_cmfd_tallies
 
@@ -216,6 +218,7 @@ contains
                   ! Bank source
                   cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + &
                        t % results(2,score_index) % sum
+                  cmfd % keff_bal = cmfd % keff_bal + t % results(2,score_index) % sum / dble(t % n_realizations)
 
                 end do INGROUP
 
@@ -707,8 +710,9 @@ contains
   subroutine compute_dhat()
 
     use constants,  only: CMFD_NOACCEL, ZERO
-    use global,     only: cmfd, cmfd_coremap, message, dhat_reset
+    use global,     only: cmfd, cmfd_coremap, message, dhat_reset, current_batch
     use output,     only: write_message
+    use string,     only: to_str
 
     integer :: nx             ! maximum number of cells in x direction
     integer :: ny             ! maximum number of cells in y direction
@@ -742,6 +746,9 @@ contains
     nxyz(1,:) = (/1,nx/)
     nxyz(2,:) = (/1,ny/)
     nxyz(3,:) = (/1,nz/)
+#ifdef CMFD_DEBUG
+    open(file='cmfd_dhat_' // trim(to_str(current_batch)) // '.dat', unit=125)
+#endif
 
     ! Geting loop over group and spatial indices
     ZLOOP:  do k = 1,nz
@@ -826,7 +833,9 @@ contains
 
               ! record dhat in cmfd object
               cmfd%dhat(l,g,i,j,k) = dhat
-
+#ifdef CMFD_DEBUG
+              write(125,*) dhat
+#endif
               ! check for dhat reset
               if (dhat_reset) then
                 cmfd%dhat(l,g,i,j,k) = ZERO
@@ -848,6 +857,10 @@ contains
       call write_message(1)
     end if
 
+#ifdef CMFD_DEBUG
+    close(unit=125)
+#endif
+
   end subroutine compute_dhat
 
 !===============================================================================
@@ -897,7 +910,7 @@ contains
   subroutine fix_neutron_balance()
 
     use constants,  only: ONE, ZERO, CMFD_NOACCEL
-    use global,     only: cmfd, keff, message
+    use global,     only: cmfd, message
     use output,     only: write_message 
 
     integer :: nx         ! number of mesh cells in x direction
@@ -932,6 +945,7 @@ contains
     real(8) :: r1         ! right hand side element 1
     real(8) :: r2         ! right hand side element 2
     real(8) :: det        ! determinant of balance matrix
+    real(8) :: keff_bal   ! keffective for balance eq.
 
     message = 'Correcting neutron balance'
     call write_message(1)
@@ -942,6 +956,10 @@ contains
     nz = cmfd % indices(3)
     ng = cmfd % indices(4)
 
+    keff_bal = cmfd % keff_bal
+
+    print *, 'KEFF BALANCE is:', keff_bal
+
     ! Return if not two groups
     if (ng /= 2) return
 
@@ -994,12 +1012,12 @@ contains
 
           ! Create matrix and solve for effective downscatter and thermal flux
           a = flux1
-          b = -ONE/keff*nsigf21
+          b = -ONE/keff_bal*nsigf21
           c = flux1
-          d = -siga2 + ONE/keff*nsigf22
+          d = -siga2 + ONE/keff_bal*nsigf22
           det = a*d - b*c
-          r1 = ONE/keff*nsigf11*flux1 - siga1*flux1 - leak1
-          r2 = leak2 - ONE/keff*nsigf12*flux1
+          r1 = ONE/keff_bal*nsigf11*flux1 - siga1*flux1 - leak1
+          r2 = leak2 - ONE/keff_bal*nsigf12*flux1
           sigs12_eff = ONE/det*(d*r1 - b*r2)
           flux2 = ONE/det*(a*r2 - c*r1)
  
diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90
index 3d7786ab8..6a2dc8ed5 100644
--- a/src/cmfd_header.F90
+++ b/src/cmfd_header.F90
@@ -83,6 +83,9 @@ module cmfd_header
     ! List of CMFD k
     real(8), allocatable :: k_cmfd(:)
 
+    ! Balance keff
+    real(8) :: keff_bal
+
   end type cmfd_type
 
 contains