irespr + MF32 (#296)

Co-authored-by: Luca Fiorito <lfiorito@sckcen.be>

sandy/core/endf6.py

@@ -6,6 +6,7 @@
 """
 import io
 import os
+import shutil
 from os.path import dirname, join
 from functools import reduce
 from tempfile import TemporaryDirectory
@@ -1486,6 +1487,55 @@ def inner(
             return method(self, groupr_kws=groupr_kws, **kwargs)
         return inner
 
+    def _handle_mf32_alone(method):
+        """
+        Decorator to handle files with section MF32 without section MF33.
+
+        Examples
+        --------
+        91-Pa-231 in JEFF-3.3 has MF32 but not MF33.
+        >>> tape = sandy.get_endf6_file("jeff_33", "xs", 912310)
+        >>> err = tape.get_errorr(chi=False, nubar=True, mubar=False, err=1, xs=True, errorr33_kws=dict(irespr=0))
+        >>> assert "errorr33" in err
+
+        91-Pa-231 in JEFF-3.3 has MF32 but not MF33.
+        >>> tape = sandy.get_endf6_file("jeff_33", "xs", 912330)
+        >>> err = tape.get_errorr(chi=False, nubar=True, mubar=False, err=1, xs=True, errorr33_kws=dict(irespr=0))
+        >>> assert "errorr33" in err
+
+        17-Cl-37 in JEFF-3.3 has MF32 but not MF33.
+        >>> tape = sandy.get_endf6_file("jeff_33", "xs", 170370)
+        >>> err = tape.get_errorr(chi=False, nubar=True, mubar=False, err=1, xs=True, errorr33_kws=dict(irespr=0))
+        >>> assert "errorr33" in err
+
+        95-Am-241 in JEFF-3.3 has MF32 but not MF33.
+        >>> tape = sandy.get_endf6_file("jeff_33", "xs", 952410)
+        >>> err = tape.get_errorr(chi=False, nubar=True, mubar=False, err=1, xs=True, errorr33_kws=dict(irespr=0))
+        >>> assert "errorr33" in err
+        """
+        def inner(
+                self,
+                **kwargs,
+                ):
+            """
+            Parameters
+            ----------
+            kwargs : `dict`, optional
+                 keyword arguments.
+            """
+            # input taken from
+            # https://www-nds.iaea.org/index-meeting-crp/TM_NDP/docs/OCabellos_2017.pdf
+            inst = self
+            if 32 in self.mf and 33 not in self.mf:
+                inp = sandy.njoy._input_mf32_nomf33 if 18 in self.mt else sandy.njoy._input_mf32_nomf33_no18
+                with TemporaryDirectory() as tmpdir:
+                    file = os.path.join(tmpdir, "tape20")
+                    self.to_file(file)
+                    outs = sandy.njoy._run_njoy(inp, file)
+                    inst = sandy.Endf6.from_text(outs["errorr33"])
+            return method(inst, **kwargs)
+        return inner
+
     def read_section(self, mat, mf, mt, raise_error=True):
         """
         Parse MAT/MF/MT section.
@@ -1802,6 +1852,7 @@ def get_gendf(self, **kwargs,):
 
     @_handle_njoy_inputs
     @_handle_groupr_inputs
+    @_handle_mf32_alone
     def get_errorr(self,
                    nubar=True,
                    mubar=True,
@@ -2054,14 +2105,17 @@ def get_perturbations(self, nsmp, njoy_kws={}, smp_kws={}, **kwargs,):
         njoy_kws["chi"] = False
         outs = self.get_errorr(**njoy_kws)
         filename = "PERT_{}_MF{}.xlsx"
+        filenamecov = "COV_{}_MF{}.tape"
 
         # -- Extract samples from MF31 covariance data
         if "errorr31" in outs:
+            outs["errorr31"].to_file(filenamecov.format(self.get_id(), 31))
             xls = filename.format(self.get_id(), 31)
             smp[31] = outs["errorr31"].get_cov().sampling(nsmp, to_excel=xls, seed=smp_kws.get("seed31"), **smp_kws)
 
         # -- Extract samples from MF33 covariance data
         if "errorr33" in outs:
+            outs["errorr33"].to_file(filenamecov.format(self.get_id(), 33))
             xls = filename.format(self.get_id(), 33)
             smp[33] = outs["errorr33"].get_cov().sampling(nsmp, to_excel=xls, seed=smp_kws.get("seed33"), **smp_kws)
 

sandy/fy.py

@@ -287,8 +287,8 @@ def get_mass_yield(self, zam, e):
         --------
         >>> tape_nfpy = sandy.get_endf6_file("jeff_33",'nfpy', 922350)
         >>> nfpy = Fy.from_endf6(tape_nfpy)
-        >>> nfpy.get_mass_yield(922350, 0.0253).loc[148]
-        0.0169029147
+        >>> out = nfpy.get_mass_yield(922350, 0.0253).loc[148]
+        >>> np.testing.assert_almost_equal(out, 0.0169029147)
         """
         # Filter FY data:
         conditions = {'ZAM': zam, "E": e, 'MT': 454}

sandy/njoy.py

@@ -106,6 +106,27 @@
 NJOY_SIG0 = [1e10]
 NJOY_THERMR_EMAX = 10
 
+# input taken from
+# https://www-nds.iaea.org/index-meeting-crp/TM_NDP/docs/OCabellos_2017.pdf
+_input_mf32_nomf33 = """errorr
+999/
+20 33/
+1/
+2/
+18/
+102/
+0/
+stop"""
+# same but no fissile
+_input_mf32_nomf33_no18 = """errorr
+999/
+20 33/
+1/
+2/
+102/
+0/
+stop"""
+
 
 def get_njoy():
     """
@@ -576,7 +597,7 @@ def _acer_input(endfin, pendfin, aceout, dirout, mat,
 def _errorr_input(endfin, pendfin, gendfin, errorrout, mat,
                   ign=2, ek=None, spectrum=None,
                   iwt=2, relative=True,
-                  mt=None, irespr=0,
+                  mt=None, irespr=1,
                   temperature=NJOY_TEMPERATURES[0], mfcov=33,
                   iprint=False,
                   **kwargs):
@@ -606,7 +627,7 @@ def _errorr_input(endfin, pendfin, gendfin, errorrout, mat,
     iprint : `bool`, optional
         print option (default is `False`)
     irespr: `int`, optional
-        processing for resonance parameter covariances (default is 0)
+        processing for resonance parameter covariances (default is 1)
         - 0: area sensitivity method
         - 1: 1% sensitivity method
     iwt : `int`, optional
@@ -644,31 +665,31 @@ def _errorr_input(endfin, pendfin, gendfin, errorrout, mat,
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    0 33 0/
+    0 33 1/
 
     Test argument `temperature`
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9440, temperature=600))
     errorr
     20 21 0 22 0 /
     9440 2 2 0 1 /
     0 600.0 /
-    0 33 0/
+    0 33 1/
 
     Test argument `iwt`
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, iwt=6))
     errorr
     20 21 0 22 0 /
     9237 2 6 0 1 /
     0 293.6 /
-    0 33 0/
+    0 33 1/
 
     Test argument `ek`
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, ek=[1e-2, 1e3, 2e5]))
     errorr
     20 21 0 22 0 /
     9237 1 2 0 1 /
     0 293.6 /
-    0 33 0/
+    0 33 1/
     2 /
     1.00000e-02 1.00000e+03 2.00000e+05 /
 
@@ -678,55 +699,55 @@ def _errorr_input(endfin, pendfin, gendfin, errorrout, mat,
     20 21 0 22 0 /
     9237 3 2 0 1 /
     0 293.6 /
-    0 33 0/
+    0 33 1/
 
     Test nubar
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, mfcov=31))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    0 31 0/
+    0 31 1/
 
     Test mubar
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, mfcov=34))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    0 34 0/
+    0 34 1/
 
     Test chi
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, mfcov=35))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    0 35 0/
+    0 35 1/
 
     Test keyword `relative`
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, relative=False))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 0 /
     0 293.6 /
-    0 33 0/
+    0 33 1/
 
     Test keyword `irespr`
-    >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, irespr=1))
+    >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, irespr=0))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    0 33 1/
+    0 33 0/
 
     Test keyword `mt` as `list`
     >>> print(sandy.njoy._errorr_input(20, 21, 0, 22, 9237, mt=[1, 2]))
     errorr
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    1 33 0/
+    1 33 1/
     2 0 /
     1 2 /    
 
@@ -736,7 +757,7 @@ def _errorr_input(endfin, pendfin, gendfin, errorrout, mat,
     20 21 0 22 0 /
     9237 2 2 0 1 /
     0 293.6 /
-    1 33 0/
+    1 33 1/
     1 0 /
     2 /    
     """
@@ -1069,7 +1090,16 @@ def _run_njoy(text, endf, pendf=None, exe=None):
         logging.debug(stderrdata)
         retrn = process.returncode
         if retrn != 0:
-            msg = f"process status={retrn}, cannot run njoy executable"
+            cwd = os.getcwd()
+            dir = join(cwd, "njoy_outputs")
+            if os.path.exists(dir):
+                shutil.rmtree(dir)
+            os.makedirs(dir)
+            for filename in os.listdir(tmpdir):
+                shutil.copy(join(tmpdir, filename), join(dir, filename))
+            with open(join(dir, "input"), "w") as f:
+                f.write(stdin.decode("utf-8"))
+            msg = f"process status={retrn} when running njoy executable '{exe}'.\nInputs/Outputs were moved to '{dir}'"
             raise ValueError(msg)
 
         # Move outputs
@@ -1436,7 +1466,7 @@ def process_proton(
         # If isotope is metatable rewrite ZA in xsdir and ace as
         # ZA = Z*1000 + 300 + A + META*100.
         if meta:
-            pattern = f'{za:d}' + '\.(?P<ext>\d{2}[ct])'
+            pattern = f'{za:d}' + r'.(?P<ext>\d{2}[ct])'
             found = re.search(pattern, text_xsd)
             ext = found.group("ext")
             text_xsd = text_xsd.replace("{:d}.{}".format(za, ext), "{:d}.{}".format(za_new, ext), 1)

sandy/sampling.py

@@ -322,7 +322,7 @@ def run(iargs):
 
     # ERRORR KEYWORDS
     nubar = bool(31 in iargs.mf) and (31 in endf6.mf)
-    xs = bool(33 in iargs.mf) and (33 in endf6.mf)
+    xs = bool(33 in iargs.mf) and (33 in endf6.mf or 32 in endf6.mf)  # this handle together MF32 and MF33
     mubar = False
     chi = False
     errorr_kws = dict(