From 0d901792925c33123af2a0f1b31dfc9675e6593a Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Fri, 7 Oct 2022 11:20:04 +0200 Subject: [PATCH] Organize in submodules --- __init__.py | 0 __pycache__/comparator.cpython-37.pyc | Bin 16812 -> 0 bytes aggregator.py | 2 +- baseprocessor.py | 16 +- fragdisk/__init__.py | 0 gui.py => fragdisk/gui.py | 2 +- galactica/__init__.py | 0 .../fragdisk_galactica.py | 0 .../ismfeed_galactica.py | 0 .../ramses_astrophysix.py | 0 galturb/__init__.py | 0 galaxy.py => galturb/galaxy.py | 0 .../sectors_extraction.py | 0 ism.py | 101 ----- ismfeed/__init__.py | 0 ismfeed.py => ismfeed/ismfeed.py | 0 pipeline.py | 2 +- plotter.py | 24 +- pspec_read.py | 374 ------------------ run_selector.py => runselector.py | 0 scripts/__init__.py | 0 pspec_new.py => scripts/pspec.py | 0 select_snapshot.py | 2 +- snapshotprocessor.py | 16 +- studyprocessor.py | 40 +- turbox/__init__.py | 0 turbox/turbox.py | 160 ++++++++ utils/__init__.py | 0 .../line_annotation.py | 0 mypool.py => utils/mypool.py | 0 params.py => utils/params.py | 2 +- units.py => utils/units.py | 0 32 files changed, 231 insertions(+), 510 deletions(-) create mode 100644 __init__.py delete mode 100644 __pycache__/comparator.cpython-37.pyc create mode 100644 fragdisk/__init__.py rename gui.py => fragdisk/gui.py (99%) create mode 100644 galactica/__init__.py rename fragdisk_galactica.py => galactica/fragdisk_galactica.py (100%) rename ismfeed_galactica.py => galactica/ismfeed_galactica.py (100%) rename ramses_astrophysix.py => galactica/ramses_astrophysix.py (100%) create mode 100644 galturb/__init__.py rename galaxy.py => galturb/galaxy.py (100%) rename sectors_extraction.py => galturb/sectors_extraction.py (100%) delete mode 100644 ism.py create mode 100644 ismfeed/__init__.py rename ismfeed.py => ismfeed/ismfeed.py (100%) delete mode 100644 pspec_read.py rename run_selector.py => runselector.py (100%) create mode 100644 scripts/__init__.py rename pspec_new.py => scripts/pspec.py (100%) create mode 100644 turbox/__init__.py create mode 100644 turbox/turbox.py create mode 100644 utils/__init__.py rename line_annotation.py => utils/line_annotation.py (100%) rename mypool.py => utils/mypool.py (100%) rename params.py => utils/params.py (94%) rename units.py => utils/units.py (100%) diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/__pycache__/comparator.cpython-37.pyc b/__pycache__/comparator.cpython-37.pyc deleted file mode 100644 index e38ea189721c34cadf48dcbd571e0bebb6b983cd..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16812 zcmb7sYmgk*bzZ+_rl)6icCfoxEFOeL3=Lu-5(~h53!*`gqDVztQ2;E`BO=4W+}_3P z&I`ADiN&ZlwzaefDuF^jV*fCbES!hqvht{kl1ftgk&>$%MOASgcIDAOQmL|%82;f@ z5~oUqRK@we(>*=23xJWaHMeiyKKI_f=braHjq_7eo`HX#|KTsb`_n6i@%y|P{mY?p z4!__RP$)y09iyIH&AKV~R^7(c>Nwq8J!kTn-N|Xz$e zspmUW-AcV8^+IR5Tdh~6?soQdXX-OjFLq|TbM-mX_?)3UReIl0CBN)@+fIGIGA|m7 z<-2$sD{;KwaV9|=53mY9jo_g(-SD${l-wT`Vp05e?(kdWU zn_-`KQ*ZWz@J-!sp^r3_E^PEJ`W?TOJg+E!rMb}w8#O-`NI&aLujdS8^X1E2grQ^M0{jP&t)<->AFNR#5KyZoR0US4HK$Z`M8a zf-0*i)Jp2pYFbrME330=M$Mu&r9PwP)PB?|>P2-xEuc274yr>Kuc{8KBe?EUkEo-# z&ZtkQ8m_bIlj>1i=hUB3kKwvs{g^t2>%4kg9mn;6`U$nDPGF=3^(plPS`Mlw)k$0r zsU>v^*TY!HY5Yo9%Mo>EC9i&5Jq7qj)H(IEdIq(lvX*Dnb6CqKJ}?}kUc)xM6niiC zyTl#PIXvyv@b8s=t(CO3gMP1e#SgFhey`Tw2nQQsP;2&7O>gu9fn9DDQ~XwOWxa%7 z@EmX*8lf2(+t!X1nUNFaqQV`=G$LcoUUSxR%3QaCr<66iKVWQmJx7^aUgRkIfDw5I z4CUZzH1bhi<jOD~CPthMyJOaxY%ielzrgT6op3U11%nOZcx@Y zHe$0e%q?*d+7)@2uzp>``d<(2Vg|;gEV7 z&B=SX02860_p@yu;q&X)n|d{fUA$g{v&A;LrwBl;ULK)>!^Sp-tn-CWL5vg|jaH`_ z1dYaT823K4bhY30mwH!#`x{HE?Qlso2TOzY0Po-P1;YsoogCcIvzU4wzu-9(Ms?n< zm{rTMW^C7-v-75Fx+wEj+2o&RRZT1XbIe1wXCAR!d&VBlj`{NB2;(iABjO*Wzk?s| zgCUB@3@VXz(_Ayw3wnkUSjpWE80L<-W})>sTM030TeCkjLMIb@E}q%17~AHy33hX~ z?8qJrLWso=%oR)JAs8nhww%bp90ic8JA6GOxeGohvlrAZ=?%ZuUa2)fqt`{6Y=q6Q z4QbL}5nAQu*4iTy&iY@{2eFtIDUkmx{J2bFG+02f3ck8r-7+KdCKu{N=FeN#@FMyI z9!dedVx2stkKrDh=mpL9e3g)jfP4z=_q?;DI@0%x@8Uk-G*FXEEy0C;k_{OO^hw-OO^()n_D{X`y=Va9 zZLq&=;G!unfVpMb8YV&$_e668#ceXB1s!p|IT-jo6`yaDvp2$=^#0ZM;bZUSGo z{Xj@DHrv7>W7(0=Ls_(7BZ2IrS5QG=fDX!#-3y=t;OW35o~Sib_)EJ`#qAeVy^YE$ zQT19#%K#!4D0u-Wakg@irOyhgaz$iYoK$xtx| zbX3e4YsMT_CsEOYz`R#jN|4v$`uUGrM%zJf@fxHW2?lBiRv@e~{NxI9QPx4@fN{s% zvF`v6p#zI42X$nxJNgpt0I_YZnjtJgSdd_pp9B0h`fn9_-st&OAuQ~;ur|#tH>(x1 znwQl|S*@Jarm|Wkt4(LM>PP10AcUpK_jtwl$h`b6_Hf_U%mddmldWVndrpgX+&+hQ zo=sP1eq`$JgkCheZH9%l(z?C1Uu<1x#fCbzJ~D%E^WHMH=C=++2b8hyqQ?MMtki-( z7tN<_SED(V2OS;gxzT)t838|f_RAx9{Y^jtdt38>IuJPvs73puxoCbpch11|z`A+P zz{<>bqbHnY)wz|Z3CrlMK+8iiX{0t$BRsKrZLrlrND&!kH9msShX$9|^} zl?Vf?n&q^VMJcONPDxp2dF{8ww&pxuR7AhT7N3r)s&d-{?N2Xja;GrBT`^qqOJqmb zxy#GTcmD+cAUT}oW;<8~9R4}y?$U0(OR;(P5{jh#-M7!(Rk&O95&`3M{U$8JAXIU_ z>#H_zx@{KC&Dgni>dez9pMUo0*uHl985GYwwOG|-p)5kmf$^>q#gz-LHV1y}&~y-% zjB`U@_XA;BNVY3MY{N2+-FC1->)O{}#LMV6QN$ih@0Q;SVO=jSBtp#*V>TB+8^HLc z;B$fjUNN_-`x}FLUij~n{wa)@riseh8~#Db>I({HiJgh4$gyzZr9;_yxpJ0gKGsSraZcBrkW*+lGczsZIXh(5AK%;vx_ zMSTn`0TBXy1yKWLfzW}=3e)nAl?c}<^jx&{AzVlhI-a79Ler#uWh3-Qil)dG@iNAe z6qX50nWT{lNCw%9GckGC`MZD!=sOuem}>c@C^fX`cUg>;!$a3mW>0vO()oMnw0kX+DsvL@ z)oj1G^lU-raFZ@!AzcFQ(V!#$$|jVo*mi2g!E-vp{r7U;;@k}WK6^g8{f zYP9@LC(ze0&7zwul*)eZLIbwR*5_nhdpFPHQdT?tE6vXD05nyQ+&o&2q~%>U)fhT6 z6u4rj#AF*O%(88SF@QC|5;Qf9f6!nJm$=aDb{gycjpS9un#)YOxWE|z@b3XB#=T(< zXC_!BecTd@4=^MmTB<@s=TOOJCy&6(TLbw}MX#F5S{0j4*|%(9+Sx&{Oq6#H3Os)Y z_FcsAMPZqd08oNjX~gCbFaH#V9ePVwA^Vo3)oKh8Qt0#t{-Tp0w*WJQA-;Gjv+5hl zALuT+K7<7F(5KRJOl|dT-bw-3^YTEhd%rXQvL`PW$r+5fUKO* zq;N|)0Se^Qe+EaQzslFEQRG2Yu(d7Fu@>-Q{fKJd{%43t^8iAztMxb z!))zd4S`5HZ=CHoCmW#o3jFSMv@FJjw_dui{QC0e7AyJ@4sP{>I5z+l_2ME`En_ig z(aqkfFKRc@v{0-e$ZgPRc|0*!wc2oV-Wuz1~y)W>R+=2t&aRr~yd7M94DZ zQdxA1qA+$_UDe<##)TGw*HF)i-lW>o+bp=uy`*`d;ctez*$M>{L1RVtyPWMm12(9k zFd#50<~&?tIvPjddd|SzfCGxZyj8V^k4_NpLk7#pn4TK>TgdoZs63R6zd0u3Z3{#k z!T6peVfbh?7RSWAqAg}*JR|snBHL{ zaOb6WkS@3@rFW1lI6(5~Z=toEwjxSC(K?0gEvMUCibMjTWh!eKNvTTK;-zwn@p4Ks zMBvkwYSvOtTlQrwQ*6n|WM+9dmypS;V1IoJQW4fSe0dBKTPlprCIpw?O>Bjv%G0c1vHbD078f8HI?a}!keX25hnz&n zPbjbQqfp*+?DMOmx6`*-`zt8ulXM5I)9LlhUcx+_5oWfYj6DAY{ewqQjL35i5-0ha z7g-73N@l=^~^z^}t= z0@;Dk6Xh6|;KCNhP$07@SSuKpHApeUqwc_B6qP!XMIL(D?3F>2JoiFp3V1X}2K}J1 zGv$0AwD*9q!CJfDKN}k?cHf(fn2hf&tIeRnYs#_j@JXRaArx%TwdS_3AV{Acoft_x)j6T;`^ zt}PGzK|d7H5j*~RyLXq?vi=nmR15MC4`^aIeix@;1;5iAz<9c6z98J$*Ee9VP#wg% zb`Pfsb&m(te$Z$KT`fZ4(|m}dW}!}j+>lSKDIj~F?FRyeCB9BgS*fteF=a^vmM+&pLA{T%Ko&WG1ZYh?;EWvv^* z4{_2AaTmE#uY3(oSEfz#$uqt!Z2Bi)(^taju)1Qedtl~ssa-z_Us&5mXxQ~nz^;EH zWt6Gq;gN(#a7=dvM{vEcrn`aY%O_v1ofh*-Xno8XiRu*1^Xu3aoTZ{`1HBS9`nzl^ zNt=)Sjnx~AGas)S#qQce&nbJziEn-X#9)o^Dh|4@_B%?y!-*;{1qey;NJC6@{TI=w zZ}ExUYxaILYwjNsyf_~Cp*9`XnDM|F@&KLiL}?xpnm>W4)DiH+f_Vsv^XNnDyUAfP zr=PY*X5alGZsZ_kE;3KGJvYI&|26$3PsC2jARU#@novfZF|eUHo#2BB6_e#;pNi!m!p_S zEXGkCbA1DI*%a7(LN{QmxcQj$JA!aQ)`zNdbeHByJj$L35I9N9Ha1TPm0{NyE!78D zP)_QT@kXcfY-R7a18y44~d{u51d!3F5NI8WlVV`I=Spo%{dXD$WtZ?2k7@AYIz`)7| z3?mONa&P{zDk4AUv9Pi>jZs`mxyW2V(w8dVvSBV%BNwvE!zlve^qz9*DMZMLiM*1W z(kvM0H)?k$+YcKsAOte)gLM&_Fwz3$&lbRy!^+M+%vIb5=9p82^Mg`UIE-MFU`EDA z8%&^4TQUkL@Fb|uQSedPNM`^QR)&vdX$OcLH?D&m`qvqRxZyW-_GoxGJ0yD$#1h$QNr^ao)JY2w` zLf&40xq_&QbsWaYoVj2PKQU>XJYdu;52JvAvmoi$ag%b+h<{D+tF_rdE4cOY6{OCP zS1BXOtGF^$BD{)#su;T?UUk5$8!BGu@B&S6P0`g<%$GWz={MrY0c5_JS;VCWv>BksHYS1h>Qdv|uW zkx)dtCDAdQPX)o(xv_HqHe_4N+}3nU5ez%JYXgOxCOc>Q$Gyw;?+duWewmx+iFfSR zb0U8c*M+MIJqb6Ir_bC&wpkNx*uBWcE&vX@z~lr*OYP<1!iZrfmHm&b+P zvxE@1Fyt8LMB^yS#S_p$|GJ;@4lS`bpV(S^F<5Ghp?{f}V;CvPT$r=K2+&r;`AM9O z!82S@S7k2z4{(EcF>G)|eGvRdjZtumQE<5P;?NThrhOm6cVIPXsT0PjPb zo#&yc5^EiHIqqF?=)pZ}`L)72ay+pMtivik2c3XfButpLdUM9QgR^P4tJFFNafk=^ z7kYy>__p;9(@FBs2RPC@3oV13Pj~{@xwpwbv`!FO$lXT%*kXRz%tZcI;6<$<{p+%X z>u|y+8EbNwFiKy$*k0{6YdAQlHIWcoTiNKfCbQe*q=>M{ogI%0q2EG!#1iGRNXQoK zTE7>_Ny!kLEKhm3)9Um$RGOK(V7(eET=@aWd-)Hn{YMt(S&+@+ysRHtnmo1kHS#aP%Ih*&sRio{9x|~dM{|2gxZDV2 z!Tlc4^lPNyc?@YdI6LN9Fj7wD&8qVh0_3Nx-B~yj{f6#RawHe%GCH?N&-}+L{mYl} z3*JP51SDin=xTDr1-b^cCy7f^A_bC`vmA0EmU_$HwjibRNQOk4Dzq)-NkStF2t|5cSub*t$>H+5Xs&y5FDU@=RIqVH@}QFEe2YZ+KhN48ow~oV>NB&p z<0pifAn%`K*e@Yx6&U9qPvYlo)^XaoeWpooAZAowUz zcGJSg7f5MsBkb`j);H5k<(p8yl9?e%UV9+~#rAT-KvY}E?Oex2&iYaRH0zJC-9}ae z6Jr{IZi$^5yQn3cA*>*lrM1%gW0Jf`3mo0gBa}6KY*JljM2>}x^MO}-X_7rOx#FB+ zWDB4D7tYS?R^i*dPbKb0l#Kk?DPh6v)yOIOT@-A zz~oW-t_gPSI_)=M{R+wu>Cho|8p(cec$K|sqmH9<<)~bU4@2;SfW2)ie0GE6G1k0W zIraivL6z{_m2@&wRc__ec7)#2zMh(%1&OID63#M+Y`%S~7Dg#?Ho%u#OMS+1zOExwCZ@YbMKKs6lXRqF})#q+k zZW@gPSmT0vO}(x@ul_7^RBjqO3+f9DgvygQP##pv2YGi;eGzWifmxg@0oF<(-c(;o z?lu?H1(>M#)Bs}~Vn&gPv!`p{dW&NoQd2O4rH|C@tbS4L%j##;{)zfL>*3+>$W49{ z^$2s7M1Eab9{a!gYhaYIt*{Sw;O5LdmdCMk3CVcF*-QNp3WP2aB|MYzg5mT;W(FH3 zjSCSL!}*b)vOAbDJn?!DL*N^{p5NaHY9l{P#^v`#_&A3jZ!q6FuI%pnV4;`znll$u zpG-1Edhj5lv(a|6N*h`5PR3AHSa`yRcbkGdo>R)H^8VW4(w{VH4B*yZdk8F{4 zcq2xZ7wlM?*SNzw%6#;lbT&raB`SIDlQ>1N^i_=>Q4$!x@W#dA!Pk3h_^4=Cd+K`q z4UY9US^O;)f1AZ`u=q_D-(vB1So}Q{knb}8pQ?{Fp)P)vBk|x@-(m4<95bJcDUW1f z^8Byz`CmgZoEz=+q~|H;-9QIXrvE34|HXo&He7hXV#V={%acP4_l=BeiGfKLxhpm= z^2v9R9JL>yk=n)FnI(O-kAQA}MgIWpDG5NcrLAaCe3i6QmMoS!Nhy52#B8m!O+1j{ z>^QM(6&tGFa~{dG2&`n}{*WBb@1anE*|-7Bj3;&)&= zNiHAkk;xPJ>?^skEi&qYL5WSAJpK<{aPla@w|~oeISZe`x6^;jrnv{5$i*}FAN1?r z1Gvy#JSB)3?Ot5H4^Q*b!tD~*d9-P9X~I#{Kjieq?)4@@K5(gHPdZ6m$P??JH5Km@ zAGL!oAG2|oj~pcAD*`}9Igf7>;nE88<1`jba*uBeU0`U4Tq9Y)|79EVAoTx1v3PXC z)A}x3e}RR^;`=D@^u1ed;;p83Pkbjf)WTd25^4yitJb7Qwt{@WZVWy*2#q!*0 z^0_5TYu8~-#zT^4dkKeA%gg$6M7bp{Jd$shkwy&esnP=VoWydReo_Eg6Bm=QTFAJcq?o!OqU^KmkdsJ_(i^p2>8S=A{|w6 Jo{;|8{|`hZzQ6zg diff --git a/aggregator.py b/aggregator.py index 1d43d2a..e3f349d 100644 --- a/aggregator.py +++ b/aggregator.py @@ -4,7 +4,7 @@ import numpy as np from functools import partial from baseprocessor import Rule import snapshotprocessor -from mypool import MyPool +from utils.mypool import MyPool try: from mpi4py.futures import MPIPoolExecutor diff --git a/baseprocessor.py b/baseprocessor.py index 1fb8e9b..1f21f50 100644 --- a/baseprocessor.py +++ b/baseprocessor.py @@ -10,11 +10,11 @@ from functools import partial import numpy as np import tables -from tables import HDF5ExtError +from tables import HDF5ExtError, NoSuchNodeError from tables.registry import class_name_dict -from params import default_params, load_params -from units import U +from utils.params import default_params, load_params +from utils.units import U import sys import traceback @@ -44,6 +44,7 @@ class Rule: return self.process_fn(arg, **kwargs) else: return self.process_fn(**kwargs) + class BaseProcessor: """ Base class for processors, should not be instanciated @@ -292,8 +293,17 @@ class HDF5Container(BaseProcessor): ) else: value = node.read() + if isinstance(unit, dict): + name = os.path.basename(node_name) + if name in unit: + unit = unit[name] + else: + unit = None if not (unit is None or unit_old is None or unit_old == U.none): value = value * unit_old.express(unit) + except NoSuchNodeError: + self.logger.error(f"The value {node_name} is node available", stack_info=True) + raise finally: if not open_before: self.close() diff --git a/fragdisk/__init__.py b/fragdisk/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gui.py b/fragdisk/gui.py similarity index 99% rename from gui.py rename to fragdisk/gui.py index 6561773..52055ec 100644 --- a/gui.py +++ b/fragdisk/gui.py @@ -17,7 +17,7 @@ from scipy.stats import linregress from skimage.draw import line from snapshotprocessor import SnapshotProcessor -from params import default_params +from utils.params import default_params class DraggablePoint: diff --git a/galactica/__init__.py b/galactica/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/fragdisk_galactica.py b/galactica/fragdisk_galactica.py similarity index 100% rename from fragdisk_galactica.py rename to galactica/fragdisk_galactica.py diff --git a/ismfeed_galactica.py b/galactica/ismfeed_galactica.py similarity index 100% rename from ismfeed_galactica.py rename to galactica/ismfeed_galactica.py diff --git a/ramses_astrophysix.py b/galactica/ramses_astrophysix.py similarity index 100% rename from ramses_astrophysix.py rename to galactica/ramses_astrophysix.py diff --git a/galturb/__init__.py b/galturb/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/galaxy.py b/galturb/galaxy.py similarity index 100% rename from galaxy.py rename to galturb/galaxy.py diff --git a/sectors_extraction.py b/galturb/sectors_extraction.py similarity index 100% rename from sectors_extraction.py rename to galturb/sectors_extraction.py diff --git a/ism.py b/ism.py deleted file mode 100644 index 2fa0ddf..0000000 --- a/ism.py +++ /dev/null @@ -1,101 +0,0 @@ -# coding: utf-8 - -import numpy as np -import pandas as pd -from plotter import U -import snapshotprocessor - -mp = 1.4 * 1.66 * 10 ** (-24) * U.g -z0 = 150 * U.pc -sink_header = [ - "Id", - "M", - "dmf", - "x", - "y", - "z", - "vx", - "vy", - "vz", - "rot_period", - "lx", - "ly", - "lz", - "acc_rate", - "acc_lum", - "age", - "int_lum", - "Teff", -] - - -def convert_coldens_s(n0): - return (np.sqrt(2 * np.pi) * mp * z0 * (n0 * U.cm ** (-3))).express(U.coldens) - - -convert_coldens = np.vectorize(convert_coldens_s) - - -def get_dispersion(dset, name): - """ - Compute dispersion from dset[name] - """ - vel = dset[name] - mass = dset["mass"] - mass_tot = np.sum(mass) - if mass_tot > 0: - mean = np.sum(mass * vel) / mass_tot - return np.sqrt(np.sum(mass * (vel - mean) ** 2) / mass_tot) - else: - return 0 - -def get_sinks(pp): - csv_name = f"{pp.path}/output_{pp.num:05}/sink_{pp.num:05}.csv" - return pd.read_csv(csv_name, header=None, names=sink_header) - - -def analyze_box(pp): - pp.cells["mass"] = snapshotprocessor.mass_func(pp.cells) - coldens = pp.get_value("/maps/coldens_z", unit=U.coldens) - sinks = get_sinks(pp) - sinks["mass"] = sinks.M - res = {} - dirs = ["x", "y", "z"] - res["time"] = pp.info["time"] * pp.info["unit_time"].express(U.Myr) - for i, dir in enumerate(dirs): - pp.cells[f"v{dir}"] = pp.cells["vel"][:, i] - res[f"sigma_{dir}"] = get_dispersion(pp.cells, f"v{dir}") * pp.info[ - "unit_velocity" - ].express(U.km_s) - res[f"sigma_sinks_{dir}"] = get_dispersion(sinks, f"v{dir}") * pp.info[ - "unit_velocity" - ].express(U.km_s) - res["coldens_mean"] = np.mean(coldens) - res["n0"] = pp.get_nml("cloud_params/dens0") - res["mass"] = np.sum( - pp.cells["mass"] - * (pp.info["unit_density"] * pp.info["unit_length"] ** 3).express(U.Msun) - ) - res["coldens_initial"] = convert_coldens_s(res["n0"]) - res["mass_initial"] = res["coldens_initial"] * 1e6 - res["coldens_mean"] = np.mean(coldens) - res["coldens_beam"] = res["mass"] / (pp.info["unit_length"].express(U.pc)) ** 2 - res["nsink"] = sinks.M.count() - res["mass_sink"] = np.sum(sinks.M) - - return res - - -def load_wrapper(pp, fun): - """ - Wrapper to load_unload data and map function - """ - pp.load_cells(keys=["pos", "vel", "dx", "rho"]) - pp.coldens("z") - res = fun(pp) - pp.unload_cells() - return res - - -def allinone(pp): - return load_wrapper(pp, analyze_box) diff --git a/ismfeed/__init__.py b/ismfeed/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ismfeed.py b/ismfeed/ismfeed.py similarity index 100% rename from ismfeed.py rename to ismfeed/ismfeed.py diff --git a/pipeline.py b/pipeline.py index 54d608d..e6e8871 100644 --- a/pipeline.py +++ b/pipeline.py @@ -9,7 +9,7 @@ import numpy as np from plotter import Plotter, plt from snapshotprocessor import SnapshotProcessor, get_time from studyprocessor import StudyProcessor -from params import default_params, load_params +from utils.params import default_params, load_params fake_pp = SnapshotProcessor() diff --git a/plotter.py b/plotter.py index 41eda37..0b387fd 100644 --- a/plotter.py +++ b/plotter.py @@ -36,12 +36,11 @@ except ModuleNotFoundError: print("WARNING: no movie support (missing module moviepy)") from mpl_toolkits.axes_grid1.inset_locator import inset_axes -import pspec_read from baseprocessor import Rule, BaseProcessor from aggregator import Aggregator from studyprocessor import StudyProcessor -from run_selector import RunSelector -from units import U, unit_str, convert_exp +from runselector import RunSelector +from utils.units import U, unit_str, convert_exp try: import lic @@ -56,7 +55,7 @@ from astrophysix.simdm.experiment import ( ParameterVisibility, Simulation, ) -from ramses_astrophysix import ramses +from galactica.ramses_astrophysix import ramses filetype_from_ext = {ext: ft for ft in FileType for ext in ft.extension_list} @@ -270,6 +269,7 @@ class Plotter(Aggregator, BaseProcessor): self.nums, self.path_out, self.params, + tag=tag, unit_time=unit_time, selector=self.selector, ) @@ -1434,15 +1434,6 @@ class Plotter(Aggregator, BaseProcessor): if subname_y: hdf5_y.close() - def _pspec(self, name, **kwargs): - """ - Plot power spectrum (wrapper around pspec_read) - """ - del kwargs["run"] - file_pspec = self.current_processor.get_value("/hdf5/pspec") - num = self.current_processor.num - getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **kwargs) - def _overlay_fit(self, x, y, yerr=None, kind="linear", label=None, **kwargs): """ Add an overlay : fit a curve, linear or powerlaw @@ -1527,7 +1518,7 @@ class Plotter(Aggregator, BaseProcessor): This is where rules are defined """ self.rules = { - "plot_arrays": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp" + "plot_comp": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="comp" ), "plot_run": PlotRule(lambda arg, **kwargs: self._plot(*arg, **kwargs), kind="run" ), @@ -1674,7 +1665,6 @@ class Plotter(Aggregator, BaseProcessor): "Density profile", dependencies=["axis", "rho_prof"], ), - "pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}), "sbeta": PlotRule( partial( self._plot, @@ -1965,13 +1955,15 @@ class Plotter(Aggregator, BaseProcessor): self._gen_from_log("fine_step_from_log", name) for name in ["time", "dt", "a", "mem_cells", "mem_parts"]: self._gen_from_log("fine_step_from_log", name_y=name, name_x="fine_step") + + self._gen_from_log("SN_momentum_from_log", name_x="time", name_y="SN_momentum") # Dict of overlays self.overlays = { "g": partial(self._overlay_vector, "g"), "B": self._overlay_B, "vel": self._overlay_speed, - "speed": self._overlay_speed, + "speed": self._overlay_speed, "levels": self._overlay_levels, "contour": self._overlay_contour, "particles": self._overlay_particles, diff --git a/pspec_read.py b/pspec_read.py deleted file mode 100644 index 7f17211..0000000 --- a/pspec_read.py +++ /dev/null @@ -1,374 +0,0 @@ -from builtins import str - -import matplotlib.pyplot as P -import numpy as np -import tables as T - -################################################################ -# 3D -################################################################ - - -def pspec_rho(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/rho") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 2) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\rho)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_rho_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_rho_" + format(num, "05") + ".jpeg") - - -def pspec_B(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/B") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(B)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_B_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_B_" + format(num, "05") + ".jpeg") - - -def pspec_v(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/v") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 4 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_v_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_v_" + format(num, "05") + ".jpeg") - - -def pspec_vz(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vz") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 4 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(vz)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vz_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vz_" + format(num, "05") + ".jpeg") - - -def pspec_cos(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/cos_vB") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\cos(\vec{v}.\vec{B}))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_cos_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_cos_" + format(num, "05") + ".jpeg") - - -def pspec_logrho(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/logrho") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 2 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^2 \mathcal{P}(\log(\rho))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_logrho_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_logrho_" + format(num, "05") + ".jpeg") - - -def pspec_kr(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/kr") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(\rho^{1/3}v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_kr_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_kr_" + format(num, "05") + ".jpeg") - - -def pspec_vc(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vc") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v_c)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vc_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vc_" + format(num, "05") + ".jpeg") - - -def pspec_vs(file, path_out, num, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d3/vs") - kbins = h5f.get_node(group, "kbins").read() - pspec = h5f.get_node(group, "pspec").read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 4) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^4 \mathcal{P}(v_s)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vs_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vs_" + format(num, "05") + ".jpeg") - - -################################################################### -# 2D -################################################################### - - -def pspec_rho_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/rho") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 1) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\rho)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_rho_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_rho_2D_" + format(num, "05") + ".jpeg") - - -def pspec_B_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/B") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 1 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(B)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_B_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_B_2D_" + format(num, "05") + ".jpeg") - - -def pspec_v_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/v") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 3 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_v_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_v_2D_" + format(num, "05") + ".jpeg") - - -def pspec_vz_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vz") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k ** 3 * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_z)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vz_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vz_2D_" + format(num, "05") + ".jpeg") - - -def pspec_cos_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/cos_vB") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\cos(\vec{v}.\vec{B}))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_cos_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_cos_2D_" + format(num, "05") + ".jpeg") - - -def pspec_logrho_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/logrho") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = k * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k \mathcal{P}(\log(\rho))$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_logrho_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_logrho_2D_" + format(num, "05") + ".jpeg") - - -def pspec_kr_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/kr") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(\rho^{1/3}v)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_kr_2D_" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_kr_2D_" + format(num, "05") + ".jpeg") - - -def pspec_vc_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vc") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_c)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vc_2D" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vc_2D" + format(num, "05") + ".jpeg") - - -def pspec_vs_2D(file, path_out, num, plan, color="r", save=False, **kwargs): - - h5f = T.open_file(file, mode="r") - group = h5f.get_node("/out_" + format(num, "05") + "/d2/vs") - kbins = h5f.get_node(group, "kbins_" + str(plan)).read() - pspec = h5f.get_node(group, "pspec_" + str(plan)).read() - h5f.close() - - k = np.sqrt(kbins[1:] * kbins[:-1]) - pspec = (k ** 3) * pspec - if save: - P.figure(figsize=(8, 8)) - P.xlabel(r"$k$", fontsize=12) - P.ylabel(r"$k^3 \mathcal{P}(v_s)$", fontsize=12) - P.loglog(k, pspec, "-", color=color, **kwargs) - if save: - P.savefig(path_out + "pspec_vs_2D" + format(num, "05") + ".ps") - P.savefig(path_out + "pspec_vs_2D" + format(num, "05") + ".jpeg") diff --git a/run_selector.py b/runselector.py similarity index 100% rename from run_selector.py rename to runselector.py diff --git a/scripts/__init__.py b/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pspec_new.py b/scripts/pspec.py similarity index 100% rename from pspec_new.py rename to scripts/pspec.py diff --git a/select_snapshot.py b/select_snapshot.py index 032a7f8..d954af3 100644 --- a/select_snapshot.py +++ b/select_snapshot.py @@ -5,7 +5,7 @@ Select snaphots with a criterion """ -from run_selector import RunSelector +from runselector import RunSelector from plotter import Plotter, U import os diff --git a/snapshotprocessor.py b/snapshotprocessor.py index fc245d9..e8e3474 100644 --- a/snapshotprocessor.py +++ b/snapshotprocessor.py @@ -34,10 +34,10 @@ from pymses.sources.hop.hop import HOP from fil_finder import FilFinder2D from scipy import fft -import pspec_new +from scripts import pspec as pspec -from units import U +from utils.units import U from baseprocessor import ( HDF5Container, Rule, @@ -48,7 +48,7 @@ from baseprocessor import ( oct_vect_getter, ) -from run_selector import RunSelector +from runselector import RunSelector # Getters @@ -212,10 +212,10 @@ def pspec(map2D): # Compute fft fmap = fft.fft2(map2D) # Compute the power map from the fft - pmap = pspec_new.pcube(fmap) + pmap = pspec.pcube(fmap) # Use the power map and the fft to compute the powerspectrum # This is typically an histogram of k weighted by the fourier transform value - pspec, kbins, pspec2, fbins = pspec_new.pspectrum(pmap, kmap, kbins, 1, 0) + pspec, kbins, pspec2, fbins = pspec.pspectrum(pmap, kmap, kbins, 1, 0) # Return bin center and power spectrum return 0.5 * (kbins[1:] + kbins[:-1]), pspec @@ -1434,10 +1434,10 @@ class SnapshotProcessor(HDF5Container): return {key: df[key].values for key in df} - def _pspec(self, overwrite_file=False, **kwargs): + def pspec(self, overwrite_file=False, **kwargs): outfile = self.pspec_filename if overwrite_file or not os.path.exists(self.pspec_filename): - pspec_new.pspec(repo=self.path, iouts=[self.num], outfile=outfile, **kwargs) + pspec.pspec(repo=self.path, iouts=[self.num], outfile=outfile, **kwargs) return np.array([self.pspec_filename]) def _write_particles(self): @@ -1754,7 +1754,7 @@ class SnapshotProcessor(HDF5Container): "vz_gas": "unit_velocity", }, ), - "pspec": Rule(self._pspec, "Power spectrum", "/hdf5"), + "pspec": Rule(self.pspec, "Power spectrum", "/hdf5"), "write_particles": Rule(self._write_particles, "Particles file", "/hdf5" ), "write_cells": Rule(self._write_cells, "Cells file", "/hdf5"), diff --git a/studyprocessor.py b/studyprocessor.py index d1f0c7b..ffa48a5 100644 --- a/studyprocessor.py +++ b/studyprocessor.py @@ -12,9 +12,9 @@ from scipy.stats import linregress from baseprocessor import Rule, HDF5Container from aggregator import Aggregator from snapshotprocessor import SnapshotProcessor -from run_selector import RunSelector -from params import default_params -from units import U +from runselector import RunSelector +from utils.params import default_params +from utils.units import U class StudyProcessor(Aggregator, HDF5Container): @@ -452,6 +452,31 @@ class StudyProcessor(Aggregator, HDF5Container): series["turb_energy"][run].append(np.nan) return series + def _extract_SN_Mom_from_log(self, series, log_filename, run): + cmd_grep = f"grep -e 'SN momentum' -e 'Fine step' {log_filename}" + content = os.popen(cmd_grep).readlines() + block_err = [] + time = 0 + for i in range(0, len(content)): + try: + if content[i][1:5] == "Fine": + data = content[i].replace("=", " ").split() + time = np.float(data[4]) + elif content[i][1:3] == "SN" : + series["time"][run].append(time) + series["SN_momentum"][run].append(np.float(content[i].split()[-1])) + else: + raise ValueError("Wrong start of line") + except (AssertionError, ValueError, IndexError): + block_err.append(i) + + if len(block_err) > 0: + self.logger.warning( + f"Error encountered in parsing {log_filename} (grepped blocks {block_err})" + ) + return series + + def get_logs(self, run): glob_str = f"{self.path}/{run}/{self.params.input.log_prefix}*" logs = glob.glob(glob_str) @@ -815,6 +840,15 @@ class StudyProcessor(Aggregator, HDF5Container): "id": U.none, }, ), + "SN_momentum_from_log": Rule( + partial(self._from_log, ["time", "SN_momentum"], self._extract_SN_Mom_from_log), + group="/series", + unit={"time": "unit_time", "SN_momentum" : {"unit_mass" : 1, "unit_velocity" : 1}}, + description={ + "time": "Time", + "SN_momentum": "Injected momentum", + }, + ), "turb_power": Rule( self._turb_power, group="/series/rms_from_log", diff --git a/turbox/__init__.py b/turbox/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/turbox/turbox.py b/turbox/turbox.py new file mode 100644 index 0000000..fdb8038 --- /dev/null +++ b/turbox/turbox.py @@ -0,0 +1,160 @@ +""" + Toolbox for the TURBOX project (simulation part). + + By Noé Brucy and Corentin Le Yuehlic, 2022 + + Compute and plot eta_d, from the slope of the logdensity power spectrum +""" + +import numpy as np +from plotter import U +from baseprocessor import Rule +from matplotlib import pyplot as plt +import pandas as pd +import tables +from scipy.stats import linregress + + +def get_pspec(pp, field:str, dim:int=3): + """Read power spectruù + + Parameters + ---------- + pp : SnapshotProcessor + field : str + field to get + dim : int, optional + dimension (2 or 3), by default 3 + + Returns + ------- + tupple (np.array, np.array) + wave number and correponding powers + """ + h5file = tables.File(pp.pspec_filename) + path = f"/out_{pp.num:05}/d{dim}/{field}" + node = h5file.get_node(path) + kbins = node.kbins.read() + pspec = node.pspec.read() + h5file.close() + k = (kbins[:-1] + kbins[1:]) / 2 + return k, pspec + + +span_resolution = { + 256: (0.8, 1.1), + 512: (0.5, 1.7), + 1024: (0.4, 1.7) +} + + +def get_slope(pp, field:str, resol:int, plotdebug:bool=False): + """Get the slope of the Power specturm using linear regression in the selected range + + Parameters + ---------- + pp : Snapshotprocessor + field : str + field name + resol : int + resolution (number of cells on 1 side of the simulation cube) + + Returns + ------- + tuple (float, float) + Slope, square value of the correlation coefficient + """ + # Trustworthy span od the power spectrum in log10(k) as a function of the resolution + + logkmin, logkmax = span_resolution[resol] + k, power = get_pspec(pp, field) + logk, logpower = np.log10(k), np.log10(power) + mask = (logk >= logkmin) & (logk < logkmax) + results = linregress(logk[mask], logpower[mask]) + if plotdebug: + plt.figure() + plt.plot(logk, logpower) + plt.plot(logk[mask], results.slope*logk[mask]+ results.intercept, lw=3, ls=":", color="k") + pp.logger.info(f"Fit results in get_slope({field}, {resol}): slope:{results.slope:.2f}, b:{results.intercept:.2f}, R2:{results.rvalue**2:.2f}") + if results.rvalue**2 < 0.8: + pp.logger.warning(f"Bad fit in get_slope({field}, {resol}) with {logkmin} <= logk < {logkmax}") + pp.logger.warning(f"log(k) is \n {logk[mask]}") + pp.logger.warning(f"log(power) is \n {logpower[mask]}") + + return results.slope, results.intercept, results.rvalue**2 + +def build_suite(pl, redo=False, cs0=0.28834810480560674): + """Compute an array of parameter for each run in the Plotter pl + + Parameters + ---------- + pl : Plotter + redo : bool, optional + redo the comptutation of the velocity dispersion, by default False + cs0 : float, optional + Sound speed in the simulations, by default 0.28834810480560674 + + Returns + ------- + pd.Dataframe + dataframe with the properties of the simulation + """ + + df = dict() + df["snapshots"] = pl.nums.values() + df["n0"] = pl.study.get_nml("galbox_params/dens0").values() + df["turbinit"] = pl.study.get_nml("galbox_params/turb").values() + df["solver"] = pl.study.get_nml("hydro_params/riemann").values() + df["slope"] = pl.study.get_nml("hydro_params/slope_type").values() + df["res"] = pl.study.get_nml("amr_params/levelmin").values() + df["frms"] = pl.study.get_nml("turb_params/turb_rms").values() + df["seed"] = pl.study.get_nml("turb_params/turb_seed").values() + + df["comp"] = pl.study.get_nml("turb_params/comp_frac").values() + df["L"] = pl.study.get_nml("amr_params/boxlen").values() + + df["T_turb"] = (np.array(list(pl.study.get_nml("turb_params/turb_T").values())) + * pl.study.info["unit_time"].express(U.Myr)) + df = pd.DataFrame(df, index=pl.runs) + + if redo: + pl.study.avg_time_sigma("x", overwrite_dep=False) + pl.study.avg_time_sigma("y", overwrite_dep=False) + pl.study.avg_time_sigma("z", overwrite_dep=False) + pl.study.time(overwrite=True) + + for ax in ["x", "y", "z"]: + df[f"sigma_{ax}"] = np.array(list(map( + lambda x : x.T[0], + [pl.study.get_value(f"/series/time_sigma_{ax}", + unit=U.km_s)[run] for run in pl.runs]))) + + df["sigma_all"] = df[f"sigma_x"]**2 + df[f"sigma_y"]**2 + df[f"sigma_z"]**2 + df["sigma_all"] = list(map(np.sqrt, df["sigma_all"].values)) + df["Mach_all"] = list(map(lambda v: v/cs0, df["sigma_all"].values)) + df["time"] = list(map(lambda x : x.T[0], + [pl.study.get_value(f"/series/time", unit=U.Myr)[run] + for run in pl.runs])) + + df["sigma"] = list(map(lambda l: np.mean(l), df["sigma_all"].values)) + df["Mach"] = df["sigma"] / cs0 + df["turnover"] = (df["L"] * U.pc.express(U.km) / (2 * df["sigma"]))* U.s.express(U.Myr) + + return df + + +def rho_pdf(pp): + pp.load_cells() + rho = pp.cells["rho"] * pp.info["unit_density"].express(U.H_cc) + rho_0 = np.mean(rho) + print(rho_0) + s = np.log(rho/rho_0) + values, edges = np.histogram(s, bins=300, range=(-15, 11), + density=True) + pp.unload_cells() + centers = 0.5 * (edges[1:] + edges[:-1]) + return (np.stack([values, centers]), {"logbins": True}) + +rule_pdf=Rule(rho_pdf, "Density PDF", name="rho_pdf", group="/hist") +def apply_rule_pdf(pp): + return pp.process(rule_pdf, pp, overwrite=True) diff --git a/utils/__init__.py b/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/line_annotation.py b/utils/line_annotation.py similarity index 100% rename from line_annotation.py rename to utils/line_annotation.py diff --git a/mypool.py b/utils/mypool.py similarity index 100% rename from mypool.py rename to utils/mypool.py diff --git a/params.py b/utils/params.py similarity index 94% rename from params.py rename to utils/params.py index 99749db..3f7ca68 100644 --- a/params.py +++ b/utils/params.py @@ -33,7 +33,7 @@ def params_from_file(filename): def default_params(): - return params_from_file(_dir_path + "/params.yml") + return params_from_file(_dir_path + "/../params.yml") def load_params(filename): diff --git a/units.py b/utils/units.py similarity index 100% rename from units.py rename to utils/units.py