diff --git a/sandy/core/endf6.py b/sandy/core/endf6.py index 92e3dbe1..677e9e90 100644 --- a/sandy/core/endf6.py +++ b/sandy/core/endf6.py @@ -2921,7 +2921,7 @@ def apply_perturbations_rdd(self, smps, processes=1, **kwargs): return outs - def apply_perturbations_fy(self, smps, processes=1, **kwargs): + def apply_perturbations_fy(self, smps, processes=1, covariance=None, **kwargs): """ Apply relative perturbations to the data contained in :obj:`~sandy.core.endf6.Endf6` instance of fission yield files. @@ -2936,6 +2936,12 @@ def apply_perturbations_fy(self, smps, processes=1, **kwargs): Number of processes used to complete the task. Creation of ENDF6 files and post-processing is done in parallel if `processes>1`. + covariance : `None` or `str`, optional + Flag to adopt fission yield covariance matrices. + The only acceptable flag is `covariance='cea'`. + This ensures that the nominal values for U-235 and Pu-239 are taken from + the CEA evaluations. It must be used if it aws used for the production of `smps`. + The default is `None`. **kwargs : `dict` Additional keyword arguments, such as: - `nfpy`: to pass directly an already processed :obj:`~sandy.fy.Fy` instance. @@ -2956,16 +2962,37 @@ def apply_perturbations_fy(self, smps, processes=1, **kwargs): Examples -------- - Default use case. + Default use case (write data to file). - >>> tape = sandy.get_endf6_file("jeff_33", "nfpy", 922350) - >>> smps = tape.get_perturbations(2, covariance=None) - >>> outs = tape.apply_perturbations_fy(smps, verbose=False, to_file=True) + >>> import sandy, pytest + >>> tape = sandy.get_endf6_file("jeff_33", "nfpy", [922350, 922380, 942390]) + >>> smps = tape.get_perturbations(2, covariance='cea') + >>> outs = tape.apply_perturbations_fy(smps, covariance='cea', verbose=False, to_file=True) - Check that data is written to file with key `to_file`. + If the samples were produced with keyword `covariance='cea'`, the same must be + used in `apply_perturbations_fy`. - >>> assert os.path.exists(outs[0]) - >>> assert outs[0] == 'fy_0' + >>> nfpy = sandy.Fy.from_endf6(tape) + >>> nfpy_u235 = sandy.Fy.from_endf6(sandy.Endf6.from_file(sandy.fy_cea_u235th)) + >>> nfpy0 = sandy.Fy.from_endf6(sandy.Endf6.from_file(outs[0])) + + >>> n = nfpy.data.query("ZAM==922350 and MT==454") + >>> n0 = nfpy0.data.query("ZAM==922350 and MT==454") + >>> nu235 = nfpy_u235.data.query("ZAM==922350 and MT==454") + + To match the perturbation values, the ratio must be taken with respect to the + CEA nominal values. + + >>> sp = n0.set_index("ZAP").FY.divide(nu235.set_index("ZAP").FY).fillna(1) + >>> p = smps.query("ZAM==922350 and SMP==0").set_index("ZAP").VALS.rename("FY") + >>> np.testing.assert_array_almost_equal(p, sp, decimal=4) + + If `covariance='cea'` was used to produce the samples, at it is not used in + `apply_perturbations_fy`, then there is a mismatch between the ZAP numbers of + the samples and of the fission yields. + + >>> with pytest.raises(Exception): + ... tape.apply_perturbations_fy(smps, verbose=False, to_file=True) """ from ..fy import Fy @@ -2974,6 +3001,20 @@ def apply_perturbations_fy(self, smps, processes=1, **kwargs): nfpy = kwargs.pop("nfpy", None) if not nfpy: nfpy = Fy.from_endf6(self, verbose=kwargs.get("verbose")) + + # Change nominal values to CEA values if asked + # this is needed to ensure that the samples are given for the same nominal values + if covariance == 'cea': + tape_u235 = sandy.Endf6.from_file(sandy.fy_cea_u235th).data if 922350 in nfpy.data.ZAM.values else {} + tape_pu239 = sandy.Endf6.from_file(sandy.fy_cea_pu239th).data if 942390 in nfpy.data.ZAM.values else {} + tape = sandy.Endf6({**self.data, **tape_u235, **tape_pu239}) + + # also the fission yields must be re-extracted + nfpy = Fy.from_endf6(tape, verbose=kwargs.get("verbose")) + + else: + tape = self + # --- PROCESSING if processes == 1: @@ -2982,7 +3023,7 @@ def apply_perturbations_fy(self, smps, processes=1, **kwargs): # assume HL, BR and DE have all the same key matching (sample ids) for ismp, smp in smps.groupby("SMP"): outs[ismp] = fy_perturb_worker( - self.data, + tape.data, nfpy.data, smp, ismp, @@ -2998,7 +3039,7 @@ def apply_perturbations_fy(self, smps, processes=1, **kwargs): outs[ismp] = pool.apply_async( rdd_perturb_worker, ( - self.data, + tape.data, nfpy.data, smp, ismp, @@ -3260,6 +3301,7 @@ def fy_perturb_worker(endf6, fy, smps, ismp, fy : `pd.DataFrame` `data` attribute of :obj:`~sandy.fy.Fy`. It contains the nominal fission yield data. + It i sassume they match all the ZAP of the samples. smps : `pd.DataFrame` It contains the perturbation coefficients for fission yields. Columns are `MAT`, `MT`, `E`, `ZAM`, `ZAP`, `SMP`, `VALS`.