From 0e72363dab734a1dd5f3d758a3c60aedda7335ff Mon Sep 17 00:00:00 2001 From: unalmis Date: Thu, 27 Jun 2024 20:10:21 -0500 Subject: [PATCH 1/9] Use max rho of surface to compute A(z) --- desc/compute/_geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index d9d736aa95..4ae72a7713 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -176,7 +176,7 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): label="A(\\zeta)", units="m^{2}", units_long="square meters", - description="Cross-sectional area as function of zeta", + description="Enclosed cross-sectional area as function of zeta", dim=1, params=[], transforms={"grid": []}, @@ -196,7 +196,7 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg transforms["grid"], data["Z"] * n[:, 2] * jnp.sqrt(data["g_tt"]), line_label="theta", - fix_surface=("rho", 1.0), + fix_surface=("rho", jnp.max(transforms["grid"].nodes[:, 0])), expand_out=True, ) ) From 6c8f81cfc1f0546f8dc4e3ab9036556a440bd933 Mon Sep 17 00:00:00 2001 From: unalmis Date: Fri, 28 Jun 2024 13:28:28 -0500 Subject: [PATCH 2/9] Genearlize A(z) for nonzero omega --- desc/compute/_basis_vectors.py | 22 ++++++++++++++++++++++ desc/compute/_geometry.py | 23 +++++++++++++++++------ tests/inputs/master_compute_data.pkl | Bin 7756094 -> 7788370 bytes tests/test_data_index.py | 15 +++++++++++++-- 4 files changed, 52 insertions(+), 8 deletions(-) diff --git a/desc/compute/_basis_vectors.py b/desc/compute/_basis_vectors.py index d260ddee31..01bd1c24f0 100644 --- a/desc/compute/_basis_vectors.py +++ b/desc/compute/_basis_vectors.py @@ -3523,3 +3523,25 @@ def _n_zeta(params, transforms, profiles, data, **kwargs): ).T, ) return data + + +@register_compute_fun( + name="e_theta|r,p", + label="\\mathbf{e}_{\\theta} |_{\\rho, \\phi}", + units="m", + units_long="meters", + description="Tangent vector along boundary of constant phi toroidal surface", + dim=3, + params=[], + transforms={}, + profiles=[], + coordinates="rtz", + data=["e_theta", "e_phi", "phi_t"], + parameterization=[ + "desc.equilibrium.equilibrium.Equilibrium", + "desc.geometry.surface.FourierRZToroidalSurface", + ], +) +def _e_sub_theta_rp(params, transforms, profiles, data, **kwargs): + data["e_theta|r,p"] = data["e_theta"] - (data["e_phi"].T * data["phi_t"]).T + return data diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 4ae72a7713..48763706a9 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -12,7 +12,7 @@ from desc.backend import jnp from .data_index import register_compute_fun -from .utils import cross, dot, line_integrals, surface_integrals +from .utils import cross, dot, line_integrals, safenorm, surface_integrals @register_compute_fun( @@ -176,25 +176,36 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): label="A(\\zeta)", units="m^{2}", units_long="square meters", - description="Enclosed cross-sectional area as function of zeta", + description="Enclosed cross-sectional (constant phi surface) area, " + "as function of zeta", dim=1, params=[], transforms={"grid": []}, profiles=[], coordinates="z", - data=["Z", "n_rho", "g_tt"], + data=["Z", "n_rho", "e_theta|r,p"], parameterization=["desc.geometry.surface.FourierRZToroidalSurface"], + # FIXME: Add source grid requirement once omega is nonzero. ) def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwargs): - # divergence theorem: integral(dA div [0, 0, Z]) = integral(ds n dot [0, 0, Z]) - # but we need only the part of n in the R,Z plane + # Denote any vector v = [vᴿ, v^ϕ, vᶻ] with a tuple of its contravariant components. + # We use a 2D divergence theorem over constant ϕ toroidal surface (i.e. R, Z plane). + # Denote the divergence operator div by [∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]. + # ∫ dS div (v ⊗ [1, 0, 1]) = ∫ dℓ n dot (v ⊗ [1, 0, 1]) + # where n is the unit normal such that n dot e_θ|ρ,ϕ = 0 and n dot e_ϕ|R,Z = 0, + # where the labels following | denote those coordinates are fixed. + # Now choose v = [0, 0, Z] and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. n = data["n_rho"] n = n.at[:, 1].set(0) n = n / jnp.linalg.norm(n, axis=-1)[:, None] data["A(z)"] = jnp.abs( line_integrals( transforms["grid"], - data["Z"] * n[:, 2] * jnp.sqrt(data["g_tt"]), + data["Z"] * n[:, 2] * safenorm(data["e_theta|r,p"], axis=-1), + # FIXME: Works currently for omega = zero, but for nonzero omega + # we need to integrate over theta at constant zeta. + # Should be simple once we have coordinate mapping and source grid + # logic from GitHub pull request #1024. line_label="theta", fix_surface=("rho", jnp.max(transforms["grid"].nodes[:, 0])), expand_out=True, diff --git a/tests/inputs/master_compute_data.pkl b/tests/inputs/master_compute_data.pkl index 836a76b534bbe7f51e7455853e46cadf08947e86..11ae1aa4cfc00e7da8856afa7914cb5b099a0c43 100644 GIT binary patch delta 18557 zcmeHvd0bRivaq+ip&Oc=W|6fOQBV*N5D;Y%H&hI$(P*UEp#_Boq~j!NNW?LTCPA^S z$)@;pGD)6^MkDo-(TFqI#w=ot+1#^Dl245qqse3vzdGkGeL;QW_~s?^{eIschrYMY zsXA44>YP)j>Revl@sal3=Z)(+sgm%h^=ufmuca;hXccZx9u|Oy&1e zT!|~DWjwKkH8#CUh+RMH2s`O67h-C5wadq*u7R$^@eof}0$-TONN;Xre~*oHrOuZ8 zmOamPu%$4A`KEUAm7ueqo7hF@V zZjK!TMLpx8x%$HaW3cvE5Re5v?R@uy;<5GhjiKjB0 zN&=Nc+}Dt_W3r~q3tL`Fiha^^wXv#Xk+o`3*^;I>8=Sv0AO;V7p#Nfw=AtK@5^7Sw zvvY(2%!xIT(DS6=iC0QTycnx>`1wKMgi!;byCG!=gw2^vP0mz@DowsH!H}_;i8#y> zxyuA6H;zbxxbiI3dpvkQnnFk|**cuIf5eJ+wGMYqatzr1Ug`m1E!hS*_jTqVOuR37 zuW9|R5yK$oznUVk^=R^2>@`Xu5Qp`zCvP@E;`!t-Ua?G@>H(#%H-)ms)kOJcU&id1 z^^ocVi67VK;GIdygYeN&-o#s>Yf3y=&!%|5xjTo4VEZWVy(Z|Ko05bDAEg{JLAEwE z54QW&=+JgHW$C*b$7#oak^PxS?r`#A=~!sB_>X~S^L-51zAyFRENFet@(@74(8K#o>50M$0O_&~; z;THh(H*Ecz=|gaRO2!HjYXjUeVbo+e=T{cVMZ$-U_^;qY=MIH2O`@FmM%btou`ur| zlOME~rN_hG59^}D<{u}GYBsSNf4uAVQL|?Vd+^Y6qrS0f9cO*l@5qS)^X$wtIGZ(^ z27S(Kol59PZRRc$UWv>6F%}NwX5I$dQ)%oIc8!R^xcJP=CMf7&NX$|mT97$xkJeGh z$1@l#TgJHId$FTU+3;xAXg~4YWBzDY(#v{h8+RO@ImSO5LaN62!GoU?PSu@broIb< zS7${*p*edbNl>PP_L6}*7&bOD4!b?Gs9Pr+^!}}BeO1$tzKeA@G&U0F#_YaE}#mFIJ2i+!H_(pM+ib)Riu26g8FiDe5y)g;&kF#`&JB5u3S=a4xu zBkzM54n8jigu~so++gV3lx0}|=h;KC+njyc1g#&Gc|hl!>Y%G}PR=gRfO7}RBKua& zU|3fPJ_C9lApRbH=vw?yl>vGRs&(+qu$rs*a~hs($_dY2 z|5=F}>U+3@xD@?JV;3Jbv9 z+j7kToNl%`+YK+jlKW|o#xWIoKFJHgGmfzW)*lN91=*6)xZl%kqz#=(nRIUzp)>LBZ0Gg@c=J9Sq&0T@qJb`XV1MeKEQ7 z0+EN0rX1--er(BHk$mUbDO*j#n3Vkz99Ix_2R zs$BWGZZ>|dyY1NnLJW4iciXbhy<-N*61&z|nig@t1Z-(@6S3)~IU}>E-P7OCc`=rL zD+bOzYNFqqw7KgU{DYRc9KL3w2(N!(?nVap`c{O$5q3L?e&5RZ-t%z2cd_krVISNR zZHa&$AB!vee6e_!@VN@`;djN)33DN7h$TjBvFlRtbDFtPusy*Z3Fj^r8({V2;uk#T zMoA3Uy7@9cEfBUFES}KtJwwNcABrn%mnwy(pz!#+^u>Dx2fZV88}qby#q$<4A+$YX&0!erRzZ;6GPp!t3{(P%j> zz!@)#yDZW?Ps56Krlpdan23;Iw7BBe{+6FPyf96IV_={~r!gwJseOc#!AZdc z7CPbXFiSkVbD`K3TSExf*s17btJrXB7{lMGfV-fmkcRW42$Fq|0__%|U6BNRU4iar zBm&rZ0@}?MlShwS)3SOnG#V_S@TS@F={3jU*<6NH55c=C^+?)ZGT8Eh zsQ9C5bsJ!0jO99JqzC`lUIQs=8o!Gr$$Kdohr~%Br+()QXZpEHzpD}nEa?~bj9_i9 zQtdY>41ASp!_uy$T3d$22i(#vTdu9&L+KWi&;wtlTSf??+95H+a`A?AOT!=BHM+NI zr3o9FWqIyuCG)xbJyppj&H;zBEq$i$)|~6<+!G$o?M0dEd=M@eYtbq?$M$g+U;gG) z_kju5Qa7BtQfxrCe9PYjS?9@9uBr2?f@|u0?JZa9{9VoKQSf+Mr5VN-S;7Tv#~MvA zY%8*a3Q7F2zNEn)9xAdJgqi&PXi0;~HHU8`P*!9K0iV(alOdV|tn+x2Dbp<`gV@B{ z#UXbSM0fxnENv(ZfWpp7gSh4dLc3RmC;C_|M}`Xf@qA}PmY!GZ6LkIh)gFt++@~5c z^zin@Qg>|I)-W*umb6zI@XapEiZBs6_mE{0LGPOHjxSwm7^!DuexhIKf7#;ayFXf9 zQ2k|sF19IB`}HBVk$TPkXgFc3j)c;gB}RBSvBc0d(B`36VIFi5SqYDoF(Ghdbr45D zLNIX}MXe$YKRhztrqjciv-92jo!tlNIT1R37q(RgL+RVK{+Li;i;!5!pIO3LZ@s-F zn~Bg}PQRf}^kf;mHrQ3mmeZgN{9K5lS_+VryMu7ujg0P*#iiOqwc`%B+2 zD6OhCK>Ki05ZDe`!m+TnBtUF~&pab`dg}R-A_+A4utd`S_gh12Ps65(&G8~P=c0!h@Rq>F{Jwl*&fB#z7wO5R zTpJCeN(scrcKQePATw;Lk4stq6ht|aeK=WPm z{c*yEMyh&tmBoO+{e7ckOz`NdtkuktNWIV>z5moGcUk*EBXQ9mHyvjnrerh>sU*uf z`Fx4HFbI5~um-Y$G+~2_Ej&PzJOs`^mogd?FE)?scL~V&@sBB=`JW$8q=VUoedysDO^%p`xIwksQT&!5mB3a@G@> z%y={}YRXfOw%(m_2|jIQ`MCQol9a?JmSooGAbFt701p?`Nc^eD*K<@P)cY9u^IJ%e zlf}sr9|Td@!}d*0nd)Be;hX_+HjW+76v-Hm`9Q}9I#;rD;`Z%LX=*%QBjXlNJ}q`B zci_Z7pV&Md=`_`=w?;lX&^oj{ho7V%XL?J2>`5<+P}A;Fc*;=@Si$meapj3;rJh%6 zB(O>shdEr6lcFiD#2!mkf7iFXL10Bn?eT_LlN$*^zkG6 zFe)?JdIUB-)f|Q0tD5zCAqSstD^Jx&9#Mv)QkmxY%F~nBn7J@c?^(Mln2AAcYgY&3 zscp?E`amb9L#g~+3+Lh%Mnh+tr3oV+D|fgC+@}!Rr95pFJ2yH>xS<)FX^Hj9obL`7 zMp>mzGO}i&zqlj|nIQb0Et@k-e9^YOJWY=!^-AStAZZ5R#uf!)RyALFq}fYvJUCP~ zyjZ2qh>>jzX&XK`6mqtg55jK6$@Y+NU6HUY#30R98_%l<;J0#cctaZPDvD zStCUQX7%oP?-RTyd~-)dwq698ULufjn_?_6xNm$tjTdj<|M7Zqhe(HO%4I^rK`7i@ zpNc)}>a$fgSi4CSoJKK_cDlyH(gdJn~fYX=O#qRNV7(bdxo zQ@@Fa&Ok6n)scbtR(grfpLm+)cC#F_hqOHI%Dq3xf17mT<6kyH59gLzqwt$!^<|PS z=YCr+xiUbHqwzpZ%Oo)z!~G==YgpPMZtZTcxvg@f)1G%K{V{Dt<={}w z@;D%GGz1GjuK7SKRAbSW%6Wr@27GRR)ASQ6+FEelzB5HWX`s$N`{fKX5$D9Ws9 zhy$PB&5sn{Je;B=y1~BJDkn1T(;@NAO7gte(=YVu%XuH}Q~>uRn7(ScFTAjjU62qZ&nnyNv$5qW`4NuLXhTHd= zhxLMBx4mkGx|*V`qX*DxFJj>j53=qW@cbJ3W#jQsh)o&~f}VPWpo5%>;7$S;nmbPX ziJxYc`k`Ri(-FwZUT_1@Y%WBAmbDj&#X!EWvG_vrfN<*38 zQ&d116jC0P{;@I>OaEG#>8cyHgQ|)tvKQ< z6=;P_=>^92D^+K;{*Jc3#=RKt(^iWEK^n(I@9K6Hk9-KsjK@eWlxA%l6mBGg&Pk*x z?WS-7o1UB>+`B{1k=2*ADwF@S5y|9Fv$6PT3{{0F$hG{a$92v9tmWu6wpuo|e7JjJ zEQPYRR7RuEwCYrDGLp6><(<4^0k^S1S2sX9%Pl1oJZ_^dRhs*j&1Bka|DnnYTM3cO zby6E?iHR>n)lu+$UFC5t(rL-Yh$n4ib&th?ZB!|x-`=VGj4ctgp|fFQ@T?7uk+{9w zdRz;1j0$2+C=UHnb*8uG>*3(@r$&y&LF0|uTI!{HG4;aNO=_{eweBjiWJ4!~np7q2 zH@(-Y#Ov_Hm(^ED-|QuWC;yNRBHuCypG$=&MhrurwN)R=5K6px;Ii+Sz#J(Rh=Kv= z++dd0a1||>w!A4!1|xMbOc(Ar($~D4OWr&XHdF;+r&rDUy~*3w_9g#Q8>JliQA6gA zi>hQIfzBIJ*K9^7g_M!%x6~P+;e2xBRT9C{s2T%Hr}S4kXX0Z{WQ|rQX$Pszj}lwC zEM1~V#*{aPXGxHUvb1 zA{^=*NpAn!uT*^?*$p4%1m0M($`sf(HT~a1BC_}lIC=-|LK4SFkVM6Rxm2fD#2Vn) z#FS7r1{7?|s(w{1nn|1&w$|!i68k7A#Ni@}k8{rRSjuwPbhsv0%H4WqMadHsX?q;c zI7nJkGhJMnRH!^U(hF5e`#L8$;ZNq-Urbmc&P;+#poDCUyezRzI-V{&)P#vUu;_?3 zh~RG02HSY2^&_?-bB0rZ!YRR`=qg3>n2M7t9mu(C=bEnY zKzVHE!pES!bJ$893cG7Yv!N5YxMRWS-Wn&>lEx2|IlCa~r@^@uu{v&GiWE#KNHG+? z_HqqT_+-#kGVqAWawI#-kqkys62U(=CWkD98APXIiLaYb0ehJbYGbO0ppeRVI>btgR}QSXR7@$8b(Ovtk!*cU3s0P> zbPwPWACAW(KK$uPpomWZc%@lsHZ@!P$%g38z&|Yp$VqLtZH6BZrPumw*C}Ie9)?44! z;__PS*LuMdD5?=CdZ@Ut$x2@BV?c8}9!h^+$PZ>}kCr^%C+?@tjm89-W+V7MyeRM* zu5Z^S>m{wmm;i6>khxT$#fbg#WmyF-V466?g)7=zPY0QG^^%pTq5+$FHJT~TB-(h9 zLRi!#se*9Rfg&lAEvn8U6izuc*7~iMMN;^K6huMWtE6OOwQ}I+$z$m-PXUaax0%kX z8_|c zvbX-2GqjG3Nls*vRr07Uoxmw%oQwxvYUHi{H3>|XVEon6m3#WRt2jcj6C%I=j?I8x?PHOf2yNB#M> zl});#)TIiSd*ajq;jUyI{KMR(3YWXw)UVU!mJ-Y?0dZCGQiaGJZR+&Y(_*LpUG`Fi z#^bEiAGw3vH#Mlzm#StwGAXi^eS`T+6*iBniOS5mfdpordbGInE#dy7e5I;O?sLlw z^Mu}Y5|%0)?rFQ4`JUq%2hwxHm2LeuFv7k(AB3*Q4)} z%}nLj$YS#O&3#(_)H{pGLD$P-^2+OFF?l8TXZh1Vlf|SI;u}i-u}><~Q}0sCGKO$B zmdCKHCN=(8Y3xjG68KXE{7dCGd*wd+WHYJGG)#7Vohi)-k|&(t5U;O@VTny@H&|)c zC%Z{?2})fPQqE}NKbPdp2{1Vp^u1(oIvO_IFXtlpb|+YoT=R+Q9A={k`Lz4WT`f$3=dpGG z9?fqUM4kl*B@HuwanHO+J;amQ1zYD`!5j6>efs$36`1o>!zkx%!vC3TuDY4RU2~%D zVv{eS(bd9uy~*+CtKHvgT{cva;_0~x42zXX-RNCI6%C$%s-Q_I)4&_LURW_&_Qfkk zn!k#jfWYukVZJkmLO4zR4QZwZnhra-nJq<*qp?gIsg@&M4h3^)%dviWHb$T;CIrVoVKDo%&&HCh`GH=c&mvH}kpS(>U z<#?>`M0cLxxtA`B@XvVZa_2YWrPF|1(@Q6sH|V9yz2Bgh&Z1AhmQT)SHTTK+Q*WP~ zgRbY3^UCY_t9fG*Lo8+Fx(+gt(Vm=yC^UJ%HaHY*}KUPeN@UxBS)( z7M*GN5zl|tvfQ9KHn4Yl18zLiGSyAvnFk3$Ml*OHZ1KSUms`fU3;g550Qe@>9tigf z_J_qU-#CCiThrK|pb-2)gFOl_YwSk^`5tJitKA2F(AhW3H$VwsfIqs}w+NG<)YbkL z{cKCH$G}~J-3%>ydz;6kJXu!yi1MI2b=Q+%*TX=AeT&gqw-l+fp9PqeY99oj`CEd~ zXteJZghsd}%^r{6dD!0*jExGZM%Exz|-+;_1Tx*5ef zzF@Ndr_ike><%D+`1Kn9XmEL;{c~A8%Yy7HWX;gC0R{xymrBVA+!buUPti$ zL#Vl!)&UOi2(>R0jx)`eMZk+9fONzGbe}McbaY%{n{>1?7)Q*#MI88ii0qCfPJ1|)Q^+Lf`5=1 zh>n4b7AY_YQwG`Z5sYS~1^*H=5Qzi{*(DzNNG1a;9b#Wa3gl-%0(7W!hwsGL1;IE} z5yxp|zn9fW7T4Rr=p#!SxoibTI(c);vc>nZ5=phZYRqui=wiSyzD48gUy*6j&mbrnZ+}@MjE6E;`d&9{g8eP_!oPi{7Xu_u zw7=pZjF(M?O!D~xGRf=-G<(dG?5`M|bxXxg{PL+2Zj$iTOj1eX6nlzBGjzOYg&cIK zfI!Bf3e`tM%+x}L$qZ1$TsNIH%Ttc;o H*u(!Xe;9Sf delta 17907 zcmeG@3s_Xu)}EPx8JIx^83cUr5D`&SR8WyOiiwZJw~`GguMt58gi~1|zSHXiMO>BM zw3pe##I%-HD&DJ=JrvDSQZtv`0d!NS~NW0d(|JVJ$@A7?{*=Oyw*IIk+ zwa;E_uYHa`e@4IWh0}VsYD1l9TjMs!=sr0+t)jGiP4ZUL{39cuFWR>@b87vy>(^;R|{(`U0F{id(H%RDR>#iG2 z95xtIb&ky12>lQOb#_bh$dtGkCP2;!LnK<`bWhr9&l*hl_<7wK8^(F-`a5vw3f;1# zodbrxZ^O6W)Sa~Lv`js0^Qar?Ewo{OJ9jVm%!XB8={~SQ^A^E<=VNY1Y+yO0H|_km znouhGz9LFa}8qgSa%qS2j$&Um%QLlMVC62CA>kalc?QxRpCfVfqO}cih`6?MWMyTyQht{ylEQEW(Qk`*Z@$ zhDTi|h0Q$yM3gHYBR07cwa?ygCyry5GvtHkSiQ26{WtC+`uNCI_gDv>%`(1c)7PEz zg1%2pXutG8sga6_iU$=>DrPDcDqd9DQt_r@rQ$=S9Ti_HHY$EpSbu*i0aOC1bfD6a zN+&9vsRU68rV>IWlu8#WVN|+O38&JHN(7beR3fSLpwg2{6qR07dQ<5`B^nQu_I+W5 zZk{)s={r6USH0ZVx&=<{?AI4w-5+9t=AQB0F>Y7C$8GS#xqd_7(v@-(*b@pnLiAd} z%hi$$m&8I7tWGHm*0dn?iRz)RYugT*3KIe#vqMZ@sELX-fHg9~0C#>?F4LL)VZZ%0 zsF@lQhWB;ufA3zX_V8*4DWQYHp|LvF0+z^xC@ekG|7a3~GobZUOaz?G>Sls!vqPq_ z!7XOL4V`^r9^MO?{ypUG)o#UB$a!#%jP&qe?-y+v!mln#F`-izdoEH}_cfd^ObFyc zA_+PFIQ)Duwi6J=Fqri~IUiUeY|F5kaN8K~`?X|7r4|Qc<9^#;WMEEMaZmPKPX_U$ zHyeB)W+x4JsuAYp1+{~k%0a!LVQFE6q^q-If=H<7pt!#9r2boDei*}P~AUHR{CJ0qOdbLgQgyb zprE*9j)R2rO_Hwedr1u>PAu*#k0Ui?Wet`XK+S#AO}ILJ@aKg01f-@99s&Dbl{JGG zPuP~&aMMGBCyy6);Ww8CUv@xE&-e+@c6@?l8ral6o>bJ(KHh?ZX2oA7rGG4OtA2lS zH_Wb!uR0E|j7w+-=Y@Dr)HxD-li^%pLO1O9RD$6+_&t@-5wq?|v?c56&a#dj;LF4W z6E54GIO;e!8>V|eO-|BCsLAePfW$M3V2$UJouuE4jpTJu)s#U<2abiR?y=KMp~^6`Y`kJ@#F>XD8BsK1(QA_Zr5 zWgEa%8`86|Pvz$W4Y6|)pYD zA5>|iB=sgjW@<5Cpe=egT9T$oT7t|sll&k#Ak7^wf05e3PhYoMT2zn{kmiP=`m}S6 z;5Rkh7wYqe7%=EW+CnnWv5*=zt^@9BH)PXMsGJ+>i=}@XvVt12q=vSkLkqrx$)#g_ zv8F@%5gPlhFc40p1!7sB^z!juRe^2VNL#;ggdN5X8W#>BU9#H3rVY`7&^#~NV`=)r zbWC43?7p@Z7p^b@uFUo8f|u6}`|g6^jTb9#dm&PY!Ht`4J8h$1ao26{yrRK34ujM? zCv{s|=^T!g&WstOgo}9U)r=2q^oz$bgz>_~rF*)Lz&%|@3fHLrEB>R-Fgg%5iqTtI zmO6^cXezf;8ABx#%TmW4{#vM78b4?pmGKxqXu|pT1@EODy(eNv?@8Xb3lp&Cos&Mb z33jxVPdXZDDGOG{>NyG<>nFW!lcs3(2a`PUNd2Tm+S2Mr=2=>SOvH@a@#~N%Z^9fUG@{^t&`C>#U7_^D47A{D++q{Kazw_YMwK z*t!?m%PZi3pQ#Jf+@&+tW`7=uSFdIzC$+?>nA-|xZQnblGCr5G?$~@)n74H0=xMle z^qtO+b-S$~WZT1_OMm-|y0h}bkNxe>!^63DPx$UysXut%Zf^^xV-(o8qwT>^A8X&m zQFLv5Bnq_y2*G-v-NS8#{_gBKa4gy$iq3)dSc7mcZi%`yso=B4-X5yc>|QuN!*0-7hN~SV`8(5E zSLBff5f1sI?9a6v1nq+(iDueprGZBBu67XN@Z;@rn^2))`u-SVYQ9%xTBT_-jiWZ* zy{9tnR<8FWw+ae*YU93?jihA)ui5Cx}smz_4?@~(G7I{ z*W8x6?hKprWQb}oCb&OefpLS$Do{E}2Gbb^7Rrzox(;+OHQZR&Ba5g_(KR}Y?LMlS zx6fwXZm4Ehi(&wN=d@oEWL?9$Jp}1FhafjmcuM6>6mEs8-ZUY-MEesr*7`Rh(>hr0 zviB1-laY0o{mb7o2|zaD4&O%G+hNzaB=0Y#BMg%7lUJs?emmW-tX0?g{Z57P>1eEK zIUQlJbDp$*JzGr2>gt5xqS+9CoerXRlargMUmmoY@y^q2jKA>P>ZUs`#S znJ_`^Dl*}HVc>l#F9c?;ws(duHTFO#nlRIs2gN$a9kn|Q0kCPcy`!*{S2w&-YBNaz zv12ud9A%I>d4eo!(#m%`4k?29`IyD?~lbF80m0iD^- zWTQl4Pl2~u5KQ|&;Y0wk7I5OiSZo! zMnM4M>5DU`W(TnzeVDk67FroePS}hz3kwSC#?I;XMr)VksNVsLl>8lc{3&vf&Kp(T^)* z;rwuGAINzvw>^&RQMT0=8fKQ6*3T+4V#GY!UQWG=~2!nH6*1x z*+|bhx6!liLO3)Z>K_A9@5P2;)2Fj~7~#H`E6D+;AP$;VQ23#5<&DM<+~ zY$WW#M(%GgMQ9r8?b8PvJ?g$_BPjq%hn03>kxm3SqX!3TX|wQ$M<=bhVtrik46a~_ zE{~Haj2Y4j#yT62f*DD=L6`Kh@fH1Dtxe=sO9ok#ub5_0XAcE0j@RI$4Kq@qeq>UF zEM5n(I8hBB3@mo;gdtt0M;gI1cuuMmK2X>kA_3hOXsBeqRoGp|OOR^Je zA~ro<+=o9J2TL-+l~_kIXNVN9$f$f~qAwQ@RKpc{X55XQOnymf9^_EeqG;A7;J#`7 zTL44)QP|weI!7Y!KD#9mS!ksGi&RlFu=#~z3EwzhD=LjjYd~oIX}ZXJx%-%1vJsXE zHy+vS-;(w5BX z5lK-$fWlbvr>*i9=a6e3}STnkJWESYIa- z({&JF>^{*3)|X3su&qyp(J0&j(@)tW*n^u9YCf6a5pq0S*^Gwh*cj-ssJmMGv#(0R z7n?qr;bt^12vE=#;*J+eVvX%xBvz0psE;{5u?l@o6 z9U-B@VuZ`rNsQ1-2(JwzX~`(RGH+<`_sI;6ta2+{3wO~M0=Rw#}c#~^fCZW3OU zAgiaU5PuB;|L(bMjppq=l%a=6ZOw=6X7qMc*tBhtnzk{eE!*PNHXrP`y0Qh$p=z5M zkKRxHblZD?<6AoRkk~Bdl9AUSOx9qIS)qcbYrv~(E6heXKRwq9FJ$F*!4JF5HXHSI zuj$w!s~`59nM=Fh7Pd*p;xkpb`D~Z7uFIwUYKyun+hFtU76aB`MF!ut+MaU+v)5{{ zQ3%Dj*K()&d2LaQVTDJ2Haz-nhYmWBnkd&pB5=64>&9~=8 z;Fqf^?#_K*(N>Uv-DbjD%H2T%qqo6F? z8H}5k6`W+0X6NO9sUKE%(yf*I=c3KGe69`deDE6I0C|l&*}OK`1CPZQ%r*=2VTZ4S zj-Io+1;D;xrTj><(@QJzr?zN9n)Lh4y7^Jyb1-r=TRkSQyk6KDbKcCKLYSBZB(K>e zM?BN)!x*w4pH5(DQttw)USX$p7EWbaX&z_79;Tn{77Gn~vjZ62V02dGXK_~KOh#e# zzWl}N6V&eyv4Hh@NGmJ`H$hH9XqeQS1Savi88+BNI0 z8k_%*xXfK$BjX+e*CBV5c#X#; zO$GO~)H62N3MyTo%SF?W zZN_k@dA_<-URRe2E$kaF6-iYZBRJ0bjJu(rt~sXxnV{k2$^f)_6l9uhb?&lnuLd_Q zOe8a^KzjS!cyLS!I|`h^Kue|{=FBa4QtzS89r%1dGe90LRxcwRkU8AFQu=kSq2CTF zZV`_h`RsDPuN>#s75s8`>G~PHa3z+~~a|=T0 zJ=ZIry_?1!WO5W@#v8H)9Ing`#+o&3rVHEqdm%57hi1fDoi)) zi%659DdqCZhapVAjU#ATkpYL!D~yHAYyAv7JxJ@d@EWOA0&l*pKYyVZkb)cq90$!} zp#dqvk=}5$%m-n^P}frotm@Z445L%?PiZMdXdx02;8BPMF%fGf7JjMco+ZugnyEy! zg)ohe7l$aXQzl~^{J5bNul3JTHjR57K6x2_OHQ}8_MBY{>whqAQP0eFp zlAK;dPZthgir=A^=0k48SIP^2kj05ZQiAf32qd)u8slbkQe;DtEyh5iAJ|Ol<>|iF z3#ASWf|M%!nUSK0IF=%uB-uLKF`ReC#lWj~4-A8pp#1)=;73-rfGZ4|9vc{e&VvP9 z;hg5n(?~5c1PKyjV6?d+ns>yRu;-JqZa6z4YeE|KHNIyt5H8!MiW`UfcBl`j(l{l< z<^xD%K~g`ozg&1o@8f!7^LV>F3Q`Z2I~27*=F!3tJeS5&Z!WW{zCvJ@cazLYLxo4r z5|!aB0jSPC5~*}PsE-^xg45f9VslGsrHq|!QwWS%#@|p zFy+y-!j$DAk#dpJq9mP^NJ?jN?mFBlid8;UlRph1o9s zW&k@FfSlYp0qPvOw9{NPN1hK^I~ty(tI?oC2;GTBG-IKvA@;B1V)3 zG8!6c6h0cQ@yNOyNo(=4u<1O;)dbU=u0~-ztXA2paWzMCF{At{dJmEHq^gc3O@BEB zJmw~f6mw@w3erf(!*6+B`7BcO-NG8vF=T?{2d#p5^sPw9`^Nu{`=QSnfClbE6X@O&}UW)m7<3gaAUrVgLix`joIhGRr6NyH> z6#a8JkN$a$NB>$n21W-;19Lkwi;y>Zp&wGze@WW)r6|9>wwXuy?J@Ha$4R~Uul<$i z=n(5gZtvUJ{V`hF^c%8;e$8*I7^_u`fM?j^f4e`d8I&h`2f|ert45&pFOGNgczuuK zGb5fmnq8|HbD3^aWHgcL`Atv9A?W7r1hu$2&}j~JD7={mEd`s@uO@otKJH3>^jzLh;34(qk&ZwXMJ&|lk9oT^ z3wf;^K=R7o89#Hhnjx~8?@^+DSC+`vQZ!$M{(BNeW3{IvJZqHXr>h$b8;fRkK+4Y1 zamH?p{-IRK?A)VM@-(=|gQ`BRV=un6=f>VF%hZXT>a6HRsU>LmloA9_cc1Bx)~`zG z0Ow{YCndMRW8k^5m=l}La~r>sX_^~L>+m;anuzW#Wtx2W`<2m`GELMkX=`?xA8b4{ zsU01e{DJuLt#%&IwiQim5;d97?9_l z7#2H5!{H%%B z_<0a%o!sfA(!u<1%R?EtCp|Zo`QyWpJZbry{#Nr))blsxp^W@^MdHU@Vgvk+M3j*` z)*32uhK6OL?xjScS(8ZYz^=ZfwN6SYa^Ey5W#n$O!h*Ze|2!+@&JX40!j~B(DK}0^ zNqqlVQp(5=CnV;a&~MF38M#ZX5OBHFQvQ~-^ebs8NyJ<^A;~j}H_ApCx!bI=P~B#U z#m&-D?%Xw=i=?4>++K~`9hP|BYBtKq{b7YQ^M{fAVcBy3zb2;ma^_AduPII{<(GJ3 z%3$=XyDq;!i+GBPqunG$#oH7oRg6i(- z`lB;k_~FLVhx{yc5%Q_v3SlJl-If=M8(t{w>IZvmc^>e_>!qXq103>UesG@`xlOm= zjayrJ<0J5At8?&;?tJ~L-uQ@wTk*z6B;0~G-pkmf?u50KRl9DP@k9L)*zls0-!l8v z+hCH9&bUfGI!Q@BI(LWv5ck~54@|f`+&v$G&Y5Lw7{9|i??^uO-`rcN=_a?q-1(8p z2<-eG=Z16{OzW=S&O$M|81`T1yyZHo9B>XlnUAisCIj2kK^1s7L zmu`!>h-iHDU%4%&hqXG_gSKb89SG$o4rBp$AOW$&VO?mvka`|69Ngp)7!NDxoA z)1fcSpCv{@a)LMv=rX`^ygo*}&mTMYM%I{^8MJptyE8ivOQuLTz^1x zhp8F#5lx5r;$H>+)!R0)GyJ$fe3pGRllrO{4o4P>8+8#w;L5d9`o8UXVoc2@*sw?f zH}lWl7&vFK_=F{5h%$IeVDSF2OT>+47r276#-EpnX9S^2dLTJY2tp=g?+|;v^e3r4 zeW%z9bt@Ux!xEOTnpzAqi2+nKl9cmyh+e2)!!Win41Z{-6^-!ITJc!{Rlp0{&y7~C z6F*jzvaC*AtSF`Hda9?wl}E%dd}qCQFIUB*50Q#W9uyZdRqzM@hr}w~URf1Q>%_Kj zP=;_V5X1zt9u^l__G-1jy9R9#tIRHNg~aR44aAGIz!}f{vp7reRL$rUrj8J3SWlz0 zZ4xU4;S(UG$dH#dvz`x2Jtsa!&4!;<&G0_oJkBt-F%0^W<3ZBQnkU6Sk!BQt7o={X zVINRGjl;if5e32Qp$hLZpY{;ARa~v8OJhU0^l6$EOOV=!eR;@iXwkDRtbu)g$iQQt z6Q5{lO+nivh|7>ZCrC>};$3$?FFx^0Lt4I__?t6T{<2A}*ddt_GbvoPQ+nVwRr4b2 zz|2Ux+$viaeGT_Z>!J_h2zRcEtF^N=uaGzk)h3t4(f!)RqN(<@n-#=Y#TDiiDvnFl z(0#9H6V_qkUa_4ZS>y+jc00rH{o-Eb0P-$@{A0qq9Iw(Ox4vC*znuaQ)q}9kf#JnaRG`qm^=bW5Px7Udc86*FyltUID z7Grce4|3PpN1Pn=KX1vPvBmVcD+gs9Vg2L}I5}i%BOw{Oq)RKc2YSBCFyt>cxs-YD XiEoi9a#kF0`}+jY=qL81*!h0}&ax6# diff --git a/tests/test_data_index.py b/tests/test_data_index.py index e3fd65a495..33ecaaafb7 100644 --- a/tests/test_data_index.py +++ b/tests/test_data_index.py @@ -8,6 +8,7 @@ import desc.compute from desc.compute import data_index from desc.compute.data_index import _class_inheritance +from desc.utils import errorif class TestDataIndex: @@ -39,6 +40,11 @@ def get_parameterization(fun, default="desc.equilibrium.equilibrium.Equilibrium" matches.discard("") return matches if matches else {default} + @staticmethod + def _is_function(func): + # JITed functions are not functions according to inspect. + return inspect.isfunction(func) or callable(func) + @pytest.mark.unit def test_data_index_deps(self): """Ensure developers do not add extra (or forget needed) dependencies. @@ -74,7 +80,7 @@ def test_data_index_deps(self): pattern_params = re.compile(r"params\[(.*?)]") for module_name, module in inspect.getmembers(desc.compute, inspect.ismodule): if module_name[0] == "_": - for _, fun in inspect.getmembers(module, inspect.isfunction): + for _, fun in inspect.getmembers(module, self._is_function): # quantities that this function computes names = self.get_matches(fun, pattern_names) # dependencies queried in source code of this function @@ -97,7 +103,6 @@ def test_data_index_deps(self): for p in data_index: for name, val in data_index[p].items(): - print(name) err_msg = f"Parameterization: {p}. Name: {name}." deps = val["dependencies"] data = set(deps["data"]) @@ -111,6 +116,12 @@ def test_data_index_deps(self): assert len(profiles) == len(deps["profiles"]), err_msg assert len(params) == len(deps["params"]), err_msg # assert correct dependencies are queried + errorif( + name not in queried_deps[p], + AssertionError, + "Did you reuse the function name (i.e. def_...) for" + f" '{name}' for some other quantity?", + ) assert queried_deps[p][name]["data"] == data | axis_limit_data, err_msg assert queried_deps[p][name]["profiles"] == profiles, err_msg assert queried_deps[p][name]["params"] == params, err_msg From a16afb7e57a0068a2a5efcad5a463620519de592 Mon Sep 17 00:00:00 2001 From: unalmis Date: Fri, 28 Jun 2024 13:33:25 -0500 Subject: [PATCH 3/9] Update comment --- desc/compute/_geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 48763706a9..8c83a51fd7 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -203,7 +203,7 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg transforms["grid"], data["Z"] * n[:, 2] * safenorm(data["e_theta|r,p"], axis=-1), # FIXME: Works currently for omega = zero, but for nonzero omega - # we need to integrate over theta at constant zeta. + # we need to integrate over theta at constant phi. # Should be simple once we have coordinate mapping and source grid # logic from GitHub pull request #1024. line_label="theta", From 95d94b31c02576450fae975ebb8951c98d15925c Mon Sep 17 00:00:00 2001 From: Kaya Unalmis Date: Fri, 28 Jun 2024 15:12:39 -0500 Subject: [PATCH 4/9] Update math comment in A(z) --- desc/compute/_geometry.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 8c83a51fd7..4f1e05735d 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -190,11 +190,11 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwargs): # Denote any vector v = [vᴿ, v^ϕ, vᶻ] with a tuple of its contravariant components. # We use a 2D divergence theorem over constant ϕ toroidal surface (i.e. R, Z plane). - # Denote the divergence operator div by [∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]. - # ∫ dS div (v ⊗ [1, 0, 1]) = ∫ dℓ n dot (v ⊗ [1, 0, 1]) + # Denote the divergence operator div by ([∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]) dot . + # Then ∫ dS div v = ∫ dℓ n dot v # where n is the unit normal such that n dot e_θ|ρ,ϕ = 0 and n dot e_ϕ|R,Z = 0, - # where the labels following | denote those coordinates are fixed. - # Now choose v = [0, 0, Z] and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. + # and the labels following | denote those coordinates are fixed. + # Now choose v = [0, 0, Z], and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. n = data["n_rho"] n = n.at[:, 1].set(0) n = n / jnp.linalg.norm(n, axis=-1)[:, None] From 44e03d1d8d7864493675f132d304281141d11d87 Mon Sep 17 00:00:00 2001 From: Kaya Unalmis Date: Fri, 28 Jun 2024 15:14:31 -0500 Subject: [PATCH 5/9] Update math comment in A(z) --- desc/compute/_geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 4f1e05735d..5522ed0157 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -191,7 +191,7 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg # Denote any vector v = [vᴿ, v^ϕ, vᶻ] with a tuple of its contravariant components. # We use a 2D divergence theorem over constant ϕ toroidal surface (i.e. R, Z plane). # Denote the divergence operator div by ([∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]) dot . - # Then ∫ dS div v = ∫ dℓ n dot v + # Then ∫ dA div v = ∫ dℓ n dot v # where n is the unit normal such that n dot e_θ|ρ,ϕ = 0 and n dot e_ϕ|R,Z = 0, # and the labels following | denote those coordinates are fixed. # Now choose v = [0, 0, Z], and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. From cc01bcecd00801850c8c34f18e735300aecaadbf Mon Sep 17 00:00:00 2001 From: unalmis Date: Mon, 1 Jul 2024 12:09:59 -0500 Subject: [PATCH 6/9] Use data["rho"] instead of getting grid --- desc/compute/_geometry.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 5522ed0157..3ba18a0711 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -183,15 +183,16 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): transforms={"grid": []}, profiles=[], coordinates="z", - data=["Z", "n_rho", "e_theta|r,p"], + data=["Z", "n_rho", "e_theta|r,p", "rho"], parameterization=["desc.geometry.surface.FourierRZToroidalSurface"], # FIXME: Add source grid requirement once omega is nonzero. ) def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwargs): # Denote any vector v = [vᴿ, v^ϕ, vᶻ] with a tuple of its contravariant components. # We use a 2D divergence theorem over constant ϕ toroidal surface (i.e. R, Z plane). - # Denote the divergence operator div by ([∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]) dot . - # Then ∫ dA div v = ∫ dℓ n dot v + # In this geometry, the divergence operator on a polar basis vector is + # div = ([∂_R, ∂_ϕ, ∂_Z] ⊗ [1, 0, 1]) dot . + # ∫ dA div v = ∫ dℓ n dot v # where n is the unit normal such that n dot e_θ|ρ,ϕ = 0 and n dot e_ϕ|R,Z = 0, # and the labels following | denote those coordinates are fixed. # Now choose v = [0, 0, Z], and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. @@ -207,7 +208,7 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg # Should be simple once we have coordinate mapping and source grid # logic from GitHub pull request #1024. line_label="theta", - fix_surface=("rho", jnp.max(transforms["grid"].nodes[:, 0])), + fix_surface=("rho", jnp.max(data["rho"])), expand_out=True, ) ) From 55313ede3689ec007fe505795cb04eec36472c9f Mon Sep 17 00:00:00 2001 From: unalmis Date: Wed, 3 Jul 2024 12:49:28 -0500 Subject: [PATCH 7/9] Scale area by 1/rho^2 to approximate area on LCFS and update description of variables --- desc/compute/_basis_vectors.py | 8 +++++++- desc/compute/_geometry.py | 29 ++++++++++++++++++--------- tests/inputs/master_compute_data.pkl | Bin 7788370 -> 7794246 bytes tests/test_data_index.py | 9 +++------ 4 files changed, 30 insertions(+), 16 deletions(-) diff --git a/desc/compute/_basis_vectors.py b/desc/compute/_basis_vectors.py index 01bd1c24f0..31cc469ab9 100644 --- a/desc/compute/_basis_vectors.py +++ b/desc/compute/_basis_vectors.py @@ -1514,6 +1514,7 @@ def _e_sup_zeta_zz(params, transforms, profiles, data, **kwargs): parameterization=[ "desc.equilibrium.equilibrium.Equilibrium", "desc.geometry.surface.FourierRZToroidalSurface", + "desc.geometry.core.Surface", ], basis="basis", ) @@ -3530,7 +3531,11 @@ def _n_zeta(params, transforms, profiles, data, **kwargs): label="\\mathbf{e}_{\\theta} |_{\\rho, \\phi}", units="m", units_long="meters", - description="Tangent vector along boundary of constant phi toroidal surface", + description=( + "Covariant poloidal basis vector in (R,ϕ,Z) coordinates. " + "ϕ increases counterclockwise when viewed from above " + "(cylindrical R,ϕ plane with Z out of page)." + ), dim=3, params=[], transforms={}, @@ -3540,6 +3545,7 @@ def _n_zeta(params, transforms, profiles, data, **kwargs): parameterization=[ "desc.equilibrium.equilibrium.Equilibrium", "desc.geometry.surface.FourierRZToroidalSurface", + "desc.geometry.core.Surface", ], ) def _e_sub_theta_rp(params, transforms, profiles, data, **kwargs): diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 3ba18a0711..2b26576974 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -176,8 +176,8 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): label="A(\\zeta)", units="m^{2}", units_long="square meters", - description="Enclosed cross-sectional (constant phi surface) area, " - "as function of zeta", + description="Area of enclosed cross-section (enclosed constant phi surface), " + "scaled by ρ⁻², as function of zeta", dim=1, params=[], transforms={"grid": []}, @@ -198,7 +198,8 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg # Now choose v = [0, 0, Z], and n in the direction (e_θ|ρ,ζ × e_ζ|ρ,θ) ⊗ [1, 0, 1]. n = data["n_rho"] n = n.at[:, 1].set(0) - n = n / jnp.linalg.norm(n, axis=-1)[:, None] + n = n / jnp.linalg.norm(n, axis=-1)[:, jnp.newaxis] + max_rho = jnp.max(data["rho"]) data["A(z)"] = jnp.abs( line_integrals( transforms["grid"], @@ -208,9 +209,12 @@ def _A_of_z_FourierRZToroidalSurface(params, transforms, profiles, data, **kwarg # Should be simple once we have coordinate mapping and source grid # logic from GitHub pull request #1024. line_label="theta", - fix_surface=("rho", jnp.max(data["rho"])), + fix_surface=("rho", max_rho), expand_out=True, ) + # To approximate area at ρ ~ 1, we scale by ρ⁻², assuming the integrand + # varies little from ρ = max_rho to ρ = 1 and a roughly circular cross-section. + / max_rho**2 ) return data @@ -413,13 +417,14 @@ def _R0_over_a(params, transforms, profiles, data, **kwargs): label="P(\\zeta)", units="m", units_long="meters", - description="Perimeter of cross section as function of zeta", + description="Perimeter of enclosed cross-section (enclosed constant phi surface), " + "scaled by ρ⁻¹, as function of zeta", dim=1, params=[], transforms={"grid": []}, profiles=[], coordinates="z", - data=["rho", "g_tt"], + data=["rho", "e_theta|r,p"], parameterization=[ "desc.equilibrium.equilibrium.Equilibrium", "desc.geometry.core.Surface", @@ -427,15 +432,21 @@ def _R0_over_a(params, transforms, profiles, data, **kwargs): ) def _perimeter_of_z(params, transforms, profiles, data, **kwargs): max_rho = jnp.max(data["rho"]) - data["perimeter(z)"] = ( # perimeter at rho ~ 1 + data["perimeter(z)"] = ( line_integrals( transforms["grid"], - jnp.sqrt(data["g_tt"]), + safenorm(data["e_theta|r,p"], axis=-1), + # FIXME: Works currently for omega = zero, but for nonzero omega + # we need to integrate over theta at constant phi. + # Should be simple once we have coordinate mapping and source grid + # logic from GitHub pull request #1024. line_label="theta", fix_surface=("rho", max_rho), expand_out=True, ) - / max_rho # to account for quadrature grid not having nodes out to rho=1 + # To approximate perimeter at ρ ~ 1, we scale by ρ⁻¹, assuming the integrand + # varies little from ρ = max_rho to ρ = 1. + / max_rho ) return data diff --git a/tests/inputs/master_compute_data.pkl b/tests/inputs/master_compute_data.pkl index 11ae1aa4cfc00e7da8856afa7914cb5b099a0c43..b356ddfceb21875df5431807cc19782caa6b0141 100644 GIT binary patch delta 5280 zcmds)c~sL^7RUW9$uE!)!2AFa1Qe0QqEx}6Ad5f&wN=!*jfNe<5&;9hwzlF@t6fmR z5N<2hrCO`iDoXWLt=l^kO&7RB7cX%9Y~I$Th9AO=uZdl!IZ%7rM?4o5S92@aTE@kI`DLSE%~d^T z4dsS10oHgn#`+@{X1&PA!ptQajrA5gh|b}vn0V`c&Qn-1&u#p@9XFtCiAJ=}V#6E) zD#~le)K)Q>^00LvV_;2_tqp7>MoVD^v!)cQnjPYh(mrbfot47-Sa)*aIAsSjf;9zL zr}A?49;#tXnzf#fcQmCb^VxB*D$yvC(XXjK7qKc$V?<&@aO2^*=S65+%cvVaRDC6a z?jWaWyk>3@8$X17!HutU|5|L!+P^-gD_>Ajo>^;B2idwhGhO6A=(&XPCW}l=qiEAI z4p6s0$(J;2i`gJT^B}bbZiyIlO8;{rLF=+e+?O!tMN)Z}#n;uB*$e1D$L$Cv@iUpz zA`uUA=qU>2w+1`d`n9JZ_lx*%MZOLblBx{#<#zLMFtaL#V9^&vvhfrCtVlj?mScoY z;4#LH+4!ogxf;oS7EGLbvxh|oo0C20n^|mK4Ef5&e=U;NuQ6+aA-6^4LcCV7TY`zx z1AbR9aj8~eyme*V8mz|IzHnqccUvS)fvPV>Wy$(z6~=h+z6v+FcZcd5k?_6Lr$ln} z2zy2(6CJoCB9!mtT}k9?Dy(weULGgjOXhGex`9K_*+J?iQR8k;uzL#kJ+9ux;5MRs zI4s1Nr^0p_Psy9Swgzo@*i)^7tW$g-*&RCUBN5CkDh*kXhVkoqsd+M^MTG?ny1`?P zbsT?HB!Nz9dHa6$B(HYEDwMKzj!l6Udv!{}MC{EiwC z5DnshI3ij^hd3d6#2IlxT#>Gb0g>0bA?}C=(hc!Mx+7kQH_`*~L3|NEq$kn~@ke?i z0Z1Pt5D7wpkr1RW@(dD+^h5e1VaT({0AwH%jzl1Xkmry{WH2%Wi9({0p-2o8i^L)E z$S@=UNko#6=aFP21xZDQBO{QJ$P37e$S7nql7@^y#vBF-wlSoR3K{S26X8{zqzKV)VX4@(TkZz&0~j4M9pMnFpqlb zWsf2sCW(4Vs*NHOKt?r{yx;~Ivn0&|m$gPD9mzm4iObropaD!1+zKuCCwkZNWW9C= z+oen4B}`m6hHSc?wZqM=hO>J%c39gW^QU>D%2?QOqbtei zn(ZX$3Lh;@>{#eZycSpJ1-2ywO!rN`N}gjg%e4Y{kIe4Mw2=NIvomhudv#<^tbp&$ zu{j&{U4LX0F`HQ%Wa|F5$*Zu1#4V`sbFW)#bRm0kau!(_+J2IIHs=rV5aXEcL7M#v zu5u7KywszvD%U{vIprs5q3YdSEg8@~_jDi5)<+GM<8nPo_xRjQ18m-yXCN>xcf}pf z=GujhK!Yal8Th0!JA8XpZnl=_#^zn)pmblR15`!iIufS|d94T8zx9J-A=7<-nPgT_ zek&&rW+v0Yd8GMLY?EW%C2X%l?@W7rseLm$>Km#RZBgtmcG|l?9B&*D3195>(~#MF z@_8PvG!7U^TvGG(JlvG>HDqyI{;h7%a46qEmQKol>k4PPp@KP+@;xErL4G3an;*wR z+l+o5BrHDn8;(dZrkg$F>7$`xbyg5@;Y_`FXx7E>@Sw!(M^?U_car%-BFvJ*A zy*hsqW4)qseq;gq*2C=(4WwXF`ltT#=>Bb`!`xWg=PXP561%CMAG0_)O|2# zo2`Plh;&cV&z6xS$nkwnm;JmS8oY`VU}o?jUveeQJWzn<&mse%G;BzxXMd=EJyQdV zlZJRHlSsqh;%xVDTFjqhX{6Rmyxq!G*{vQ9c5AR=k@y4Cb)kLBUk-59LS!%`$b+&Bvld!`JuKl=@c z(~wXX%W!*4y2M|eJ>pT|W=kx((#(YdtbCB~N7_~w2cw$vS*ngiSH5%Q$;zhULU*Pf zT2sw&&~PVK1Eq8PJV@KL;t)Y-RyHY@nxNE?u|D6MG`wui#T(#8zK&#MTb{SqLD%bu z*Tk$O0dUjtfThy9a2IIuV4P*IMDWqkuE{-#lS661c~dyuqD3GJM`xX-4Z+0Oq`IHQ zEOZhhrM9h)1sG(~K&3~7FUfE-V*uyPma}@Mp5U$)(eh}sSE}MQWbZ-CESJl}D#`gOWIJ%Y0JoiE-0(SF#_pTRsM<|E@$Eog)M!iy_^ht-{ zPK8fPT~zoiq_m)HV-Zg7ak{Lnbs!y~!bB}9Ow^*nL@g>z)S|*3>{N9OQGc9eCZ^|%z?~CDRle0(kX_WgI|d8xOUbS$%(KFahP42 z0AHVw+Fe|Yv_!+UwkN()E9v3D`%SuAf#*H#v4l0@NfMeV}rdnX~_~Icj>W~)b5a78xmmgLVMF|c479<5}8_e}O;SG1*n`-@$p@qZj4Svaly1_2rT78Rf$&4p0ok*>9zB0Q ofz9mmrDWyVAxTPqsq5R{OW~~SqWz_=)UR+BGnQOBAyx4I08$=9KmY&$ delta 4819 zcmd6pd015C8o)Wr%sDv9dIm(<+`=K*1WZL?6mYjtOU-Nmhh-Sq2am-~#3h&92FKyE zFqcQNO${Y|uaQb_SADz%rD>vqmUS)jqN#OTZudRkVNlDSr$6p`lrD$kc%Yq){3HI3m5BH$%rg7 z)ltO;2uIeP36?%faO`3w#|l1GnQ~m<;&hR7%GOO@x00KdUaZ(T?l3(aKe5qWjjEoX z937QOCzfy$V{|NFhmme=MS{b@B;wGANntHU$7p^y4!W>|SWBeCn@?g!IVSNTj-PnN zWr}8mx}q`cW2|MX!%H*LeXAOsnL7}Km>zte? zYKJL_KjQdW1JT$uUa&HPzbwIZ9}P~v+@yrit!4HZfa$G8{o(Q26BsC74ODiO*==x0 z-*`dTUKTBrJlV4nzm|1#gP%{bmpjC}Gr2D$C|}8*lOW+5kLPluR)c@TyJ;~;(oqf< zlRuI_j&TJwZpt(sucZ8Q4SM+56xc5ii6W@j-kM#aEB` zA^u2rBmn7w^h5%YAS4(GK|+yUNN=PM(iiE6Jc5KF{SgBajy#GCKq8PxBnpW}VvvDI zEE0#rBZH8~kip0hWGIq=3`2$^iAWM+L`EPZkz^zVNkv8>qmePlSY#YB9+`klL>@<; zK+=#&$Yf*+@+9&UG8LJIOh=wZW+2ZXGhyScXPH;Yn${v;676FaSWAo}lj%;PXXp8n zAiGIqo}}X+TZ<%7Z5Eags!6nCg#_j23Cv77j%>SUv6pJ5qYwX5KR36f4R#Zfj$|O2 z&|sf6Z#2_Fyc>%8fU%)KBY+T9Hd4@Kac*wbi|Aq<3p>A@wNI~K)~VUkDQ$EQ6XGMzvICzyTEE^2}~<2wG@2WL^grlhRx9zNe|ls4q%Hg8>j zOONH2p4Mi96Y<1&a5@hK8|Qqbuij?zf&*Dut8HY-8dCsdpUe9D2&gq=pLv9s>T-O+ znUf-^K)`;9A&7MBlGoa+k}&*d&_t9~Xoka+!)nLrws7b|<iLziipxaZQT`&}*8s#d2AvUK7>OYogM2l6lVLM?$l*y`gZ2rHNyzs0L8;DMv#0 zSwz5F`&JbR-j%;P1vqt=V@`wre*D!~FUDeT^ZdYao6{*G?H|@Q=(XQ+o#SwpOp3`0 zQC!b*ERL!MyGHk!-y&!p25I{&N4e^tToJCEvuMK@9bEg>QsT+#qR9MiE;nW59qt;T z&3OO3?`i1fq^u{EM3*^uJ|ypkOBYAK0D@mYV8`snu;0vfoF)m|F~_yby4bv{p2WRx zCRLopV+hXf0t0(w|G*LVi?Tn#M3eXalKni+e3jFc)Fk9}fmNgPX8S82Odv_yn~l|z zp4YBZ>Q^H9twmj!&hnLXO6k$Fviuin7;Dy)s^^@Kuq5Qg!T_kbQ#5TLzak!AWqZ5{ z8%DEf5o~@pFLNflA|5j0)6e%|R=}k=YfB$Nv9fj-$j$Q7{_v0KR=u9ooGNzTd!-mK z;yQQoLGkPW?^pq?R8uXLAn}d#ie=tZ^}4D;eF=?YfXM0+Nl%I!Gu@%iS+Z6Y35w&$ z=7SOV#cv4c386DHrWkP@q0h}u##C-@)rHmoIJqk$UVsf^Y3XakzNhpt(o`AmMbey> zK&Tv=zjQzxP5Fh=Z>aPX?~<}LE-A>}B@Iy}lD6JvBz5Bx@guucIvWdAw27o{TW%mk z*XEZBq%vI*3#|r0Soj+k5sVf21F-ntcTZQok{dDBgQ{q8&6leV5KD)LeF4yydvBs=c$IX&n?h&?q;)Qga}RF*e3gN&V^lx z_s-JZbZvl%3i>_d5%Hss2)($momaTcmt0OM3xl>VXLV7sZ9Y*jR#_HYNGXdE)b>SP zX#1kl1FuYy@ZbaEA@I%N`0YpRu>`5l#|$5{ub{D>M=cZum=gptK`PS0?oe~Q@aJD3 zNqMC3@ktcGZ)@3e9;7wXEW&HU3+Ij|n@<(`l8OqmJFR?zl`lc%hh-^x68K`F5AiIt zc7aBZITg=8_-?)JL92xM-CLy=Prk9J-xw+?#pA2QjwO-DLPUtzk!8f9N#$tKR36Hz z*2DSwqS253kJ`;dJ=pIQ%^Jv?hdrdrn1?A{W?>-bJghblSRY|c6I(g|?j3V44bxL@7T2h4z8Jw07B_c{RXrrl~= zEvRv^&HP~lP<>O?`bRg+{C>xBjO1ZO-{Ek7&6|gj`)ZzU#YmI_Cs^P-G;WDL5u^R}rstT|B>r6(6A%h5~=beke)2C}zCLUEB?9>{h}CQ<&nP-34ccZKdN zkV2WODF5m+Jf1<+UPnNy!>YetAdsAb+is zE7(gN`B5X1^j(-5H8eG~@!2g`Fqfcew>*|(g%mYw+Pj!l+$)zcii&2Xs<%M)<(N=7 z@SYsNv*F34YP2bsL=ahoUk}J>9Fs)S4#}tZG0E73nqy5t>aJ4{)9vBH`V*fca_xg> z(%PZ9cO2V~D)xZK!LZ7g;mO|}&f-(3;$Ko05wcGXQcrBjhl?2cq&H<|j9z|*qC?17lm@@*~4 zq>_vpwg)5qb3-S?$slF&*XBhbQqTGnbmZ~1I*I7K|mMd}`%Y=j1 ORh(jyq4}CThyM>?RfAOk diff --git a/tests/test_data_index.py b/tests/test_data_index.py index 33ecaaafb7..1eb3ea6639 100644 --- a/tests/test_data_index.py +++ b/tests/test_data_index.py @@ -40,11 +40,6 @@ def get_parameterization(fun, default="desc.equilibrium.equilibrium.Equilibrium" matches.discard("") return matches if matches else {default} - @staticmethod - def _is_function(func): - # JITed functions are not functions according to inspect. - return inspect.isfunction(func) or callable(func) - @pytest.mark.unit def test_data_index_deps(self): """Ensure developers do not add extra (or forget needed) dependencies. @@ -80,7 +75,9 @@ def test_data_index_deps(self): pattern_params = re.compile(r"params\[(.*?)]") for module_name, module in inspect.getmembers(desc.compute, inspect.ismodule): if module_name[0] == "_": - for _, fun in inspect.getmembers(module, self._is_function): + # JITed functions are not functions according to inspect, + # so just check if callable. + for _, fun in inspect.getmembers(module, callable): # quantities that this function computes names = self.get_matches(fun, pattern_names) # dependencies queried in source code of this function From 9266c10f6501de1cf56750b32bc39bf678c05c56 Mon Sep 17 00:00:00 2001 From: unalmis Date: Wed, 3 Jul 2024 12:53:26 -0500 Subject: [PATCH 8/9] Fix label in basis vector description --- desc/compute/_basis_vectors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/compute/_basis_vectors.py b/desc/compute/_basis_vectors.py index 31cc469ab9..44433252cb 100644 --- a/desc/compute/_basis_vectors.py +++ b/desc/compute/_basis_vectors.py @@ -3532,7 +3532,7 @@ def _n_zeta(params, transforms, profiles, data, **kwargs): units="m", units_long="meters", description=( - "Covariant poloidal basis vector in (R,ϕ,Z) coordinates. " + "Covariant poloidal basis vector in (ρ,θ,ϕ) coordinates. " "ϕ increases counterclockwise when viewed from above " "(cylindrical R,ϕ plane with Z out of page)." ), From 5c3d5f06965b7c48adf0091bffdc8131341e6ac0 Mon Sep 17 00:00:00 2001 From: unalmis Date: Wed, 3 Jul 2024 13:10:00 -0500 Subject: [PATCH 9/9] Add clarification to description of perimter(z) and area(z) --- desc/compute/_geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/desc/compute/_geometry.py b/desc/compute/_geometry.py index 2b26576974..36bb1ac32a 100644 --- a/desc/compute/_geometry.py +++ b/desc/compute/_geometry.py @@ -177,7 +177,7 @@ def _A_of_z(params, transforms, profiles, data, **kwargs): units="m^{2}", units_long="square meters", description="Area of enclosed cross-section (enclosed constant phi surface), " - "scaled by ρ⁻², as function of zeta", + "scaled by max(ρ)⁻², as function of zeta", dim=1, params=[], transforms={"grid": []}, @@ -418,7 +418,7 @@ def _R0_over_a(params, transforms, profiles, data, **kwargs): units="m", units_long="meters", description="Perimeter of enclosed cross-section (enclosed constant phi surface), " - "scaled by ρ⁻¹, as function of zeta", + "scaled by max(ρ)⁻¹, as function of zeta", dim=1, params=[], transforms={"grid": []},