From 057f09444aa49f25f0f301c15fc77604094ba327 Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Thu, 22 Mar 2012 13:58:48 -0400
Subject: [PATCH] Added estimator for diffusion coefficient.

---
 src/constants.F90 | 19 ++++++++++---------
 src/input_xml.F90 |  5 +++++
 src/tally.F90     | 15 +++++++++++++++
 3 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/src/constants.F90 b/src/constants.F90
index be0f712fc7..9cb5c99d55 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -243,7 +243,7 @@ module constants
        EVENT_FISSION =  3 
 
   ! Tally score type
-  integer, parameter :: N_SCORE_TYPES = 15
+  integer, parameter :: N_SCORE_TYPES = 16
   integer, parameter :: &
        SCORE_FLUX       = -1,  & ! flux
        SCORE_TOTAL      = -2,  & ! total reaction rate
@@ -252,14 +252,15 @@ module constants
        SCORE_SCATTER_1  = -5,  & ! first scattering moment
        SCORE_SCATTER_2  = -6,  & ! second scattering moment
        SCORE_SCATTER_3  = -7,  & ! third scattering moment
-       SCORE_N_1N       = -8,  & ! (n,1n) rate
-       SCORE_N_2N       = -9,  & ! (n,2n) rate
-       SCORE_N_3N       = -10, & ! (n,3n) rate
-       SCORE_N_4N       = -11, & ! (n,4n) rate
-       SCORE_ABSORPTION = -12, & ! absorption rate
-       SCORE_FISSION    = -13, & ! fission rate
-       SCORE_NU_FISSION = -14, & ! neutron production rate
-       SCORE_CURRENT    = -15    ! partial current
+       SCORE_DIFFUSION  = -8,  & ! diffusion coefficient
+       SCORE_N_1N       = -9,  & ! (n,1n) rate
+       SCORE_N_2N       = -10, & ! (n,2n) rate
+       SCORE_N_3N       = -11, & ! (n,3n) rate
+       SCORE_N_4N       = -12, & ! (n,4n) rate
+       SCORE_ABSORPTION = -13, & ! absorption rate
+       SCORE_FISSION    = -14, & ! fission rate
+       SCORE_NU_FISSION = -15, & ! neutron production rate
+       SCORE_CURRENT    = -16    ! partial current
 
   ! Tally map bin finding
   integer, parameter :: NO_BIN_FOUND = -1
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index b5e29a5331..bb9841a14e 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -1006,6 +1006,11 @@ contains
              case ('scatter-3')
                 t % score_bins(j) % scalar = SCORE_SCATTER_3
 
+                ! Set tally estimator to analog
+                t % estimator = ESTIMATOR_ANALOG
+             case ('diffusion')
+                t % score_bins(j) % scalar = SCORE_DIFFUSION
+
                 ! Set tally estimator to analog
                 t % estimator = ESTIMATOR_ANALOG
              case ('n1n')
diff --git a/src/tally.F90 b/src/tally.F90
index 4702ed9baf..db57a56802 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -345,6 +345,20 @@ contains
 
              score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu)
 
+          case (SCORE_DIFFUSION)
+             ! Skip any event where the particle didn't scatter
+             if (p % event /= EVENT_SCATTER) cycle
+             
+             ! Temporarily store the scattering cross section
+             score = material_xs % total - material_xs % absorption
+
+             ! Since this only gets tallied at every scattering event, the flux
+             ! estimator is 1/Sigma_s. Therefore, the diffusion coefficient
+             ! times flux is 1/(3*Sigma_s*(Sigma_t - mu*Sigma_s)).
+
+             score = last_wgt / (3.0_8 * score * (material_xs % total - &
+                  mu * score))
+
           case (SCORE_N_1N)
              ! Skip any event where the particle didn't scatter
              if (p % event /= EVENT_SCATTER) cycle
@@ -1186,6 +1200,7 @@ contains
     score_name(abs(SCORE_SCATTER_1))  = "First Scattering Moment"
     score_name(abs(SCORE_SCATTER_2))  = "Second Scattering Moment"
     score_name(abs(SCORE_SCATTER_3))  = "Third Scattering Moment"
+    score_name(abs(SCORE_DIFFUSION))  = "Diffusion Coefficient"
     score_name(abs(SCORE_N_1N))       = "(n,1n) Rate"
     score_name(abs(SCORE_N_2N))       = "(n,2n) Rate"
     score_name(abs(SCORE_N_3N))       = "(n,3n) Rate"