Merge branch 'develop' of github.com:pyne/pyne into develop

docs/theorymanual/variance_reduction.rst

@@ -4,7 +4,7 @@
 Variance Reduction
 =============================
 
-:Author: Elliott Biondo
+:Author: Elliott Biondo, Kalin Kiesling
 
 .. currentmodule::variancereduction
 
@@ -175,11 +175,118 @@ The expected results are:
 
 Notice that the value in the :math:`ww` vector that is a division by 0 has been replaced with 0.
 
+************
+MAGIC Method
+************
+
+The Method of Automatic Generation of Importances by Calculation (MAGIC) is 
+a global variance reduction technique in which an initial particle distribution, 
+in the form of fluxes, populations, or weights is obtained and then used to 
+generate mesh-based weight windows or importances. This method recognizes 
+the initial particle distribution will be poor in some highly attenuated regions
+but upon iteration of the MAGIC method, the solution will improve. Below are the
+steps for the MAGIC method. [2]
+
+1. Run MCNP, in analogue mode, to set up a flux meshtally. Multigroup cross 
+   section data and a high energy cut-off, corresponding to a mean-free path no 
+   greater than the mesh voxel size, should be used.
+
+2. Process the resulting meshtally data by normalizing the fluz to have a value 
+   of 0.5 in the source (or highest) region. Use the normilized flux to create
+   a new weight window file to be used for MCNP.
+
+3. Modify the original MCNP input to use the generated weight window file and 
+   run again.
+
+4. If results are are sufficient, no further iterations are necessary. Else, 
+   repeat starting from step 2 until desired flux results are obtained.
+
+5. If a high energy cut-off was used, reduce the cut-off energy and repeat 
+   iterations until the particle distribution is acceptable. The final iteration
+   should be performed with the appropriate energy cut-off and cross section data.
+
+...................
+PyNE implementation
+...................
+
+The implementation of MAGIC in PyNE uses a PyNE meshtally object, which is the 
+result of a meshtal file processed by PyNE's mcnp.Meshtally. Using the results 
+of the meshtal file and a specified tolerance for relative error :math:`t` and 
+null value :math:`\phi_o`, the flux will be normalized for each energy bin and 
+then be used to generate a wwinp file to be used in a subsequent iteration. The 
+steps are as follows:
+
+1. Read meshtally and determine maximum flux :math:`\phi_m^k` for each enery bin :math:`k`.
+
+2. Normalize flux :math:`\phi_i^k` in every mesh voxel :math:`i` for each energy 
+   bin :math:`k` according to :math:`\phi_m^k` to obtain a new :math:`\phi_i^{'k}`. 
+   If the relative error :math:`e_i^k` 
+   for voxel :math:`i` in energy bin :math:`k` is larger than the tolerance 
+   value :math:`t`, then set flux :math:`\phi_i^{'k}` to the null value :math:`\phi_o` 
+   instead.
+
+
+ .. math::  
+        
+        \text{If } e_i^k < t \text{ then, } \phi_i^{'k} = \frac{\phi_i^{k}}{2 \, \phi_m^k}
+        
+        \text{If } e_i^k > t \text{ then, } \phi_i^{'k} = \phi_o
+        
+        
+3. Use new flux values to create a weight window tag on the provide meshtally 
+   and use PyNE's Wwinp class to create a weight window mesh.
+
+...................
+Sample Calculations
+...................
+
+In this section, the expected results of the test_variancereduction.py unit test
+"test_magic_multi_bins" are shown. In this test, a 3D 2x2 mesh is given. Each
+voxel contains flux data corresponding to 2 energy bins. The mesh is described by
+the following flux and relative error data.
+
++--------------------------+----------------------+--------------------------+-----------------------+
+| :math:`\phi_1^{1} = 1.2` | :math:`e_1^1 = 0.11` | :math:`\phi_1^{2} = 3.3` | :math:`e_1^2 = 0.013` |
++--------------------------+----------------------+--------------------------+-----------------------+
+| :math:`\phi_2^{1} = 1.6` | :math:`e_2^1 = 0.14` | :math:`\phi_2^{2} = 1.7` | :math:`e_2^2 = 0.19`  |
++--------------------------+----------------------+--------------------------+-----------------------+
+| :math:`\phi_3^{1} = 1.5` | :math:`e_3^1 = 0.02` | :math:`\phi_3^{2} = 1.4` | :math:`e_3^2 = 0.16`  |
++--------------------------+----------------------+--------------------------+-----------------------+
+| :math:`\phi_4^{1} = 2.6` | :math:`e_4^1 = 0.04` | :math:`\phi_4^{2} = 1.0` | :math:`e_4^2 = 0.09`  |
++--------------------------+----------------------+--------------------------+-----------------------+
+
+First, the maximum flux for each energy bin is found. In this case the maximum 
+for energy bin :math:`k = 1` occurs in voxel 4 :math:`\phi_4^{1} = 2.6` and in 
+voxel 1 :math:`\phi_1^{2} = 3.3` for energy bin :math:`k = 2`. In the first 
+energy bin, the flux values are normalized by :math:`\phi_m^1 = 2.6` and in the 
+second :math:`\phi_m^2 = 3.3`. If the error tolerance is set :math:`t = 0.15` and the
+null value set to :math:`\phi_o = 0.001`, then voxels 
+2 and 3 in the second energy bin have errors larger than the tolerance and are 
+therefore set to the null value while everything else is normalized. The following
+is the result.
+
++------------------------------------------------------------------------------------+------------------------------------------------------------------------------------+
+| :math:`\phi_1^{'1} = \frac{\phi_1^1}{2 \, \phi_m^1} = \frac{1.2}{2*2.6} = 0.23077` | :math:`\phi_1^{'2} = \frac{\phi_1^2}{2 \, \phi_m^2} = \frac{3.3}{2*3.3} = 0.5`     |
++------------------------------------------------------------------------------------+------------------------------------------------------------------------------------+
+| :math:`\phi_2^{'1} = \frac{\phi_2^1}{2 \, \phi_m^1} = \frac{1.6}{2*2.6} = 0.30769` | :math:`\phi_2^{'2} = \phi_o = 0.001`                                               |
++------------------------------------------------------------------------------------+------------------------------------------------------------------------------------+
+| :math:`\phi_3^{'1} = \frac{\phi_3^1}{2 \, \phi_m^1} = \frac{1.5}{2*2.6} = 0.28846` | :math:`\phi_3^{'2} = \phi_o = 0.001`                                               |
++------------------------------------------------------------------------------------+------------------------------------------------------------------------------------+
+| :math:`\phi_4^{'1} = \frac{\phi_4^1}{2 \, \phi_m^1} = \frac{2.6}{2*2.6} = 0.5`     | :math:`\phi_4^{'2} = \frac{\phi_4^2}{2 \, \phi_m^2} = \frac{1.0}{2*3.3} = 0.12122` |
++------------------------------------------------------------------------------------+------------------------------------------------------------------------------------+
+
+The values :math:`\phi_i^{'k}` are then set as the new weight window values.
+
+
 **********
 References
 **********
 
 [1] Haghighat, A. and Wagner, J. C., "Monte Carlo Variance Reduction with 
     Deterministic Importance Functions," Progress in Nuclear Energy, Vol. 42,
     No. 1, pp. 25-53, 2003.
+
+[2] Davis, A. and Turner, A., "Comparison of global variance reduction 
+    techniques for Monte Carlo radiation transport simulations of ITER," Fusion 
+    Engineering and Design, Vol. 86, Issues 9-11, pp. 2698-2700, 2011.
 

pyne/mcnp.py

@@ -2032,7 +2032,7 @@ def __init__(self, f, tally_number, tag_names=None, mesh_has_mats=False):
         ----------
         f : filestream
             Open to the neutron/photon line.
-        tally number : int
+        tally_number : int
             The MCNP fmesh4 tally number (e.g. 4, 14, 24).
         tag_names : iterable, optional
             Four strs that specify the tag names for the results, relative

pyne/variancereduction.py

@@ -19,7 +19,9 @@
 
 from mesh import Mesh
 from mesh import MeshError
-
+from mesh import IMeshTag
+from mcnp import Wwinp
+from pyne.particle import mcnp
 
 def cadis(adj_flux_mesh, adj_flux_tag, q_mesh, q_tag,
           ww_mesh, ww_tag, q_bias_mesh, q_bias_tag, beta=5):
@@ -126,3 +128,102 @@ def cadis(adj_flux_mesh, adj_flux_tag, q_mesh, q_tag,
         tag_ww[ww_ve] = [R/(adj_flux[i]*(beta + 1.)/2.)
                          if adj_flux[i] != 0.0 else 0.0
                          for i in range(num_e_groups)]
+
+
+def magic(meshtally, tag_name, tag_name_error, **kwargs):
+    """This function reads a PyNE mcnp.MeshTally object and preforms the MAGIC 
+    algorithm and returns the resulting weight window mesh.
+
+    Parameters:
+    -----------
+    meshtally : a single PyNE mcnp.MeshTally obj
+    tag_name : string
+        The meshtally tag_name (example: n_result or n_total_result).
+    tag_name_error : string
+        The meshtally tag_name for the error associated with provided tag_name.
+        Example: n_rel_error
+    tolerance : float, optional
+        The maximum relative error allowable for the MAGIC algorithm to create 
+        a weight window lower bound for for a given mesh volume element for the
+        intial weight window lower bound generation, or overwrite preexisting 
+        weight window lower bounds for subsequent iterations. 
+    null_value : float, optional
+        The weight window lower bound value that is assigned to mesh volume
+        elements where the relative error on flux exceeds the tolerance.
+    """
+
+    tolerance = kwargs.get('tolerance',0.5)
+    null_value = kwargs.get('null_value',0.0)
+
+    # Convert particle name to the recognized abbreviation
+    particle = (meshtally.particle.capitalize())
+    if  particle == ("Neutron" or "Photon" or "Electron"):
+        meshtally.particle = mcnp(particle).lower()
+
+    # Create tags for values and errors
+    meshtally.vals = IMeshTag(mesh=meshtally, name=tag_name)
+    meshtally.errors = IMeshTag(mesh=meshtally, name=tag_name_error)
+
+    # Create weight window tags
+    tag_size = meshtally.vals[0].size
+    meshtally.ww_x = IMeshTag(tag_size, float, 
+                              name="ww_{0}".format(meshtally.particle))
+    root_tag = meshtally.mesh.createTag(
+                        "{0}_e_upper_bounds".format(meshtally.particle), 
+                        tag_size, float)
+
+    # Determine if total energy or single energy bin or multiple energy bins
+    if tag_size == 1 and len(meshtally.e_bounds) > 1:
+        total = True
+    elif tag_size == 1 and len(meshtally.e_bounds) == 1:
+        total = True
+    elif tag_size > 1 and len(meshtally.e_bounds) > 1:
+        total = False
+
+    # Reassign arrays for total and not total case
+    if total:
+        root_tag[meshtally.mesh.rootSet] = np.max(meshtally.e_bounds[:])
+        max_val = np.max(meshtally.vals[:])
+
+        vals = []
+        errors = []
+        for ve, flux in enumerate(meshtally.vals[:]):
+            vals.append(np.atleast_1d(flux))
+            errors.append(np.atleast_1d(meshtally.errors[ve]))
+        vals = np.array(vals)
+        errors = np.array(errors)
+
+    else:
+        root_tag[meshtally.mesh.rootSet] = meshtally.e_bounds[1:]
+        vals = meshtally.vals[:]
+        errors = meshtally.errors[:]
+
+    # Determine the max values for each energy bin
+    max_val = []
+    for i in range(tag_size):
+        vals_in_e = []
+        for ve, flux in enumerate(vals[:]):
+            vals_in_e.append(vals[ve][i])
+        max_val.append(np.max(vals_in_e))
+
+    # Apply normalization to create weight windows
+    ww = []
+    for ve, flux_list in enumerate(vals[:]):
+        tally_list = errors[ve]
+        flux = []
+        for i, value in enumerate(flux_list):
+            if tally_list[i] > tolerance:
+                flux.append(null_value)
+            else:    
+                flux.append(value/(2.0*max_val[i]))
+        ww.append(flux)
+
+    # Resassign weight windows to meshtally
+    if total:
+        meshtally.ww_x[:] = np.reshape(ww, len(ww))
+    else:
+        meshtally.ww_x[:] = ww
+
+    # Create wwinp mesh
+    wwinp = Wwinp()
+    wwinp.read_mesh(meshtally.mesh)

tests/test_variancereduction.py

@@ -13,9 +13,11 @@
 
 from pyne.utils import QAWarning
 warnings.simplefilter("ignore", QAWarning)
-from pyne.variancereduction import cadis
+from pyne.variancereduction import cadis, magic
+
 from pyne.mesh import Mesh, IMeshTag
 from pyne.mesh import MeshError
+from pyne import mcnp
 
 def test_cadis_single_e():
     """Test single energy group cadis case"""
@@ -100,3 +102,111 @@ def test_cadis_multiple_e():
 
     assert_array_almost_equal(ww_mesh.ww[:], expected_ww[:])
     assert_array_almost_equal(q_bias_mesh.q_bias[:], expected_q_bias[:])
+
+
+def test_magic_below_tolerance():
+    """Test MAGIC case when all flux errors are below the default tolerance"""
+
+    # create mesh
+    coords = [[0, 1, 2], [-1, 3, 4], [10, 12]]
+    flux_data = [1.2, 3.3, 1.6, 1.7]
+    flux_error = [0.11, 0.013, 0.14, 0.19]
+    tally = Mesh(structured=True, structured_coords=coords)
+
+    tally.particle = "neutron"
+    tally.e_bounds = [0.0, 0.5, 1.0]
+    tally.n_total_flux = IMeshTag(1, float)
+    tally.n_total_flux[:] = flux_data
+
+    tally.n_rel_error = IMeshTag(1, float)
+    tally.n_rel_error[:] = flux_error
+
+    magic(tally, "n_total_flux", "n_rel_error")
+
+    expected_ww = [0.181818182, 0.5, 0.2424242, 0.2575757576]
+
+    assert_array_almost_equal(tally.ww_x[:], expected_ww[:])
+
+
+def test_magic_multi_bins():
+    """Test multiple energy bins MAGIC case"""
+
+    # create a basic meshtally
+    coords = [[0, 1, 2], [-1, 3, 4], [10, 12]]
+    flux_data = [[1.2, 3.3], [1.6, 1.7], [1.5, 1.4], [2.6, 1.0]]
+    flux_error = [[0.11, 0.013], [0.14, 0.19], [0.02, 0.16], [0.04, 0.09]]
+    tally = Mesh(structured=True, structured_coords=coords)
+
+    tally.particle = "neutron"
+    tally.e_bounds = [0.0, 0.5, 1.0]
+    tally.n_flux = IMeshTag(2, float)
+    tally.n_flux[:] = flux_data
+
+    tally.n_rel_error = IMeshTag(2, float)
+    tally.n_rel_error[:] = flux_error
+
+    tolerance = 0.15
+    null_value = 0.001
+
+    magic(tally, "n_flux", "n_rel_error", tolerance=tolerance, null_value=null_value)
+
+    expected_ww = [[0.2307692308, 0.5], 
+                   [0.3076923077, 0.001], 
+                   [0.2884615385, 0.001], 
+                   [0.5, 0.15151515]]
+
+    assert_array_almost_equal(tally.ww_x[:], expected_ww[:])
+
+
+def test_magic_e_total():
+    """Test total energy group MAGIC case"""
+
+    # create mesh
+    coords = [[0, 1, 2], [-1, 3, 4], [10, 12]]
+    flux_data = [1.2, 3.3, 1.6, 1.7]
+    flux_error = [0.11, 0.013, 0.14, 0.19]
+    tally = Mesh(structured=True, structured_coords=coords)
+
+    tally.particle = "neutron"
+    tally.e_bounds = [0.0, 0.5, 1.0]
+    tally.n_total_flux = IMeshTag(1, float)
+    tally.n_total_flux[:] = flux_data
+
+    tally.n_rel_error = IMeshTag(1, float)
+    tally.n_rel_error[:] = flux_error
+
+    tolerance = 0.15
+    null_value = 0.001
+
+    magic(tally, "n_total_flux", "n_rel_error", tolerance=tolerance, null_value=null_value)
+
+    expected_ww = [0.181818182, 0.5, 0.2424242, 0.001]
+
+    assert_array_almost_equal(tally.ww_x[:], expected_ww[:])
+
+
+def test_magic_single_e():
+    """Test a single energy group MAGIC case"""
+
+    # create mesh
+    coords = [[0, 1, 2], [-1, 3, 4], [10, 12]]
+    flux_data = [1.2, 3.3, 1.6, 1.7]
+    flux_error = [0.11, 0.013, 0.14, 0.19]
+    tally = Mesh(structured=True, structured_coords=coords)
+
+    tally.particle = "neutron"
+    tally.e_bounds = [0.0, 1.0]
+    tally.n_flux = IMeshTag(1, float)
+    tally.n_flux[:] = flux_data
+
+    tally.n_rel_error = IMeshTag(1, float)
+    tally.n_rel_error[:] = flux_error
+
+    tolerance = 0.15
+    null_value = 0.001
+
+    magic(tally, "n_flux", "n_rel_error", tolerance=tolerance, null_value=null_value)
+
+    expected_ww = [0.181818182, 0.5, 0.2424242, 0.001]
+
+    assert_array_almost_equal(tally.ww_x[:], expected_ww[:])