-
Notifications
You must be signed in to change notification settings - Fork 28
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* update * update --------- Co-authored-by: Luca Fiorito <[email protected]>
- Loading branch information
1 parent
593a57d
commit 7876500
Showing
4 changed files
with
4,399 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,265 @@ | ||
{ | ||
"cells": [ | ||
{ | ||
"cell_type": "markdown", | ||
"id": "46d48b51-b649-4bc7-875b-648dcd6d2f56", | ||
"metadata": {}, | ||
"source": [ | ||
"# Create perturbed $^{239}$Pu thermal FY files using CEA evaluation" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "80083589-cecb-4332-8aa9-02278c5a9663", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"import sandy\n", | ||
"import matplotlib.pyplot as plt\n", | ||
"import seaborn as sns\n", | ||
"from os.path import join, dirname\n", | ||
"import pandas as pd\n", | ||
"import numpy as np\n", | ||
"import random, sys" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "a7f2d039-5285-4b0e-9f84-55294a747dce", | ||
"metadata": {}, | ||
"source": [ | ||
"## Extract FYs and covariance data for Pu239 thermal fission" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "d572351b-4450-4c30-8d81-b3cd1a0cb310", | ||
"metadata": {}, | ||
"source": [ | ||
"If you don't find these files, download them from the [NEA GITLAB repository](https://git.oecd-nea.org/databank/nds/jeff/jeff-lab/fission-yields-pu-239-th)." | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "4e7b4c6d-5c73-4e80-b90a-1ea90137cc87", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"path = join(dirname(sandy.__file__), 'appendix', 'fission_yields')\n", | ||
"\n", | ||
"file_eval = r\"jeff-4t3_cea_pu9_cons_28-09-2023.stn\"\n", | ||
"file_cov = r\"jeff-4t3_cea_pu9th_cons_28-09-2023_ind_corr\"\n", | ||
"tape = sandy.Endf6.from_file(join(path, file_eval))\n", | ||
"nfpy = sandy.Fy.from_endf6(tape)" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "cb8852d0-e541-4580-a4ea-e4139267f513", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"idx = nfpy.data.query(f\"E==0.0253 & MT==454\").index\n", | ||
"ifyth = nfpy.data.loc[idx]\n", | ||
"\n", | ||
"corr = pd.read_csv(join(path, file_cov), sep=r\"\\s+\", header=None)\n", | ||
"u = np.triu(corr, k=1)\n", | ||
"corr = u + u.T + np.diag(np.diag(corr)) # must do this because the matrix is not symmetric\n", | ||
"acov = sandy.corr2cov(corr, ifyth.DFY.values) # absolute covariance matrix\n", | ||
"cov = np.divide(acov, ifyth.FY.values.reshape(-1, 1) @ ifyth.FY.values.reshape(1, -1)) # convert to relative terms\n", | ||
"cov = sandy.CategoryCov(pd.DataFrame(cov, index=ifyth.ZAP.values, columns=ifyth.ZAP.values))" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "ce955aa8-843b-479c-a15c-69ee6f978e60", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"fig, ax = plt.subplots()\n", | ||
"sns.heatmap(corr, vmin=-1, vmax=1, cmap=\"bwr\", ax=ax)\n", | ||
"fig.tight_layout()" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "3c338b9d-f74b-4f20-922e-37f04c0326e2", | ||
"metadata": {}, | ||
"source": [ | ||
"## Generate perturbation coefficients and write them to file" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "57d4f1b7-2762-4153-a814-9fe09c3fe994", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"seed = random.randrange(2**32 - 1) # create a seed\n", | ||
"seed = 1444271357 # or use always the same\n", | ||
"nsmp = 1000 # sample size\n", | ||
"smp = cov.sampling(nsmp, seed=seed)" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "a0aaa60b-0c9f-4993-91ea-e66db301e0f5", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"with pd.ExcelWriter('PERT_94239TH_MF8_MT454.xlsx') as writer:\n", | ||
" pd.DataFrame([[seed]]).to_excel(writer, sheet_name='SEED', columns=None, index=False)\n", | ||
" smp.data.to_excel(writer, sheet_name='IFY_TH')" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "ea504633-a1fc-462b-8808-9a45f8361d52", | ||
"metadata": {}, | ||
"source": [ | ||
"## Read coefficients from perturbation file and generate random FY ENDF-6 files" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "17d1ac55-f384-4484-a0bb-7a4cc2dee76b", | ||
"metadata": {}, | ||
"source": [ | ||
"Skip the part above if you already have the file of perturbations." | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "632bbf34-9030-4b70-bc5b-51501b3dff77", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"smp = pd.read_excel(\"PERT_94239TH_MF8_MT454.xlsx\", sheet_name=\"IFY_TH\", index_col=0)\n", | ||
"smp = sandy.Samples(smp)" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "ee52ce9c-f451-409c-b686-9bf83a446728", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"idx_ify = nfpy.data.query(f\"E==0.0253 & MT==454\").index\n", | ||
"idx_cfy = nfpy.data.query(f\"E==0.0253 & MT==459\").index" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "b52f4587-052f-4bcf-96ee-aeea456829b6", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"tape_rdd = sandy.get_endf6_file(\"jeff_33\", \"decay\", \"all\")\n", | ||
"rdd = sandy.DecayData.from_endf6(tape_rdd) # this can take a while" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "ced64de5-69c9-4571-90f1-4b6d2394f2ff", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"smp_min = 0 # write ENDF-6 file only in the sample range [smp_min, smp_max]\n", | ||
"smp_max = 99\n", | ||
"file_template = \"pu239_fy_cea_cons_{}.jeff4t3\"\n", | ||
"for ismp in range(smp_min, smp_max+1):\n", | ||
" file = file_template.format(ismp)\n", | ||
" f = sandy.Fy(nfpy.data.copy())\n", | ||
" f.data.loc[idx_ify, \"DFY\"] = f.data.loc[idx_ify, \"FY\"] # just for me, i copy the original IFYs where uncertainties should be, so i can compare them to the perturbed ones (anyways I don't use uncertainties)\n", | ||
" f.data.loc[idx_cfy, \"DFY\"] = f.data.loc[idx_cfy, \"FY\"] # same but for CFYs\n", | ||
" f.data.loc[idx_ify, \"FY\"] *= smp.data[ismp].values # IMPORTANT, this does not update the CFYs, which in random ENDF-6 file are inconsistent with the perturbed IFYs\n", | ||
" # f = f.apply_qmatrix(942390, 0.0253, rdd, keep_fy_index=True) # Run this if you want to update the CFYs (slower), or else comment it out\n", | ||
" print(f\"writing file '{file}'...\")\n", | ||
" f.to_endf6(tape).to_file(file)" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"id": "337a51c3-8c9c-46d8-9ab4-46b7d6d1958d", | ||
"metadata": {}, | ||
"source": [ | ||
"## Compare sample estimates of mass yields with deterministic uncertainty propagation (sandwich formula) " | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "b6e97e07-dedb-48f1-b59f-4d0ae97e31c7", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"mfy = {}\n", | ||
"for ismp in range(smp_min, smp_max+1):\n", | ||
" mfy[ismp] = sandy.Fy.from_endf6(sandy.Endf6.from_file(f\"pu239_fy_cea_cons_{ismp}.jeff4t3\")).get_mass_yield(zam=942390, e=0.0253)\n", | ||
"mfy = pd.DataFrame(mfy).rename_axis(\"SMP\", axis=1)\n", | ||
"\n", | ||
"S = sandy.Fy(nfpy.data.query(\"MT==454 and E==0.0253\")).get_mass_yield_sensitivity()\n", | ||
"C = pd.DataFrame(acov, index=ifyth.ZAP.values, columns=ifyth.ZAP.values)\n", | ||
"cov_mfy = sandy.CategoryCov(S @ C @ S.T)\n", | ||
"\n", | ||
"mu = nfpy.get_mass_yield(zam=942390, e=0.0253)\n", | ||
"sigma = cov_mfy.get_std()" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": null, | ||
"id": "261f5b3a-645b-4966-a06b-2d97862e1bf9", | ||
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"fig, axs = plt.subplots(1, 2, figsize=(10, 5), sharex=True)\n", | ||
"\n", | ||
"ax = axs[0]\n", | ||
"ax.errorbar(x=mfy.index, y=mfy.T.mean(), yerr=mfy.T.std(), marker=\"s\", ms=4, linestyle=\"none\", capsize=2, ecolor=\"blue\", color=\"dodgerblue\", label=\"sample estimate\")\n", | ||
"ax.errorbar(x=mu.index, y=mu.values, yerr=sigma.values, marker=\"s\", ms=4, linestyle=\"none\", capsize=2, ecolor=\"k\", color=\"k\", alpha=.4, label=\"original data\")\n", | ||
"ax.set(ylabel=\"mass yields\", xlabel=\"A\")\n", | ||
"\n", | ||
"ax = axs[1]\n", | ||
"diff_mean = (mfy.T.mean() / mu - 1) * 100\n", | ||
"diff_std = (mfy.T.std() / sigma - 1) * 100\n", | ||
"ax.errorbar(x=diff_mean.index, y=diff_mean.values, marker=\"s\", ms=4, linestyle=\"none\", color=\"dodgerblue\", label=\"$\\\\frac{m- \\\\mu}{\\\\mu} \\\\times 100$\")\n", | ||
"ax.errorbar(x=diff_std.index, y=diff_std.values, marker=\"s\", ms=4, linestyle=\"none\", color=\"tomato\", label=\"$\\\\frac{s- \\\\sigma}{\\\\sigma} \\\\times 100$\")\n", | ||
"ax.set(ylim=[-100, 100], ylabel=\"relative difference / $\\\\%$\", xlabel=\"A\")\n", | ||
"ax.legend()\n", | ||
"fig.tight_layout()" | ||
] | ||
} | ||
], | ||
"metadata": { | ||
"kernelspec": { | ||
"display_name": "sandy-devel", | ||
"language": "python", | ||
"name": "sandy-devel" | ||
}, | ||
"language_info": { | ||
"codemirror_mode": { | ||
"name": "ipython", | ||
"version": 3 | ||
}, | ||
"file_extension": ".py", | ||
"mimetype": "text/x-python", | ||
"name": "python", | ||
"nbconvert_exporter": "python", | ||
"pygments_lexer": "ipython3", | ||
"version": "3.12.0" | ||
} | ||
}, | ||
"nbformat": 4, | ||
"nbformat_minor": 5 | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.