From 6e5849e1ed8e1cea65e8ad1f8d713a24f1313aa9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jul 2018 10:02:40 -0500 Subject: [PATCH 1/5] Support tabulated fission energy release. Remove compact library --- openmc/data/fission_Q_data_endfb71.h5 | Bin 67543 -> 0 bytes openmc/data/fission_energy.py | 404 +++++++++----------------- 2 files changed, 135 insertions(+), 269 deletions(-) delete mode 100644 openmc/data/fission_Q_data_endfb71.h5 diff --git a/openmc/data/fission_Q_data_endfb71.h5 b/openmc/data/fission_Q_data_endfb71.h5 deleted file mode 100644 index 6e74e90aa3664e3e1f924ff54f6491a2007cdf96..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 67543 zcmeHw2VB(3_J06TY=8wtY->TWAt230G6)Ewf;16SSb7I3_OdD}RU`+joS`{4aQFZa61OioT_&Ua3jbLPyiAzaA z!s8ct&znUFq%<;JE_7qp!@-3-Up4EkG~0KECJ`Kp$4-7<1jcoc+gh^prWX zl{x0#?tyOOF}|k$?tv>D-JE?qeaCZVDszCp1J{R(7S0U2c69gh;^No3{?5zD=TU~b z{=Uml)+o&kuAgt9GyWEo^$f1RR{-`>^!4Ugp;4xWXi1#?{oS3>Y-454%peDh5Be&s zB3v)@S(MVy*TKs@ki7Kr@?FM7>tlgFsmz(>>*j;=M3ctcWoV`Ff@5`ZcXmPx65#2M zzs1A$84@Y{y19g+0IRu)u% zX#M6m&Ik&0L=nSpqj__7dCiGD}9T@0~)f;7o7I}`Nv5&i} z8&MwvG#%vX#&y8{;T1>Uj=p2AuY*6z4Hbkb*UQ`8H!u*b2AYqig1jBjN}_zxDx12y z`k*(gP|2YQe_x*!T<4|ui%d`%b5~%(M4v+cH}~~f>5CJi6-Ucu9)vbrtev6nv*3F9 z294(!q2fdnZob}rxJ=NeEc`uv34W-6tbARGLZgo_<@#`(oc+NGF;mpEwY4UwO`fdDQBqS; zRUM<;!xTa_UKZ8(zPQHwp-salkmKu*N&wehe2Ko0ST7EY${#xrNWD9x{tJCb68%a`bd&0h{)(xPWjRND z_#r3^C+<)`efLc_Z5zNVx{wtT6~J0k$89$-H`j+IlbDL3P3Rj2h#2~+FPJdJiXli5 zYgkXJ#y4WGZx49Odtq;eppFc0a9y=T+u*epYJjNV6>H>NAq4m5zwH5USRbN%9U0y? zJ>#m|2Cr6k5(QDgTk!Mp=|2aSplmI?!tbfDX$n8i!u&t^I3sezyIa1=Ef=jV{Xa*4 zEG6np+u-fjlPsO6;GG&A|Dd%MZY{i`_}iOB!rPIxxBtHJDQ$zdL0ZJ{Dt-t*{7=Ix z)Y7(SFRwpQ$`++~ zzMnH9;x^!)D3Xn^ER3Iy=a=)4;HxE%*7%Qv-^*lfQ;~|q@_D;5Z9v2i!(IcdcDKCu z!@Mg)PAkp%<-%%OWl8^DZxN0_-gdfQ*)>7$tLuSok z(BElLS{PcI&PFe|f7TovK7<`t8T?oEJu~yg%3vygw{H#6U`%z-TPh1`_~GM5ZrI|3 z%SXVibOJVtKOAh6a$#S;ZnmbI`(X01U$Hj^+}<2F|6yG%l^z+j%NacXb~a;2VF8ss z>*2jfpz_mrhu#MoUYr4&?7{G<=F-9Lbi9sUHwLMbN6cEakq}urb zUCudwx83mGnk*mrY-v1SkW-Ef(>eE&kC^gDST2O7!8ZxK1bYy}RP*%km74Iwn(F}< z-XvrBg6UFByfC>8PN-qU>v%L1uV}yN-7crH;svP&GVxl70&FP$fDJE5<}nj5A%7OU z5U-EXK>V1|k;GNC{Mq*}_noPV%ZE2E#|890qnR?C^BnGTYi#5dopG3a(jV`L0Y^1o zs_V(sQ|XbrTU|kH){32rBFm`sSNh+I1RpIe_MCl5!^_@mtv#4H`qRcnE*-C|DkFe3 zpF?BQVe>fi^E)2CyotUSf)>9+vY}G(uG>Okt*52_-fz0&cpqlG5KpPWW@`C z`&cma3&R)V&w^K#<=bcw^yJkMZ=V`|RpPo6t6#jp+=yNT{n-p_e8bVGWEjQ7N zU_%?)ed$6v{Y5`B+Woon7f{po)z@OE>8!X-vD9?BH5>(&Nq76y-}+9=_XmBgy=;v= z80K@#GhcE&_4@-^%rXYMI@k0wG|ht7L|DYj^i;*$2~_+JDP6K{*Q>fv*x$#!AT_uL zksq3eOxbOpcqm1W-|O9k^@k3eC)y8hGjgbISiHi`v9PTBgpFZa;{HziL9Pv)$N1i2 zxb||(V{1Qxlmkq9to`&af5J*XJ$62${RsK9v>)rahhhNn9joEHW#odmA-H@HhXNN+ zv+u6q>*x9KI-8U6dwz_@@BVRg!Rd%wV#&vuwAsR_g61Z z@)T9)KIYBOO15JDyS(=3+^cYi5_ROxz1K2>6>WnzTY;3jsNjw7V>SQZLCJQnJ>V@L z*74zWQy4vsfEPyn$kcy->KIp15;44=D?DNh%nVH-%pog|jt#UNLg#!o=_nC8R!Q$Q z1a&4!$@AQnC;$BgI@l<5yh*(PKPuYB}4QE!5Sum~8+_ zzQ}kg2Hiu;iAtmIlU|T7`k)MyHSt?t$Pne;2by|5T@F%2OD78n-d*D0qZG1`wvto}H-)-NSR&Jby45vP4TD3%zyBn+(0i32Hx0nF0mj z1jR44VQ3S6X?gD}fRa$Yve4GVLAua9yg2qFde4B4;JI(gZD&HSaaz@wh}n=O&g_uJ z<=K!uPFr=+&lJ+YbLFNTH-j4R%Vw!(3m_$&gST4!0>~4;{Gis$77E7+R!ElFLP z_=|7tp^^CI+*LeBNFOKI+!XHuS>i9w+kV0oDntonp^f|x9?*82*~|X3J)v-%pfEqo z2Wr4C!8OeQC)7glnkOodLQrUD3q=>)mlwaqykTafZ zy<^sTr~!ZMSWeD*C>*~u+AFaElEZ1`wIuhUd96L#rku0Kkna=Kw^ypfkHcC8>kkv7 zKnMdniN*$jA?%`(RP59CS%_HB}O(FV`OB)f{m6}O9 z^N+w!uGp38OqR#<1sf`n%OBq)NMBvRx77P`a?OhSME}xsACmhnDM>s<8;*Y(zfk>N z%HQc#^rez}gXb7STHPiJuXAWs`yzyY@+ncLpQaL7ZOq|p7iJ$WTvH5PEXep}jm-e)R3 zwEMa%7`?`GN&fRvDm_ojJPLF@VfIkgi-s4nQ;Z8%U%8X$(VuP)^w2Q@LE?v{1|Gi! zJC5EwKIHjS92efjMr88E^rg2>9EQtpYwcNgyFZ?fMAstTL+kb@4fW(d_n7|fsdp)X zmv$6l+dVle!`%Qzri|1q5WkD(ACqRc|5Z=3+yBzT^%?TV{3|O?v&&!oI+NZ7Jz@iH zv*NW*a}EP9A%7OUtmn>+0h!Tra;IOY8JEJKIYP53Z-u3qqE;0_z=igT1VBsr2b7)1ttpOGk~DSkm#jRc#NFZWr&|AVJ4# zTD}n&Y~wVw%;X82<}K)&Q9Oab>r4}JL$Y^YIPEBWLnk+M&dLFJK9bQDD46Njj`du{ z_r2{`W8V0dz$^9&vNGQLj*`@3c=jCb>*G6%u>4m0k6t*;WavBQZ~YXJ!K%NuJq#K2 z`B8q;XL;UbqlZ3@F=gnxkUtAvA$*+}V6u0B@3xybUVAp>gl_ExR zY0J<_*a(72%ZN$K1#1!8B{S*S(`&7j4cgtrkKbh@%>R=?Va`Z`M@M!ntd{zw)V9v; zx3Q$CMAf^uLMI=&RKj>`p!QN>>dB3o(pfA8J%fO;>@X*bQ%;0H_HcEng};r&R+>;qyQukmE+ z+}&_l=8oJ=`-t&m=*m50XwTVwY~M`gpK$)*x~AeZ(ZA&QBUu%BAKfPOgTFf&tlW5P z@89WPB1v-0@uEg_IDfc3>-ZA13N8Qd@d6nAhFtdXrIKUJ@g*UDmi}c#j&>~Aq+^}7 z`C}dbo7bR6`<-z42yEB5fTF~a7i83O;jPE+7mOP-43k^dorYL2w4m0G6Y`lxf6KrX zr0p7JTkBRprAHQBi2~wDg;$5h((vMi+HgUn;H-(1CLJ$(F*7iHSJ0AyJ=5Xl=C0YR zUXpm#zd`J}WDJ;Bv>n!Q9IbaSNs7SnOeUf@+rj3wIpAwtxIK6B`U?bJyj6(Jh|yPm z9N@rn8!L_oERW%MC9kz**l$lBKYo1oAZxsa;BExd|MFsG-!onA=d$(h5&Z#7ys(CC z-<1_FME^U}UnAttf>*&QwOC-ADZ6XIr8<7%b?fdw8VI}s0$o7ACE*($ROP~)Z59O{ zUq)c^kqvzp3mV)KN7tKtrqQ3!b_Jnpit37o6;kOBhMtQ8Z=%7#hlgo+3C^2yf%>)* zMd|T$yvox}!L8K(qs%2A!S(?!o;jrB_6a(YeV&Z3af6k1$CpgJSw-S<32|G{6uPO3 zL)Bl@0(C>IUhnWo&;8CtgQ1(jYZKRybS?KpgPPbdv9=*?^{*!c0sGv(oQydu|`0*%GjjCzPa zf{muj>~o33?SBaJtw)aWMxmiS2jJ4=%8FwReF^=NFGuz)UomC=K3)EgZ$G*`=i}q* z5LC5^ZMG-BnN2RB=cVDRGUCdq^n%2m(YSrg3m$@>mCl$YjxYFYFAF>3zPWeg35Z!Ls|N5%x9jMR$xYBJ2?1U*kn~#nOV+ z?@51Ounx}0f=84`prb(W4M|2&I#fSH@w4{Pg?K(MSec-An2(?BUyj1h!TCg5SU;iv zkM#FJYfFjpd|f*UqZdpmCCbCe8-YKMz3Q)y_UP7yO(IIuk!@Vl(-S`CS@)v zqv^=4U#HpwE<~f_!{zg&Ph#8P+BH?ga5-Aeif#|MRFIAjSDssV)QKhYieV<;Ool$@gXySwI3Lrp#q9UsxbH7H5Or#;{bo6+&LtM7-mSPa9T z9$~%@qQH;)&2Or2(fYfBaCa`S(L4C4*qGkm4KX(XYtCG9SXcN2HrRMuS_>IVz`G#_ z!Oj1t>)}D>-WQ+j5+m~GrXhb?7>ur!_6gKUnc3g%h7K6rRAKG!3iP6w^jNs@cZ9XS zD~Mji?C%Qsv-Ecb7ly?GzeTDeZv3d_Q{*Fv`RNGeHt;X&ZB3)f$L~Z=Ea=ws@P_4q zdMZ5mv8N}ZhbhdiOu%NEH+gMPa5TNU~}quP7CoKsxzb;YIY_a24O>^+38{n`*39B#cmbjrW3gAsz3-o&SrSA?A}*`R0$&g)s-4y4r4G;+_U7Ju))b73gWcuZ;hhN2NbkiP|TJ>sE7jO4IN{@?y9k*t5%K z&OJI_+S%qn?`!pY!Q@Bq8-?+SjkVeYE{DqzpC~bO2saqc3^+f>{xPXvYPS%1*_8CB z%hdQp98C;?7vh1)7kz@OLS^7zds4#HE4E|#g09OKGW<1|Jgv3sSm$dvdd_9VEBHNY z{6Ca8lNpbLjpy5hK4im-$L(y%iWkP81uuc_7Svy6JoH+`+Iqf;uZ>>ecU(S#1(7ab z(SB#0yxTdj%>8(k4M*|ucFxj?ap20UE9O(=8))=rEnUG*XUj5+U->k8zcJAuLSbUb z$j(ox`UN3XgTdoP#_3(AAED{j(Lgg2uV=J;s6QG<6J}(Par-zuzv>2p;=nwtcO%S9 z2r`|p0DMt>b4Kx1GK~*r52g;9zSj5II$D0TqBR3AoUca}oBW~pN3Ek@dd{AK7tSB8 zFAHAhW=O|@Yo1-dM7T!rsrBQTh=YZv;#1Xo_)^O+YdjzxTpumdJta1kMi1Tmh3tXl zmw%kTgGw(r-$b{!3#y6nb0is$cY&agI|%+>a@@VDfQYXV+S~GT3Hw`UZ)fER?Qz1N zAidBL@1&9THfIsqo-lrT6?nW2G{iZc=)@86JP{?Dgg*jnFBsjcGSZ(Q^hfM_ieUMW zd<=g6TxE^ZMW9QDVA;{j$@m^ro>)3!ygk7`q#jRizDmUVK*e9Nw7}#jQ6Bj@eK0ys zRTEJjo=I=wIW+UHj<)r}?--JNMP=h|+Z{KgwcR=~vJAfu3<~rAWZ0cInqblriO_QK zSAzS`=f8S5a&+6^@}EVPT2ydt2%Hnt9&puYc6_*gel>aBHn@h(7BO5|`jdm(11_Fc z$A>HG(B0t#TnHKfMg8}uHg47&5yNGdzQK}!i&`A)`mf`&L=h!{ke`Q+?2;&wduh`) zxYn777_L1DZtq(Qm+)^?TmE*{7<7E?+Hp4{zin_S&l53Ra6k8#t%XYz+BM6Js7OcF zt_+2ZtI<+cRZNTq_d`zD6e!aySHaCYT%bF-aOn0Ibjww$+zg~XJgI-F{wW+jIP2BD zaMGjTSBdOY!sfJZh0DH%mCaw&1-EE7V>AYgu?P2 zV?p?(D`~F&b$qIPntu2Iqlb^LM_eqY$|r37n^-W>)XZM2$7d=%f?V?jhdH5}b{?pp z((@{&M}rk_N2yAk%Hdy=+$`8GtA@)NvqZQSW}_pvM!;ddd5}TnciclT-fNs||Ts}_Kn{$FI?Ugp#IV3y8gC2nlB z9o@5tSw6o<=HRSQONeD958N_hrAI9#HhRo*`h6u2oIh$wvEaq2-x~+eFsJ>)cE!-- ziCJ0?e@N@E*I7c9CvW?wcra)+SItu8HjN&$#O(KKyc?#llS*&;mJGtHEC@0 zh>?;R11}+e7QCdQ7RCZa+eLG`O4afy@)1CKPGDS$(vJ*5DOEm#*VVCr|GBX;@?;&2 z{ztYGK<;wy9-8ukN{@u7MFVVpo2@DhFQiu*7u1?4771R`@p^v36pSeTg06{E0I!k} zUy#vnJmJYh@gi2rv8M8e{%O44$TB_Om}l;AH+E0C`6ed)5)^|+xEnxtYJ`y(+OC9W zhV%?Y+;Af!K?LdmWyewRB|SqPN?2O<`!S+C*t@Vp$B;IY1A?rtBlr#SC)S?uL<#dk zmeE!Yfix}YGe@RehmH?dzShfKXes53;-kU&GtN76K>?o>KX&4& zC_8}pE5x~R*5Nc0vR|Yn7_9B8=kA*Ye=(7A>#7IHu5~3sn7}(=-TJ_U%Al@9{=yqb zDKx1p)uac1#rg|z4GT6BgXE9}(kp+})Mw#ebY`8q{`qw2QapdvH|8KYHi=-_4)#HE zw-LGw(*~n=O})=%+VF~Ene>>cug>fffvsf*(Y^;$J9>{{H|xl zt|3)xuzW#r4s%iiwx~z1@vN>Pp5fpb47{*O=800Qixl#n#4+iIp-X2jIVWJnOVMi< z11}+e7QD{wxE>2~7Zplwd|%H8VI%w#7ZUPO`sD;n#|ToDdRM@zJ10iHS001oHBB)N z_%!*S7__N@O3#zFMFZ+cY2e4u9+c*zHI$kyQ z=3uJWoZy}&nXu&GBY{u7NxUeNhRii$6QRT8q@gt#i1kjzQ#;Sk;0NlvJQRC=nZT=j zCGt`)u{3KI2ev)>ZpniEN3nb)j>oQF#nV~!i#KCFQ(rLGkGUl4qOFkfb0$4ziSLkL z$IH`(sb50=EO==zJB?0iD2w{JQvWkw$&?qwUqQ%cCpt-Ck5ttFtJ%5mwd2-5Uq4gC z@ha~f2fD17ez0#@1C^e)+|3n0r@z!5*k3@U59t~c4T3*#b8M#5@G8$4>j1p#(qt}K z(D9N@umB^ke&6!BI1?`3b=bS~Il1U#$V;T|oT=yGWjkO+!%jhmZSX}gA>{rWMEZ+s z=i2p)_=@hM75htOk^03&0*eE$RUDiOpD6z+{?Yl$-_b9~b1clzFDz*8m>a7r6&leQ zVbEiXyo}Pk&Zb{Hk8S!4{SxwL!Hbt25ep(>#uXS$`^@K!pZKtD5bhck$cH(BC*sx8 zV>;)<=SJD*Jc}KI;br>0F&0dL?=30MZJ^O#ed`30c7#28U|B?^M=CGUT@9QogsTAh zB1f3y=;AJ?wimlW&{hQH4=L>xPU;{geY{E~=!HQj1qfm43tbI0pJ=WI`}y|?)4o7a z72|Jjl}84B^yz`WFRU^emnTm!5=$dT3rW`p8bpJp(V!UGBotoF$VE7PgtQ_~4|!Mv zEGP^jCcPl(0~XNIvtE!vDMZjqNNpeD{aGJmvFDXK4^P0$!3&ng&-gq4;-_uKkfCt0 ztEhU{CEuHWwMTc&*@=)(NA6ww?GFoR8(d9ZWGfIAT;{R|$WZX#JN*CbFu9yd$A>F^ zTVbcR!4>NxVz@4j+ML%~?Ghu)&=yxiw*>@~77s(q1vzS%I{yEKXo;9@)t2&sdv`IwH7W>L^a%VB`VU9|8{x4H@n$3 z?Me<1FB{c&P&bE{ZAcnTKee zq2~V9z_<5C!NZeJ;Cf&8jY#++f?`9e0!2HAL#Nssl! zr(Mrtqlc!hv0+$3LjEk45Mqc93cmNBh0MO)z;7PWecy#cxO@a7ZCrrC&qVJWi$Zwz zg5oJ==T$KIAZo|rL2!>(g|BU!sPs^m39i6NeC^JWsn4nORj%?e;Kh*eQ#osBcpXbdVpx*j&@B`1|6^NYYG z{yehvdTz#cA0__${5$C;6C+5xtdL>S>3tnuebXUuT#n@zdwtJwZie*7@T&1Fjt3VeNMxARHB#vzpF&q~$i8IottGis z`iO0^F(5JYsd1h^4KIQ3RtNA>`QfGs@921C46z2UxhM0~KNXPMNDLzM=*w#)wwr-s zkoY!OBJF7DjGm-tNHz*tB>z^wdaWGaEXFKlW#T#lFUT4(sgU<^R&LN4b5`clvzkNG z(DJu-`XA3@6B92?UwZ6j#S5wE&BP094~x39;>GL#ovB|!{w#R$Tq{WZ(yiwYywh*l zsU3uTlooq{)!VZ7nJ>wMFJ0i=4>;8ulbf~SHB!HdKU3*>o2@*7T*iYplF5Zsdi#ad z(ZB%{zcvjo2=b5z1J^HK#J7*5eGl#J=L6hPUXx49htMFey=Ga}iaSrz^B2Y;7!HgB zUH*rWhh~UhUI4r+za`XW--P#14$42RKKK?!CpbHZ^hj+XgVZW<+;grl@xu9{%Wbd) zZS4R|y#5`1{4}P1;q<5{g9Wd0u}g6vD|tb~)uP*cYWYK=RKUfdV}{=uyNgq|wVusDb{6->YYv9iq~od;gdo^ih*c588kbK7yBz{P8(Kzu&1FJ$+*{Jt%~; zfS8^{KF%>7PYd$klyC`8jPME^3rqyE^2YG>wngfSGJ<#!%kDBA|1JB_)z}&BeaxK zJG;ezrDk^fmwM0)SKioQ2e2(Qv43U--O*qQTLD{LovOl9rSSPd@PcRW@cnYoq3dVJ z{cBftcCJ|pi)+t%ojJXpa405iPq??3vERg5FE|pxH_f@CP z`^?Sg^LU6?YiE-oj_5j{Y_p9_kFd|4sBL7@VpX2);L&&sGNVf^^f1a21<{J%n|2A2s zg20h?1W{d8xU@(Czy~(>(1OPBvHaGCcp;B0n2rWRG-ht`JXXgL^s<~ukIfc)Th8iu z<-NGVbPNgkv)~1NN{s`(88aW)th(;+;)wn0k zZSW@k6ork&;znNyygno&R<#CicEi@N;UZ17JBxDv7GC?hF&#bFY)XH)h|MvCJUzms z#~er}*ReYCpm7(Ojv*m`7QAHlXT$;D8&lP+j(p~?mMK{7pNq>!a7^6|yt#6w{-J3$ zd`Z(%D$jN}CZ8~^Q}IB`D)-SR^+qZ^lE-le*z(ei(YaK5rDa=UfU9`u)zv{S_~`S( zqr6|812E0H{@(Z+ecd3p2pcffFC!$oGzWfo?uPz?_oRMNLMSEH)xX#pLxxbg+(2&M zI;I!APLe<$X5)9hO)E8`a z9_I+Fe(@eFF!92_hrLXDt(R1ZY2X{=d8-H4;PT-E*PXs^uyB20;eC%6w0(mi>sND#a+3u69B0$gizCT(i-`NO zBLK~@2x*qm)57<#o!i#C9=Dt1DJnu!e}lexw%?57@8!8?NGdZF!7c#+6Re?k%pgEZu;M5#s8mfR!nSGJ=mga zC?tvcC1jyHoJ*n5PP{2-e%!nl(#0>!M<3Y>rF6mw;-|hn2nC^SLl)Y-()bW0i4w>{ z%Pr)NL6|R07TWZ1ek7EO(`IBV#zUL%%M0nMXP_AT^2n_|7odAMZQ(A7WavHqV$J1~ zlOcV4XXYybA(xSdu=xs7)q zNlEN-=z;b3A!j@{W%F(`TZ#*yw|LQy51sr1dVpU>md>nzUBllhHkA4}8 zmQql*EC$RnpL;Lgjb>USi(ff_C99))bdRQ+)?Nm-ASHIi{V`MW;7=3$RYNc0E*J>9 z^ARa4?~ZPAyc~do9XSXBqcz$b_G*5^nVk{<2N^WrtV6lsew$E8zAz?GD7_Ldo73EWz%}W2IdJ!WA zO~{|cv^Evr91k)hPOjTv+`#`N)8p-juK3otQ0{wAuwG~PXruYL@chpn1uy(?*It2N zPy%qDGIhY9*Nrs#<(<93R0kc8MHdUG^#3?sP|joIi{U3hjFdiq;QNR!2UcGt@Di*> z3=8j{y&KPgHE(ZURnr;Y(yX;HXo7=rOuVoyi0u}#;w9Lu$HWVxhpSle;>mwv;w9wI zf|pYBrg(5^!b7hxqXzz6-}QDOm3VvRSuAx2T8}rs?~;@YZ}dF-4xaKayiD&qb{_k( zkxCCWukZkS#wDGo?ovplKc~Mr23$I<&~sBDUB9~3I)JLrT!({K==wFpZUO0PETG4l z{na3XicBV|86Y zQq@e?u#i6sUd=1G@!-fbXN^UfpIh`Dt+)Nl{sDm4QsE`FzE@soCjj-lpBasAjWl}I zy8)nh-qEOC8Lz4I|HBIv=6?}hyvT66Yfr)Och}wtgsKE$t~q`wbpi5Vogs+bZBX$jvkT!a38|7Lpnl{PKL zkvA-#R-TBecv^XcY+283>#ps`I{4Tjs_xpoU|p&9=&mInC+g9WyK4uUf1Ga{T>5cD zEky>GV}%P5x$_@3hO9%0inQ23Ef-`3JHY8D=MFI!@NdTqevTJ0Tt<7>xwVIOIYtl_ z>B!nO((1$QwrN*XqKM(j&Pl!jZXW5B4#GMA4}dqwSvBUQ$ZpwH2-NL>bfdQIfTg+TStSo!drr|^>Uer0n~ zX5!Hm(1U8kf7>o})#W{Kuw4C7`!3}4#P5-a)ep!f!Ag!lGH>BeJ(Dej=?cw7RD-_X zm>XKBb9TVfU)Br#u>4j>$@JdNG+nVMCkZW9)0L+>nrXUX3sKF_WHnu(#QRLsRmh*k z6yh!Kk^r!vrkxL)_&RAW&u(AB>Oq!|>wC)-GmDPsHT2>IrMq1=Q^>b2l_>_vQ4x1#GS! zp44onsf$HjY;s_YYJ$>hn7gM5`Lp2F9Mc>RF54$4mpeD|PkZl}H_8{6k6=Td8<>ZT z+&M742fVB^X*!uZ_c@Ee`!<@8AXJALhz-k1NXufO7mR5-tR?xmy$}5 zaU)F(hpjW=_l91}zD>AJe1DK_$S|FAFZqZH>^gSJjj-s`SU!>w%#0$!TpzZtSf|?x z3jCP6k7Ls-cjU9~E{{zA#H7d8n7$UszNUZWKx>ApQ^=pCJug>ohzB88qdl&7ZR9`E ze4F@f10kOYTJFH0$)WTRzX<-=xiTo(0$(_9dCyA;V7R$XC(EcNDm@Y)=>fW5-#E!! zqm)Wt_lxJ%BB~1SwMO0Krf8RqzQ|kL5of5k1;kE*;yb-IX)5>$stTU}VrWc&( zZO@+3~9>$~1m`tvvD(e#Vc*^gG9@?k-=^0bJ7?V@enHL*)1Pf@uV zT=bn5wnuku$XTKy9oZOq^=16Zw!wAovWVfjyz^yfYvB?j%g~l_*Q);lSHmK&zHi$G zmmStS#|}|xS6EHq)b@buP83m(7Kc#F1(Dwoe&?xdYNI1h`Z`$;QNgt-=0r++z-1oY z@!?W)Ju7INcHO)oVz>^u7%I01Tp8y&K3vQBXVoR&2-kZWZA-jx4zC4?1sL=2aUUB7mjo)CVQ zs7OcFt{1OHPew~wJ{OH3i2WKIIdKxrO~tFybp*MKPwb0Mqr0i3CfI`W{13;cjW2+k zX7tEB#la^_q7KA*r2J^Fn475w;0G_)2CMP03pOD#iYI>ec=f};DZengKfsqW;O@+BRM;>BzQ}G^uW4ci>w|D>XtdkRY$9J$t z3|;eQjTjQ#IKiCkCgjiJrb0a5#{;)>3H9;Xjr>`MF03=sCFC>J)*YC+luE|TD1;a9 z4BuE$hP%m^yQ6Mb@mVu=>Sj05=$Fp)0E4z)E_v5qL*-wE@Qw{-VRx*ERs48H^@wY3rN_9&18}O;vNNVd1YVrw zO-T9y)6^G|OZbmE2S|v2!QF=tbS?@>x;J}f%o;;@h}?)9sTOHCUaj2UL2wpzV_|M| zY6K;^DBqP^&ZTT_azRB2lOFXduFRQ{!%Ba?n?1v=DCEzAmr|K{0@&Z>?HK>`2ENR+ z7rg||xO@Z?r@Mnxi%FhiIJxj3-b=}^wU(IMRb`t*B2YQG;ZAplCK|oK*aNs46`McV zQb46QO*V-I-zTnqyQg0zZ5YZp(h;P^nva@0FPa(=2JM??2Yy{xg>DUy4`|a9)*|SpDm9xP?1N(t*i5=?fv-{rLAuFEzuHPgS;1)j6=^qtqS7G@7u4Zvk2GvB z38)(`0beX{jt+du`&)Qj`Not#=B6FwUdASWXiI`M17Bw#<6tsx<+_%g^s@BDguQd$*1L zU#Rpvg=`P7gInM4qQ-M7yEO_&NXNtlEI$nV(j=<0;_rNe`I$nLl?7*<0pS`Ts zl)~Bzw!--iT7-MvqX5}#^X$gUzH4BuZ5m6H*M1>zxo{mBG4SVB-j}ibTh%jbBVQaK z@_8=E_8arE_IP!JDZ`1p>U!od6t;?`X(sD1lwglP^LvQ3|HaGMBPcigV$$Q2>%X!N zLkYT9Fe501{8{jdun$iF3$3+E2M=iEPfF5{&ve7(gN!nA2R|hmD{?m$!pF{^mhasg zpZpuY#48ar)MhQae4&X(pFG9`yf+*1c~Zl3Dt(wdxwRwS=Fy)WW9Y*+&~I+(tpjOp zW=MYRZ*FNocl#Tv+fm)lnU`%dhLoO=97W|BiaS)irM0dGq0!S8V@N-Z zs7Ocdt~uBr-PAU?@-svXm$Kp8b_wxnxZm;Nn*VAcxs~MKHHJ24i5M>ZL&gW%L%U)g zbbPq{-vzYMcGbuhFr+AWqFAu(+QIH(inRW&VD}kEFf~@C5?V~}?=CG{1eTg$o81?~rw=rBYYP0= z{&R2O=(1bVJK!!>?`G`HCs!HQzJm-J;tDryo5VjYHXyKROeE3YMf{M>X??)aEq!1I z=RSS5Oy7;=w>nI28MM{^>=riJxn3tkaPsC|-otBx~WvXS5Uz@~0n<8k>Qcc;07rqxR0XHR(q zUmLwn&U+6&Q!hSlbs})H-f=DOcoUVL7dPGm92bmza@6T5mA*!8ZYnY_i!4C&tIj4ID&- zf$_TJAT|b}!NWHg@KraMX)PUeiNH%`J(4eRbf5N%k?_i$ySzsQ;Ck0u{Sr9lG4%zT z;mDs+&We}Awgn7*#O%q0v#g8p@h*<#GVl`eXTi&`@mK;twq_o)hnx7~(A)2QALH`j z9lGxhZaq6=-Z&cF|1az0rjuv!eFQlt&n1Gw;qF>r3!16)(A!E6@M3O(dQ4IYmA-tq zT`UN{mwaveC^}x{XC1*vYYEkbc67X~q!)uakm-_JSp=_FiR@{YPkaw~Y4u2T&s!>o z#vUZ`dO_eauLQaGX<>%#I8*-MUFx=+D>sOIG;B}XW+*2ib0SR5Z;)Wuuj?A@`X!gf z3b>dE|mV&7B8if3%}(2!Mj47`N=S@5zSb0GoT?c9B|`_D%HgNK=e zufN9S!@ED$1C$&Nv>o2R7{1r(Zh_7lJYP!uOd?nu5iJlu*i5BIerbAwoS!!zh8-%U z(*N~ew*TGvKPAu}{Y(8kqF5c-*zWJt1!@o3G~{=D+05uX=VaSDy7pyc!xWX@$U%CX zTzkM3Q`qt0I;k|jf7{@ac_CuBOz);GYY(`>Uv+%AqF)ZoZyQ{jRlUb&Q`O=zN3y$@)Y|z zfk)mN57&04?^<>7i#-^HuNqwpyPhBC{KT9b|3{}?eO>EuBXDmBTroq#u$oU!yYh=f zOddvW)`SM~>pz~kemXUb=7>)C7s2boS>M7&J&9z;+pHt$o51K+OIw)_)!G zly#b%Ks|yO-nf~sKl@9hhy^xE>Nv0&+e zm5&?y)A5q)>;&>;X0)9cBm_7KYYo0 z*II{F#NyvjXLQ#^sf6Of`O1950i9eeZ;<{PXffjLn$S7vn;4u>Y4>4cKE7LAYw<$V z&6)NFHVw&j2fM#!`gbNhwqnTPMXdfBp8QazzedQP1uwxlg+#FT?jHZ9FU|aj0f8F5 zf8p{G=v2CcAEOqZ{aT$1KOAJ${Xm`~CO7NIK}lf3;_*Fp@BK=phXU(8fPK9^FpVyx z(#Lx|js>~}W$ASzX?O{G$U1?6t4D8nBu3XSWd}P@g+*nQl)-nW3~=ptkDS(s`fChL zzsm>YZ-D)}4}3dxE$Od$c^i3`mp19nNj3f(i;YS8hsj+RAs^&%dFgu_)!s0%yI&}# zcdf;XcQcA<|Ks-Fa8~;tnPbJY|FOlP{Wh}N|K%6DFztUKe-^wTSUC}#wK-H{JiUqU z`~Blq$0x+J#uJa-LHDA=juRgi!NFgr+}m~ppVqkinobfh`mDC?bzU=-p11Y22N;O1 zlQh_*luD0$xDgBNOD7af=tIMc=PvF9YULYNDb>;S%Qes*WHy<@gYfE4|l@;upGimJb@#%pUKQ~uaA zJpXETyh``kF#I*x3S}R6vrh99JS}13CFIY7SF^WHBIxINC4G8HBcCFl2p6{iFnvo^ zTH@u`FhxFovlb?SbDhi_@|QMK>5+Ok5Y*fn+kY&-l1h&@>Hi&nR{})0j@oU(uK#+a zZ87pizYn$x^Zz9BSotEsq$Ps2<>K$wLiEoKTL?L|ty`z}kt~&{Y}~CIs&}>)E>Xlr zG*l85>B#PyF?#*FwhgZOPa=j(uKz)s*1{#!qPA$)t^We;vQk)pPFvSA6X*Q?<@r;i zDY90?a2O~A!r@*m<**)y{7WCJDyF>w3_~(uf z*QL$IfE4$iocs(@VnJ_1v&u=#8do&R9N@F{1iG4OY}4-{t>7^^no0|vk^tvt9 z=OZp3ztp8k;QAV?Lt`d6uU-`~p4&<-y$ne()`Lp09+h~{w?re0)ziHXX56FMCbMefHc)N|) z^#IKmqDBwAmj}PMH=ozleJF1KPhXt`77m7%%>UFxrRRNH?+K*7Z&se8SwyA(6CL~? z&Hp6uKMDL#0{?Ldw5?mG^n)z0sBGLDoA$qHk6yTbBTXz9tXXlV_VC#5a-XSk@TYpyP*rM1VAqxx8R(EF2%UQ+LCsGP5+rB&@ohZaGv zp(aCmOvv=4TXi&g_3-a4Y5f2-4Rt;8V}a5ep!L_N|3%#Y|M-`VKe`{)uS5RYX#MB< zsZG-AOKY!POG{6lz@w(sRi3Y>rB%0Y zq@~rSY@($nz1>Vps~;r4iD>!J@~CUhBzJy9)4*FzU8Qkv3q6)rjR~iv)hx)(k;EV= z?7o_vx2i-8jbze=SXx@0OD-W!KBA63B1K8Kpr{4kfA)H-5HO={1wp(ZDPv85!{RpJ*>ldALWgh^b{)jM-(X^j*QT3XY| Ji%65d{|^$F9uNQk diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index c602ba2c6f..1a3ddbf901 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -7,172 +7,24 @@ import h5py import numpy as np from .data import ATOMIC_SYMBOL, EV_PER_MEV -from .endf import get_cont_record, get_list_record, Evaluation +from .endf import get_cont_record, get_list_record, get_tab1_record, Evaluation from .function import Function1D, Tabulated1D, Polynomial, Sum import openmc.checkvalue as cv from openmc.mixin import EqualityMixin -def _extract_458_data(ev): - """Read an ENDF file and extract the MF=1, MT=458 values. - - Parameters - ---------- - ev : openmc.data.Evaluation - ENDF evaluation - - Returns - ------- - value : dict of str to list of float - Dictionary that gives lists of coefficients for each energy component. - The keys are the 2-3 letter strings used in ENDF-102, e.g. 'EFR' and - 'ET'. The list will have a length of 1 for Sher-Beck data, more for - polynomial data. - uncertainty : dict of str to list of float - A dictionary with the same format as above. This is probably a - one-standard deviation value, but that is not specified explicitly in - ENDF-102. Also, some evaluations will give zero uncertainty. Use with - caution. - - """ - cv.check_type('evaluation', ev, Evaluation) - - if not ev.target['fissionable']: - # This nuclide isn't fissionable. - return None - - if (1, 458) not in ev.section: - # No 458 data here. - return None - - file_obj = StringIO(ev.section[1, 458]) - - # Read the number of coefficients in this LIST record. - items = get_cont_record(file_obj) - NPL = items[3] - - # Parse the ENDF LIST into an array. - items, data = get_list_record(file_obj) - - # Declare the coefficient names and the order they are given in. The LIST - # contains a value followed immediately by an uncertainty for each of these - # components, times the polynomial order + 1. - labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') - - # Associate each set of values and uncertainties with its label. - value = {} - uncertainty = {} - for i, label in enumerate(labels): - value[label] = data[2*i::18] - uncertainty[label] = data[2*i + 1::18] - - # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not - # converted from MeV to eV. Check for this error and fix it if present. - n_coeffs = len(value['EFR']) - if n_coeffs == 3: # Only check 2nd-order data. - # Check each energy component for the error. If a 1 MeV neutron - # causes a change of more than 100 MeV, we know something is wrong. - error_present = False - for coeffs in value.values(): - second_order = coeffs[2] - if abs(second_order) * 1e12 > 1e8: - error_present = True - break - - # If we found the error, reduce all 2nd-order coeffs by 10**6. - if error_present: - for coeffs in value.values(): - coeffs[2] /= EV_PER_MEV - for coeffs in uncertainty.values(): - coeffs[2] /= EV_PER_MEV - - return value, uncertainty - - -def write_compact_458_library(endf_files, output_name='fission_Q_data.h5', - comment=None, verbose=False): - """Read ENDF files, strip the MF=1 MT=458 data and write to small HDF5. - - Parameters - ---------- - endf_files : Collection of str - Strings giving the paths to the ENDF files that will be parsed for data. - output_name : str - Name of the output HDF5 file. Default is 'fission_Q_data.h5'. - comment : str - Comment to write in the output HDF5 file. Defaults to no comment. - verbose : bool - If True, print the name of each isomer as it is read. Defaults to - False. - - """ - # Open the output file. - out = h5py.File(output_name, 'w', libver='earliest') - - # Write comments, if given. This commented out comment is the one used for - # the library distributed with OpenMC. - #comment = ('This data is extracted from ENDF/B-VII.1 library. Thanks ' - # 'evaluators, for all your hard work :) Citation: ' - # 'M. B. Chadwick, M. Herman, P. Oblozinsky, ' - # 'M. E. Dunn, Y. Danon, A. C. Kahler, D. L. Smith, ' - # 'B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, ' - # 'D. A. Brown, R. Capote, A. D. Carlson, Y. S. Cho, H. Derrien, ' - # 'K. Guber, G. M. Hale, S. Hoblit, S. Holloway, T. D. Johnson, ' - # 'T. Kawano, B. C. Kiedrowski, H. Kim, S. Kunieda, ' - # 'N. M. Larson, L. Leal, J. P. Lestone, R. C. Little, ' - # 'E. A. McCutchan, R. E. MacFarlane, M. MacInnes, ' - # 'C. M. Mattoon, R. D. McKnight, S. F. Mughabghab, ' - # 'G. P. A. Nobre, G. Palmiotti, A. Palumbo, M. T. Pigni, ' - # 'V. G. Pronyaev, R. O. Sayer, A. A. Sonzogni, N. C. Summers, ' - # 'P. Talou, I. J. Thompson, A. Trkov, R. L. Vogt, ' - # 'S. C. van der Marck, A. Wallner, M. C. White, D. Wiarda, ' - # 'and P. G. Young. ENDF/B-VII.1 nuclear data for science and ' - # 'technology: Cross sections, covariances, fission product ' - # 'yields and decay data", Nuclear Data Sheets, ' - # '112(12):2887-2996 (2011).') - if comment is not None: - out.attrs['comment'] = np.string_(comment) - - # Declare the order of the components. Use fixed-length numpy strings - # because they work well with h5py. - labels = np.array(('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', - 'ET'), dtype='S3') - out.attrs['component order'] = labels - - # Iterate over the given files. - if verbose: print('Reading ENDF files:') - for fname in endf_files: - if verbose: print(fname) - - ev = Evaluation(fname) - - # Skip non-fissionable nuclides. - if not ev.target['fissionable']: - continue - - # Get the important bits. - data = _extract_458_data(ev) - if data is None: continue - value, uncertainty = data - - # Make a group for this isomer. - name = ATOMIC_SYMBOL[ev.target['atomic_number']] + \ - str(ev.target['mass_number']) - if ev.target['isomeric_state'] != 0: - name += '_m' + str(ev.target['isomeric_state']) - nuclide_group = out.create_group(name) - - # Write all the coefficients into one array. The first dimension gives - # the component (e.g. fragments or prompt neutrons); the second switches - # between value and uncertainty; the third gives the polynomial order. - n_coeffs = len(value['EFR']) - data_out = np.zeros((len(labels), 2, n_coeffs)) - for i, label in enumerate(labels): - data_out[i, 0, :] = value[label.decode()] - data_out[i, 1, :] = uncertainty[label.decode()] - nuclide_group.create_dataset('data', data=data_out) - - out.close() +_LABELS = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') +_NAMES = { + 'EFR': 'fragments', + 'ENP': 'prompt_neutrons', + 'END': 'delayed_neutrons', + 'EGP': 'prompt_photons', + 'EGD': 'delayed_photons', + 'EB': 'betas', + 'ENU': 'neutrinos', + 'ER': 'recoverable', + 'ET': 'total' +} class FissionEnergyRelease(EqualityMixin): @@ -249,14 +101,15 @@ class FissionEnergyRelease(EqualityMixin): - incident neutron energy). """ - def __init__(self): - self._fragments = None - self._prompt_neutrons = None - self._delayed_neutrons = None - self._prompt_photons = None - self._delayed_photons = None - self._betas = None - self._neutrinos = None + def __init__(self, fragments, prompt_neutrons, delayed_neutrons, + prompt_photons, delayed_photons, betas, neutrinos): + self.fragments = fragments + self.prompt_neutrons = prompt_neutrons + self.delayed_neutrons = delayed_neutrons + self.prompt_photons = prompt_photons + self.delayed_photons = delayed_photons + self.betas = betas + self.neutrinos = neutrinos @property def fragments(self): @@ -345,92 +198,6 @@ class FissionEnergyRelease(EqualityMixin): cv.check_type('neutrinos', energy_release, Callable) self._neutrinos = energy_release - @classmethod - def _from_dictionary(cls, energy_release, incident_neutron): - """Generate fission energy release data from a dictionary. - - Parameters - ---------- - energy_release : dict of str to list of float - Dictionary that gives lists of coefficients for each energy - component. The keys are the 2-3 letter strings used in ENDF-102, - e.g. 'EFR' and 'ET'. The list will have a length of 1 for Sher-Beck - data, more for polynomial data. - incident_neutron : openmc.data.IncidentNeutron - Corresponding incident neutron dataset - - Returns - ------- - openmc.data.FissionEnergyRelease - Fission energy release data - - """ - out = cls() - - # How many coefficients are given for each component? If we only find - # one value for each, then we need to use the Sher-Beck formula for - # energy dependence. Otherwise, it is a polynomial. - n_coeffs = len(energy_release['EFR']) - if n_coeffs > 1: - out.fragments = Polynomial(energy_release['EFR']) - out.prompt_neutrons = Polynomial(energy_release['ENP']) - out.delayed_neutrons = Polynomial(energy_release['END']) - out.prompt_photons = Polynomial(energy_release['EGP']) - out.delayed_photons = Polynomial(energy_release['EGD']) - out.betas = Polynomial(energy_release['EB']) - out.neutrinos = Polynomial(energy_release['ENU']) - else: - # EFR and ENP are energy independent. Use 0-order polynomials to - # make a constant function. The energy-dependence of END is - # unspecified in ENDF-102 so assume it is independent. - out.fragments = Polynomial((energy_release['EFR'][0])) - out.prompt_photons = Polynomial((energy_release['EGP'][0])) - out.delayed_neutrons = Polynomial((energy_release['END'][0])) - - # EDP, EB, and ENU are linear. - out.delayed_photons = Polynomial((energy_release['EGD'][0], -0.075)) - out.betas = Polynomial((energy_release['EB'][0], -0.075)) - out.neutrinos = Polynomial((energy_release['ENU'][0], -0.105)) - - # Prompt neutrons require nu-data. It is not clear from ENDF-102 - # whether prompt or total nu value should be used, but the delayed - # neutron fraction is so small that the difference is negligible. - # MT=18 (n, fission) might not be available so try MT=19 (n, f) as - # well. - if 18 in incident_neutron.reactions: - nu = [p.yield_ for p in incident_neutron[18].products - if p.particle == 'neutron' - and p.emission_mode in ('prompt', 'total')] - elif 19 in incident_neutron.reactions: - nu = [p.yield_ for p in incident_neutron[19].products - if p.particle == 'neutron' - and p.emission_mode in ('prompt', 'total')] - else: - raise ValueError('IncidentNeutron data has no fission ' - 'reaction.') - if len(nu) == 0: - raise ValueError('Nu data is needed to compute fission energy ' - 'release with the Sher-Beck format.') - if len(nu) > 1: - raise ValueError('Ambiguous prompt/total nu value.') - - nu = nu[0] - if isinstance(nu, Tabulated1D): - ENP = deepcopy(nu) - ENP.y = (energy_release['ENP'] + 1.307 * nu.x - - 8.07e6 * (nu.y - nu.y[0])) - elif isinstance(nu, Polynomial): - if len(nu) == 1: - ENP = Polynomial([energy_release['ENP'][0], 1.307]) - else: - ENP = Polynomial( - [energy_release['ENP'][0], 1.307 - 8.07e6*nu.coef[1]] - + [-8.07e6*c for c in nu.coef[2:]]) - - out.prompt_neutrons = ENP - - return out - @classmethod def from_endf(cls, ev, incident_neutron): """Generate fission energy release data from an ENDF file. @@ -463,11 +230,111 @@ class FissionEnergyRelease(EqualityMixin): if not ev.target['fissionable']: raise ValueError('The ENDF evaluation is not fissionable.') - # Read the 458 data from the ENDF file. - value, uncertainty = _extract_458_data(ev) + if (1, 458) not in ev.section: + raise ValueError('ENDF evaluation does not have MF=1, MT=458.') - # Build the object. - return cls._from_dictionary(value, incident_neutron) + file_obj = StringIO(ev.section[1, 458]) + + # Read first record and check whether any components appear as tabulated + # functions + items = get_cont_record(file_obj) + lfc = items[3] + nfc = items[5] + + # Parse the ENDF LIST into an array. + items, data = get_list_record(file_obj) + npoly = items[3] + + # Associate each set of values and uncertainties with its label. + functions = {} + for i, label in enumerate(_LABELS): + coeffs = data[2*i::18] + name = _NAMES[label] + + # Ignore recoverable and total since we recalculate those directly + if name in ('recoverable', 'total'): + continue + + # In ENDF/B-VII.1, data for 2nd-order coefficients were mistakenly not + # converted from MeV to eV. Check for this error and fix it if present. + if npoly == 2: # Only check 2nd-order data. + # If a 1 MeV neutron causes a change of more than 100 MeV, we know + # something is wrong. + second_order = coeffs[2] + if abs(second_order) * 1e12 > 1e8: + # If we found the error, reduce 2nd-order coeff by 10**6. + coeffs[2] /= EV_PER_MEV + + # If multiple coefficients were given, we can create the polynomial + # and move on to the next component + if npoly > 0: + functions[name] = Polynomial(coeffs) + continue + + # If a single coefficient was given, we need to use the Sher-Beck + # formula for energy dependence + zeroth_order = coeffs[0] + if name in ('delayed_photons', 'betas'): + func = Polynomial((zeroth_order, -0.075)) + elif name == 'neutrinos': + func = Polynomial((zeroth_order, -0.105)) + elif name == 'prompt_neutrons': + # Prompt neutrons require nu-data. It is not clear from ENDF-102 + # whether prompt or total nu value should be used, but the delayed + # neutron fraction is so small that the difference is negligible. + # MT=18 (n, fission) might not be available so try MT=19 (n, f) as + # well. + if 18 in incident_neutron.reactions: + nu = [p.yield_ for p in incident_neutron[18].products + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] + elif 19 in incident_neutron.reactions: + nu = [p.yield_ for p in incident_neutron[19].products + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] + else: + raise ValueError('IncidentNeutron data has no fission ' + 'reaction.') + if len(nu) == 0: + raise ValueError('Nu data is needed to compute fission energy ' + 'release with the Sher-Beck format.') + if len(nu) > 1: + raise ValueError('Ambiguous prompt/total nu value.') + + nu = nu[0] + if isinstance(nu, Tabulated1D): + # Evaluate Sher-Beck polynomial form at each tabulated value + func = deepcopy(nu) + func.y = (zeroth_order + 1.307*nu.x - 8.07e6*(nu.y - nu.y[0])) + elif isinstance(nu, Polynomial): + # Combine polynomials + if len(nu) == 1: + func = Polynomial([zeroth_order, 1.307]) + else: + func = Polynomial( + [zeroth_order, 1.307 - 8.07e6*nu.coef[1]] + + [-8.07e6*c for c in nu.coef[2:]]) + else: + func = Polynomial(coeffs) + + functions[name] = func + + # Check for tabulated data + if lfc == 1: + for _ in range(nfc): + # Get tabulated function + items, eifc = get_tab1_record(file_obj) + + # Determine which component it is + ifc = items[3] + name = _NAMES[_LABELS[ifc - 1]] + + # Replace value in dictionary + functions[name] = eifc + + # Build the object + print(functions) + return cls(**functions) @classmethod def from_hdf5(cls, group): @@ -485,17 +352,16 @@ class FissionEnergyRelease(EqualityMixin): """ - obj = cls() + fragments = Function1D.from_hdf5(group['fragments']) + prompt_neutrons = Function1D.from_hdf5(group['prompt_neutrons']) + delayed_neutrons = Function1D.from_hdf5(group['delayed_neutrons']) + prompt_photons = Function1D.from_hdf5(group['prompt_photons']) + delayed_photons = Function1D.from_hdf5(group['delayed_photons']) + betas = Function1D.from_hdf5(group['betas']) + neutrinos = Function1D.from_hdf5(group['neutrinos']) - obj.fragments = Function1D.from_hdf5(group['fragments']) - obj.prompt_neutrons = Function1D.from_hdf5(group['prompt_neutrons']) - obj.delayed_neutrons = Function1D.from_hdf5(group['delayed_neutrons']) - obj.prompt_photons = Function1D.from_hdf5(group['prompt_photons']) - obj.delayed_photons = Function1D.from_hdf5(group['delayed_photons']) - obj.betas = Function1D.from_hdf5(group['betas']) - obj.neutrinos = Function1D.from_hdf5(group['neutrinos']) - - return obj + return cls(fragments, prompt_neutrons, delayed_neutrons, prompt_photons, + prompt_photons, delayed_photons, betas, neutrinos) @classmethod def from_compact_hdf5(cls, fname, incident_neutron): From 9d2f804ca09ca6591c378c85a66315c099c75101 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jul 2018 10:57:19 -0500 Subject: [PATCH 2/5] Refactor how recoverable/prompt/total are calculated using sum_functions --- openmc/data/fission_energy.py | 122 ++++++++-------------------------- openmc/data/function.py | 42 ++++++++++++ 2 files changed, 69 insertions(+), 95 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 1a3ddbf901..0f30ecfea0 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -6,25 +6,18 @@ import sys import h5py import numpy as np -from .data import ATOMIC_SYMBOL, EV_PER_MEV +from .data import EV_PER_MEV from .endf import get_cont_record, get_list_record, get_tab1_record, Evaluation -from .function import Function1D, Tabulated1D, Polynomial, Sum +from .function import Function1D, Tabulated1D, Polynomial, sum_functions import openmc.checkvalue as cv from openmc.mixin import EqualityMixin -_LABELS = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') -_NAMES = { - 'EFR': 'fragments', - 'ENP': 'prompt_neutrons', - 'END': 'delayed_neutrons', - 'EGP': 'prompt_photons', - 'EGD': 'delayed_photons', - 'EB': 'betas', - 'ENU': 'neutrinos', - 'ER': 'recoverable', - 'ET': 'total' -} +_NAMES = ( + 'fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas', + 'neutrinos', 'recoverable', 'total' +) class FissionEnergyRelease(EqualityMixin): @@ -141,27 +134,33 @@ class FissionEnergyRelease(EqualityMixin): @property def recoverable(self): - return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, - self.prompt_photons, self.delayed_photons, self.betas]) + components = ['fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas'] + return sum_functions(getattr(self, c) for c in components) @property def total(self): - return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, - self.prompt_photons, self.delayed_photons, self.betas, - self.neutrinos]) + components = ['fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas', + 'neutrinos'] + return sum_functions(getattr(self, c) for c in components) @property def q_prompt(self): - return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons, - lambda E: -E]) + # Use a polynomial to subtract incident energy. + funcs = [self.fragments, self.prompt_neutrons, self.prompt_photons, + Polynomial((0.0, -1.0))] + return sum_functions(funcs) @property def q_recoverable(self): - return Sum([self.recoverable, lambda E: -E]) + # Use a polynomial to subtract incident energy. + return sum_functions([self.recoverable, Polynomial((0.0, -1.0))]) @property def q_total(self): - return Sum([self.total, lambda E: -E]) + # Use a polynomial to subtract incident energy. + return sum_functions([self.total, Polynomial((0.0, -1.0))]) @fragments.setter def fragments(self, energy_release): @@ -247,9 +246,8 @@ class FissionEnergyRelease(EqualityMixin): # Associate each set of values and uncertainties with its label. functions = {} - for i, label in enumerate(_LABELS): + for i, name in enumerate(_NAMES): coeffs = data[2*i::18] - name = _NAMES[label] # Ignore recoverable and total since we recalculate those directly if name in ('recoverable', 'total'): @@ -327,13 +325,12 @@ class FissionEnergyRelease(EqualityMixin): # Determine which component it is ifc = items[3] - name = _NAMES[_LABELS[ifc - 1]] + name = _NAMES[ifc - 1] # Replace value in dictionary functions[name] = eifc # Build the object - print(functions) return cls(**functions) @classmethod @@ -361,44 +358,7 @@ class FissionEnergyRelease(EqualityMixin): neutrinos = Function1D.from_hdf5(group['neutrinos']) return cls(fragments, prompt_neutrons, delayed_neutrons, prompt_photons, - prompt_photons, delayed_photons, betas, neutrinos) - - @classmethod - def from_compact_hdf5(cls, fname, incident_neutron): - """Generate fission energy release data from a small HDF5 library. - - Parameters - ---------- - fname : str - Path to an HDF5 file containing fission energy release data. This - file should have been generated form the - :func:`openmc.data.write_compact_458_library` function. - incident_neutron : openmc.data.IncidentNeutron - Corresponding incident neutron dataset - - Returns - ------- - openmc.data.FissionEnergyRelease or None - Fission energy release data for the given nuclide if it is present - in the data file - - """ - - fin = h5py.File(fname, 'r') - - components = [s.decode() for s in fin.attrs['component order']] - - nuclide_name = ATOMIC_SYMBOL[incident_neutron.atomic_number] - nuclide_name += str(incident_neutron.mass_number) - if incident_neutron.metastable != 0: - nuclide_name += '_m' + str(incident_neutron.metastable) - - if nuclide_name not in fin: return None - - data = {c: fin[nuclide_name + '/data'][i, 0, :] - for i, c in enumerate(components)} - - return cls._from_dictionary(data, incident_neutron) + delayed_photons, betas, neutrinos) def to_hdf5(self, group): """Write energy release data to an HDF5 group @@ -417,33 +377,5 @@ class FissionEnergyRelease(EqualityMixin): self.delayed_photons.to_hdf5(group, 'delayed_photons') self.betas.to_hdf5(group, 'betas') self.neutrinos.to_hdf5(group, 'neutrinos') - - if isinstance(self.prompt_neutrons, Polynomial): - # Add the polynomials for the relevant components together. Use a - # Polynomial((0.0, -1.0)) to subtract incident energy. - q_prompt = (self.fragments + self.prompt_neutrons + - self.prompt_photons + Polynomial((0.0, -1.0))) - q_prompt.to_hdf5(group, 'q_prompt') - q_recoverable = (self.fragments + self.prompt_neutrons + - self.delayed_neutrons + self.prompt_photons + - self.delayed_photons + self.betas + - Polynomial((0.0, -1.0))) - q_recoverable.to_hdf5(group, 'q_recoverable') - - elif isinstance(self.prompt_neutrons, Tabulated1D): - # Make a Tabulated1D and evaluate the polynomial components at the - # table x points to get new y points. Subtract x from y to remove - # incident energy. - q_prompt = deepcopy(self.prompt_neutrons) - q_prompt.y += self.fragments(q_prompt.x) - q_prompt.y += self.prompt_photons(q_prompt.x) - q_prompt.y -= q_prompt.x - q_prompt.to_hdf5(group, 'q_prompt') - q_recoverable = q_prompt - q_recoverable.y += self.delayed_neutrons(q_recoverable.x) - q_recoverable.y += self.delayed_photons(q_recoverable.x) - q_recoverable.y += self.betas(q_recoverable.x) - q_recoverable.to_hdf5(group, 'q_recoverable') - - else: - raise ValueError('Unrecognized energy release format') + self.q_prompt.to_hdf5(group, 'q_prompt') + self.q_recoverable.to_hdf5(group, 'q_recoverable') diff --git a/openmc/data/function.py b/openmc/data/function.py index 3d09e44fc8..f7636f25cb 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,5 +1,7 @@ from abc import ABCMeta, abstractmethod from collections.abc import Iterable, Callable +from functools import reduce +from itertools import zip_longest from numbers import Real, Integral import numpy as np @@ -13,6 +15,46 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} +def sum_functions(funcs): + """Add tabulated/polynomials functions together + + Parameters + ---------- + funcs : list of Function1D + Functions to add + + Returns + ------- + Function1D + Sum of polynomial/tabulated functions + + """ + # Copy so we can iterate multiple times + funcs = list(funcs) + + # Get x values for all tabulated components + xs = [] + for f in funcs: + if isinstance(f, Tabulated1D): + xs.append(f.x) + if not np.all(f.interpolation == 2): + raise ValueError('Only linear-linear tabulated functions ' + 'can be combined') + + if xs: + # Take the union of all energies (sorted) + x = reduce(np.union1d, xs) + + # Evaluate each function and add together + y = sum(f(x) for f in funcs) + return Tabulated1D(x, y) + else: + # If no tabulated functions are present, we need to combine the + # polynomials by adding their coefficients + coeffs = [sum(x) for x in zip_longest(*funcs, fillvalue=0.0)] + return Polynomial(coeffs) + + class Function1D(EqualityMixin, metaclass=ABCMeta): """A function of one independent variable with HDF5 support.""" @abstractmethod From ef9070ce76dc664e7112ff764944f733e8ee6474 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jul 2018 11:43:53 -0500 Subject: [PATCH 3/5] Update tests/scripts --- scripts/openmc-ace-to-hdf5 | 16 ---------------- tests/unit_tests/test_data_neutron.py | 9 --------- 2 files changed, 25 deletions(-) diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 index 29c987e07c..1933a6e4a7 100755 --- a/scripts/openmc-ace-to-hdf5 +++ b/scripts/openmc-ace-to-hdf5 @@ -25,13 +25,6 @@ follows the NNDC data convention (1000*Z + A + 300 + 100*m), or the MCNP data convention (essentially the same as NNDC, except that the first metastable state of Am242 is 95242 and the ground state is 95642). -The optional --fission_energy_release argument will accept an HDF5 file -containing a library of fission energy release (ENDF MF=1 MT=458) data. A -library built from ENDF/B-VII.1 data is released with OpenMC and can be found at -openmc/data/fission_Q_data_endb71.h5. This data is necessary for -'fission-q-prompt' and 'fission-q-recoverable' tallies, but is not needed -otherwise. - """ class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, @@ -55,8 +48,6 @@ parser.add_argument('--xsdir', help='MCNP xsdir file that lists ' 'ACE libraries') parser.add_argument('--xsdata', help='Serpent xsdata file that lists ' 'ACE libraries') -parser.add_argument('--fission_energy_release', help='HDF5 file containing ' - 'fission energy release data') parser.add_argument('--libver', choices=['earliest', 'latest'], default='earliest', help="Output HDF5 versioning. Use " "'earliest' for backwards compatibility or 'latest' for " @@ -142,13 +133,6 @@ for filename in ace_libraries: print('Failed to convert {}: {}'.format(table.name, e)) continue - # Fission energy release data, if available - if args.fission_energy_release is not None: - fer = openmc.data.FissionEnergyRelease.from_compact_hdf5( - args.fission_energy_release, neutron) - if fer is not None: - neutron.fission_energy = fer - print('Converting {} (ACE) to {} (HDF5)'.format(table.name, neutron.name)) diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 03746430d9..9682b2e8f3 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -115,15 +115,6 @@ def test_fission_energy(pu239): assert isinstance(getattr(fer, c), Callable) -def test_compact_fission_energy(tmpdir): - files = [os.path.join(_ENDF_DATA, 'neutrons', 'n-090_Th_232.endf'), - os.path.join(_ENDF_DATA, 'neutrons', 'n-094_Pu_240.endf'), - os.path.join(_ENDF_DATA, 'neutrons', 'n-094_Pu_241.endf')] - output = str(tmpdir.join('compact_lib.h5')) - openmc.data.write_compact_458_library(files, output) - assert os.path.exists(output) - - def test_energy_grid(pu239): assert isinstance(pu239.energy, Mapping) for temp, grid in pu239.energy.items(): From d9afa804b86848af02542bc85dbf009412ae509c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 20 Jul 2018 14:27:29 -0500 Subject: [PATCH 4/5] Updated data format documentation --- docs/source/io_formats/nuclear_data.rst | 36 ++++++++++++------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 2e553a4edf..cfa5173909 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -85,33 +85,31 @@ temperature-dependent data set. For example, the data set corresponding to **//fission_energy_release/** -:Datasets: - **fragments** (:ref:`polynomial <1d_polynomial>`) -- Energy +:Datasets: - **fragments** (:ref:`function <1d_functions>`) -- Energy released in the form of fragments as a function of incident neutron energy. - - **prompt_neutrons** (:ref:`polynomial <1d_polynomial>` or - :ref:`tabulated <1d_tabulated>`) -- Energy released in the form of - prompt neutrons as a function of incident neutron energy. - - **delayed_neutrons** (:ref:`polynomial <1d_polynomial>`) -- Energy + - **prompt_neutrons** (:ref:`function <1d_functions>`) -- Energy + released in the form of prompt neutrons as a function of incident + neutron energy. + - **delayed_neutrons** (:ref:`function <1d_functions>`) -- Energy released in the form of delayed neutrons as a function of incident neutron energy. - - **prompt_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy + - **prompt_photons** (:ref:`function <1d_functions>`) -- Energy released in the form of prompt photons as a function of incident neutron energy. - - **delayed_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy + - **delayed_photons** (:ref:`function <1d_functions>`) -- Energy released in the form of delayed photons as a function of incident neutron energy. - - **betas** (:ref:`polynomial <1d_polynomial>`) -- Energy - released in the form of betas as a function of incident - neutron energy. - - **neutrinos** (:ref:`polynomial <1d_polynomial>`) -- Energy - released in the form of neutrinos as a function of incident - neutron energy. - - **q_prompt** (:ref:`polynomial <1d_polynomial>` or - :ref:`tabulated <1d_tabulated>`) -- The prompt fission Q-value - (fragments + prompt neutrons + prompt photons - incident energy) - - **q_recoverable** (:ref:`polynomial <1d_polynomial>` or - :ref:`tabulated <1d_tabulated>`) -- The recoverable fission Q-value - (Q_prompt + delayed neutrons + delayed photons + betas) + - **betas** (:ref:`function <1d_functions>`) -- Energy released in + the form of betas as a function of incident neutron energy. + - **neutrinos** (:ref:`function <1d_functions>`) -- Energy released + in the form of neutrinos as a function of incident neutron energy. + - **q_prompt** (:ref:`function <1d_functions>`) -- The prompt fission + Q-value (fragments + prompt neutrons + prompt photons - incident + energy) + - **q_recoverable** (:ref:`function <1d_functions>`) -- The + recoverable fission Q-value (Q_prompt + delayed neutrons + delayed + photons + betas) ------------------------------- Thermal Neutron Scattering Data From 95108ce311cb3f01af92b6cdbde8faccd3761805 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Aug 2018 22:01:54 -0500 Subject: [PATCH 5/5] Bump up test from 1 MeV -> 5 MeV for 458 bug --- openmc/data/fission_energy.py | 45 +++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 0f30ecfea0..fed4b21a39 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -224,8 +224,8 @@ class FissionEnergyRelease(EqualityMixin): raise ValueError('The atomic mass of the ENDF evaluation does ' 'not match the given IncidentNeutron.') if ev.target['isomeric_state'] != incident_neutron.metastable: - raise ValueError('The metastable state of the ENDF evaluation does ' - 'not match the given IncidentNeutron.') + raise ValueError('The metastable state of the ENDF evaluation ' + 'does not match the given IncidentNeutron.') if not ev.target['fissionable']: raise ValueError('The ENDF evaluation is not fissionable.') @@ -234,8 +234,8 @@ class FissionEnergyRelease(EqualityMixin): file_obj = StringIO(ev.section[1, 458]) - # Read first record and check whether any components appear as tabulated - # functions + # Read first record and check whether any components appear as + # tabulated functions items = get_cont_record(file_obj) lfc = items[3] nfc = items[5] @@ -253,13 +253,14 @@ class FissionEnergyRelease(EqualityMixin): if name in ('recoverable', 'total'): continue - # In ENDF/B-VII.1, data for 2nd-order coefficients were mistakenly not - # converted from MeV to eV. Check for this error and fix it if present. + # In ENDF/B-VII.1, data for 2nd-order coefficients were mistakenly + # not converted from MeV to eV. Check for this error and fix it if + # present. if npoly == 2: # Only check 2nd-order data. - # If a 1 MeV neutron causes a change of more than 100 MeV, we know - # something is wrong. + # If a 5 MeV neutron causes a change of more than 100 MeV, we + # know something is wrong. second_order = coeffs[2] - if abs(second_order) * 1e12 > 1e8: + if abs(second_order) * (5e6)**2 > 1e8: # If we found the error, reduce 2nd-order coeff by 10**6. coeffs[2] /= EV_PER_MEV @@ -277,25 +278,27 @@ class FissionEnergyRelease(EqualityMixin): elif name == 'neutrinos': func = Polynomial((zeroth_order, -0.105)) elif name == 'prompt_neutrons': - # Prompt neutrons require nu-data. It is not clear from ENDF-102 - # whether prompt or total nu value should be used, but the delayed - # neutron fraction is so small that the difference is negligible. - # MT=18 (n, fission) might not be available so try MT=19 (n, f) as - # well. + # Prompt neutrons require nu-data. It is not clear from + # ENDF-102 whether prompt or total nu value should be used, but + # the delayed neutron fraction is so small that the difference + # is negligible. MT=18 (n, fission) might not be available so + # try MT=19 (n, f) as well. if 18 in incident_neutron.reactions: nu = [p.yield_ for p in incident_neutron[18].products - if p.particle == 'neutron' - and p.emission_mode in ('prompt', 'total')] + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] elif 19 in incident_neutron.reactions: nu = [p.yield_ for p in incident_neutron[19].products - if p.particle == 'neutron' - and p.emission_mode in ('prompt', 'total')] + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] else: raise ValueError('IncidentNeutron data has no fission ' - 'reaction.') + 'reaction.') if len(nu) == 0: - raise ValueError('Nu data is needed to compute fission energy ' - 'release with the Sher-Beck format.') + raise ValueError( + 'Nu data is needed to compute fission energy ' + 'release with the Sher-Beck format.' + ) if len(nu) > 1: raise ValueError('Ambiguous prompt/total nu value.')