From 9a2b2803057666173050a10c17c726ba4a590ff0 Mon Sep 17 00:00:00 2001 From: Joanna Piper Morgan Date: Mon, 1 Jul 2024 18:34:04 -0700 Subject: [PATCH 01/25] removing spyder autodoc --- mcdc/iqmc/iqmc_kernel.py | 7 ------- mcdc/iqmc/iqmc_loop.py | 7 ------- 2 files changed, 14 deletions(-) diff --git a/mcdc/iqmc/iqmc_kernel.py b/mcdc/iqmc/iqmc_kernel.py index 0e914dac..9b0bf296 100644 --- a/mcdc/iqmc/iqmc_kernel.py +++ b/mcdc/iqmc/iqmc_kernel.py @@ -1,10 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Jun 10 21:52:17 2024 - -@author: sam pasmann -""" from mpi4py import MPI import mcdc.adapt as adapt diff --git a/mcdc/iqmc/iqmc_loop.py b/mcdc/iqmc/iqmc_loop.py index 353f66a7..5da490e5 100644 --- a/mcdc/iqmc/iqmc_loop.py +++ b/mcdc/iqmc/iqmc_loop.py @@ -1,10 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Jun 10 21:52:29 2024 - -@author: sam pasmann -""" import mcdc.kernel as kernel import mcdc.iqmc.iqmc_kernel as iqmc_kernel import mcdc.adapt as adapt From c5f3b022cec3a5ea02ceb8d363bc2e54bfe405b2 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 17 Jul 2024 15:02:57 +0700 Subject: [PATCH 02/25] Update index.rst --- docs/source/index.rst | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/docs/source/index.rst b/docs/source/index.rst index bc0be10b..df2da497 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -114,8 +114,7 @@ To build the docs To Cite MC/DC ============= -If you use MC/DC and would like to provide proper attribution -please cite our article in the Journal of Open Source software +If you use MC/DC and would like to provide proper attribution, please consider citing the following articles (as appropriate): .. code-block:: bibtex @@ -138,6 +137,19 @@ please cite our article in the Journal of Open Source software doi = {10.21105/joss.06415}, } +.. code-block:: bibtex + + @misc{variansyah2023developmentmcdcperformantscalable, + title={Development of MC/DC: a performant, scalable, and portable Python-based Monte Carlo neutron transport code}, + author={Ilham Variansyah and J. P. Morgan and Jordan Northrop and Kyle E. Niemeyer and Ryan G. McClarren}, + year={2023}, + eprint={2305.07636}, + archivePrefix={arXiv}, + primaryClass={physics.comp-ph}, + url={https://arxiv.org/abs/2305.07636}, + } + + If you are developing or working with specific numerical methods please take greater care to cite the specific publications where that work is presented. An exhaustive list can be found on our :ref:`pubs` page. From acadc09ce7d268a07e10c3c645f5507f9fcae06e Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 17 Jul 2024 16:39:58 +0700 Subject: [PATCH 03/25] revert changes --- docs/source/index.rst | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/docs/source/index.rst b/docs/source/index.rst index df2da497..bc0be10b 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -114,7 +114,8 @@ To build the docs To Cite MC/DC ============= -If you use MC/DC and would like to provide proper attribution, please consider citing the following articles (as appropriate): +If you use MC/DC and would like to provide proper attribution +please cite our article in the Journal of Open Source software .. code-block:: bibtex @@ -137,19 +138,6 @@ If you use MC/DC and would like to provide proper attribution, please consider c doi = {10.21105/joss.06415}, } -.. code-block:: bibtex - - @misc{variansyah2023developmentmcdcperformantscalable, - title={Development of MC/DC: a performant, scalable, and portable Python-based Monte Carlo neutron transport code}, - author={Ilham Variansyah and J. P. Morgan and Jordan Northrop and Kyle E. Niemeyer and Ryan G. McClarren}, - year={2023}, - eprint={2305.07636}, - archivePrefix={arXiv}, - primaryClass={physics.comp-ph}, - url={https://arxiv.org/abs/2305.07636}, - } - - If you are developing or working with specific numerical methods please take greater care to cite the specific publications where that work is presented. An exhaustive list can be found on our :ref:`pubs` page. From ec93ef9725f14947f0936d0ebe80bfac15c6c24f Mon Sep 17 00:00:00 2001 From: spasmann <46267220+spasmann@users.noreply.github.com> Date: Sat, 20 Jul 2024 23:34:11 -0400 Subject: [PATCH 04/25] Bug fix: particles sourced on cell boundaries --- mcdc/kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcdc/kernel.py b/mcdc/kernel.py index ca31c9d9..4e689ee7 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1724,7 +1724,7 @@ def region_check(P, region, trans, mcdc): if positive_side: if side > 0.0: return True - elif side < 0.0: + elif side <= 0.0: return True return False From 99eaa06fd19cbebbf0ad67f8bef32377d6f51dbb Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 25 Jul 2024 10:38:30 +0700 Subject: [PATCH 05/25] remove unnecessary file --- examples/fixed_source/slab_ce/build_xml.py | 88 ---------------------- 1 file changed, 88 deletions(-) delete mode 100644 examples/fixed_source/slab_ce/build_xml.py diff --git a/examples/fixed_source/slab_ce/build_xml.py b/examples/fixed_source/slab_ce/build_xml.py deleted file mode 100644 index 5a04582a..00000000 --- a/examples/fixed_source/slab_ce/build_xml.py +++ /dev/null @@ -1,88 +0,0 @@ -import openmc - -rho_uo2 = 10.97 -rho_h2o = 0.997 -rho_b4c = 2.52 - -a_u235 = 0.05 -a_u238 = 1.0 - a_u235 - -A_uo2 = 270.03 -A_h2o = 18.01528 -A_b4c = 55.255 - -N_avo = 6.023e23 - -N_uo2 = rho_uo2 * N_avo / A_uo2 * 1e-24 -N_h2o = rho_h2o * N_avo / A_h2o * 1e-24 -N_b4c = rho_b4c * N_avo / A_b4c * 1e-24 - -N_u235 = a_u235 * N_uo2 -N_u238 = a_u238 * N_uo2 -N_o16_uo2 = 2.0 * N_uo2 - -N_mix = N_u235 + N_u238 + N_o16_uo2 -N_u235 = N_u235 / N_mix -N_u238 = N_u238 / N_mix -N_o16_uo2 = N_o16_uo2 / N_mix - -uo2 = openmc.Material() -uo2.set_density("atom/b-cm", N_mix) -uo2.add_nuclide("U235", N_u235) -uo2.add_nuclide("U238", N_u238) -uo2.add_nuclide("O16", N_o16_uo2) - -N_h1 = 2.0 * N_h2o -N_o16_h2o = N_h2o - -N_mix = N_h1 + N_o16_h2o -N_h1 = N_h1 / N_mix -N_o16_h2o = N_o16_h2o / N_mix - -h2o = openmc.Material() -h2o.set_density("atom/b-cm", N_mix) -h20.add_nuclide("H1", N_h1) -h20.add_nuclide("O16", N_o16_h2o) - -N_b10 = 4.0 * N_b4c -N_c12 = N_b4c - -N_mix = N_b10 + N_c12 -N_b10 = N_b10 / N_mix -N_c12 = N_c12 / N_mix - -b4c = openmc.Material() -b4c.set_density("atom/b-cm", N_mix) -b4c.add_nuclide("B10", N_b10) -b4c.add_nuclide("C12", N_c12) - -mats = openmc.Materials([h20, b4c, uo2]) -mats.export_to_xml() - -# Create a 5 cm x 5 cm box filled with iron -box = openmc.model.rectangular_prism(10.0, 10.0, boundary_type="vacuum") -cell = openmc.Cell(fill=iron, region=box) -geometry = openmc.Geometry([cell]) -geometry.export_to_xml() - -# Tell OpenMC we're going to use our custom source -settings = openmc.Settings() -settings.run_mode = "fixed source" -settings.batches = 10 -settings.particles = 1000 -source = openmc.CompiledSource() -source.library = "build/libsource.so" -settings.source = source -settings.export_to_xml() - -# Finally, define a mesh tally so that we can see the resulting flux -mesh = openmc.RegularMesh() -mesh.lower_left = (-5.0, -5.0) -mesh.upper_right = (5.0, 5.0) -mesh.dimension = (50, 50) - -tally = openmc.Tally() -tally.filters = [openmc.MeshFilter(mesh)] -tally.scores = ["flux"] -tallies = openmc.Tallies([tally]) -tallies.export_to_xml() From 61fc93f537e2c55ea8b1c4c8801398d9f641e3bc Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 25 Jul 2024 11:51:39 +0700 Subject: [PATCH 06/25] bug fixed --- examples/fixed_source/slab_ce/input.py | 6 +++--- mcdc/adapt.py | 4 ++-- mcdc/card.py | 1 + mcdc/input_.py | 25 ++++++++++++++++++++++--- mcdc/main.py | 2 +- mcdc/type_.py | 4 ++-- 6 files changed, 31 insertions(+), 11 deletions(-) diff --git a/examples/fixed_source/slab_ce/input.py b/examples/fixed_source/slab_ce/input.py index 4067209e..97246daa 100644 --- a/examples/fixed_source/slab_ce/input.py +++ b/examples/fixed_source/slab_ce/input.py @@ -61,9 +61,9 @@ s4 = mcdc.surface("plane-x", x=2.0, bc="reflective") # Set cells -mcdc.cell([+s1, -s2], uo2) -mcdc.cell([+s2, -s3], h2o) -mcdc.cell([+s3, -s4], b4c) +mcdc.cell(+s1 & -s2, uo2) +mcdc.cell(+s2 & -s3, h2o) +mcdc.cell(+s3 & -s4, b4c) # ============================================================================= # Set source diff --git a/mcdc/adapt.py b/mcdc/adapt.py index 4cb2ac92..3904f1fd 100644 --- a/mcdc/adapt.py +++ b/mcdc/adapt.py @@ -408,12 +408,12 @@ def local_group_array(): @for_cpu() def local_j_array(): - return np.zeros(1, dtype=type_.j_array)[0] + return np.zeros(6, dtype=type_.j_array) @for_gpu() def local_j_array(): - return cuda.local.array(1, type_.j_array)[0] + return cuda.local.array(6, type_.j_array) @for_cpu() diff --git a/mcdc/card.py b/mcdc/card.py index cb1d5d50..ecf66167 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -50,6 +50,7 @@ def __init__(self, G=1, J=0): self.uq = False self.flags = [] self.distribution = "" + self.name = "" class MaterialCard(InputCard): diff --git a/mcdc/input_.py b/mcdc/input_.py index d1eb0d35..3c4903b0 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -326,22 +326,41 @@ def material( for i in range(N_nuclide): nuc_name = nuclides[i][0] density = nuclides[i][1] + + # Create nuclide card if not defined yet if not nuclide_registered(nuc_name): nuc_card = NuclideCard() + nuc_card.name = nuc_name + + # Set ID + nuc_card.ID = len(global_.input_deck.nuclides) + # Check if the nuclide is available in the nuclear data library dir_name = os.getenv("MCDC_XSLIB") if dir_name == None: print_error( - "Continuous energy data directory not configured \n see https://cement-psaapgithubio.readthedocs.io/en/latest/install.html#configuring-continuous-energy-library \n" + "Continuous energy data directory not configured \n " + "see https://cement-psaapgithubio.readthedocs.io/en/latest" + "/install.html#configuring-continuous-energy-library \n" ) + + # Fissionable flag with h5py.File(dir_name + "/" + nuc_name + ".h5", "r") as f: if max(f["fission"][:]) > 0.0: nuc_card.fissionable = True + card.fissionable = True + + # Add to deck + global_.input_deck.nuclides.append(nuc_card) else: nuc_card = get_nuclide(nuc_name) + card.nuclide_IDs[i] = nuc_card.ID card.nuclide_densities[i] = density + # Add to deck + global_.input_deck.materials.append(card) + return card # Nuclide and group sizes @@ -1888,14 +1907,14 @@ def append_card(delta_card, global_tag): def nuclide_registered(name): for card in global_.input_deck.nuclides: - if name == card["name"]: + if name == card.name: return True return False def get_nuclide(name): for card in global_.input_deck.nuclides: - if name == card["name"]: + if name == card.name: return card diff --git a/mcdc/main.py b/mcdc/main.py index 7cf1b8d9..25f0e604 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -438,7 +438,7 @@ def prepare(): # CE data (load data from XS library) dir_name = os.getenv("MCDC_XSLIB") if mode_CE: - nuc_name = input_deck.nuclides[i]["name"] + nuc_name = input_deck.nuclides[i].name with h5py.File(dir_name + "/" + nuc_name + ".h5", "r") as f: # Atomic weight ratio mcdc["nuclides"][i]["A"] = f["A"][()] diff --git a/mcdc/type_.py b/mcdc/type_.py index 2714c697..3dfb9b2c 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -315,7 +315,7 @@ def make_type_nuclide(input_deck): dir_name = os.getenv("MCDC_XSLIB") for nuc in input_deck.nuclides: - with h5py.File(dir_name + "/" + nuc["name"] + ".h5", "r") as f: + with h5py.File(dir_name + "/" + nuc.name + ".h5", "r") as f: NE_xs = max(NE_xs, len(f["E_xs"][:])) NE_nu_p = max(NE_nu_p, len(f["E_nu_p"][:])) NE_nu_d = max(NE_nu_d, len(f["E_nu_d"][:])) @@ -624,7 +624,7 @@ def make_type_source(input_deck): if mode_CE: G = 1 # Maximum number of data point in energy pdf - Nmax_E = max([source["energy"].shape[1] for source in input_deck.sources]) + Nmax_E = max([source.energy.shape[1] for source in input_deck.sources]) if mode_MG: G = input_deck.materials[0].G Nmax_E = 2 From 71315dd975ffeaf63e4b18bf4408857164d994ce Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 25 Jul 2024 13:35:14 +0700 Subject: [PATCH 07/25] fix local j issue --- mcdc/adapt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcdc/adapt.py b/mcdc/adapt.py index 3904f1fd..4cb2ac92 100644 --- a/mcdc/adapt.py +++ b/mcdc/adapt.py @@ -408,12 +408,12 @@ def local_group_array(): @for_cpu() def local_j_array(): - return np.zeros(6, dtype=type_.j_array) + return np.zeros(1, dtype=type_.j_array)[0] @for_gpu() def local_j_array(): - return cuda.local.array(6, type_.j_array) + return cuda.local.array(1, type_.j_array)[0] @for_cpu() From 78bf3959b30d787e1b6c9b689ff321a78516d49d Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 25 Jul 2024 16:20:49 +0700 Subject: [PATCH 08/25] add regression test for ce physics --- .gitignore | 2 + mcdc/input_.py | 6 +++ test/regression/slab_ce/answer.h5 | Bin 0 -> 128632 bytes test/regression/slab_ce/input.py | 72 ++++++++++++++++++++++++++++++ 4 files changed, 80 insertions(+) create mode 100644 test/regression/slab_ce/answer.h5 create mode 100644 test/regression/slab_ce/input.py diff --git a/.gitignore b/.gitignore index fc02ca5f..284efac0 100644 --- a/.gitignore +++ b/.gitignore @@ -55,3 +55,5 @@ docs/build docs/source/pythonapi/generated/ *.csv + +dummy_nuclide.h5 diff --git a/mcdc/input_.py b/mcdc/input_.py index 3c4903b0..4aa4db7d 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -322,6 +322,9 @@ def material( # Set ID card.ID = len(global_.input_deck.materials) + # Default values + card.J = 6 + # Set the nuclides for i in range(N_nuclide): nuc_name = nuclides[i][0] @@ -335,6 +338,9 @@ def material( # Set ID nuc_card.ID = len(global_.input_deck.nuclides) + # Default values + nuc_card.J = 6 + # Check if the nuclide is available in the nuclear data library dir_name = os.getenv("MCDC_XSLIB") if dir_name == None: diff --git a/test/regression/slab_ce/answer.h5 b/test/regression/slab_ce/answer.h5 new file mode 100644 index 0000000000000000000000000000000000000000..993ed042a94bff4c1eda0b5da18736fbc73a56eb GIT binary patch literal 128632 zcmeHw3w#^Jwf;JeLA;6+ponJ$Ev-s#{&fko=D7+9U?>Eaz(u9BvLY*%h%E^{5+!hp zAwX3Mt$9|Vw5Ts|TW;%?LLDIB(3We;qb9rzwA8myYG~-+{CLzn3U_wre3G@&&PIwX z?>eje;xjYn+u50K&YYRqo!$L<+maJ^+j-xeCFSRENOPrTHHbeS&@cT~gPMdm^c@u7 zT1m=KynhZ4G)T%nNpe#9Jj!2-JM)YaPLLQ9`9soS(pVD@As?jMa{FMQ?Sw_`N=>~~ zT{p>LvNMWG;Yci%^mPUUD|t%N-`%~21>>nmGTakXZBfy*aJAuIJ^O8IMRp5+t+itp;k++$1aQ3+&%MM`;D^(xVx0SrJl^&gjzmeUE%n-`eMere z81-+Xacp8Yjz@b~n-#nCl`;T7fN>1{sE_Oi`9)HJ?r>+&hj9wM$pb)zlKX)ssg;f! zJdFpC@93SR=;gAVq60L(0jKpdD1Vmh8n~8|c?O}wxrDvDslNcH?LQbFv-3CtxFL#< z?ST|3V4a1YAMozN3z8-txF2ZJoe#7TaV{!<^I<&Tp?DiLREXju$MQfQ#Rr%0z!=4C zVJ!yW2Rg`{$Eo;k4-e?snYLf);o3jX^NCSEr|c!}2O9G8fflk`^eU|#<2YLuGteM4 z@wEgI2S&kTBU|Sx;b!LexbO3*88oq&l0Ik792M%Im2qQ3TeAwtw+xR-QVs!AOVk^M z+yMcBY(AJP3&&2#fHx?96adW{4;Us32TVy+f>RWQ>^TNS)23^3P@^aT&6?FTo#UpQ zQH5#a<+86V%`nGa^g;RXh05=7 zia$doAEWdavlsf0}$C51E3d-l}5IQoG^^B+wL9O-8{urf@Q)d4c~B3r{(j28v_mdYGb8 zdbN0#g^=ztcwQW>qSAPECoW-y7lu(XwX5vTQkB6Sd@5lz?O+a{#{_&TByzs$)jM+k zdSKu@EeegtOqIK zc`iA>4mo@m>F^I|K+<22bm;JuF_pEH+J{>$UjtuW<;w~_kbrpg{Alk`BIh&$Ga`-n zK)v;3^?ZbM>}7l9BW4jkiv5SyXsfpez5Zc#|G~#e^bgd3A2tnA8dUkubsSglpFFv6 zJts}=Pa;ojHXkx$;H1T;pQ_f6iVYr6R;-4xhUiT7^sMpG-3kVpdZ9!)23H;U1`~E4 zKDeSg7)XX!1$mB`w_-&e>-I;23(^wiv5>3$fmEtTl7NE*p0@_BwMkNsKN*aN{oS~& zg7ch7`vFiMpFf&PvJC~2lwU;gVoRt$MI~j8r2KHcGlPLpBz%5~ zs;SQR6+Oz4h%S$IvpB2+G7(M113_PymENBWM1-h9gD6G#tUk$WPjUHfP9J_ zuo-|8 z52+_UBF33rTY4R~{B@zzXyEU$(W87#%g%3iY zyo4X{E-bG{Y6k_JM&)B{J0ac0-nub*wOgk7P%PNKlc^!(8&6c6}g$y7Wj1wvt8CreXgjAVjI_C&T{ zc@DcXoJhvQ%Tvm|O%|0p6Ft7A47ehkNGM4p{mZ-gozoTF{x0SCh@U;Yt;m#`y&_jo zTse&rS71^i;7_u9r&1yqNrc&L=5SIuD@!lLeM;FAv0zXsXEMs35N9%VzJkJ>Rg%zu zc_@QIO2j)oa#KQnO@P10qEYbl`@&L)K@(_XGH$=Mr#R3phjX-@O0_Ei>B{#zn zmxoDXl&V~KgVhrFu zNC&yO2`c}fvKV1Mz55ko*Kz&WMO^CPwD0lYXsP76{~_AI3U90Tz5pC=OYmOn=l);N znFr)7K{>8BVZ3WOU0iIctd#kJ7{_N_l9W9|jB}B?f0-C(c5UfAoHL64S^4`U{$Xjm zsy!dC+kCK2P$HyKP;X^yVRC`}}KdLWAcC z!yZZ+{m)X})RZ=gjGsL<)q3I5exUUWwYt8ykoo~E1YDV`AC&S@>=#gK$m+*UYwYVP z<%4?$?YshfXj`}T{G@TJ8Ih%Y6x#>t4_STv54Xg)0YB`e)|f;ovdhJ|EUEhgVw~Bv^|PyUEn!!sd^pr5EX*Hkwf^cP zK9Dmk1(u5V%*s+eoSYAO4JAqsS^fMvXi@u^Ekt!C_ z$0(X%Qul_%xKUE~ohQbPkvg|hj2kEQBzw12J%6SB%%PsX(64IccT$${J4XAuhG^o( zdtW_U%08DF7%7v5+_40s+;?B!?po8MMZ}32d;+boB0Dg#)ZDA1xaQ_7R zx<;tnPMYBB-RGDb_ip#n{Tldj?R^yZuRhwR3Y>O- z8S>v=r15xI?R%70>CCILrluC_T|O&d`$kgY;x?1IW3?Cu`-`Q;I2WmV`Quw+u!*(P zXIEdfq+LxPAD7#Fuzpb3-}OzKk0L&9u=yzB<9jwAMSR?1^HIdd&uu=6__))OeCVx9 znt#PA$9~Xzc{Kk%a#t@UkO&03yOm8}STD2{t7{g+ywgoN_N^6X&5ES{aul1Cdg#BzIJ0Yw=UJ8BAHG+|{&QxhLM>k@AKb#U`31K2s@)!%1$fOx|S33M4O zfB>5&S;_}r>7m}0KX3BXEVBmc>yIx)`dTH}6J?RBMp>hw2r+m$gNwYu;07V!g*2205xUN456xuUTpz z)$;LyCHa6K4w1DwE)$zAW)&S@6XSYG9W&?-Q1^aaoS)gX)7Q^-oX_j-&`@CeVTF*! zhTVxDtNYW`dik?`IvtUDpAUw=Bu4z=q~c=nlQf}fv{b2~_ApS~-^<$`w| zzfz3Lk$UnfF>Zj=j&F%^L!?e$EyfL#dgNL$uGZR8yW-;#yf5si^7lu(b?gW&QR|x2 zCO@NnJIz#l7o9Hi10L{Dyt$tTIw+p~ZyrGJ>HRg2qrC}c3yUzYwf&{?_x88yjQ{5D zFIaDIufZDe+v0}AsN!8O#yLqHx>1aCkvjbyG0yDT>9Z@3CGBeZ_-M2FV4b0`@1x!3 zqlgc$%|{U*9X20De1vR1iuj1xd=&AKwj>|e*MW6yn$=19Ie*V3^j;p7Fh6tB`iGLG zR3yAA7*FJ1_=9z%2Y?z$%sbtbU*;xp#x11IepifZC3XMz#5fPB2X7YR%&s+_XH|L~ zvscIdb7rSPEng`g+`^`R-$#ac!wLy>p(i3g6yv<49{Z6P*FoxuAB%Bj*BbY;YVTtY==3*p+m$gN zwc7VFNc_N~!BR2^m_i!!QLTM!n5FhnEgzdK$p`#!2-bQl0_7!U72cne*CVx;0->KG zrMyIMN&l?8-YTin*Uuibq+J>FQ7b?E*evnGm9CdI>)1ukDqTh`=OQZY@g(hiAGK+| z(*EaEGp_xQlc;&tmjGALrdXN6u2Ex!ZWF zGhk#Nv5mZjA*)x^E`L=F@Nuaj|6ZLCz5fgTRQukZs!p;U9ArE$RXV?Z$I|m_ZdZ=H zU3J_eHUiAv2JaQ)T%;bpPmFVudUQyPYbEvgCNZv!)RXs%ab8k8ekaDk-wxA#NJS z>LW-h65q!{K81Kf)hFZrNTS=HWCy8n%r9Hj*w+6sYQILxsvxVqpK=tnBe#!2 z2{VPD2|`BeJI}+sklNAS1Vy)tP)7#n`y9M8vh()%lfih{-;D(7uZT*@dytd}_QgTL zx^R*`&xoj$h^FF7H0kf|RvukU1ere4B_)G_P$YbQN_o~%`JunNc~tkWkv(T0<%eo} zrEZ3{VpcJr-OmYo8h@y`niz+c`V>oM-@#Gd>xGZiFevP?rT8%2nD#nJrinw5Pyd5F zKCCH+Zid^bqICS`BNKS6e7<$Y9-K<;2$!=2!`!YK^Vdl-kBEzn)!5h{#W?uKi4iew zwrFc_SLayTu1fjf7J+pi`Ws~B{cgDqKjvKL@OgzKVYLc5TkPMi$F?z`omXHq@jOwS zFZxS+adOo$z3(6l(Rs(7;#nKbi=EGM9LimFf$Re<$$-{=G5$+W7gvX`*Yck=&(DKB z#jvf8>Dg~Doj3Ok&x+gcTwVVg`xD2x=*;6U@qmZou`N6hqWI)%JkU>ZTUeU`*l!0p z*a50e$D2H$XJ;e-;(GVr#JHgyc*bpnqwMLsGKKr!Ub0*Gd95AeI9nDoP&(h_BN2G4 ze13VJ&ipc0d}5Qv__}UxeqGn`oVcLPq#pdU7ze*T{Ja?FBK7EB#5l8SjqR#d`?%v| zZ?IS}lR>}~vXqbFbsg?oL1I5-;GG-Kc5E_YH4==!OL5yf33!shlmg43`Wew z2j=0Je}|};!*hALq4!d{hh#`IuZES=>r5ql-URO}?V<1N`PYWk9J(kUu){|+N)(7v zd)BsJwfFxY*0JBLpfoeD!o1T>MHwFxm#Kx+lm9EmwUXNLf*9u^wev+W&g|M!yW$p! z>*w0>qIy611lgI@=2^AsXD?c^pOxB2t@zkt^TC`$VV?PhCHa6KZYM`Ksu>@({?AJZ zE#q2EalC_~o|lW$wXqN-w7;U&SSxPa&uTS4dyni4GZ65ULBJH!n2&1bJO8$1Kbt;2 zn)c)k&i#6!@y$K|n?62vvm_tf4?Cn7HDJd-3N2YKl%~YNSIcP&X}uI2{%0YrTnOG{ z^#7OB7Sj4CI6hHGD;I+Gd0sE4Eu`fr*#1T#tz4*|5{KR@r!AySpIz-`NxPapJ`S)X zAFvecyIdaSJ&8%55>ljcEkKdVwSZtG81G^UUGZp2>6_t1G#QV^!T}amezACTxgwL{ zo}eOA)GrqA9Ev8286{RDG# z_Y+9(@yP z@+Wmua_If?dZe~dpnr}yKQG}18^ySMfBgVwFqbs8t6Hr`k0m?P@`V|tIe4ibAII5z zuqH2D7n5y1iumx@d=&Aq$mXMnkK=7ViuhP;^HIdd5=-#Gimd*UetyTCJy#!C*&U2@ zB|~cZif+DA?hHl};p7_j+hm%SM zJ;8*U5oR5j->I2b@dxux7ai;Uw>XmyQjb5w11N`f|E_~@6K_*MN3>FUCco3YJP@<_ zw)e{GXFApB`HzX6I8)AG05p#5fn>yt|2U@MB{RF%JHzbr&%X{l9HjF|IrdrKv`C zRl9vKuaaG%rYMWK?5ciz{M+V(nLuH^yBoPK+dh<@(b#Z+%}0@aESMF16y95+YU;|KEVS(X`kbF6K?GAVw@w&AY=mR`H`F#S+kSLv3In*AMXI#k6 zP))b)v{=7%5CPfGfr;|Rhe(|}M2s7v!(E>g<5~#k-cO8klDhdYF>ZwNOMgm?ga2-J zi*b~lp)j_qTHUV;lAUSGZLSfdetfK$C45x+-u5o4H}_u!1I-nJCJq_fN2Pp3Z9Z6g z7uK<>XE`6$u2(M7;lXw_mU6vknE>p?$muYVEDKI;|9K-CvAL6E1eJ7&k_^>;f^4 zV5SxI^9xx^_=WoMvEJr`^@_rL4){MS81>7=;6jr=l(N8X>zp6z@Q^d@7&WzEwddn!mgWQVASaRJ*hf4c!7OibZ!r$N zqVHp3T!_x=bcu0O+t7^R+S}FbmbNQnK5DgIyqkE^+RyfSJ@Ou#59S~W^TLg@gpW$c z>3c2hKlr$`gXAJ-@7!10zl%qKG^s;Ji*bXb9+Snm45`z{h;cbmXO0!)Mya@&{}AJf zYqjk$V?SJ}T|HoFyDH^_&wQ|d6s-e^`>~bpOLr|orIdTN&V;cD0CPgDIdgGS%R zOc`-=N$!W6B-&mgwH%I~QNEmQr1l&k#$gqm+e?gVqtmn0D_`pAPc8L3 zoB@=!QL5E^`0tkJXJ8!9_w+7vfPOsY!|OrlJA^*=nNnU5`p$(CFI18p zd-7tAA(_`PUgJUba8x3y(jzVRM^^fHReZ}+D^}pVYC4^KYC8LFs`|uLc$MFtE z?q6Ci(E1sEUgQ2Qjy0kJP-3m&aNCDeNOH%s{E%;g*2=+tjTj$m;#-spN;yV?dg@)+OBU z6yA@N^@{szD^;Wai(;IcaGA5kI2Ym4|0%}7FG^>KabVkfrWl9*A395nqv9Eg`q|aF zmawb(@zHJb!K|~e-o4P~qlk|yZ9azPW$_F12F?!W*f6P6^2Xba)sXsH$ zTjokbR61!-US5=1`#zsB!nw`^0tXh$35uq>#W>s#>l+s1igUM4taSX=+FPa9BR7&= z0_)0w<>LI%|Azu%oX%-VlIv$zPg=sR>c_|5Y(AJZ7S@$- z+I$r8(X==3tnduzx)}saA?w#)_Okh4Enes^hgg&kJ|bfDYGm)>Rl?0InQJ_%T{p3q zl0Ih+|KMAp-bYZqIa&p)e{zMOi9@CRMaxI6?pq&AJb`JDl0m=}67x6*-4N;cy|%o? z#cY8(jzF;O#c&T4%*X{6xD1o^o_x}!Qn2P22S}|YR(zZ}S^BOU3oN%pI zh;e+j%D`zc&O!P0Yrc;`$$8OFq0x2X{00a&eyJGOLAWuU{n?~uYt|-p-OB$c&u8-z zBX{igs{H+uAxqlT^zo6NrF>L+JtRXspl84ma|oDP!vA}zl&vJ&HT?|Vn@xthlMB+4 zq)PTshAP=Z8T9?x1^j&sh4&)3bZW0Ru7tkLdunDqdTen9#`auk`}JA6{c6ugc9!!| z?K(TB!%Mcy4=lSS!0e)d6g_5DPh;U0D*-i?$9 z#QiHYpwH63U(1<>YH(w;i zX%5@wsIecbble}*>94tME^9q9s68J;HXp3v3jTc9=A(#@5u1-9K1OFbAGNwpK1Mug z$BRro+~iiDqtu^z)ngl)c)3yfT1i}os{FmW|Eptn7{Rx?%0gf%A3Xn=f1dqio%-|k z5>{#(l-MVcr951{u!QoN;y#{>_A3Cbeb0oG6#aC)1vu?>kWo5~`|CiIKko=7J2PUO zgVedJ#JGGf{i+yOzBHw&rFLc3e)gJ%};N)m@m^t{RA9v>d3J zUb`|cA3NE6FcT=O%jemA6!EdojPPMJZ)yGx_j~Q=zP8-17{Gny0jhPHWH!g`Lht_X z5bYBhdt8k3KF{^(F)l%FsUPR#2vD_qj{YF(cgPtf$RJ<}iG4+GDoyCS+VU0`6Z0VN zO=6sv4)gv%j0@qsKNRD5D>3kmVw{KaOMg#{Ya?9dW-$)_sqOn>T>b3oU`yCl{rLE( z%?GpP!u;oBHXlWN9AXJR$Q?=8M-^PLg53gVpVL{82>6xT+{#CD*sqXIJF4hFqkUJg z=xX&=x!yjkRuo_>W&Qg1$1Tx*_2c6cHXqDJ3;xb+2|kS0L#}>4zrZ|fE4!f$n;eOct{S?^Ao71ItZ``Z99_Eq|mfY~c`7CK7VhpXKW&k?nMAL=LR z)L22F2pI%SA#oqCl}a;sCzvZ=JgiHaZx!QO=Gq;Oz z@H^RGi*e;yC{3-OU8R4@OAO{v7IWDZE0xQan{z+O(lfz<$Slc$0jjJ*h^|+$V3?&v&8w zqI={Y%zx$!XM69Gw>-J`RiD~;k32TK?ki2#Ju3gLSn`XCkNW-MHbK^8Qz?I_q=eSIFy5nX}}fO_#}c zHg7s2IdPfnfBHjD?flhs@|S+mwBn+pFP1NQ_S)Ni`urvGeu-OGod1=Jt9_>Ou2fy+Gdo%XbGCt-MTLa>o&$TpPJge*fw5!7IOit$bV4F7N+x&JFS@sl|7X zy?ecU`#}qXKXzX$zxnhTH(xXV8hPy8_Q}4Ru8|Mvy#L7c*WDoBeMI2+!!El)K5WTu zlfk3DEid~`dtbuYC-0tYyy1jf`{aFgIrsXGZ}rJrCcKaA`|+#gFWtU=+4yBw%j3J= z|K6cryjni{#!Gg+`Fr1#FF3{>*nEAT{P10;wePw2)$$(Ke&fJXMy`W>^X_H;wHdG( zm=z4vZ$6Z<#C)iJd}M7tSPv>(FX^}WDB@$#=A(#@VTNLeg5evsyELEwU(q@A!y=I?fDp`eql{M&{Ro3zk;}$J*XKf)x zJ(Nj+a8Il|90(_U0e>u+QXh3xHL++c)$LcF-t@8OIFoUD98-z~_?5_5z_-G$JhMsp z@cezN!;#MDYL>QI`J+mf)|m|Pn3C@62_`~PXSBy3j`%v+Q=&bwXhMBXRKaw1@+UBv z-Ev3L~uGW_aje5Ha3D~VK|cr?`=WVh3#)p1`k#2zLM zMY}tdV^}rFA_{KCRFOwxbl!|DxOeC zU~$s2Iz?E0=mITE&!eVlO7TAS=&G7e37%_|QmOjDnm~8Zmq^Bg$v~dh<#B%`5b6#l z60E(t*>kGY#GQ#AUo_r{$7ZF}`8*%SS31Uu%&bC9P~wr|_!^CVjYhv#qhC8!5AW&W zr>pp{aeuD%IrTA}ex0qjjP*QV`uG^P`CtuMm@iM*d=&99X-PitoPCZQ(dd2fNacGi z9*-KnSBx7Z_4tq&H$m$5--&U5dNYoAAj-|W}mw9mI9|LsK@jl(pJkRYpHA2!o<1+*9#HrM(P`_Wsd zc&&eeMWMXKvX!E3e-`81gmXV9#=%Z=Pl|D1aq?*~4*kaQj2P#l@;RRs;|OLbjO=PJ z$^o)^yK<3T^=ryxF1s?~1N{rKdOn6N(a-9~$EeK*>luanYvVQ_MSM)o3O*Xi`xRuh z_vanA@eu$c2i&T;{Hf8pMEWbQ1;$C#!j{_y1KM>7tT{VYJiqA)`qlkJiO%3Eg&Al) zRnk;CuDR8m8u({>^YiFEDSwRjh*v-EyQtrKPz5?ey2_f=<}j?&9F*OJhJWo@&|lkL z&{RVU+|J@NXICBn&0E4;cBS{es@@lQX>%`eX*5o_$2ku6xt*ESc3)lXIMGIo43YhS zpky@ZSH6hCHiyRi>LypdgEtZGp?$g>$-a>iOR Date: Thu, 25 Jul 2024 16:22:04 +0700 Subject: [PATCH 09/25] back in black --- test/regression/slab_ce/input.py | 62 ++++++++++++++++---------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/test/regression/slab_ce/input.py b/test/regression/slab_ce/input.py index f020b873..8f1208d1 100644 --- a/test/regression/slab_ce/input.py +++ b/test/regression/slab_ce/input.py @@ -8,37 +8,37 @@ os.environ["MCDC_XSLIB"] = os.getcwd() # Create the dummy nuclide -with h5py.File("dummy_nuclide.h5", 'w') as f: - f['A'] = 1.0 - - f['E_xs'] = np.array([0.0, 1.0 - 1E-6, 1.0 + 1E-6, 2E7]) - f['capture'] = np.array([0.01344, 0.01344, 0.00384, 0.00384]) - f['fission'] = np.array([0.06912, 0.06912, 0.00619, 0.00619]) - f['scatter'] = np.array([0.26304, 0.26304, 0.15024, 0.15024]) - - f['E_nu_p'] = np.array([0.0, 1.0 - 1E-6, 1.0 + 1E-6, 2E7]) - f['nu_p'] = np.array([2.5, 2.5, 2.7, 2.7]) - - f['E_chi_p'] = np.array([0.0, 1E5, 2E7]) - f['chi_p'] = np.array([0.0, 0.0, 1.0]) - - f['decay_rate'] = np.zeros(6) - - f['E_nu_d'] = np.array([0.0, 2E7]) - f['nu_d'] = np.zeros((6,2)) - - f['E_chi_d1'] = np.zeros(0) - f['E_chi_d2'] = np.zeros(0) - f['E_chi_d3'] = np.zeros(0) - f['E_chi_d4'] = np.zeros(0) - f['E_chi_d5'] = np.zeros(0) - f['E_chi_d6'] = np.zeros(0) - f['chi_d1'] = np.zeros(0) - f['chi_d2'] = np.zeros(0) - f['chi_d3'] = np.zeros(0) - f['chi_d4'] = np.zeros(0) - f['chi_d5'] = np.zeros(0) - f['chi_d6'] = np.zeros(0) +with h5py.File("dummy_nuclide.h5", "w") as f: + f["A"] = 1.0 + + f["E_xs"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) + f["capture"] = np.array([0.01344, 0.01344, 0.00384, 0.00384]) + f["fission"] = np.array([0.06912, 0.06912, 0.00619, 0.00619]) + f["scatter"] = np.array([0.26304, 0.26304, 0.15024, 0.15024]) + + f["E_nu_p"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) + f["nu_p"] = np.array([2.5, 2.5, 2.7, 2.7]) + + f["E_chi_p"] = np.array([0.0, 1e5, 2e7]) + f["chi_p"] = np.array([0.0, 0.0, 1.0]) + + f["decay_rate"] = np.zeros(6) + + f["E_nu_d"] = np.array([0.0, 2e7]) + f["nu_d"] = np.zeros((6, 2)) + + f["E_chi_d1"] = np.zeros(0) + f["E_chi_d2"] = np.zeros(0) + f["E_chi_d3"] = np.zeros(0) + f["E_chi_d4"] = np.zeros(0) + f["E_chi_d5"] = np.zeros(0) + f["E_chi_d6"] = np.zeros(0) + f["chi_d1"] = np.zeros(0) + f["chi_d2"] = np.zeros(0) + f["chi_d3"] = np.zeros(0) + f["chi_d4"] = np.zeros(0) + f["chi_d5"] = np.zeros(0) + f["chi_d6"] = np.zeros(0) # Create the material dummy_material = mcdc.material( From 62d969829ad77e2eb7320079acdf6533ac70e338 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 25 Jul 2024 17:31:16 +0700 Subject: [PATCH 10/25] add test for filed source --- .gitignore | 2 + mcdc/__init__.py | 2 + mcdc/input_.py | 35 ++++++++++++++ test/regression/slab_isobeam_td/input.py | 1 - test/regression/source_file/answer.h5 | Bin 0 -> 162064 bytes test/regression/source_file/input.py | 58 +++++++++++++++++++++++ 6 files changed, 97 insertions(+), 1 deletion(-) create mode 100644 test/regression/source_file/answer.h5 create mode 100644 test/regression/source_file/input.py diff --git a/.gitignore b/.gitignore index 284efac0..6b6cd7ed 100644 --- a/.gitignore +++ b/.gitignore @@ -56,4 +56,6 @@ docs/source/pythonapi/generated/ *.csv +# Regression tests dummy_nuclide.h5 +source_particles.h5 diff --git a/mcdc/__init__.py b/mcdc/__init__.py index 03f94118..853af00a 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -22,6 +22,8 @@ uq, reset, domain_decomposition, + make_particle_bank, + save_particle_bank, ) from mcdc.main import run, prepare from mcdc.visualizer import visualize diff --git a/mcdc/input_.py b/mcdc/input_.py index 4aa4db7d..cdf3207d 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -1958,6 +1958,41 @@ def check_requirement(label, kw, required): print_error("Parameters " + missing + " are required for" + label) +def make_particle_bank(size): + struct = [ + ("x", np.float64), + ("y", np.float64), + ("z", np.float64), + ("t", np.float64), + ("ux", np.float64), + ("uy", np.float64), + ("uz", np.float64), + ("g", np.uint64), + ("E", np.float64), + ("w", np.float64), + ("sensitivity_ID", np.int64), + ("rng_seed", np.uint64), + ] + iqmc_struct = [("w", np.float64, (1,))] + struct += [("iqmc", iqmc_struct)] + + bank = np.zeros(size, dtype=np.dtype(struct)) + + # Set default values + for i in range(size): + bank[i]['ux'] = 1.0 + bank[i]['w'] = 1.0 + bank[i]['rng_seed'] = 1 + + return bank + + +def save_particle_bank(bank, name): + with h5py.File(name+".h5", "w") as f: + f.create_dataset("particles", data=bank[:]) + f.create_dataset("particles_size", data=len(bank[:])) + + # ============================================================================== # Reset # ============================================================================== diff --git a/test/regression/slab_isobeam_td/input.py b/test/regression/slab_isobeam_td/input.py index f93d87ab..f566a232 100644 --- a/test/regression/slab_isobeam_td/input.py +++ b/test/regression/slab_isobeam_td/input.py @@ -1,5 +1,4 @@ import numpy as np -import sys import mcdc diff --git a/test/regression/source_file/answer.h5 b/test/regression/source_file/answer.h5 new file mode 100644 index 0000000000000000000000000000000000000000..0a4543de8b39a6de4c288f4209dc7e71ff872cfc GIT binary patch literal 162064 zcmeEP1zZ$c7hg)T2^BC=R6s?#<1D>&hp4~;OG!xxn5a)dG4NT~h>D1XjfjekjR`i2 zfsNhz&F~!VZdVD(5~OD+N(+mMNr?9gVg-&QQ{tHs5y?a_HX$lLERscN760#SJvUV6rD zSG;h}c4Rp*$j(vmG|~?U)r<$!kltcDPIFShA^nX)DG-|F}K%E%AAWp1T zGYr5I1OK&b@Wd2c6uz@oL;~he6lW6DpnQq~_wFmZdVAdUeZ)!fep$k1>|Ps}ghm84 z<^B6W{R`rnl4W!*#%1C3f8GH*OdkL59iTA?r!a2aVjC;IQ+Aj;AJ>em3Gct$I5umz zT=rN`hzo$Vkd!AtOOt(E1TuaV>30Q6uY{0b_VEDaKjN~3YLxyG;a{WtH|V@~sJssd z_YvVfBit8+`-bg`PTCobPS!_=PF`P_PS#U|PBI~iaH0q&hH$MAt~J7mBb+$GwL!SH z2-gnb+9OPU8^1PFBUj{d1o-#YobL^g8YOPtV3)1s9^dL6(UP z6ialbnb1F-5_|Ym1E}xch+NK51NjGbb)&|~ zqwM_d)Hr_6ZT@kka|<6=HupH&`w!!Tv=AXQp1`(`r3hg`?69Of)XYP)+R?PWBVWRM zKH|@Vr!>AYq5Afs#HpVUIlZc>E9L_y;CR#d;eJAiOHRUp-y$301K&3<*>pZadFac3 zm5-Pf!be^I!L77S*9Ttq!SDWqJ?X0*L*e_EBm677LD~OF+aYna{*z3uUC&8E?T-Yz zP|lHw1AHm1&7JHy_lK&D{RFNo_jI_*kfS!fv@9^p8JysOPFa%nVYT!ESt zK|sKOfB^vm0tN&O2pAABAYeeifPeu30|EvF4De_gm#l-z&ow;x&K=6}wT!J|FXxU2MA zdv6)*=(yit;2Z9P13#A74|r?-I3K6_qkt^A?241?$xxJ{q*E`eMI|(BXXNwT#y^i( zjO+@D8dfrELi`x0&xaDHn)vxHx0m3KH^1E_{Z^Xl^HJCD;QPn2>G!|Bh4xE7i;t}> zf)DJ6v6kxTQmbnU%FfWI#^K&K%YYh}fwJ=rsd49a zS0&PG<9_VA$HlN%?0UvW6K{wUGGP=Og~zNk8nz!Y_2>Wnb^3O&Md(|7KEyeVP}{FH z*L>py;sckA+Zua39RIztVLlrBzR^h@dXV*dZ}eYYo$TktX_*LLY(#xx7Rr{zFEsFg z=c4p9Qyu`6oR4HC`^5l?qZsAonp1#?tWun9g#$m9TMI zXz0Cvx-&_IJYu8(tv_QoP4_1=Vo<(P2pPcrTt0iSN*uT!AY+2RErkL69K;{Fxe8SA z5^}mjteLa_)xP(Hzn4}zj|7&Y?{o0Naev1H6<_bY|6;U(6?ONVgC^NPWh7fZMB5E&-1_@pTH^e~iPa zdMHrhfm41g=_t9}^T+u(m4O0PUO1&?%7*NK`>wE`_*3^W)#rn>1nzzm^Mdsq$D5y@ zBwX}wk@flbS$j4pN|qS$<9wWox(S-R zFHUKh^4IB=d5h4i#`DqGbvHX6JY?{t0yUnGrhfOwjR!Aje8tPJ<6~Hh-~%@f-WI^U zKg-b%eNh|%|G8g_?L*o7wb&b|xc6t(?90OYuo~EZg}+;FW5L;uMeNHWzhfRm?!RJ@ zsTBLP8rYYWg6jXP1@7nYXiUwIx(2wOX-sMywyMkkYFsAD&JCo-@q2E4yA|hD8IPlz zZ$CjGvRhm-W(G5kFx67&WdGWoLv_YnyhW0TD*L-I>1wEg@qde~C2H1aPpnaJb$Gx8lm;Yap z^~YJ_oLk5IG<_XAfycNildq}8^E;9S8tns~sPd?_1WI=v%>#gvvHtldQsd-Nc1#jA zj^A^Co?eY_k$Uy>_?Rl-gD{5L{YY8)M9=Z~Yt@q2E=d|*EuinR5^6l%4_qYCF~)VLIs?Vn1G!`2dm1JrQ7 z&HZ_H)=G-JX>qQ57x-Z|AtFK@5kD>64>xwc)QN{)WO3;Ue@9JG^~{jOpKT&O8OJ#+ zk;I>A!uu^5DBT82bOQc9o}*t}L|7n;77(2f6~v59_KQ!BVSn-gfB!BbEQ-Y>_FP0E zhZcp_0iw`6{~zD4<6Pk%anASbh@URw9H`iub&Cv?{<}qeextP;#<im_W1;2L7$RH2(cgU1Tq${ngUJ-}Zv%8@Sfh3uO6mffQFd`QHO?Dlm#(J9`J?QLHPkqM&u!Swn))1j3lDqaZ@p@m zkLKFxwcbe@M^WbwPH)Ty5-iV608)S>17 zj?}k5Sx%YFzwhHJ(nId~C-&_B_pLNb0c1Z{yYC}u19h=@lv}!q8kdE#D>hT(@=rp?Fz*)Ivzd{@UimMk3mQ||BfG~=cJsG=I~ml`|1 z{)qY`ZqJ?7B6uV9O1ws|ocB^!8;{;{_fg~IQFg(8YMc_vE;&Gr(?i+i#ndTzQK{!c^gCxVKql}&w~qL12<)Q{R4rqu@dAk@(OPVE#~N$%dB52`Mx!Y^b1{T>5p zjKp;#nenXHFlGc!K=z8GobMp9BY0mN&LB)Uo_L=TM`>}<39;;GJToGK{pw;Ii^z`( zp~bTTL!-h*C9vOhWFOexmB=X1zYXbg+8OdxZmZ;3pdjOqVt{);C)TI(lXY(s_d||; z>@Ru|zjMO=t{47s5(0|#QLXN2JY(GJBxzzK6ql1(LdFZY>cO+X|7b@2{+sNXFprHN zZ@Ea5RFdV?(wbEpaQ?Q|0{cq#NqLH#yv8Q z$2(+9@b^z;0PAl)y1`H+D=NtVFF#p+pUnU8kQ%4-lFVnkCMEt`f1JcUA8G3IlQ1+s z;@03yd}+RvFub{4`63R}wNsVKm?BYArIB%3)dY2kk<0o!W<9=6e!+OG$y_(J|VKE??6AX>ckI~ijIeAMA%oPduye2f?HQHPHS0zT^SF;T!r z9X?W81Rq3|IY&G*gp{$&z>yKGsF3(jPCPgwA&DIzf1f>>c)lMO$9^^+A5A>TBhmvC zVq;lR@ibOY5b-!K9zC_E(O5}goKKEMvf?-e;ZX-3-*F0W!XG^Dlt<^r+@LPRA7$sG z`6aFg_k2Jc;l!^}m#2ravu;!4@MzxuIyLTh&EH02RlGuhu&LBVx@D+73*?Tud)k6I~BW{2ou^ zzehO<`R7dR5>Ux$Nb-ju{atGsp2v#6pstr5%9gFB#wnrf56`G^@+iCNDK!rJ&loQK z;64lc#gfO|Id$V$FUp_PO>b(fA-%%;_^@ocUdeSK>xl0Z*Oeun1{73MC{SB5) z*AIQf2Y>sWU#B1Ey9(^(pX)Gp1bh$*QM)envPJMgj;CAkes1-S)bf1P9jHl=!;W$k zr}ZM^xhS65hKz@zcm~>uR)OLLUC4a?>QzwX1fm`TcsyN<1{ft|^U)+^$pEkCR#}#0 zzJxM0F3E??&m2xlTmd{^_D2b67s&udaqk~cLbwX8y0{OKWL(^blvpnrSL))ppXVl# z*)-W!xc+2%jy{YJ*=2$#b3UTMiuGNURB0e+LH-ye!|YJz#hU13xD_14ckj|^)H%;AO7xNxX0lsXguz}gseBt z_eMUDj-Vv4-2Psr$~fPlcTCRQ+(nE<^pMR8F zn}*x7T!>qLYccFl1&(&0#$g!0%7*o-x%`4fi|`A-jt@%#AB0iVjyJ6S5I&mfhv&CQ zKg?^rihEs$7mlGy`p;^#{?oAiHFex(jrirB?_~0;Yv$CeVLlo=&vD|xLl$52fGgMh ze0a5JKFIYq(VE{Wiur(<`mw~LJpcDU&c~_PDn5u%k8eWJ2`OCTVq8v3kCvH9audYXi3iu$LM(wyT><{6ivHf&-i}oL6ztoB*k2-*qkJLRu5;FhN zZrpyL?s)|$o6(yZmxi({x>Mt_PYFs(0Ms+)CT;1CK|FMSsaAWlxswJ_X%B9zkQ<#v^4M=P!x(edr;B$dtyG z@#A_jP;CxcpPyRa?~n}O+KSMcYdk!! zMc5f;9KSCV86Cv(bFijyq`4){(JzJ>OMLk+0)5cW(Jz1*9~erDjSBILW3hs2&JT-Z z`2|ELLJ>F8A%ot#~jV_3yWf>1+fB0(unVb z29Es4wK&&jg@v%9Ml&ObuL+JM!0b1W90|z`j3<7ym%Os0-#=sQ{J>;(WdoUkq3rJw za<21_TL?*x<$R3PFP0e;mJmndiu>VDowC0l7ZpbU&_Txh;)&bjq>$e~^dmlxS}$IE zlLXSnoWvkmlSl9x%OYO3ii;*UNChyXMxtBs3rGkK#@FSfql^DDo%miWQNJ*Ph_HAz z-&lV?;1mBNvSE8}?6_(TkMTlQi(s_*`Ple}^U>6G=p8(GX$j*9QJNBR7->@2bNpDn zt5LEI?7gCM2*AtD!-f>$D`UmJ-K1MYOeX* zC8THEJn)yt`E`6;7Vtr6ckT0yD}M+djh!!*q5CG+QECl7t2W5TAzu4K)zJG{xavtL z&|LdNdF`W-9ZmX|%m8jbLsp<+`)=xd<~EOZ=5as#Z7=wFF80<*eX%C|n0{F_A1hx=jl=W*YIkZ} z47yC|Kx!Nw&6oD0#^L#&gcmjL$9LEC)P}ydl%zq{N4Dvg5S-b98L01fNKL@cw~<-V zuxz?t{eb)swirwrf6Yrszf=3&BN~2?iL1*c^-x^Gl-jP95UyH>8Yho%724D|>=#4T zsd1R?Vohosw*OKsY8y%Iu#*~bHv{|L^UW*<*$xF^ebiSVyc{u^}OJB0gy&i{yTpAqg0!hJ(H8l5bk zMkmV^qLbGbrjzv)p_4Zvig2O`Cx&pX5Uw@Ci6fji!nHxTwg}e_;o2iydxYzNa2*k@ z6T(TLbP0r)pmXTcltkxABAgV$Ng~p=YyOP;|{O+_Q$FrK5)sn1^H`!Nv>On(uz^_LSK^sTx;&{`BflX znhQpt1TImH;uUVxIJ_U`!#ZjlSMe<-+OWPg_Iac_(kr}@gQ>{E0lto4&5B-y-@e?sa8vM`}N*hxiF~qQ)s99K(ScCy#J? z_S87+jf$equekZ``wT z;KZLS9N=qNzklK8qq+92ZACm`rZFvPIKUSYkK@GAz_nxzcYW(tQw?Q@4yDE^q0=h* zQ{(VC@o){a*5J4J1>eYoG66*8uv+S&}Mc7%xd^B}Gour3E!<4hHahIR41&@%Sz?luY3DzDz|rI$<_7E@Ly9 zpP!HLsL~DFX;a@zJ&o)VujpZxuqVWE{*t8s7pIk=YQ;p=<%_$vYA7KIHyk$r#jB9a z8bw_mIa;OBLj9<5;;1}{5!AR`grkq7#^ocNa|AWcAK@5*)VQCgS2tRuUi~~ip8lbH zGcXRktj zpQ-k{Nb;BtF3GE3DV;#tcXi5-?Yv=qZfyIB8jx3NneEs7e6;_=`Dp4qyE6}7GPyLH zzoHH6M^pLeA>e~(>e_Yjel3y@tgXZc|5(Hq`p{SL*#10;737D$wig#2F`7kwFOTzS zK+Z=3|B*-(!v3mWbW{-WVZ3OfAiuao7AuDR5x>Q-D_RT!c?aPVetgykl$;++jVnet|43>a{@z#t`aN5W=AOSXP^1BEuIsNTp!xaG67WGZT&+Jh67W%n4+{Ywb@;IV z!}(~gb#iya6L)`+_AA`vzrIF&f9lOSR|x&iHm-eh$)@iQ_28j*xCj5&w@eAF&j(rl z@Bg072p;!Owio41{U4clpF}1qBP|4L!jB8Z&+*dGeg%x?{?3#v%8)?oni$8u4pNRT zlM#avKfZ1$N=}KU#)+eB@dRodw$e}*HSWi?sh?V3ulUu@SUmKsdF`yB_YN4zB=7R| zsT8FA2N-%wX1b|V^_a6NMWOuk^T!S0jLSty@hz4l15Pjl-k7;_1}5P*i;B3~C(NO2q!7Nz^zxDz9QHHI9LB z)zheP*gxeXy+o(g6TeQc7Pbhz`gMFP67WH2W9|4SQ@}?ZJ{GqKKCnZ<`>0sK!NenQ z;^*{K;sTlM$K326%_EN5bhM)io!HR6tC;9S&ZBZ(`>^zI1^%V{b@siaMYP|q<726S z4?;(4{oS$_!AC>$p|m;VI0TQw{-qfTV1Eb$cs_*3&2)5M^632ue4Ct3fQX=NM>ITd z;{N@U$OQ5X=9n-}0FQEIaY0Dd_gf^Q`M+DU9N8(JLEws{;Q(Jq{2WgY)kc0HX70z; z!}F3Pu6E)^tDyvm1=QtXc1km-akw2T=2GM2(RHi;rpDFyo%z%_?)%~u|8;t`sYU43 zuj6B%fDb|`YsUv?1bo!t`X@o6j{U0XCCxc%ap(-i62qc`q7#X5 zXTRi<5F5uPfrz7A;+~?2`#}>lBE1In&z$;rKjPIDPC51PWDdJloczG#zzCLKTzo7m zKCp(@fLLZ!U}yv@E{kG&_Iyj-LIHtgS;`kq=Z9(JAC=sNz&6Mh~a0|a~!4Ou&X7$V@K z4j(!#k`MfzT^4dg^p)5n{n%>p<6gBD)VKna?TicladEV2gezE0U0w+~KgpRISAcLC z7SuQ%)oF6_&+|JbEmE%<=A*gZld(elaJP3_liTuVoY*iQjU9Je^U#Z|pHcm0FSB7j znoB>NTBv@IE7|yY19t;7{XMAf=-@7q-|a=f2gk!9Hz!0B@FC!%)`si*)Bm=9;Bl=y z`aW}24)^$=mOFOH8R2B}sBv_J6JJk_i$S=QP1HDlgj3o;jmtvSl{cftAy~5sfnL?6 zvZ4Jd#c2Pvc-_+kvHu7I-1}8XFHEDA=GV2`kk~NJ`YPVnh;!*plro8oPfA@4LeKi}f9sNA|urOh`VkhTxXsFFX5nd}IjtAXK3C`C68Mk2-wh3;3wRM{$ebqoH|8=)oU*LF^3$ z^3P^~dtQRaLA2t!FE{HXn%?Qq=H<6zD{Tey)6C!`MkEVgjn`tDx9UlCqrW(Y=V^Fc>&?YKcXjZnrBnisgt!s^?YM%z%A`}_5@7%7OmU2Z1n^?x zV316e-Mt<>^SRGUU8wzb{?r{N-09r$)X{e*+P`W8vLad|p2x`6lz-~t_XirHy`lG1 zN2?IW2PVBQ=}GRL*9{8n`+ZIMDqMr(E{{7tx<)MtB1`u3wvvlm9=G2;dDiw_}=~J^tZDYhK_D0 z?j;2Gfwrh!5~RJuwP$w#V~MqMJTCTy)MHayGuA7@MCWz0dWtDQQ3LJr>hF_VrP51H z=05AN+cM^viJ8vH*|+=}t~ahfxBTKT^nlo8Y4Ew=UEz~B5DM~g%`dK&0c{Jl^7Hon zAgh~-Sd`jJll0{SWpg9cVfv)p%ZabcOp;sQi0i*p8Du(4d}lk=5VUvwoe+Fkj3*!X zy2j#dibGcQgB;W@{$ffSW|>LC!cH9%+>4ZJ&fjRSa(iu;esDBu&8JR-)oRjJGzN&w z$rOi44vK$uJl(S<9kbiQ!kQvW=Py>W;Gw8<`(rU3=Is5uoO#&{w(X>Dm#haDMW37) z1e#@)LuO4@fy1R@V}lNALnvyOSzPV147JM%)Goe-54}zoT{5{ZGCoPM2OXljM$a!1 zR)Xm-FWYLbZ37mSuO^p1l7~0Lwgj&4sRy5)kIHtAp~3aIc1O~8>D9DLTu@ZzMF&l| z8E@b>qjV4mqjvcgGeUaF>%O3Ay7<%RTN)r9kXUwWnGu|CtLSn0gCeN%XqPuV$#zLc z?Q-+<&E?9e+Hmy3n1cgPsW-Y^@Nm~OLwXP z$W)Fz5#HGVu)WM~pk2(Is}5$j?Fx~WcO!1Zn?kFPl2LD74~B0#hjpx@w4to~UV6;O zfpGpp*=VINx=?`HrIM>%rmHNynRiJa&Y*Vr;Pq^PP^%%ZLN4-VsSu%OXQd;rY}5gX z(SiF*cN&9mrDN*LWBvbiyL|Q7IU_Gu8`jvYzO;UQ{dT!b8+AMQlLF87%18DRj@xB^ znS*kYHaKp3RK7J>6~fDv`W}3v3bx%oINj1y2Zh*zk;}UoL&1^bQsrriP{`FT>S>Q= zj~!N3ufLV(w!wSAP&L>%CGVmy$&T4p|oILH|}0S?@|^;NgZY&K>s}f`qK2iE@q}Oj>*Gm}>M8*tOfh z=ba%P3Q)TkYn&XMcGDP63BAgjoY)_n8?cv+$X<*`bu0<(qYwJMBzJpXFos?2-p4LG zssN%_-rFlIZv`uQf9<+tPCrnM46^84Dp{|6aJxjKcDX%l*prL>%-{}vC{5v_I_yC0 zqJr9`ApA(#M4AP(J+P|diN(rr5VcDygWUl|Z{^|PHZl2#d%AG@Yv7<8wwfUI+HXw9 zj1HhJn^n61uvSgGNVyklDLu4@72Uh*Ee|$@^yiO4lBV^6yvx(hL|H4s=J$COYbW)E zoQLs4KAhKq?#bIqmf7n=D~Fub{u;`#QFX)vVGRZxH|%h@W4SH_TXlcQnyCwBUFnlk zQ;a}HDI$Az8XYophqUn?Cj?Xb8W!!^t`7yf&Yl<&DGimVUBpqlsG@c`jM~Mz-|+hq zD|A6-uFVdIr$#mH!po24olg!l`K}Cw>x=JItm+A!rwG~Wk1+tb7~5mthFSovK>Kv% z(*ZT@lCX7Od{;X)n7TSyHOs&LI1~G^ji_Dr8^;9(ytIZdU#nsw_p880)Gm)^u01rq zod)zOoG#RBH+C))-0U*f1;cDkvo9@N?t9y)XwcFDsKZ0X_Kvh9XIdXEPm<)n?_ zvPBQ`$WEpp*?VI9AvVsySZc=f3ATWvV=^+nO&kQrt!>WrFd7VoW3#%sNGL$I8|(1d zL_^r|!C}|4o4Rn>o@^I!?-A-ZPdkBn(0hyA&z9BsJ3AwgS!4!76*j}3Y z!~`~9nN5Fa-WG!PZt?ROWB`NOIV|lepl}2r@r9d)xl)O6CHTyT9Ugy&=ho0yVyQbe7jVI4*B}M zPfzgIhW1%YOVS3L)YR*hu&%qBzX;%VasIq&zzJWAns)IwR+m20Umd2PcG;-$w~qS8 z`t33+JJk2?Ru%Zd-2V6MXmuF$F1A`xUjr@`DHN=DHW)50vdbNNSQntv>e%7$)xq%m zZbb)Q7wD5_nwt}21P-nHiEauqg|5b-_ZVtc5OijL;Y5}>jGuW(g5{|Q=UgsbStn`) zkzH(#*%+z8Sk1H96Lc&C{po~03wItKp{oWySG2RcKt zq=wSC*+alyq<7v$xq+bhIdIITKnv*aIDLCNdqWt`qg}3?ELq2x>H?EK3w4?5zyMXl z(dq6Tj6wNBrlv=^BG@B)xzoxx#d4_%s3y(LUVg+D*4%!hc>bOWxSIrpGeT^@U`5Yp zv0eZt)Qc7`P_+b&1xjO1&$NbUXBFF*lNj*E#_HhQBbsnH?m`=@(ZTmh#SqR5TO5kW)P!XcskS0!mQWD5BVJk05-#TJi+oqph2+G0<-Pl9z`{l4`>qrH z_|ox_-9nBN<_=bbo`OO@5Dt4IGSj-}GW~7wBsFeb}jr;(yfy-bLBWoh#Ro=zG651ns!M=_Ob)n%XH5!eWG`}L(+~#sr@_K zLZN-9*RE5{Ag;8SH7ch+=(nyMa`u%0cn6%zE`P5JMgE52-M?5tp?d2ViaCAZOIQI* zBFF)zj);AIVxJWVTg}Z9a;HO@eGucxRy(jb_eFe1hy!eP`71rI!W6WH`x#bb>4LZV z=C+vwG+-98mn9PmLpGk(2KVF@OlvzMSk}wXZ()E9ygSnWakULGPB6MZg5_Zb(^sk4 z9ezrOi>uQf6m&C(>T})fr8S0tw2Hx-J2s|ZR$^^^VTC5#uW&Dt%kBXK%;PosIGTb& z$%{q8Ii_$vZN#-(a{xkI2*}d|Z)7jw+i#AaLg>u}SDl#G85R&VjX818 za94PvzE5ZBV;dNdpTAtBw>2mWOEL;C7{W5&r$Iq9MNm(=m?O799n>pKj%E&FfV%F2 z(^W&Xp!2lQ)WZvGVA?_B`CZ3Yz{P!^_Yd4khf`AFU8^$P;QJBVU5N`F;KiF6(Z|d5 z;M^j8X7M6pXm?s!F?@>&bW?2a5iybuQfR(Ebb1%7d(TWDH9*W?&)W&A89EVX3oJnN zasTdLp6EfG#n6_P0eZ_DSEAT#-Xx65T4)U2p{b(ai;A*P}M`>*=U{K=wV^h71;mZM;Gr8gB(05(z z=UBPHHT}a?$FrrAB^99Ekk=m0h5B$kRo3!sccT5!k?Lr5%&uc6a6OW+Y{x}eNI>lp z9MDzE)k_6tij68$QZj;b=aYBLS*{3coaV);tulv9tx~bvL7K46%J+nQ8zbnPEIswR z3j<=`EmTZcYy}g$o*(*VlM}R&`LrO|%N5i=DOqTZvVi+35+~Q2^#QYKM>H73xJ$xO zxkq}w2`udH)o$)6E0{EL_QDR``$PY6ujbjdv4Oq9h8+fF*~6asTXRlNwt^R1{%VuA z#ThKL7mOKjg|L@lFE*N-Fb7rrfkpl|jKLn+%RA)93Z9=Wc>8KF4FOhg^{=rfFW#U-NA<2xvs3M1?ssN4n<5*K2%XvPL1!(Pe%Uhb`$q$~ zdGO4XK9fx0=L2y`sfS^@h6%>T{)T^@pN((vS6a9;B=FO{=CoaJ%dmJH5J#k21_XXfU(i z*a47vD7>HD4I9`sa_v0rzWu>t{GKV@2AV^(zE-)_i@qSGV`w$Y(h~CCzuhC8r3V8u zvl3!hgx)M2*k$71W)PrdzIlHy2WTa6W8KAL&NcS3IDFKZeM%OvuyX&tcaPO!>A({^ zUx%xK(#NyGBdxVS#vx&q*Ck6(ocis+n;`>1b^Kpe2iG~li^R&v=#efUGRkI4G%+8F z?td+K+B-KeJe&Q7rsM)&-bESZMVNxULCk7{GlXA#sB13qsz3O(Q*atyVGh2An_RLc z_JOI^#m9~lf9gPV8Inl@J@@?tvF#O_{K4IyW zkYhJe%p%7Ma?BUc9MPNTH@0i_eX_L%;j3 zHw2Z^@2|?T)nV3?t{YuLY$5-fr^D$Bim+e(>aNw_8DKFyMEGzURfy^OX7>hTTZr6# z+iKTqHOR8~(pkQv6|_w(GJS3B2(uR7nm1{qC2Z>>Yf(7J5zZ;~NH{;+9&VV)7_VNz zfJys;x$z>!A*XlSbMYID;b~&w9Mc&N&}om;WcuG)pxq(3-IHN<5MXw! z=l(?=Fyq_eySZcRVPf3=Lji$q(C__o(Y+sB!N;)Oca4?A_#db* zhuwS?)DwD7T6;C*y$a}+Z-0F3uo*}_K5KQa$N&o7x4C&M#0o0fb=>CBSqF+r3`#1E zoj~Ho$(yX<#5^>&XkqDl6KFeBq@Bk~Gf?cAw_pk3Pn|{l^!k-KLGl4}x7RZnur(v) zW0wk7xOo2SWFHYL=>1fF-iH~sKr0PjId*(+Xrn*;+@}SG@X*wvf4+_>JTIEJ_@I>) z+<7E+hu#`Dl2+&x%o4gkYz-Z9R)-Z#5>tjZYnff&1Opr{(7B%fRtYAK?Wh|>=<}LeJH}kvVh{fN zlvitiw1vY*{K_wc(V_hBm^8UccX0CUHNE>DH%JzeEwhf(gl!%dgtTWEKu_m4x}CG_ zK!K>u#fG~S!RnV1gv5lH@ukD*aQvN}Th>bl=wH1OwiVez>#JUI{WM)+VPCnc*)HCY^R=hJ z%P%%C|8Ae5<{f+hy0ngUzG(`DLaSqUo73U&stYn-C2gQwz16~s{+h%(6KO92q9^Z} z%oq%t^~*oZQ(%Dj{T22VJzQbYBNczKW3Djj_|VujZCpU}^}AKWhB|`##0l~v(k$Ub z*OCRsr32w`m-&VlFPMY)o4FcKcG-cL?F&}(? z;I>?@7Myw?Jb!s_bFkX-bg|7uEBLbi{;~0sT)|PdlkVfru8>v4FdjVB4#t{wHl7h@ z134pJm~DtMh5aw)jL#;XPYs&zT{fbC4my|ThDQ_q_sj77G8)|#9t19(D>Bv&R#;s5 zl7GMrq8|<1GPAos426~Hd!|}K`^DW(-f8U!a`&EX*fWFh{{tVF=@?o=W)j`)ZV}P1 zBYVMq?DCVGvy+pJ;L?84K}VGwL3n+!Q@qV^cuc>ldZ4u{lx$sp$kKxjJ}YcaFCyl{ z&lm2E?DX0J#_m*Gu&IL@yxu!){=%8WdXLU@`s2f{(1x|`Nf@!dZkVh1Vbn7^sJ@lh zCzk;5ZJ~kW^c7ApG50;}JjwvDzNql_h&ALK4seW~)DI3@GNaF4ZUs{{ZPpH$HvrC_ zZ>KGCMw%Fp47?e0$qclEw<@jm>w577}^>h%nmiRfK%Z&GkrP%w9f7~uak)>+@zn{Da6u*o%h6UnB=&?^*MK| zF3UMX!bsZ~gGwiupSxzxllyj%7<$`6si!#%mn$mK%&`Rjwe!ZO zn6#%vl+dej6HW$o^?;5My$*eR>?|Ow%yVBnTQ2Uc{*(6z5!G`M~N$L;2t~0czUF#3JmXB9n`=AF0V*&?MSDC{1 zq9?cVwsnV8i3RI6TvdlLHi1zGA{gMRJ*MxHIAa(k#++;~mLVZ~HWUJ4J zF3#}KLr0voRT~V4E{}{~Zwp~zT@J-6(?K%T)U4#3If!oW_$bkqnBRB5x~Mb577poj zuwHP*16t?5Gc2594WfFfWy`LbgY&ALHcymI;Y4-!!Uu+q@G7RC=ilK(yUsI|*B|K# zPm{L4@4dhbVSZ=FThN8o~^ zZ8I%EdHl*z)!T@EI^1}Lfwu*uYiyofdDOEQ$V8M?yw!HjuIX9&CO z9#rtAhZQ7Lbbb5!tQDMsjly#>ZD2{ar5)$JwI|l0D_@ID=?Bv_|MET?V+Lb)rCpM_ z>IN*!)+Qn6tzclg;7eLW{~w?hQhuuo;RhUqI^0dSg+AfwTh=>CK&v~c#h%Nxq0gde zeI`|C!i4+XHg(9-0PUz&eJmEKL&btwk+c=WI^SNaZn>#yaB$w#l+trbFu_+hMp@h# z&K+I0$6%@vnA=t=gxeD9aPNk1nwM-3r{?TAeCQ}Z#{=^WthZZ&)a+YJ6WS5|tlNsK+k#AtjdZ=pE$z6+`I#P{yRE0p7&C`t;o@dmrdOp^NFtNBXd` zGGb)ETb6KRd`8m9YBT5*bNtSp%U&>E+o&RJybXk=_IPOg&;~x)DXp0EkXRR8mUz%9 z&J&_8-p*lE4u$EiU&d*yu!K2x&vYr=$bf|YDNphYyxf@9&jpA0P^w$u6% z^J3^Q$Y@ph`%%O=y~o$Cv)&H^XC=oz55{Q0JI$rDUT@QZg!MC8jmjmSt5&~KO_^*4 z;i3l~F4gM{_D-eK7dhxaWa6>@Zuj)y(8#a0{f%{Da97n63XgQ5t;2NL_9etTy0_jb zT|)+Fw%fu??q~>u7f1S~67z?Py5UEA63dCbmfS6Q``Q>7L*HI@a4>}m7*H#-9ggY%H=^1VjecsG_R`A9Om9!$L!)~ z1AVq)2LFnR8TcSQfN?zG8ng=lw}lRs4UysQp5kA0um_pvi*uQ>H$^a?%b`R;|$s3XL> z_-vc*tNK|%yWOX@WTiVnyKB2YJ+SbC(!HB+4e99&AC8s@eIwS<42KO>YDf6j!!C*p zRT|MxPMiHL-;xf0TgO20baRN?;ViaNf$&drg;wag*g)j;$q%gd5zk**b-0o&DFNMg zcD*C?lF;{<4#g)4|K!={O<{Ly3s6y#x!wPmCisSnR)}rZfcwSKN=vGg!C2UBamWG_ z7&K$y(0S2jV4)ebb=q@H_`cP%+k6XauyQy%Yhx>In7L(anAL1Ym=?T5(pSYCmh@vB z9JkIGqIVsVc)gJh+a6lT_^xL_pMt8A>0`}c;l-S8Lkb-rq|fEbuj%%1?_Ae2j(NjC z%lGU1xeG?X`HLbBGNlf{dTTRIZnh=N3ZF4U+Q%9;l|M9l_<+y{(QjkrH`>AIQn%rQ z*V%%?5cdI5D;*#`VMJLsCecq#uRQ&HfD7223B7Ui+;FhpTmX9MLt)yn8zD|z+@S2( z)Z>r5T;RN6=FYMxEy!jD%|E-=9P$bx`@Vi_0iIVXzlSM0LWkY$io|8?VB&>SswWHx zf4Mn-=0I&1=(x>xLzgsj*z@;YtqbeTVA9+L2ZTSnK+Yb+`EBxu^#sFd<9*}%5$n_s zH1xY#K-NjK;OM_yAQ;(8!TzdELv>|f?#NjKB6NDfhh1M?3@;f$?=Bzb-`MN`W)tp; ze)+5kw{JgK+y9dqeCvDW@j@kCSeCH`X09-SJFoI17KoSv(`r}gSuF-Uyen)yGE)N* zUrjnVXSXBVI~^wHY-|RA{Md=Rr_wVI0#sIe%bg?EZ=b34mK<@?5`0%beBkrk4&GhZ zdVQ6HI~WGil-!E#A$W`Y$g;n@Ao-BPE+HQ;$l3NlWy4utI5Wmcr;C*{XhwE4`pZup zuIr!q7_f0{d?2K+L!24L&kD{(>iDhyAtkPOb%E-|}tpzG#BBL3T#bGh)5@ zQtJ`YkKEv$ybMkCi6vO}d5}=i&la@GuKCWm>kJ$67jCdN04LOP0@zt|{PLFnySX@6fF0ymC~ zshBsJ4plP~Og4#{!xz@XBYoUW;Y;V0Dqf;qaAW^j`Hs1EU=w}bS!SFs1SM>oz1Gkf zx(g}L%>swR?ZUu!RboL<6PmZgY|(u2F}nkS>tYDD?iY-G9A!5+7oPdt7acqfCc4T&`e1kPiM*RClRg-N zQ{slbMAt9ro*}3sd-2H{t)(&893%_7}n#cT z?Z}{9(*E$vd`iad5HpdvSL#zpS)a zju`)aQCB!}hIk&j<<>>!AWIwAG5pB8X=1KW_|?j9@ilu;Ik@(-%6ek`<$a3Er}kEm z=HjQm?3x46V!cj_p3s8iu9Y7(iT-bws*P)9uqSlR>GHaqIKO-T+|C73dNuRE)JqT6 zOh2gwg$EuDnV6?V>;p(jUH+~WjKBG;=`f%{*vGa<-cK|8(K0B@yxk1?D_OhK;xj{nrtNp%vxWcsVw-e6zGC)D=MU>1= zHOQxTTEDE_a9HCXv)XH?8$3Oxcqu~81}<5;q|5jB1)-P4^J2qoAyv~Y@sX?-q!vB> zC~4*jMy-!AQir$$t49x)7w6nzxIvetF0q~=?pR7S12duzqsp!75URDrtccaXm{hp9HQnGyY zY7hA2Jk~7dv>)sopJM%`lP{#tk2o6?=?3?+&ppl>?**6BtEEZ{jp1m{!UuDwxe@cF za--2(JP3bVxxZ^WCt}`bDY+nlc%E`!M`4YW8?oL}I9-Wh2UGGQSC^`J!-m^2w{vIu z!1J7Qwjk^d%RMrrzuvcktfWFCRTWq0x7^lg9^X~NQP3T9Qilpp)SyPBSW;#~n9)p-? zkAKnH-v!RC9DZc!B?i#$#NOKQ+5{e&mlmzrO3dG!=LT-r=~dGnxLtapcJWg`q}ut7 z1-x6cNwn>mQ+e>Sf4yJ2~Z0bkeqo_dvG3YD8e zD{dBf!eyx;qV3a}@aiv}_fN0e!ozdN^7enV1Ie9TGgqm*!`vNlvU$pCn6UtKEA#cm09eW5pj5(mDa^|xSoZm8{qm+y@m za-*{&M0Q`>DRGV?lx>mhJ;2El?kMGr(3E!vu`JaD2C+^%_U4rNpNMhh@cxzNd)k=Q zYmbEPir&K>DS?U6{%5C%=RzsEY18J15bOAHc|EsjT0l;0M)9}9#Pj8V*^|0j1H9Lh znK5}C11|Y(ekQcT1-vdu#I5{5j0+z;8o%4u8s>#+%x2tBgm?2tJd4_C1?%$-+h|!j zfcTK5sxNlhz~*~a?+g!eRLC}+rV+7;j&tY$4+U&T}s3ue<3E1#V`WC3@Jx1(8>spDNrO`i!b%s3XVtxK&|b<{ed9Q5$bIx&^GO#ADEXMV(>K`* zg7+B@)=MGwb*ynbYZ^!FFXB3sIqQnV{8Dz4@G}h!c(DJxYCA{bb3hwHIt>*~$-miN=ziAi3mt1#-p@pkmdP=$)J}_aQ)5u*Q&CxtHf<@0W&xSjC63;hA z&Xm4o4cnKBT6?zgflq@rZ#j6I=%3#FZS$_g1LkEaz4L7A4%v63U(0THhoz$q9-g_? z1HM&cuqKXm1Id{>4(Y>OpfI{%=$JPyz%1MNsl6G|pJ_c?IjPVYo_SAOOVjs;wA}bM z+5_F8=Wnf#wEB{pVM%b8>`*TgD8GJ4}q@)(qb}O{*U~ojz~gR!KduoiC(R zGPy4(pKt#lQ^F-2}KV*v3S!Ls$jG&w5}U(|8g zfCu()-*wD;xZ(*z4mjlK6ncS7X7SVsA;kPbq|?h6Bm5!##MGtBef&Z9c-}zgGmaqE z%|W)^)1lyOtF`F)bvMxOeNw0Jk~Q=n@FgfR+Y8)x%vPvyWxyf5Z>PJl+~C5RfuAQX z9R`jk*Sl+8as-QMvFqJl`M|Ns5x0aoy20pUivmK%66?`VKO`|_ec*)93YVESO~ zEZVT~L8YZ1#CQ1q<;Vd#xVcJ}#?5pAALXfug-cvuUXa}`?cFvo!9-?&o3R5dyWP8R zEb-jV`$?8kX|6lmEu-xo+hquhUvK~IjEw9iC5QycIZP%n$-G1*=bWnC_uJi1TUSWKmX>6akYb6>EQXc+*Uiy9 z7<{qsl@V;4bRDLfjqtAb)}F2vS~wAGm$lKw4pJW{ie|RBVTFCkCuc`*Eb8n0P`k(h z9K&Je*|8SLo7LB4;dO_8mc$LQtCnbMo=t9OFu;J+<>s8d&frP&i9O|I3)_cMc2|f# z=ix36@hKaq^>uO82swfB3F4>g$on7RB(sLo8Y@PVc^M4t5pl?1kA?&Bw_bh8m0V89 z-eGt-lIDOtd8WrY=15=Iw%#Lw$1D-;m%00+oE=h#a*I^$sT$E+!qd!jSk$WOGxivb1*H+a;DAgKEe*ARzAo;C6U80H0Hv`z#X@w7+ zwuWN!*|G;DADh)HHuSDFMX$8Z{Z4TO%$L6q+%-nU_Mk71S6GTc!E`rUUcw5D@^s#o z_14F6s(eQgpDH#D(#3vk(7>~P3vnSBz>dnA<1VHPRo|oL+ll@@bjFH)IfED)eA8x_ zT~yI)dd<9w#O-SwFNjAyu_EXDreHn0G3Gb4Z3^Qhx~R*w!4N@f7=YKc6 zPlX@+dV+4~I*{L=IOro;W{mC6&m@gSyMQHW&?}tL9TKukxLPAQbgA>|G&{tsyqv?$VulH7 z`k=CpCFq;e?YZyTpfpLy-%{NHLLzrq8Oitlsd}lmTEra%87;cKLDnRG-l0{{ZwJ+2 z`A;%JmPl@@kryKM%Hgcpw5UN+9t{_K2?V}^>$~Yz99djoB zv#b^}@-%wb4=P|HUF73A5|@ot$TtdJr{QgvipR-ka=3Od%XAOP!*@2eBtK2i!Rlsh zk-`aG{8H8xeJLc2KJz?-RWCHqccR&$i_|4zgxgL|?68Du?vX87yS2a|Hj$fr(SyWm zFPi)r+@Q3<(%al-D+F^?Ue_IV#H#q8*XbfGfr{8mONbphi&L^=NFDpQRDH5A(E&|N zj#3pD>mc5w>zl|0bKD6$dNk;RBbMDiw|MS{3-%UWDB;()gnR6vC%=^KaAMjrdy>l* zc4FiFTF&}7lb5O~8DNDsn_*D=*%d<7*0VcvHp5hue(z6bCy2>+sy`2MhDy2QN@Eol z+}$1*;_7M#{}hv5W{nnjx#q$CkA)0nc-h&kg8d8obD^SYa(3H_`0t}<}Rl1BF} z(VyMrAo#GU*Us7;C&F(&8zOPs0mhfRjEO#PUY^svvRf2yKdWhcFVMz#{u!Ao_lVwE zQ70n1!wzo+ERJXmd18uhSLZu(Z|vJ_;@@`Jfz+#e65g(JLkUI7lcnAnZl8`C9N%IK zuWTDzPbQ*Uh3cg*ueCuwBD(9?EWyl`nAm$)8!T&F_Rw@3;BB}5(Lk;z)~~btve?lK z3?1365Bte^+2yo^>46ivjW)IaG}VVkRsQ|D-xiSH!kWBy(h+fthkog3Sd%2hXRAun&(bR)?>w|<(? zJszj}H@@m`*(uQ&G^!EfZ2Ek!JL%&(-uv@aiz>L7ZhlDHqltOn?Hw;pD&e!kDOeq!z3dbLVbH>H^!9yJ$L6 zIB0SDMQ4`*HoZ(xuv_8^4=->rne%Mp$kN_(@U`G!{PamuCyWr$ zxzOm2C0@dY!AETnRj5^I|gM7EKbi;sP855z}@fP z&%Ut3{QRad5p_@eVtQrp^}H3-W>&X6rnuwtsrakbjLxV%BCCIS)DSOpr_@-Zobc|G z+?K2@7Fc>!??lHzZ8-G2v3*H&h;Q@M$RTShn2azNh>*I?8mDkAb44!%OJs9DnRN#H zi8C?FNdJ_Jc;M2OIo@DQU+$ZyjY>njq)pkzYOYy{IHyu}FCS&Z3&14z; zW(y&Aul^!qYiK9$6VIlYBEi?q;6k!G;*Dqa$B}&X=U7-+;2R6vY?fLc(CGq^@hQ#g zXWj8#XCzQ}vlT=NzP#VR{Dj3btAAD79J_eMwUqN|@f+@a(B z2^s}9DDGaW|MRsYetoS!Tub6zjZuvhJ6C4}F!}NwBU2>$m`t4L>Klr zy!7dq4a}$}6^^ZrP%=|gsLCXLz;5<}UHk0d!npL-u9M#AOne*DdC?gH8HXY>Wh}AE z;<*uPl08lpPcII?XN|NJNxK1}pB!O~t{fqCLn*pz5j|ZtXkU7EaAe*cp4-l>n;Egi zY|oeX6muu+J(K)Om&8v8K1_3Rogwv&Glw<1NgT&_(?EtU%n(wm;&*7zlf-%O-Q>Wi z0hCwEtg}t$#LgV|cZxrb-X;g(+H+aqN;>?RTavu~G@N1u>ybLFa+jiyBCpap>y`ziY z=8G-{1`hl(d9Q}-Qp>`n_GB(5_oIqTj1D*=#^f1oY;nIsPO3f10!xBfW76+fgOBwU zS8|yH;x`$N4Tk$d=JNo9%)B###(&K55WkyzyeVq^TNhYWs!~5aCUx6D>8op6EHH9f z=X~lFqAMA8tZU7%1f9;@X-{ECaGT9m{#faN2ZtZU4?nPh*wm`E?@zg7n7JqG*k%W` z2)qniPU2&sLbc~(^G;YiHx+!f!WvdG0ZNNg-SI$DF5#xS4cgAE4wHOsgzf=uYvvdm zDE9fcy z$-eYMt?+pt`~Cz^XSB94rd2YZIXBHa&W}{BVuz| zvl94v-uRVBIGF=B+LiIkOOn)=e=vO(mxXhPk7h_773~c%SteOh__%zGv+~jw@O*RG z*+TjxyVpFgE+X~v>azOV$0|tQ{58Y;T*l3iq18V5Zbk)X&+V%($efy-zTX?y^IQ zBlB@XqQhTvv2AWYZwh&4o*LPub{I{VvP(2`fZNXWORI<-CRI1Tc=^Nu1}`^Q-%)c% z_yw0W5|YG^6!zY^a>^dRjvb50yk(2yHr4*6Mef)>FB%bH<%F%H!rGqscF3;T%{oQ% zKt}Xxqm}mz(Cne%FUesAF|8HOrx)3hzWT5^)g}k1opbNC@3BLAZE+MwfEx@ajXMu) zAo_30&1Zu>77!|5ymP)JJ_7E^*w| zLUri$sb^M-2#U_&p8Cd#k_Ed2-97n!oL&uU>OLDL7!lnfG5gZZ7Co@p|JY_h>M(V0 z^dC*BkviYpo8_F$j!2@LQOP6q(gBGdcjKfCvGUY`!DY+0;X%8V*2m?pI5rqNG)n3o zw(oD!)%4n8Q$W&Y4a+TvVr2Q*uV#YW){wKmiac>??&NOsBX-z#h*ejE*$(%LE0?Xi z;fu_z8JDeW+)?w1b>vwl@%tvWyU&g~fy3f0y%KrPd(79=bhvnsI<4u`0UZl)o;+@7 zVCRPBQL)hdt!`-JOnP-;xenYOFmS08eWabkxlsRz1-eJVgP3Amz&@&dlGn%+#cG>$ zI7t3$s#1D+=#({Bj^_RL3pR)JmsvUUGbApYZk`Qcv&G)a+MW4_ZSlP6&vucj!%l@0PM?(mv2@u2*96#Y<`HePya^12-` z0L$F+qzY0WZm4$=zE`RU%geNtFLux%kzRQr@}eNl+;y-n|4hSZ@PJ`6$@|pHbV?NH zOu%`pac)T!(bumpFTAV?uq30jO4iKIFde+y65ku_+&``z~|9bx5RxR7*WQx zXpa}JeU!fAkm!o)H|G{EpMv&X3r@#*4G3sm2DvgoFQ=}0)EwN5d5cv&2Iekc z_-ZrUrsjfG>|9Sj{dPe5F7-1eaU?Ezp_+L#jKsMz%}?oV?D4hMW7`J?J9OG-mvMa~ z@!_krHk1jyyPahW|{Ksu<#pcKE-n z5rNlRo=__jd1RZj@Jn&)V_3v?OWiGHF#EpjJGq1O+uSIly%{0()}i<#f>a62-!>a{ zC-W2!`tmOCiYI+~lkXa2NS~PAq-M~D#irPhC{oj5Wscsm$&a6EoiXH=;`&t69%nKs z=9L0oIFZn=!t>GwFZgO6eNXd5?urW^cXN2cCS~KuG(%ql73WE_2jB$pZ_sG`1q@tU@HlXbe|+| z(DCMiuFJ;)wXalM>oAvge-3O~uDRxdv@q7su$Z?1WUe4` z?UUck>iGF3)Lli8)Sb+EJ$7!F!P4XBEeT8J!tTzNEGV2ewSVK4SiqDI4?JEckk41?hzCy&c0C)l zLWl0C^X12$pgjw8lfUDLV$&UIkv=X^m0k40J<$`Xe7BDZS~(-!i>2Ah${sJLH{E|u zo+nq}!#H!Q6LvQqbh<0=hy*{?xYyjuh|fvOJ#x<$CtDSbgB+>OO|xqqjI2JoAx-sN%Sc9{lO_WWpQ-4r7Jy74hJ@`ANJj+itmaA zFMUW~s#NgB*g7>i47%#6-DuZ?WvU2s|5ZKoh19QW?a{!_h(n$a$XxfexJl_v2aWN7 z^2&eTT2s;|Lia9?ywAf2cXsF#dq}r?O2w1!|83XkKGryQ(9{GYHXpJD^NXUr@AIuN zd8_zD`m`ei%Pc*Og&eT6)JI8))SXUsybr%d`a~XZ#U-(Ru)yAnYS{_OM(_@O_&p=h z9n81m1~_6#-$E^?cRK0+$e(OjImFi%5@=wIi$WLm-sZR-Tm zSMum##El`XWTg?AmuGO7qVvSaR;d|@<3yheb9vJ6$r6!Y&rDCfcZBV9d5UnpE7+1x zhMCElgYkrD=FvV=qj%*bqzO!Ta=_>}; zrZv0W&r)zR=L>)MGBQ7Wz>`b#C=GnhMT^pHl(BMVZ}1jZIXsX0ebxKB5uQi<_&Ktl z)N>;QTvA1W0=DeB#Y7MIagh1Ptc?|PtzJHRWaS9sW8dFq%)3BLb!FVPG-nKkRo;A~ z>j9Rn?hmv`Uw`U}lv_6}9I$oY$E(-FNFB~!t*|o00ZJ?VC3?j?P=C6^r2mH@I-+=a zOz-I<@%UD25dka2C)L;&kl%Ta{j1US(;iqmd+il(mMvU#OHupG0vBt3H9v^=#2e!o z*Tf2U@Q6Mvn6e^$*c*P`2!81a^{>s3Z+Cg%*MQKoce%PqtZ!BNZfS?jNeVbB?t!Q) z$uoiFdQjy>tM05d_FQ-W{4_xuo6kAh`F*qlUjvmRrqKpkw~rgYRwVU8-Ze$_7Ume_ z+>|zZ)eTzvGg5}8>~TEraO9K4j_{KzKeBf=)g?UMVDvrJud1oiD43Sn= zRe>{_Sah{^<|?}hnOEuQdgDR*@3)ugF)Wiu!Wivd;U{x2+r1sz7DM8a+#DvZy-uJG zb$vAaNaByZ(g$9~IHPdId})Tc15zfbehaM9l6z}UUw+CvPfLjO@}=n)7~)Q!SmQPlpVBW9*ru3hg{s-gwYT`m z`O~<$@g1pG<#4o|-RkZPgZ*6L-^P5n};iD0SuEcN2i59LQb^CQ*+q*8Y zS`eMREt|^iN_ZWEc`CK$IQf(L@bDAjR~7Uk^+=uf(S_#K<*nAxW4)#7cTgRS%neuW zwcDb=b*0{VRV#2k{90I<;z;6-w=d@dtWl%Mk>D*y>fQWo#Vvq$LNOnhbwk!7$V{TstvvNf(){JA;*!h{#>mrk@Zd3d4fan9LHYj2#& zbFMHV{@(q{kdFY#*Q|ofA7Avf1g-ehwB<9>$Kc|1i`BsdJ9SqLTV~lns#exhht#bT ztn26d$MkV~@tIX09IOy}Pf*X~hAl2z#_Vr-N_3rP4FmmFMu^$Pm!#6Wmf>A^si?2I`l=E!D$hilqg2b-~ogKJR-opmox~_XI`BG3`IlNBH zLJJYWUvwy67wN( zSY!$1=Ii<%kI9_Nw_n|+2kg;Un#);D;;6a3pAJ8;*8|HE|9l%#w^Nac6S*`?`Y^*z z$8bBkA#cQqS!R_d#L5pHRwsHDUmByuDRLh;m(T9tDzb%*Zc16{z zn&$Hd9bnJRK;QS0^iALA9>|m+eF%Hz$G=y3pml#^%Pw0Nyk;7Szb{4lfGQ;XidJ}n z!*+efd2&BKUB|bc&Cm}sW)F9iU3WtHQC6G8RubRt*Ix3HkJOD7WlW_vlDZJfN`)FL z2QV^r+}{vx33EQL`+L%@;hSS5nPY8)6XFSWigz9Gt0Xt>5ThmVJS%_~)57E%0G^7WlA+J4&4*`k9=6YN$#21wg1W=`12h9eE#8({wF>x&jKHIkI&JsSDU#0 zBOg}a4#vjL38CzNlMk!){}3N`dr-RG0|g=M$odlc`wlZ`|HOw4V*1ovd695G7x=JK zxf^r@+_LMm67WlAcvAXfgv^HSFKk;GD%Fg`EYv%eFKJ0}ZEj^O21mW~U zA^3!*BULWaMYjaNmKR>BhI9Ble<--b2 z?tJW&!-n)YVf`wK2uh36p0(FeaM*urbNEl6 zf4GXl3wBwyV3#iwBd(FdYya|LD{93WV`_LIHfW)FR9Y7EYXobp_Hm-A>+2K4ld|wx zu*-kq!}fE(Robd9gcWyK_s-5qAk-wA-)NB}q!-T1f3ge1BZX-7SSFa>xSpooE{o$` zBOIQMe2|KaTHEX=4bBx&9H|PT_(n|>&r~4KbHOej{>d(dtT}bQG~Lf=X;ltz9ZXGm+`FyWB?JKp6#b9W`Ly?X_>Ies+mL=sy{C@Gw_7EUm3Zprk6b>`S9U7T*vjI} zf?c*BR;%An>Iwh5eAt~l`-*vq&|3dG`@xAfGGHjKD?av#_b-R^zsrYJTd+%`Z)pDw zOVVE$xXrP(N)U4kb~(FXmr#z%suOW|z-buEeQAv5+=&*ZYAR~f%XY?Ko-j;^v#q%q_G*EyOE$77& zsjz%mk@<*uUgLGO_L*NIyR=(luP0=!F%In>*azYXQ2UZR3 zSV_TV5mm-GAwCqSx3{g-6#mPH4SS`XQsF?wDgB|M#kmwXFWBW1Pu#`HFEqqe+HNr1 zB7tpZGmgL6Dub6YiJia2sAwGF9hd(>ILqcO=Ua;wLx3vzqKj|0mzoM#%v33V`LLSzly9hxaQurO_qr6U&k`a1 zmk;ZcYAxXOQJTzQ$<>tplKmIEY+JC)zvIImvAw;)`!E-7-KDJEmm&etJ4@OQouk5U z!7h83x3VlHT#S7;&dS_x6^7#QS<&F~fAAK^|A`N)R9CEP=O_2qE`=#mzRwqn;L3ts zN>`)~uJWM#i(M|9Tg`NcOG2tenb_^=6wDz~%j5XR{RKCF9^5@oX@6&&RTWvZOA zkW$>>buLVl%<4N`9KE{=X$y9db}E-9Lqz|U4;!{{UVO$|wX$lcpzOX}(kUrTxNTZ- zSAWZ45aMN~S5{ z?J{^BUU(4=OJBvpRLp%44!{OXWia1<=bGZ*8aN z1MB0Is}+Q+dzxPIelFo$=^7lVm0m=_Iw9dCsdmC|6qQMRLh3)iKlUl{I>|wK`(X2X z^1EU!bV_P+#h@3PcIeFr;lb*TluK{QhW<3s@2P4JduZEG^u;W|uS??IX2TcTxyH`bp?`xgQl8~$zou&l};o-u0 z(Ue(v+Uug6FlwLk zwuJ)U7TTV#YvtjTwwQj1^dmlLxNXWy;;4V-!?L|$KHnqC3|Fg7tAg%|;)w5e)03k# z$b{S(WiS-{dyoF%!|G{fddkMjVabABE-u(b=a`32Ww`{tFEtq-BmU>&bhNRBgA{Ho zoR{=i-OA(piC^h_>$1W@7Aj@sJRf_N(ez#7G6xqGV!mViu@Q2Jv(_uR>?s3Uo)KJj z3Sj#tUY%%8QSeN<95vmdh?uxnygMDmuruj;Bh^;{5&CZCIm@WHlz4$5^dyPvCR*>) z6W-@77pa#zeS}wRE^TJaNA`cb{WaS{Cyi~UTV_9!=b_h^zi)b27ZlBZMJL))FyJOrB{I%M`eV7DrHaXeV?20F-%AWDzgkzv=gQ!H zl(L6uFcn-MciHTJCkr8eLDP=oG}QZT7t}bVh?@_GX}+3NBpI8FbCz#FHm~xM3-p4R zR2Qjf5tc>iOm4u5lQev8d!Dl)pObKzY~Ed}mIs&P)174j@|d|Oet+jd3as4Lk57y% z;mr7=ts8o1FnepvwvJ5(z4G?c8+u88vv6Lt7Vfe6kAm!MA^f;sZTmVjfZT)m>lL~@ z$$Pn5V_y=Z9NutLYq7^sa6BP1w0e++*{rBv0aSTd8O+v8M++ciz(Y2T+}A>Tw>aJK zqM$1LsrBAM;$QhQgttT~;#R-6{N;po_~tlsXlNS+_vTK-UMBTfHFjZd`8Fz!E}WOU z1{YhK=h*-9Vapcyu<0C(N0;sR*Zikqfe-t{TPG$lp9cr+nAc55Dx!Mx-q52f*wCsh zWT0Lt2cgGrL)VqF!d=SnMnj$g4(E;aq<)jY>jSLg?_bEFvFlUbQ6U+`&T0;^a}Yj9 z^7fS)c8Y(`%jmx7#80=WIB3ptOOfzCQdMGKmS~iI*=?k;lV5k7W7CdD;InhZ>%vfCj~^$phX3 z(7jdEUu407KX&Ck@+Q16BqjMtO{Sbtp~mRl~Ef^v|CI6yOmZ)?R;51sQ!0PI`Y(#VZa=mn*-;ao8~R z3okqI8-{)F7VGgNrbndnhqWYR4EEn?5ah3Z&j%{8DY1d*+~H!T{Ch*R0^2m zHaLx}R)B1FWqj#z0m$^pCoYMjlKv3C4Byo>oM>CjTu>o_tqbSH%Fx`i_Phkndk5Ri zkl**ucKE}GW%oZ5`tX4ew)~o6k|pj=pkt=1H1PTuOha~p^e*KSG19evvJltl;mj$EUer&Axgy+b9S>- z3&I;>KFP%9Nx|>eSqB41-|+nTolS?j*I_ZMp=TrE)AlnhDiHf3g%*pJo=|Ujtl_bp zEJzi@(Ds?CpI?;_v7XXsNv{rPrU`a0G6%k4?o{Rg+2c$M)WgO#0LvX}}Z&D%VoB8p&X>07R6t$RzEqm!;5^6 zP^oTZ!o$h6{PJsTj{w@@9_VU2E1+1k|Kx*lD(1&)1-8ACz%@C$jG!h3Y;=rX+r=dD z*B<|t538VeNb$ZDghNj&oV-fjhaII|)86EMcR9vWK;{pR4ZqyhOecz?Ta|mF$bA{& zx|;3$Wx_vx#JTR^c8R~|CGbn!L!l-r5^lZbKY2zCd(`XrvPgc{!{~T7%}NO?_8-mN zwpSKe<$`-kO9}V!KwjOGOaU0A`LvFc2E5l2b-PV`q_B>wkn!X;0;Dr_r#zS?JPQHA z=F9#{SdyunJU69|ZPn-5V#&FF;j(7cokN<0bGA-!zKja<-zB@|$(|tvotX0usqzpM zN=dz|EQLJm(V!by@R z=bWSreEc|rE(wbP%8g9TgEaWEc$MxG6U2qY z;J1`k5p1PX*~-~2hiUa}3r=4-WboRyYm>ND;+&fFcSa@5oMH4Ct|s?NiA_X$xIB8p z9+Za>eqCmNw8yCfYjIQbap|Q%DolsVpJv~6TEW&mJfT3vt{WAvOlew_rB(OS$R~6F&`DOWWmwm=Nv!G ziGQo|iJB>(pww1-W9us_q;7_Y9K5W6mE%wD);e(GJ*zX<>~VPvuHHzQ;Fm(B@<;uf zPoyECx&A%fZZRClnoMT=sfvAK+e*bvX~=)8eZ%Q04XX~;Sen!+z`<_Eo1r5znEt%( zcbmIBl79t;Ga5_5jeT!?_-0v*xic%RjZ{F#u@W^`68|pKIJPWLM-iKi?6f2OR1hP_ z#XCBo3fkI>HfAC6xZ6_}kp5E{F*~iJl2|8=~9QOR9QU%Rn*6QY+6$0nw*c59N^aB0jKjQNt#A%)31DpL<4u zC6C48&@;sEKM!LMAY8IdJF5>CtMFoo@q zfD+-$=4n`cCi@qL+H{V4z7m7lHj`Hz@n_st zip&N%Z2BIUVK}OWz&s7ju39Qugj_f7;vu?J$qLctJrtADZzi* z&_Mi%8v62B)%EX5!^z9(4Z8yk<2kbPO?AYd=&Co!dQsss!uNUjGz|rO1tE{??hdz8fO_CpdznEBsYWIxfHs?nA-3OMQW0(##OKIPa?Ill0fc=Bocic#X{Z-+1Y^|@LU z{i&&I($9+^Fd|Z1o5VX*8l8A=8ku)<6iGF(C421*yH{Vh&4-Wl*?Tt<9e1*1CNbWU z=;1pRB+imN<=s74n71jSNbA((mbFrZw^|YX!+^vEna2k@1ZcRlSUH@bif{^?!r8w+ zl11J9u#?n1M2Gs~%_$iwkB$lDiGmI_Xhhn_=-6q(UE|doA$?&~Wd#5BTdRQCtHHKP z6;wE#;5quW3+VN*qkp6-hd`F*k;hT;`0CHe_B~SyVhOu93-GHV{O!CmSD!l0H9T&e zyi3JhbDQim3mQ%_xPM98rHQN(Pv#B_q9?{ITD65o1y+YO89g4+@P;$r#88TYj|1su z1?)2LR(1=h^HM=owA=jJ=L%q68o?l+tP1v*c5jtQ92`jv{heMS0qRq;n!XN_CooE< zUr=O2aCzKplN;$%cW&d(IZq?}m~x6I`5uQXuT@yH$zn%x=JTy(GAKNwz2m#CIC7=s zgO*4T-hN7bUuc39_Q{sniCBnY>y~3H5=g%BO?Zm|`_Njn2^TGoFp;!4KgU@ z84}vvEdj4dvpcsM#8J-c`yks-24@3C3s!|H;FcoOx%1l6sHPL1&afu=`k^iI3WTG% zi#EXY*o=b6tNizrNxqw2o43^>k&3Kv&W&0`FZLY^h`;tu1Di6xOr@%8;6t9UnLmRn z+Rcj>=Q$JoXVGJMD;`o$$-H*QyHA1CFI_61c&p;s8Mi&g56Qi*=6b)qOAZqIN|N@H zcqc^PNNy3CJ2HF7MHiL;JYwx)ijkyZ+rhmd_bmxWOJSLi0-4i4=XXErEcrcSbl-P< zu2IEIi+)gmv?K;kqxdd~D{zkFDiz1!(^mJ6dSHYF;xQ3~f)v$EIkk;ABM#j3BRqWJMT*y8!1IND}?m8rgXBxE7M@&*+7*)BvW0>sU`k~up6CnuKVfJU-?#g&vqkm~+ zzbwSA->dQdECCU|kkOZg6ikh>(4^QDVb-$YERQbPtCO2D!g+&+Z+GtL$tKGq=n?-j z_S+!!Krz>}< zAlj@5ix=nbWGX4cX=M26vjG*9sO~9!>#2x$yVdhWBQ&8)XDnf+1UM*}e$C)hL9}`Q z%Fm{%@O-|gVl~NcV=WcEQu5`nd1HF0t_OLJ-;#3kf)!Dh{_^BiQU_YLr_%k2C5g8l zox00)N(pOX<~kFIez5PwXV2@Uvfz<=cKXC;6|xU^*yWml3bt2sRW2K11%F(u{GEPT zTu6F7xzUx}y9?)~bb$|>vcQLBUM|(?`A!0n$4_}3DWu_+o7e~0kAmpZ8vXp7T@-oS zd4^MFrJ()osfK4K(FxslS!wnWoyb0V@AsoL%ne&zKDJB(H+$(aHrpv-YrOME)r(T_ z8{~S-aZ?7@J6Zb1f+!fR43y!Rl0{15r~?CulSQd9=UzQm!qKc0OM6yzRIXBtGVdi^ z6{#2HM9u@F!GjMkB*{LXS5tRBlYMQ!-X*nM*Ha?8K6|-WARl&K2-Xo4m&CKppK=#3 zC3{7}wk;28Qp8vO+7D+Ul`vIVeZMS09bdz)6&*UPjH65KAFq8%1N~d!Eo@)35r6gb zVM*dY){h?-`Zz?+$-$E0B;B>$HWJG#J!W!&g(7E~{cW3Klay2zZg;py%(yN@Vh zO_OWnMoS5Nco&WIqp~nHF?i+@BaT>Z1Bvooawwa7K_QR1;-y)8e37|h0beRLlaFI2R*tE z`aDy_vXT^*8!P2ef1~Wm6B0j3-4AO}awU6wVmYTI>J{+jS_!kgk}L{!y&8&19a(kB zC4QxiDrgk_SaCpE9drfnZV5kQhe~CR(5i4^H{pcHgW61Fo^+#Y!y92_F7RPp-c+tp zzbgz^YwwbRRx#AaIh-lnpp4t&^Cp5-GN`MM3f3n&v(N1_-uHG5O+hdQ{ciA(J#{X! z0XbkDRpyB+@ zoAkqk``Oc~w`}ES5uA46Stcx}j2_9R_g9IIURikcqlD>k_qalTPreu?*Q?>AOLjV^E#!3z!P2|4R179iN=+@F3@tJ#};5O*MSBEg5zv z+zXrdl+BH#&Q%?^@6+{{WS_uAH`eNO61RU|ulSqXJ7G>Nzl{%*@8#2U_#u_N_nl|g z{U+aM?WrZHECNE1RDZYM{5}=E%0GveK2b-b&+uCsk1DF@&m|8NUdA82M86<(KVQ8# z3KJ%-Gq6hj<-6~ z?M@5X3t`ZG>#zq1w|Cp4V{cWEecz6e@4F^go#-@IH4z^3*?g*Yno=GqEkbC0vqJ1F~1Avw8>bKXuoEY0_l%VX z{}(N?AFwX)Q)Q|Q7+tfD9Vhygw9xlykoVa3hWHmr5hbkr@OU%lQ=&gD&qhtf_K+fk{<`H@Dr&WSHXp|hkspIuYqow zlxmcw2J+pK!W7gsv10{Ccc_*ct`2BlzG6ajF~L5It7LDGw?ezKWWO9a$MlLCWG_d9 zkj=%q6{KGF`Gdw+qTjDgPf66WC38;BPGU(P6c9q^y6mJTg`kM8x~L>6;GD1iuMB7W zJs*G0NzH;?k~*u4;)Y0l>`3GRgWYnlo!s9r8YGP2r4JTGDv&*M*Clf0@5(`)_3Ecp zzA7+LHZ&^!E)M>jBz;yAZ%LYURYz8}P9ONVL1)|?OUa+RzCiC$fr$c2*qST?i>G78g(tofezMPeHN(=RR_xq-Zzq$P6 z^-{koSuE*OobstsN3S!3N5clPpN)xqRa=JwZtS}_a*0(L^nn{p6${lN;OfdwjBJidI_74cCX+2M_DWlR450Q3gu#khv&+z;!>#vnn4aAM7Lixx#6&j(rp| z&(-*;DP94tZXMe9=t&;sGTO4ilZt|i0@fD^zjjTXPlk&MH|qL0PZa%>N8mGVKQ
oan8P2HDgT_8ApDoB>P#69n{vu9MeHhl|~&W*^3!|P4Rxb zj4Gn5&i0>6Rm9}OA2NG638&7d;#?}xvmUYp>TM?ewr%O<^ZrJfU|G7uHHpMytu=E$ zYvxoC;mp7i!mWf|r%Ks2h7etS$yC~jdL7ixhILtAR|98T;)aS0S$v2J?Fy{{F8(<8 zHa>yWjZ0pXz6>G&>GM~O-fd995#5f$3nU(PS(88MQ>qERXJ&H+xrAfrGy8bcFFA<6 zD?GpLfE+3cJSpK{lp%FNM&GGX1p#)cVk|_*C|fr6_Spy(Ci0budpBy5J=v1&&nozF zDQD4VA4P?KIUm2=5*?2nq{7ij+VGQ<1j_W6U9KnpuJ0)S>as1Q&*H?A=d~B*FfYIA zVH4Sh&i?gcmrVXT+-w@qG3Hi5_lXxQi%Hy&AN<4%ucT3TS@uzmE(Kc*)TAd#UVqRc zBdv+(E;Lz|kT9wm^4_`7Z#b(8hwl{8a5fc`j|IC5U)6%-_6z4_PHZN5=B||XHg$NI zG@gA<-qYsNV-=*B?Y?i!%`2S-b6b^Gz)@ zEfOhzZ77dg?Zf-&ThwrAX@{k{v>LfbtJ3;*o>v6dg$CO31Azg`ur)nzBsd_Z{ciALwi-^J1opUNcp==hM$P)V#N zxNJ5x4GBoWgGJonPp;UrxKRSpf+qX{p`?$% zFmrvX42dfaOG?yi5kvU+i6bZcq_A4;$-}OAvS0Fgm};7?BEp;_G9R-lLDKd1wP97F zORZlX=x~Pux6jLCE~(1{$;?;I&=k=*k*X1>3K)sRD>cRe--lhM2kvO0f7zt}ie9o` zPeWjunyih|l*hN8U#7wPj2&me7fm!&Y-uueR|DI}F6tp`3PJ^9UzRGWgIil@_=v4M z{3atB6-(4`v}eHoQnMEPLTs3_Nj+rto7MSbvZr;`CEEIFdvy%1=ksn3r^1&b&f6&m zq`tIkyy&hb>SQjd?0!$`gY~=?5!-=JzV^AXk|d5h?MYFK2a2q_`6I?jTqbmE2$w0C zTkE}2*;@$}9JKU&D{WNjze=AOk zyrkbE>X4R}zy>8GwVAE0UM!C@;ZH-%uE>D?e0}tDD=Ez0uk2epApjAvKEo^a%D6il zm0f13fMkt`wo9%G(41MB8tkKh{qp=fUYwVK^aQ0kcC9Ais~ojDyi^+*+r_*WE!9KC z((3}1*U5gjD{CIT%2z^^dfC>+EZVquFQy>dSOrJbP#^8Pa5E4yAdJ3TcobQyNJmr4o%QB+)<`s5HAkN8dJ^?@@gP=ofkK+vBPhPUL-_ z+`pFZB^l<`5I4Ngra_1pTd|ZB#^L+#f7ypPW%`Jj`_Wo}RQ;wDlVp9UFL6FV!{?lI zc8~qB8N82HMn*&0Bsd=7@kW_Q27WseJIMel`d1P!hoDc^*R*t*o1e)bL2tJaW_jCB>VRJ2688i;4R(wU2UpW^oSiAmKRtjrW5=MSHz(?og)`@(1W6#k zqq^<^U|nFh3R$5R^Dc>sDQitN=rbs=8j4ndhHH0hqFyS%uG`s@O?LBO*K{PQJ&*wF zVy9();85|Wee+)?A@2R8aA)RI9q2i!EM~n&8*&_j65R3-e;1k*eC#3)j9rrLPtU0X zgDh8m3H_HAC)ujnj}f3^Tuiodk1mXg@Y!17d}Og&T{SV11gYK5zK>83L9hCqUG*Vj zxTKcgLz1Ne=~O$}@SGmJ@D=1$LtR@a(J}J#rU0S&f;0!48jSnIY1p`sptI@Q(UJR9 zXqrh2)yB9W$Jy6%Gv>jZO7F5Lstn-F+tF{vw52wZ2C*|(Sm)MxqG326kF%+Xt&CyjpXN`|9)DP)kj^yH)^7Y&Gh9!gVhxxXl(9aBzE<0}FfIx#L7m@*hn75H=tCZPbw;IB7sP?= zZa5+ZALS=dh({ABP4-{)wo+lQxSD1J;^czThdmOSSAtuw`&>#X>YYBl61-z69d*ED z2bSs@1L3Tjc@Qq|=XHC?^5c1~vVCo>eE!On=)f#nkei5kPkL8& z0(%`1;}GS;k9leEvV_B-PJ<4UmAWVW`A9$}OD#qjF$nQ3{Wvs4fimwi!W(YtfVNx{ z|I>4b8=B>Fj&U;J{wYDt>r7q9B#InTK|g_yib`CH&HxpO-Qm8nba;0lSzRTT3e@ic z{Dn>=u+5?u?XknUY~a?mDvAN(DLsQ9H4R|nzF3Zz{W@^&MR@%A<77BtdTT><7|Iop zjXtGCQ^D>~Ax(7&#$8LDx}V$P{g_AO80E!r@fUsA8Nui43y_CrudiYkM*p{j~ms_$D)7R7K34F%mayav<7FK*d#wJNIPG_M18^SK9^Xj;5*J~#E1H~ zhA+-U02}I5UUuDT^&b1VPEl@%X%cBev6N68Tprgr{&xH1YNJU)XtPlbqh-^5Nhn8 zhb+>7{K!^VguD=rjBJB*RaiHXQQXpl`Jz{(ozbC6RUy)mN#u7tKd4q{i1Uc`0{?4`IPOcZ1e#?G{o9?^ zxUvD*7{zJt--PX}yhDEB5;_dOK0v9;p~K@9b{a7_kDha`J}S#{XRGCL_r2jeQY#kqV$SO`fvfwq|c?--don2Pn;XNidtR5Bb- z4V3mR)r7bOX>+`CO!(aCT)hW%cL(9RZc(9~9AeU?iI@P<1Ro>#?TNu~Acz^)s zo))=@+%Z6(h#LN=dHOK5abMImbAWd)4gEK@2;klp>%Llz2Ie7Gn=XbL!ZV|sL3zXt z-Hv(eNZ$*v-nXP};X@(_rlu&eW1jJ9Yt3kDI_9G*Y}hOu8Iacd;<CC+oQ80WFO$E7Ip zr~}o-J=D9&5K^y$sU7y`M3F4g!8rpEc=*-U0d-O;whfPQx1bL2{?CmJc{=R8zU+WS zKc1J}k*m8-8EjNrMDthCU|6pBAnP3RehA`{FQ3vNN_Fy~>vjU%id}dYRlLQwO=LYh6so*xx>=|z>;(Fy{U$ZZep=?k5qupjSSY$nN zpbC9dy`o8{FJ_@XY?Q}!*NYU$dRz(&eXL_0|E9KEoB~1uqOZBb5QkQ>t6b7dLZ3e4 ziD+C;F+3YC=~tx(k#b7La;}I&ZdA&1p)tW`UD&2XS%7QjoO~RlXrK@taH&H>7xM28 z7QF~nf%>JEI-IYxU`3qgx%KKo|4kqEMyEq6ijRShe*2^Y@}Dxqv@OQapF6TXIc1~~ zb!oi0?A<4HU;$TNwAUv!;1`uln($Hq-;giMS=2S5(0$dRntl>&^xAYiP!9V8|9sP4 z8w#9~UNsVQ8h|@|*Tct_#?bQq(ZnXyC0`H_kuQh1XR6${s%?lXR{957@OINcwR_jD zQWuQV#0-w#`h@e)K>JbFrLcMnH zb_f8}<_|qo#JufsLON&NG!1Gq6Q3;Igt+o$-nVJiObDK8I(Gx(sJ-psOLCC^QqB3y zGwVDR-Y|4;hh0WrHIAB&gUExdnof-rQ9$3T&#Xy-sFU1zD022yD*C`R!*bPK$Sc;` z6(W0<2yO-!r=xPnAZtBpZ#}027bhPtGPEPZ;mbRE%u%PWlq|R5!5HEJtjmrceoTS1 z4)UFcmt?_+wqbU3wKPmrZWf*|Acs0sc`u{}WFbw?`+1a~Ds=Q_&R(06f{@P+A8thI zLY&c!6!E~7z<+P8=bb5a$ZQsV7d%e}l15|ZRUtp}Li-QlIocSTe|;(P+O#^_?&~SQ8(w{>E8NIC4e#`( zWMX~aX>{M?+o<23)KL4_K>+y59x{wisl%K|)6QHg|X+iw3*;HD-0YTO$?o;RKz=qV@WG(c8>ywj>eLI15z&oNlBRc4CCw!A$(+CsR z^Vk=2S0gXa<&>zN61H#G(Isi=1SoOyIqQ0y1k!;ln3fW zFDuAgA*}~(_w&V@C8+SUf1z&I1oxl!HIucA*`2pOsqH*#{VqW%KMAc`X5p2Dq=icJM1cQ~PT_(dA zUww?2_=Y^r_TuvJR>WDPr1plGpTzo${(TdwD+Sg&_HYjRAbvilIrbXoq3@Oh7v_Sg zAa!59^<9A>oDk9+oNPt@ko7(@$*&l9TJN`Y;$uMACyt$Yh{wMd=bik9^`Xc1@dxD+ z4MFD7=SG8kGHfd=Qo9jBfa8a|g_{w_sfpfr#S7OJ#xLKSZhI8{Y#qH{ZKy(A@nRu+ z{9BX}e=w07jpNoUd*x&0X#*&2UA~FmjEd`QQU>fTDR6*4%D+7X^Pfu#vS$1ZAbsM> zFf$u{0{5o)5>RjNw?3?qBVGPM?*XExh0i>0Q%chW3XPkxsTNy@v7p4tNnqG^dYYH z5X%teJ=vynQ#FWxM2M>g*PGKpBU5Fdyo3tef;l}C!%S%1=@e3wOo4cveO5c0nOMi_ za^2#`1jEsLcY_-A!P{o&$cgQ`V5+D7=BO|U(zkqg6|)BOy=ndHy*;P{c|h(3y^{$y z+jAAKrqJNb*ltl%%p*=4i%A@kG=!UyrJ^&rR9I-Wc0hSC8Db;PsZXK)+ML+a{+Tj5 zn0^tjb01_vfZDvqd$m-EIurJo+Mx?V1-CqP(AQUec%Wc%iVElBrY|QJkznG0=5|kZ z8a$Auc5DE|Tk0MMsUzq=E{f9|`)9cS|xbEwhKI{=TCwH=!D!iQkC7Sgr^0_%p;xcZD z!njjoY%B6W3{D)8xY()*8y*rf$`)z?t4;9bH~Dz{mYOe8b$W2GCY2s&D-X<`xNMEj zY9P<)pDekl#E6!|I1 z&trfcf5Y`4+<(cq`LyQ^+HBSb(dr>{N4)OOE6aCH zBoV-#Cqu1P2>n2eisNit$$#p@S`K!P^F0xO&`6)Cd$tx9QhV+4^WEy){Mjv#q$d#TGrh){|r%m%wADgr3o<{d42CSP+kM+fT zT|seh?~Pc1rM&^6yVWRQt$BIp_ro+8k>OY@+)0J??`2G4FmAA0bGRiy3Vj$oYCKkn zkibJ)!ZANh7v!(sKIn{jIKiW8MA`)L#<|rlSI|%Vy~*hAFXl9`J7rPH%*A}MxZ>uP zbUIYpl{(~_(vbHkJtBJ57~JJY1LFkrLHg9TuU~MSNSvpg3-bb~AqD0O7}7yx9ov~` ziXK=MNt^JfP@pP0_=tdtE?7Ug)c5j>HXQhtc%>KPf_M^}yk;*6iY-=2spB|3zCP1- zi~$g6+a1l{i#X~oC+b+0KCt9j6!K{xZ}yiy>>Ay(VLpdMVB^ixdjtg#=a@<@cV2;d zu)NMi7b{l)Uu~KGJ`qLuu;z10leRkSJ1xM)g6qQ)W}J)DUDcpKH?qGL?@LecOZPsk zyKH-18MvEI4@j;b`qv!SfTID68@8do?EP0g;vz9*#5d_yX^4Xiza1bMB5owEUYh_$ zbSU@upzer%!n!i)VjGSc!1U}V|2}UTaE4X{Jg>#LM8EHr5(OYX#r{Y&=Iz35VWNXI zsQbJj-`G9?*MVpSuhYiyRb-uTNLz&nVp+-KmKV_fhf`2N|;6u92!;<4kD0bF5w zo&6z*0P4P9b02-ccyG~gqNy+imf3u?Z9si#4r=(woDK5ce(A#wk2(guTfP$NZHHEN zX$e5(FMZhV*BS4J{E>hCaM?DuXsn<4Y3^O-A_qF$*RICOs=(qpZKIjH8X&i+o;#vU z8JI*)v+bCdwk0N~eL$ZVaz+0~e_W53QNVF>>nRnOI6v3H@`VT<{PC14I6t(Uw%9`t zr-8+E%Hb!|7^gR{cUmY$huAJA2lE;QM17yBoi{gt7l*_Tm+Ya!K6xjnBe-7Fu$yUErGSt6a218Ma9*2`WH;ia;-g)VqQNNV;~cH77+2v?%u`#}UVI<8^ywtzZOi zoa^1zA`bUbvOHIT!2nhV(|M{C0{ELk#d_pjXk>lZE#O9j@H^$bPRKLS*=`x}5&hP? z%u)*akvH(FzhAS0hYCq=hog(M=*WX<*0Axw`0?!Pk1^QZDyQ1M^bZkX#(uVM+c8`R z8mN5X8Ro|(s~*s27vtaU;^!RUMto7B?S{b+1D+1ktu}WuAmbKSlj~O|97rbf_~;X1 zI)HKE_)QuV7~Za<;5=0`lj5Sd&j3O;jj{^1>H+g*>%aoU2h&QnB#rJRg30bDq8TZ; zE@F>kwoo?i&-B!9&cQgM#B;E8&;S_k7H+Qc)qzV^!Ht6P6foi|I9xS|^Xo5tSRa<* zbE%gUpf>B~9tUNSKlNej&G%Q!O)Y^j>QjB%1w{zkDzqa+2KgrZ^tqr^74Vr^{i1{& z{dNXFm%qqBKAt@1Kv=Cdw21;QSG6p#EoOHwu~x@*Do+ny#PuytzLruGGYN3?_4l=g zz5s0$;qhixBOo+x5isr1fk%>_{m*$wuted!6=yC1Om)mA>dfe%>+;q$q?Q3@5{gDN z10s-BV@|j%rvP8*S8_Y%e=}zMclYsOzHJzsKuSWK{iLIiZWh4%@X)f?G;I)Z9NF?} zs}3w$?mIuikqLY1M?>qDFhF>DJO_&p#%0eg^7ddJ5*-&KRZAvAc~i{%c8W2~2um=g zkO%hibt_x^5;}bEUMIO*h5~1KTiv;kkIBDr&ZrFQGM7}xx1SIq!S`LypRCQNz~&XA zA>YpEfNhc3cJmGzY%eM?tFXtsq1$$8S{n_{#`z!I_DBT;?o=mCVf{vQkIUq$5`EN% z{qc{5Mfew?fuDu%kKcs1{#0)zc>a7{zs|^r1&^Hj%g=8f#4WF%$Kj;%&(HtE7BDrU zulaG=_P?&(^w+EXum1CElO@@Q#J^g(Dcq_=$C&c?-n=F1I!Tuq^ zEG(`rzT2IATzp-&diZ$w2LI7lf1L2Mm;UIH^MCI9(TVYg^;b`3W&OG1*X}>hW5@S@ z?BkaLmY;o<#-hW+qBF<-$LaWKTz|aJ#>V!?sa(H4kDc|M@BiGx5Bxmw|H%K=R;vwv zv#N4E9q)Gh`{RDMH<$OnKThh`=L!4#`{RB;CyM{lKmPjt{HLyuEyD5U^al8W11$ zV>_E*BA)YSy!^Mn&VMf8pBDJ11^#J)e_G(57Wk(H{%nE&TwKTXw}|WhbM3+Rw`h<5 E0yh&5>Hq)$ literal 0 HcmV?d00001 diff --git a/test/regression/source_file/input.py b/test/regression/source_file/input.py new file mode 100644 index 00000000..34fc0cd9 --- /dev/null +++ b/test/regression/source_file/input.py @@ -0,0 +1,58 @@ +import numpy as np + +import mcdc + +# ============================================================================= +# Set source particles and create the source particle file +# ============================================================================= + +rng = np.random.default_rng(seed=7) + +N = 500 +bank = mcdc.make_particle_bank(N) + +for i in range(N): + particle = bank[i] + particle['x'] = rng.random() * 5.0 + particle['y'] = rng.random() + particle['z'] = rng.random() + particle['t'] = rng.random() * 5.0 + particle['ux'] = 1.0 # All going right + particle['uy'] = 0.0 + particle['uz'] = 0.0 + particle['E'] = 1E6 # Arbitrary + +mcdc.save_particle_bank(bank, 'source_particles') + + +# ============================================================================= +# Set model +# ============================================================================= + +# Set materials +m = mcdc.material(capture=np.array([0.1]), scatter=np.array([[0.9]])) + +# Set surfaces +s1 = mcdc.surface("plane-x", x=0.0, bc="vacuum") +s2 = mcdc.surface("plane-x", x=5.0, bc="vacuum") + +# Set cells +mcdc.cell(+s1 & -s2, m) + + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= + +# Tally +mcdc.tally( + scores=["flux"], + x=np.linspace(0.0, 5.0, 51), + t=np.linspace(0.0, 5.0, 51), +) + +# Setting +mcdc.setting(source_file='source_particles.h5') + +# Run +mcdc.run() From e591b2b7d29df7a62cff28225c96c8de56c34ae5 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 26 Jul 2024 13:50:42 +0700 Subject: [PATCH 11/25] mpi fix --- test/regression/slab_ce/input.py | 67 +++++++++++++++++--------------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/test/regression/slab_ce/input.py b/test/regression/slab_ce/input.py index 8f1208d1..2220758a 100644 --- a/test/regression/slab_ce/input.py +++ b/test/regression/slab_ce/input.py @@ -1,5 +1,6 @@ import numpy as np import os, h5py +from mpi4py import MPI import mcdc @@ -7,38 +8,40 @@ # Set the XS library directory os.environ["MCDC_XSLIB"] = os.getcwd() -# Create the dummy nuclide -with h5py.File("dummy_nuclide.h5", "w") as f: - f["A"] = 1.0 - - f["E_xs"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) - f["capture"] = np.array([0.01344, 0.01344, 0.00384, 0.00384]) - f["fission"] = np.array([0.06912, 0.06912, 0.00619, 0.00619]) - f["scatter"] = np.array([0.26304, 0.26304, 0.15024, 0.15024]) - - f["E_nu_p"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) - f["nu_p"] = np.array([2.5, 2.5, 2.7, 2.7]) - - f["E_chi_p"] = np.array([0.0, 1e5, 2e7]) - f["chi_p"] = np.array([0.0, 0.0, 1.0]) - - f["decay_rate"] = np.zeros(6) - - f["E_nu_d"] = np.array([0.0, 2e7]) - f["nu_d"] = np.zeros((6, 2)) - - f["E_chi_d1"] = np.zeros(0) - f["E_chi_d2"] = np.zeros(0) - f["E_chi_d3"] = np.zeros(0) - f["E_chi_d4"] = np.zeros(0) - f["E_chi_d5"] = np.zeros(0) - f["E_chi_d6"] = np.zeros(0) - f["chi_d1"] = np.zeros(0) - f["chi_d2"] = np.zeros(0) - f["chi_d3"] = np.zeros(0) - f["chi_d4"] = np.zeros(0) - f["chi_d5"] = np.zeros(0) - f["chi_d6"] = np.zeros(0) +# Create the dummy nuclide (only master, rank 0 +if MPI.COMM_WORLD.Get_rank() == 0: + with h5py.File("dummy_nuclide.h5", "w") as f: + f["A"] = 1.0 + + f["E_xs"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) + f["capture"] = np.array([0.01344, 0.01344, 0.00384, 0.00384]) + f["fission"] = np.array([0.06912, 0.06912, 0.00619, 0.00619]) + f["scatter"] = np.array([0.26304, 0.26304, 0.15024, 0.15024]) + + f["E_nu_p"] = np.array([0.0, 1.0 - 1e-6, 1.0 + 1e-6, 2e7]) + f["nu_p"] = np.array([2.5, 2.5, 2.7, 2.7]) + + f["E_chi_p"] = np.array([0.0, 1e5, 2e7]) + f["chi_p"] = np.array([0.0, 0.0, 1.0]) + + f["decay_rate"] = np.zeros(6) + + f["E_nu_d"] = np.array([0.0, 2e7]) + f["nu_d"] = np.zeros((6, 2)) + + f["E_chi_d1"] = np.zeros(0) + f["E_chi_d2"] = np.zeros(0) + f["E_chi_d3"] = np.zeros(0) + f["E_chi_d4"] = np.zeros(0) + f["E_chi_d5"] = np.zeros(0) + f["E_chi_d6"] = np.zeros(0) + f["chi_d1"] = np.zeros(0) + f["chi_d2"] = np.zeros(0) + f["chi_d3"] = np.zeros(0) + f["chi_d4"] = np.zeros(0) + f["chi_d5"] = np.zeros(0) + f["chi_d6"] = np.zeros(0) +MPI.COMM_WORLD.Barrier() # Create the material dummy_material = mcdc.material( From fe4e73090243afdf848f8b55b8e414bf804d726c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 26 Jul 2024 13:55:37 +0700 Subject: [PATCH 12/25] back in black --- test/regression/source_file/input.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/regression/source_file/input.py b/test/regression/source_file/input.py index 34fc0cd9..77fe6743 100644 --- a/test/regression/source_file/input.py +++ b/test/regression/source_file/input.py @@ -13,16 +13,16 @@ for i in range(N): particle = bank[i] - particle['x'] = rng.random() * 5.0 - particle['y'] = rng.random() - particle['z'] = rng.random() - particle['t'] = rng.random() * 5.0 - particle['ux'] = 1.0 # All going right - particle['uy'] = 0.0 - particle['uz'] = 0.0 - particle['E'] = 1E6 # Arbitrary + particle["x"] = rng.random() * 5.0 + particle["y"] = rng.random() + particle["z"] = rng.random() + particle["t"] = rng.random() * 5.0 + particle["ux"] = 1.0 # All going right + particle["uy"] = 0.0 + particle["uz"] = 0.0 + particle["E"] = 1e6 # Arbitrary -mcdc.save_particle_bank(bank, 'source_particles') +mcdc.save_particle_bank(bank, "source_particles") # ============================================================================= @@ -52,7 +52,7 @@ ) # Setting -mcdc.setting(source_file='source_particles.h5') +mcdc.setting(source_file="source_particles.h5") # Run mcdc.run() From 47615f34ccdde2ba33536017b5d9cb27717c71d6 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 26 Jul 2024 13:57:31 +0700 Subject: [PATCH 13/25] back in black --- mcdc/input_.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/mcdc/input_.py b/mcdc/input_.py index cdf3207d..f0a9306f 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -1960,19 +1960,19 @@ def check_requirement(label, kw, required): def make_particle_bank(size): struct = [ - ("x", np.float64), - ("y", np.float64), - ("z", np.float64), - ("t", np.float64), - ("ux", np.float64), - ("uy", np.float64), - ("uz", np.float64), - ("g", np.uint64), - ("E", np.float64), - ("w", np.float64), - ("sensitivity_ID", np.int64), - ("rng_seed", np.uint64), - ] + ("x", np.float64), + ("y", np.float64), + ("z", np.float64), + ("t", np.float64), + ("ux", np.float64), + ("uy", np.float64), + ("uz", np.float64), + ("g", np.uint64), + ("E", np.float64), + ("w", np.float64), + ("sensitivity_ID", np.int64), + ("rng_seed", np.uint64), + ] iqmc_struct = [("w", np.float64, (1,))] struct += [("iqmc", iqmc_struct)] @@ -1980,15 +1980,15 @@ def make_particle_bank(size): # Set default values for i in range(size): - bank[i]['ux'] = 1.0 - bank[i]['w'] = 1.0 - bank[i]['rng_seed'] = 1 + bank[i]["ux"] = 1.0 + bank[i]["w"] = 1.0 + bank[i]["rng_seed"] = 1 return bank def save_particle_bank(bank, name): - with h5py.File(name+".h5", "w") as f: + with h5py.File(name + ".h5", "w") as f: f.create_dataset("particles", data=bank[:]) f.create_dataset("particles_size", data=len(bank[:])) From ceda27b18a30876d7866da480ca1edc859b45ba2 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 26 Jul 2024 16:44:00 +0700 Subject: [PATCH 14/25] put iqmc tests back in --- .github/workflows/mpi_numba_reg.yml | 2 +- .github/workflows/mpi_reg.yml | 2 +- .github/workflows/numba_gpu.yml | 2 +- .github/workflows/numba_reg.yml | 2 +- .github/workflows/python_reg.yml | 2 +- .github/workflows/verification_man_mpi_numba.yml | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/mpi_numba_reg.yml b/.github/workflows/mpi_numba_reg.yml index d65daacb..fe291d34 100644 --- a/.github/workflows/mpi_numba_reg.yml +++ b/.github/workflows/mpi_numba_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Numba and MPI run: | cd test/regression - python run.py --mode=numba --mpiexec=4 --skip=iqmc* + python run.py --mode=numba --mpiexec=4 diff --git a/.github/workflows/mpi_reg.yml b/.github/workflows/mpi_reg.yml index 66ec4f85..792f40e1 100644 --- a/.github/workflows/mpi_reg.yml +++ b/.github/workflows/mpi_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - MPI run: | cd test/regression - python run.py --mpiexec=4 --skip=iqmc* + python run.py --mpiexec=4 diff --git a/.github/workflows/numba_gpu.yml b/.github/workflows/numba_gpu.yml index bc50312e..ba0f427c 100644 --- a/.github/workflows/numba_gpu.yml +++ b/.github/workflows/numba_gpu.yml @@ -34,4 +34,4 @@ jobs: - name: Regression Test - Numba run: | cd test/regression - python run.py --mode=numba --skip=iqmc* + python run.py --mode=numba diff --git a/.github/workflows/numba_reg.yml b/.github/workflows/numba_reg.yml index 9d5124be..493c0c9d 100644 --- a/.github/workflows/numba_reg.yml +++ b/.github/workflows/numba_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Numba run: | cd test/regression - python run.py --mode=numba --skip=iqmc* + python run.py --mode=numba diff --git a/.github/workflows/python_reg.yml b/.github/workflows/python_reg.yml index 04718ea7..5c811c4a 100644 --- a/.github/workflows/python_reg.yml +++ b/.github/workflows/python_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Python run: | cd test/regression - python run.py --skip=iqmc* + python run.py diff --git a/.github/workflows/verification_man_mpi_numba.yml b/.github/workflows/verification_man_mpi_numba.yml index e3b5864a..9508d38f 100644 --- a/.github/workflows/verification_man_mpi_numba.yml +++ b/.github/workflows/verification_man_mpi_numba.yml @@ -32,4 +32,4 @@ jobs: - name: Verification Tests - Numba and MPI run: | cd test/verification/analytic - python run.py --mode=numba --mpiexec=2 --skip=iqmc* + python run.py --mode=numba --mpiexec=2 From e3dd034177a62d3c5930e41919de18b95c89f9ea Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sat, 27 Jul 2024 04:46:31 +0700 Subject: [PATCH 15/25] fix iqmc BC bug --- mcdc/iqmc/iqmc_loop.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcdc/iqmc/iqmc_loop.py b/mcdc/iqmc/iqmc_loop.py index 353f66a7..4ac23e6b 100644 --- a/mcdc/iqmc/iqmc_loop.py +++ b/mcdc/iqmc/iqmc_loop.py @@ -79,8 +79,8 @@ def iqmc_step_particle(P, prog): kernel.surface_crossing(P, prog) if event & EVENT_DOMAIN: if not ( - mcdc["surfaces"][P["surface_ID"]]["reflective"] - or mcdc["surfaces"][P["surface_ID"]]["vacuum"] + mcdc["surfaces"][P["surface_ID"]]["BC"] == BC_REFLECTIVE + or mcdc["surfaces"][P["surface_ID"]]["BC"] == BC_VACUUM ): kernel.domain_crossing(P, mcdc) From 4da65fe3a349a1aa28cb4d92d20bffd687ea1ba4 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sat, 27 Jul 2024 06:02:33 +0700 Subject: [PATCH 16/25] iqmc_kornreich_pi fails. Rehide iqmc tests temporarily --- .github/workflows/mpi_numba_reg.yml | 2 +- .github/workflows/mpi_reg.yml | 2 +- .github/workflows/numba_reg.yml | 2 +- .github/workflows/python_reg.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/mpi_numba_reg.yml b/.github/workflows/mpi_numba_reg.yml index fe291d34..d65daacb 100644 --- a/.github/workflows/mpi_numba_reg.yml +++ b/.github/workflows/mpi_numba_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Numba and MPI run: | cd test/regression - python run.py --mode=numba --mpiexec=4 + python run.py --mode=numba --mpiexec=4 --skip=iqmc* diff --git a/.github/workflows/mpi_reg.yml b/.github/workflows/mpi_reg.yml index 792f40e1..66ec4f85 100644 --- a/.github/workflows/mpi_reg.yml +++ b/.github/workflows/mpi_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - MPI run: | cd test/regression - python run.py --mpiexec=4 + python run.py --mpiexec=4 --skip=iqmc* diff --git a/.github/workflows/numba_reg.yml b/.github/workflows/numba_reg.yml index 493c0c9d..9d5124be 100644 --- a/.github/workflows/numba_reg.yml +++ b/.github/workflows/numba_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Numba run: | cd test/regression - python run.py --mode=numba + python run.py --mode=numba --skip=iqmc* diff --git a/.github/workflows/python_reg.yml b/.github/workflows/python_reg.yml index 5c811c4a..04718ea7 100644 --- a/.github/workflows/python_reg.yml +++ b/.github/workflows/python_reg.yml @@ -32,4 +32,4 @@ jobs: - name: Regression Test - Python run: | cd test/regression - python run.py + python run.py --skip=iqmc* From de1a8f14cbf520ecb39a4194cee1cae7bdb86773 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 29 Jul 2024 14:34:15 +0700 Subject: [PATCH 17/25] fix source_file test input mpi --- mcdc/main.py | 35 +++++++++++++++---------- test/regression/slab_ce/input.py | 2 +- test/regression/source_file/input.py | 38 +++++++++++++++------------- 3 files changed, 43 insertions(+), 32 deletions(-) diff --git a/mcdc/main.py b/mcdc/main.py index 25f0e604..ffa9d264 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -885,22 +885,29 @@ def prepare(): # ========================================================================= # Source file + # TODO: Use parallel h5py # ========================================================================= - if mcdc["setting"]["source_file"]: - with h5py.File(mcdc["setting"]["source_file_name"], "r") as f: - # Get source particle size - N_particle = f["particles_size"][()] - - # Redistribute work - kernel.distribute_work(N_particle, mcdc) - N_local = mcdc["mpi_work_size"] - start = mcdc["mpi_work_start"] - end = start + N_local - - # Add particles to source bank - mcdc["bank_source"]["particles"][:N_local] = f["particles"][start:end] - mcdc["bank_source"]["size"] = N_local + # All ranks, take turn + for i in range(mcdc["mpi_size"]): + if mcdc["mpi_rank"] == i: + if mcdc["setting"]["source_file"]: + with h5py.File(mcdc["setting"]["source_file_name"], "r") as f: + # Get source particle size + N_particle = f["particles_size"][()] + + # Redistribute work + kernel.distribute_work(N_particle, mcdc) + N_local = mcdc["mpi_work_size"] + start = mcdc["mpi_work_start"] + end = start + N_local + + # Add particles to source bank + mcdc["bank_source"]["particles"][:N_local] = f["particles"][ + start:end + ] + mcdc["bank_source"]["size"] = N_local + MPI.COMM_WORLD.Barrier() # ========================================================================= # IC file diff --git a/test/regression/slab_ce/input.py b/test/regression/slab_ce/input.py index 2220758a..ad519db2 100644 --- a/test/regression/slab_ce/input.py +++ b/test/regression/slab_ce/input.py @@ -8,7 +8,7 @@ # Set the XS library directory os.environ["MCDC_XSLIB"] = os.getcwd() -# Create the dummy nuclide (only master, rank 0 +# Create the dummy nuclide (only master, rank 0) if MPI.COMM_WORLD.Get_rank() == 0: with h5py.File("dummy_nuclide.h5", "w") as f: f["A"] = 1.0 diff --git a/test/regression/source_file/input.py b/test/regression/source_file/input.py index 77fe6743..2ffc2e26 100644 --- a/test/regression/source_file/input.py +++ b/test/regression/source_file/input.py @@ -1,28 +1,32 @@ import numpy as np +from mpi4py import MPI import mcdc # ============================================================================= # Set source particles and create the source particle file +# (only master, rank 0) # ============================================================================= -rng = np.random.default_rng(seed=7) - -N = 500 -bank = mcdc.make_particle_bank(N) - -for i in range(N): - particle = bank[i] - particle["x"] = rng.random() * 5.0 - particle["y"] = rng.random() - particle["z"] = rng.random() - particle["t"] = rng.random() * 5.0 - particle["ux"] = 1.0 # All going right - particle["uy"] = 0.0 - particle["uz"] = 0.0 - particle["E"] = 1e6 # Arbitrary - -mcdc.save_particle_bank(bank, "source_particles") +if MPI.COMM_WORLD.Get_rank() == 0: + rng = np.random.default_rng(seed=7) + + N = 500 + bank = mcdc.make_particle_bank(N) + + for i in range(N): + particle = bank[i] + particle["x"] = rng.random() * 5.0 + particle["y"] = rng.random() + particle["z"] = rng.random() + particle["t"] = rng.random() * 5.0 + particle["ux"] = 1.0 # All going right + particle["uy"] = 0.0 + particle["uz"] = 0.0 + particle["E"] = 1e6 # Arbitrary + + mcdc.save_particle_bank(bank, "source_particles") +MPI.COMM_WORLD.Barrier() # ============================================================================= From b4c6c409222b132c8a6434ed5ee98a866fe4c66f Mon Sep 17 00:00:00 2001 From: Kayla Clements Date: Tue, 6 Aug 2024 14:19:04 -0700 Subject: [PATCH 18/25] Added mpi4py.util.dtlib to docs MOCK_MODULES --- docs/source/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 68308523..168fc4d9 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -18,7 +18,7 @@ # On Read the Docs, need to mock any python packages that would require c from unittest.mock import MagicMock -MOCK_MODULES = ["mpi4py", "colorama"] +MOCK_MODULES = ["mpi4py", "colorama", "mpi4py.util.dtlib"] sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES) From 4914dbb9f1668ccebfe6c826c36a2e355c916b98 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Tue, 6 Aug 2024 15:53:56 -0700 Subject: [PATCH 19/25] moving visulizer to a sub package for mac support --- examples/visulization/godiva-mockup/input.py | 3 ++- mcdc/__init__.py | 1 - mcdc/visualizer/__init__.py | 1 + mcdc/{ => visualizer}/visualizer.py | 0 4 files changed, 3 insertions(+), 2 deletions(-) create mode 100644 mcdc/visualizer/__init__.py rename mcdc/{ => visualizer}/visualizer.py (100%) diff --git a/examples/visulization/godiva-mockup/input.py b/examples/visulization/godiva-mockup/input.py index acbd361b..729ca58a 100644 --- a/examples/visulization/godiva-mockup/input.py +++ b/examples/visulization/godiva-mockup/input.py @@ -1,5 +1,6 @@ import numpy as np import mcdc, h5py +from mcdc import visualizer # ============================================================================= # Materials @@ -73,7 +74,7 @@ mcdc.source(x=[-22.0, 22.0], time=[0.0, 5.0], isotropic=True) -mcdc.visualize( +visualizer.visualize( start_time=0, end_time=1, tick_interval=0.1, material_colors={"water": [1, 0, 0]} ) # ============================================================================= diff --git a/mcdc/__init__.py b/mcdc/__init__.py index 853af00a..c902b276 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -26,4 +26,3 @@ save_particle_bank, ) from mcdc.main import run, prepare -from mcdc.visualizer import visualize diff --git a/mcdc/visualizer/__init__.py b/mcdc/visualizer/__init__.py new file mode 100644 index 00000000..3f800a36 --- /dev/null +++ b/mcdc/visualizer/__init__.py @@ -0,0 +1 @@ +from mcdc.visualizer.visualizer import visualize diff --git a/mcdc/visualizer.py b/mcdc/visualizer/visualizer.py similarity index 100% rename from mcdc/visualizer.py rename to mcdc/visualizer/visualizer.py From 85fb12a5a5f8292083fd150e1ad0fecfce99718c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 08:03:24 +0700 Subject: [PATCH 20/25] set up the front end --- mcdc/card.py | 129 +++++++++++++++++++++++++++++++++++++++++------ mcdc/constant.py | 7 ++- mcdc/input_.py | 5 ++ 3 files changed, 124 insertions(+), 17 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index ecf66167..299fd375 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -1,6 +1,13 @@ import numpy as np +import sympy -from mcdc.constant import INF, SHIFT +from mcdc.constant import ( + BOOL_AND, + BOOL_OR, + BOOL_NOT, + INF, + SHIFT, +) # Get the global variable container import mcdc.global_ as global_ @@ -16,7 +23,7 @@ def __str__(self): for name in [ a for a in dir(self) - if not a.startswith("__") and not callable(getattr(self, a)) and a != "tag" + if not a.startswith("__") and not callable(getattr(self, a)) and a != "tag" and not a.startswith("_") ]: text += " %s : %s\n" % (name, str(getattr(self, name))) return text @@ -88,8 +95,8 @@ def __init__(self, type_): # Set card data self.ID = None self.type = type_ - self.A = -1 - self.B = -1 + self.A = None + self.B = None def __and__(self, other): region = RegionCard("intersection") @@ -117,6 +124,21 @@ def __invert__(self): global_.input_deck.regions.append(region) return region + def __str__(self): + if self.type == 'halfspace': + if self.B > 0.0: + return "+s%i"%self.A + else: + return "-s%i"%self.A + elif self.type == 'intersection': + return "r%i & r%i"%(self.A, self.B) + elif self.type == 'union': + return "r%i | r%i"%(self.A, self.B) + elif self.type == 'complement': + return "~r%i"%(self.A) + elif self.type == 'all': + return "all" + class SurfaceCard(InputCard): def __init__(self): @@ -146,23 +168,32 @@ def __init__(self): self.sensitivity_ID = 0 self.dsm_Np = 1.0 - def __pos__(self): + + def _create_halfspace(self, positive): region = RegionCard("halfspace") region.A = self.ID - region.B = 1 + if positive: + region.B = 1 + else: + region.B = -1 + + # Check if an identical halfspace region already existed + for idx, existing_region in enumerate(global_.input_deck.regions): + if region.A == existing_region.A and region.B == existing_region.B: + return global_.input_deck.regions[idx] + # Set ID and push to deck region.ID = len(global_.input_deck.regions) global_.input_deck.regions.append(region) return region + + def __pos__(self): + return self._create_halfspace(True) + + def __neg__(self): - region = RegionCard("halfspace") - region.A = self.ID - region.B = 0 - # Set ID and push to deck - region.ID = len(global_.input_deck.regions) - global_.input_deck.regions.append(region) - return region + return self._create_halfspace(False) class CellCard(InputCard): @@ -171,12 +202,78 @@ def __init__(self): # Set card data self.ID = None - self.region_ID = -1 + self.region_ID = None + self.region = "all" self.fill_type = "material" - self.fill_ID = -1 + self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) self.N_surface = 0 - self.surface_ID = np.zeros(0, dtype=int) + self.surface_IDs = np.zeros(0, dtype=int) + self._region_RPN = [] # Reverse Polish Notation + + def set_region_RPN(self): + # Make alias and reset + rpn = self._region_RPN + rpn.clear() + + # Build RPN based on the assigned region + region = global_.input_deck.regions[self.region_ID] + stack = [region] + while len(stack) > 0: + token = stack.pop() + if isinstance(token, RegionCard): + if token.type == 'halfspace': + rpn.append(token.ID) + elif token.type == 'intersection': + region_A = global_.input_deck.regions[token.A] + region_B = global_.input_deck.regions[token.B] + stack += ['&', region_A, region_B] + elif token.type == 'union': + region_A = global_.input_deck.regions[token.A] + region_B = global_.input_deck.regions[token.B] + stack += ['|', region_A, region_B] + elif token.type == 'complement': + region = global_.input_deck.regions[token.A] + stack += ['~', region] + else: + if token == '&': + rpn.append(BOOL_AND) + elif token == '|': + rpn.append(BOOL_OR) + elif token == '~': + rpn.append(BOOL_NOT) + else: + print_error("Something is wrong with cell RPN creation.") + + def set_region(self): + stack = [] + + for token in self._region_RPN: + if token >= 0: + stack.append(token) + else: + if token == BOOL_AND or token == BOOL_OR: + item_1 = stack.pop() + if isinstance(item_1, int): + item_1 = sympy.symbols(str(global_.input_deck.regions[item_1])) + + item_2 = stack.pop() + if isinstance(item_2, int): + item_2 = sympy.symbols(str(global_.input_deck.regions[item_2])) + + if token == BOOL_AND: + stack.append(item_1 & item_2) + else: + stack.append(item_1 | item_2) + + elif token == BOOL_NOT: + item = stack.pop() + if isinstance(item, int): + item = sympy.symbols(str(global_.input_deck.regions[item])) + stack.append(~item) + + self.region = sympy.logic.boolalg.simplify_logic(stack[0]) + class UniverseCard(InputCard): diff --git a/mcdc/constant.py b/mcdc/constant.py index ab603f52..da60c12a 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -20,7 +20,12 @@ REGION_COMPLEMENT = 3 REGION_ALL = 4 -# UNIVERSE +# Boolean operator +BOOL_AND = -1 +BOOL_OR = -2 +BOOL_NOT = -3 + +# Universe UNIVERSE_ROOT = 0 # Events diff --git a/mcdc/input_.py b/mcdc/input_.py index f0a9306f..02079bed 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -732,6 +732,11 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): # Assign region card.region_ID = region.ID + # Set region Reverse Polish Notation and region description + if region.type != 'all': + card.set_region_RPN() + card.set_region() + # Assign fill type and ID if fill.tag == "Material": card.fill_type = "material" From 6cba33b4c04f6b1a2c8f69268550708d16e84c91 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 08:16:00 +0700 Subject: [PATCH 21/25] add set_surface_IDs to CellCard --- mcdc/card.py | 11 +++++++++++ mcdc/input_.py | 27 +-------------------------- 2 files changed, 12 insertions(+), 26 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index 299fd375..797619a5 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -274,6 +274,17 @@ def set_region(self): self.region = sympy.logic.boolalg.simplify_logic(stack[0]) + def set_surface_IDs(self): + surface_IDs = [] + + for token in self._region_RPN: + if token >=0: + ID = global_.input_deck.regions[token].A + if not ID in surface_IDs: + surface_IDs.append(ID) + + self.surface_IDs = np.sort(np.array(surface_IDs)) + self.N_surface = len(surface_IDs) class UniverseCard(InputCard): diff --git a/mcdc/input_.py b/mcdc/input_.py index 02079bed..848577f2 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -750,14 +750,7 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): card.translation[:] = translation # Get all surface IDs - surface_IDs = [] - - region = global_.input_deck.regions[card.region_ID] - get_all_surface_IDs(region, surface_IDs) - surface_IDs = list(set(surface_IDs)) - - card.N_surface = len(surface_IDs) - card.surface_IDs = np.array(surface_IDs) + card.set_surface_IDs() # Add to deck global_.input_deck.cells.append(card) @@ -765,24 +758,6 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): return card -def get_all_surface_IDs(region, surface_IDs): - if region.type == "halfspace": - surface = global_.input_deck.surfaces[region.A] - surface_IDs.append(surface.ID) - - elif region.type in ["intersection", "union"]: - region1 = global_.input_deck.regions[region.A] - region2 = global_.input_deck.regions[region.B] - - get_all_surface_IDs(region1, surface_IDs) - get_all_surface_IDs(region2, surface_IDs) - - elif region.type in ["complement"]: - region1 = global_.input_deck.regions[region.A] - - get_all_surface_IDs(region1, surface_IDs) - - def universe(cells, root=False): """ Define a list of cells as a universe. From bc68a49303a0935fda9fa739ced4b331cf6ccd3e Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 12:18:52 +0700 Subject: [PATCH 22/25] working on kobayashi-new --- mcdc/adapt.py | 10 +++++++ mcdc/card.py | 4 +-- mcdc/kernel.py | 51 ++++++++++++++++++++++++++++++----- mcdc/main.py | 73 ++++++++++++++++++++++++++++---------------------- mcdc/type_.py | 68 ++++++++++++++++++++-------------------------- 5 files changed, 126 insertions(+), 80 deletions(-) diff --git a/mcdc/adapt.py b/mcdc/adapt.py index 4cb2ac92..d80e6def 100644 --- a/mcdc/adapt.py +++ b/mcdc/adapt.py @@ -416,6 +416,16 @@ def local_j_array(): return cuda.local.array(1, type_.j_array)[0] +@for_cpu() +def local_RPN_array(): + return np.zeros(1, dtype=type_.RPN_array)[0] + + +@for_gpu() +def local_RPN_array(): + return cuda.local.array(1, type_.RPN_array)[0] + + @for_cpu() def local_particle(): return np.zeros(1, dtype=type_.particle)[0] diff --git a/mcdc/card.py b/mcdc/card.py index 797619a5..d107cd9f 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -126,7 +126,7 @@ def __invert__(self): def __str__(self): if self.type == 'halfspace': - if self.B > 0.0: + if self.B > 0: return "+s%i"%self.A else: return "-s%i"%self.A @@ -207,7 +207,6 @@ def __init__(self): self.fill_type = "material" self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) - self.N_surface = 0 self.surface_IDs = np.zeros(0, dtype=int) self._region_RPN = [] # Reverse Polish Notation @@ -284,7 +283,6 @@ def set_surface_IDs(self): surface_IDs.append(ID) self.surface_IDs = np.sort(np.array(surface_IDs)) - self.N_surface = len(surface_IDs) class UniverseCard(InputCard): diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 4e689ee7..72461eb3 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1708,8 +1708,40 @@ def split_particle(P): @njit def cell_check(P, cell, trans, mcdc): - region = mcdc["regions"][cell["region_ID"]] - return region_check(P, region, trans, mcdc) + # Access RPN data + idx = cell['region_data_idx'] + N_token = mcdc['cell_region_data'][idx] + + # Create local value array + value_struct = adapt.local_RPN_array() + value = value_struct['values'] + N_value = 0 + + # March forward through RPN tokens + idx += 1 + idx_end = idx + N_token + while idx < idx_end: + token = mcdc['cell_region_data'][idx] + + if token >= 0: + surface = mcdc['surfaces'][token] + value[N_value] = surface_evaluate(P, surface, trans) > 0.0 + N_value += 1 + + elif token == BOOL_NOT: + value[N_value - 1] = not value[N_value - 1] + + elif token == BOOL_AND: + value[N_value - 2] = value[N_value - 2] & value[N_value - 1] + N_value -= 1 + + elif token == BOOL_OR: + value[N_value - 2] = value[N_value - 2] | value[N_value - 1] + N_value -= 1 + + idx += 1 + + return value[0] @njit @@ -1864,7 +1896,6 @@ def surface_shift(P, surface, trans, mcdc): dot = ux * nx + uy * ny + uz * nz + J1 / v else: dot = ux * nx + uy * ny + uz * nz - if dot > 0.0: P["x"] += shift_x P["y"] += shift_y @@ -2531,7 +2562,6 @@ def move_to_event(P, mcdc): # Also set particle material and speed d_boundary, event = distance_to_boundary(P, mcdc) - # print(P["material_ID"]) # Distance to tally mesh d_mesh = INF if mcdc["cycle_active"]: @@ -2693,13 +2723,22 @@ def distance_to_nearest_surface(P, cell, trans, mcdc): surface_ID = -1 surface_move = False - for i in range(cell["N_surface"]): - surface = mcdc["surfaces"][cell["surface_IDs"][i]] + # Access cell surface data + idx = cell['surface_data_idx'] + N_surface = mcdc['cell_surface_data'][idx] + + # Iterate over all surfaces + idx += 1 + idx_end = idx + N_surface + while idx < idx_end: + candidate_surface_ID = mcdc['cell_surface_data'][idx] + surface = mcdc["surfaces"][candidate_surface_ID] d, sm = surface_distance(P, surface, trans, mcdc) if d < distance: distance = d surface_ID = surface["ID"] surface_move = sm + idx += 1 return distance, surface_ID, surface_move diff --git a/mcdc/main.py b/mcdc/main.py index ffa9d264..2803540f 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -337,6 +337,24 @@ def prepare(): root_universe.cell_IDs[i] = cell.ID input_deck.universes[0] = root_universe + # ========================================================================= + # Prepare cell region RPN (Reverse Polish Notation) + # - Replace halfspace region ID with its surface and insert + # complement operator if the sense is negative. + # ========================================================================= + + for cell in input_deck.cells: + i = 0 + while i < len(cell._region_RPN): + token = cell._region_RPN[i] + if token >= 0: + surface_ID = input_deck.regions[token].A + sense = input_deck.regions[token].B + cell._region_RPN[i] = surface_ID + if sense < 0: + cell._region_RPN.insert(i + 1, BOOL_NOT) + i += 1 + # ========================================================================= # Adapt kernels # ========================================================================= @@ -352,8 +370,6 @@ def prepare(): type_.make_type_nuclide(input_deck) type_.make_type_material(input_deck) type_.make_type_surface(input_deck) - type_.make_type_region() - type_.make_type_cell(input_deck) type_.make_type_universe(input_deck) type_.make_type_lattice(input_deck) type_.make_type_source(input_deck) @@ -370,6 +386,7 @@ def prepare(): type_.make_type_translate(input_deck) type_.make_type_group_array(input_deck) type_.make_type_j_array(input_deck) + type_.make_type_RPN_array(input_deck) # ========================================================================= # Create the global variable container @@ -519,37 +536,16 @@ def prepare(): N = len(getattr(input_deck.surfaces[i], name)) mcdc["surfaces"][i][name][:N] = getattr(input_deck.surfaces[i], name) - # ========================================================================= - # Regions - # ========================================================================= - - N_region = len(input_deck.regions) - for i in range(N_region): - for name in type_.region.names: - if name not in ["type"]: - copy_field(mcdc["regions"][i], input_deck.regions[i], name) - - # Type - if input_deck.regions[i].type == "halfspace": - mcdc["regions"][i]["type"] = REGION_HALFSPACE - elif input_deck.regions[i].type == "intersection": - mcdc["regions"][i]["type"] = REGION_INTERSECTION - elif input_deck.regions[i].type == "union": - mcdc["regions"][i]["type"] = REGION_UNION - elif input_deck.regions[i].type == "complement": - mcdc["regions"][i]["type"] = REGION_COMPLEMENT - elif input_deck.regions[i].type == "all": - mcdc["regions"][i]["type"] = REGION_ALL - # ========================================================================= # Cells # ========================================================================= N_cell = len(input_deck.cells) + surface_data_idx = 0 + region_data_idx = 0 for i in range(N_cell): - for name in type_.cell.names: - if name not in ["fill_type", "surface_IDs"]: - copy_field(mcdc["cells"][i], input_deck.cells[i], name) + for name in ['ID', 'fill_ID', 'translation']: + copy_field(mcdc["cells"][i], input_deck.cells[i], name) # Fill type if input_deck.cells[i].fill_type == "material": @@ -559,10 +555,19 @@ def prepare(): elif input_deck.cells[i].fill_type == "lattice": mcdc["cells"][i]["fill_type"] = FILL_LATTICE - # Variables with possible different sizes - for name in ["surface_IDs"]: - N = mcdc["cells"][i]["N_surface"] - mcdc["cells"][i][name][:N] = getattr(input_deck.cells[i], name) + # Surface data + mcdc["cells"][i]['surface_data_idx'] = surface_data_idx + N_surface = len(input_deck.cells[i].surface_IDs) + mcdc['cell_surface_data'][surface_data_idx] = N_surface + mcdc['cell_surface_data'][surface_data_idx + 1: surface_data_idx + N_surface + 1] = input_deck.cells[i].surface_IDs + surface_data_idx += N_surface + 1 + + # Region data + mcdc["cells"][i]['region_data_idx'] = region_data_idx + N_RPN = len(input_deck.cells[i]._region_RPN) + mcdc['cell_region_data'][region_data_idx] = N_RPN + mcdc['cell_region_data'][region_data_idx + 1: region_data_idx + N_RPN + 1] = input_deck.cells[i]._region_RPN + region_data_idx += N_RPN + 1 # ========================================================================= # Universes @@ -981,7 +986,11 @@ def card_to_h5group(card, group): elif value is None: next else: - group[name] = value + if name not in ['region']: + group[name] = value + + elif name == 'region': + group[name] = str(value) def dictlist_to_h5group(dictlist, input_group, name): diff --git a/mcdc/type_.py b/mcdc/type_.py index 3dfb9b2c..b17c38d6 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -34,8 +34,6 @@ nuclide = None material = None surface = None -region = None -cell = None universe = None lattice = None source = None @@ -47,6 +45,8 @@ translate = None group_array = None j_array = None +RPN_array = None + global_ = None @@ -513,46 +513,21 @@ def make_type_surface(input_deck): ) -# ============================================================================== -# Region -# ============================================================================== - - -def make_type_region(): - global region - - region = into_dtype( - [ - ("ID", int64), - ("type", int64), - ("A", int64), - ("B", int64), - ] - ) - - # ============================================================================== # Cell # ============================================================================== -def make_type_cell(input_deck): - global cell - - # Maximum number of surfaces per cell - Nmax_surface = max([cell.N_surface for cell in input_deck.cells]) - - cell = into_dtype( - [ - ("ID", int64), - ("region_ID", int64), - ("fill_type", int64), - ("fill_ID", int64), - ("translation", float64, (3,)), - ("N_surface", int64), - ("surface_IDs", int64, (Nmax_surface,)), - ] - ) +cell = into_dtype( + [ + ("ID", int64), + ("fill_type", int64), + ("fill_ID", int64), + ("translation", float64, (3,)), + ("surface_data_idx", int64), + ("region_data_idx", int64), + ] +) # ============================================================================== @@ -1233,12 +1208,18 @@ def make_type_global(input_deck): N_nuclide = len(input_deck.nuclides) N_material = len(input_deck.materials) N_surface = len(input_deck.surfaces) - N_region = len(input_deck.regions) N_cell = len(input_deck.cells) N_source = len(input_deck.sources) N_universe = len(input_deck.universes) N_lattice = len(input_deck.lattices) + # Cell data sizes + N_cell_surface = 0 + N_cell_region = 0 + for cell_ in input_deck.cells: + N_cell_surface += 1 + len(cell_.surface_IDs) + N_cell_region += 1 + len(cell_._region_RPN) + # Simulation parameters N_particle = input_deck.setting["N_particle"] N_precursor = input_deck.setting["N_precursor"] @@ -1298,8 +1279,9 @@ def make_type_global(input_deck): ("nuclides", nuclide, (N_nuclide,)), ("materials", material, (N_material,)), ("surfaces", surface, (N_surface,)), - ("regions", region, (N_region,)), ("cells", cell, (N_cell,)), + ("cell_surface_data", int64, (N_cell_surface,)), + ("cell_region_data", int64, (N_cell_region,)), ("universes", universe, (N_universe,)), ("lattices", lattice, (N_lattice,)), ("sources", source, (N_source,)), @@ -1385,6 +1367,14 @@ def make_type_j_array(input_deck): j_array = into_dtype([("values", float64, (J,))]) +def make_type_RPN_array(input_deck): + global RPN_array + N_max = 0 + for cell_ in input_deck.cells: + N = np.sum(np.array(cell_._region_RPN) >= 0.0) + N_max = max(N_max, N) + RPN_array = into_dtype([("values", bool_, (N_max,))]) + def make_type_mesh(card): Nx = len(card["x"]) - 1 Ny = len(card["y"]) - 1 From b01a039759b6ef2d7554fb3be81b72a15ce741c7 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:48:32 +0700 Subject: [PATCH 23/25] debug halfspace checking --- mcdc/card.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcdc/card.py b/mcdc/card.py index d107cd9f..a5670127 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -179,7 +179,7 @@ def _create_halfspace(self, positive): # Check if an identical halfspace region already existed for idx, existing_region in enumerate(global_.input_deck.regions): - if region.A == existing_region.A and region.B == existing_region.B: + if existing_region.type == 'halfspace' and region.A == existing_region.A and region.B == existing_region.B: return global_.input_deck.regions[idx] # Set ID and push to deck From bc81119611f867df38bd198954f9e26e020f53ff Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:51:12 +0700 Subject: [PATCH 24/25] add sympy to pyproject.toml --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 3dc378b2..4fed70cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ dependencies = [ "h5py", "colorama", "black", + "sympy", ] [project.optional-dependencies] docs = ["sphinx==7.2.6", "furo", "sphinx_toolbox"] From f793529b6b391f10d5e0a86cf4563d3e4b7be50a Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:53:40 +0700 Subject: [PATCH 25/25] back in black --- mcdc/card.py | 58 +++++++++++++++++++++++++++----------------------- mcdc/input_.py | 2 +- mcdc/kernel.py | 16 +++++++------- mcdc/main.py | 22 +++++++++++-------- mcdc/type_.py | 1 + 5 files changed, 54 insertions(+), 45 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index a5670127..fa53a1fa 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -23,7 +23,10 @@ def __str__(self): for name in [ a for a in dir(self) - if not a.startswith("__") and not callable(getattr(self, a)) and a != "tag" and not a.startswith("_") + if not a.startswith("__") + and not callable(getattr(self, a)) + and a != "tag" + and not a.startswith("_") ]: text += " %s : %s\n" % (name, str(getattr(self, name))) return text @@ -125,18 +128,18 @@ def __invert__(self): return region def __str__(self): - if self.type == 'halfspace': + if self.type == "halfspace": if self.B > 0: - return "+s%i"%self.A + return "+s%i" % self.A else: - return "-s%i"%self.A - elif self.type == 'intersection': - return "r%i & r%i"%(self.A, self.B) - elif self.type == 'union': - return "r%i | r%i"%(self.A, self.B) - elif self.type == 'complement': - return "~r%i"%(self.A) - elif self.type == 'all': + return "-s%i" % self.A + elif self.type == "intersection": + return "r%i & r%i" % (self.A, self.B) + elif self.type == "union": + return "r%i | r%i" % (self.A, self.B) + elif self.type == "complement": + return "~r%i" % (self.A) + elif self.type == "all": return "all" @@ -168,7 +171,6 @@ def __init__(self): self.sensitivity_ID = 0 self.dsm_Np = 1.0 - def _create_halfspace(self, positive): region = RegionCard("halfspace") region.A = self.ID @@ -179,7 +181,11 @@ def _create_halfspace(self, positive): # Check if an identical halfspace region already existed for idx, existing_region in enumerate(global_.input_deck.regions): - if existing_region.type == 'halfspace' and region.A == existing_region.A and region.B == existing_region.B: + if ( + existing_region.type == "halfspace" + and region.A == existing_region.A + and region.B == existing_region.B + ): return global_.input_deck.regions[idx] # Set ID and push to deck @@ -187,11 +193,9 @@ def _create_halfspace(self, positive): global_.input_deck.regions.append(region) return region - def __pos__(self): return self._create_halfspace(True) - def __neg__(self): return self._create_halfspace(False) @@ -208,7 +212,7 @@ def __init__(self): self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) self.surface_IDs = np.zeros(0, dtype=int) - self._region_RPN = [] # Reverse Polish Notation + self._region_RPN = [] # Reverse Polish Notation def set_region_RPN(self): # Make alias and reset @@ -221,25 +225,25 @@ def set_region_RPN(self): while len(stack) > 0: token = stack.pop() if isinstance(token, RegionCard): - if token.type == 'halfspace': + if token.type == "halfspace": rpn.append(token.ID) - elif token.type == 'intersection': + elif token.type == "intersection": region_A = global_.input_deck.regions[token.A] region_B = global_.input_deck.regions[token.B] - stack += ['&', region_A, region_B] - elif token.type == 'union': + stack += ["&", region_A, region_B] + elif token.type == "union": region_A = global_.input_deck.regions[token.A] region_B = global_.input_deck.regions[token.B] - stack += ['|', region_A, region_B] - elif token.type == 'complement': + stack += ["|", region_A, region_B] + elif token.type == "complement": region = global_.input_deck.regions[token.A] - stack += ['~', region] + stack += ["~", region] else: - if token == '&': + if token == "&": rpn.append(BOOL_AND) - elif token == '|': + elif token == "|": rpn.append(BOOL_OR) - elif token == '~': + elif token == "~": rpn.append(BOOL_NOT) else: print_error("Something is wrong with cell RPN creation.") @@ -277,7 +281,7 @@ def set_surface_IDs(self): surface_IDs = [] for token in self._region_RPN: - if token >=0: + if token >= 0: ID = global_.input_deck.regions[token].A if not ID in surface_IDs: surface_IDs.append(ID) diff --git a/mcdc/input_.py b/mcdc/input_.py index 848577f2..02c14c77 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -733,7 +733,7 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): card.region_ID = region.ID # Set region Reverse Polish Notation and region description - if region.type != 'all': + if region.type != "all": card.set_region_RPN() card.set_region() diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 72461eb3..a02f321f 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1709,22 +1709,22 @@ def split_particle(P): @njit def cell_check(P, cell, trans, mcdc): # Access RPN data - idx = cell['region_data_idx'] - N_token = mcdc['cell_region_data'][idx] + idx = cell["region_data_idx"] + N_token = mcdc["cell_region_data"][idx] # Create local value array value_struct = adapt.local_RPN_array() - value = value_struct['values'] + value = value_struct["values"] N_value = 0 # March forward through RPN tokens idx += 1 idx_end = idx + N_token while idx < idx_end: - token = mcdc['cell_region_data'][idx] + token = mcdc["cell_region_data"][idx] if token >= 0: - surface = mcdc['surfaces'][token] + surface = mcdc["surfaces"][token] value[N_value] = surface_evaluate(P, surface, trans) > 0.0 N_value += 1 @@ -2724,14 +2724,14 @@ def distance_to_nearest_surface(P, cell, trans, mcdc): surface_move = False # Access cell surface data - idx = cell['surface_data_idx'] - N_surface = mcdc['cell_surface_data'][idx] + idx = cell["surface_data_idx"] + N_surface = mcdc["cell_surface_data"][idx] # Iterate over all surfaces idx += 1 idx_end = idx + N_surface while idx < idx_end: - candidate_surface_ID = mcdc['cell_surface_data'][idx] + candidate_surface_ID = mcdc["cell_surface_data"][idx] surface = mcdc["surfaces"][candidate_surface_ID] d, sm = surface_distance(P, surface, trans, mcdc) if d < distance: diff --git a/mcdc/main.py b/mcdc/main.py index 2803540f..d47f7e91 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -544,7 +544,7 @@ def prepare(): surface_data_idx = 0 region_data_idx = 0 for i in range(N_cell): - for name in ['ID', 'fill_ID', 'translation']: + for name in ["ID", "fill_ID", "translation"]: copy_field(mcdc["cells"][i], input_deck.cells[i], name) # Fill type @@ -556,17 +556,21 @@ def prepare(): mcdc["cells"][i]["fill_type"] = FILL_LATTICE # Surface data - mcdc["cells"][i]['surface_data_idx'] = surface_data_idx + mcdc["cells"][i]["surface_data_idx"] = surface_data_idx N_surface = len(input_deck.cells[i].surface_IDs) - mcdc['cell_surface_data'][surface_data_idx] = N_surface - mcdc['cell_surface_data'][surface_data_idx + 1: surface_data_idx + N_surface + 1] = input_deck.cells[i].surface_IDs + mcdc["cell_surface_data"][surface_data_idx] = N_surface + mcdc["cell_surface_data"][ + surface_data_idx + 1 : surface_data_idx + N_surface + 1 + ] = input_deck.cells[i].surface_IDs surface_data_idx += N_surface + 1 # Region data - mcdc["cells"][i]['region_data_idx'] = region_data_idx + mcdc["cells"][i]["region_data_idx"] = region_data_idx N_RPN = len(input_deck.cells[i]._region_RPN) - mcdc['cell_region_data'][region_data_idx] = N_RPN - mcdc['cell_region_data'][region_data_idx + 1: region_data_idx + N_RPN + 1] = input_deck.cells[i]._region_RPN + mcdc["cell_region_data"][region_data_idx] = N_RPN + mcdc["cell_region_data"][region_data_idx + 1 : region_data_idx + N_RPN + 1] = ( + input_deck.cells[i]._region_RPN + ) region_data_idx += N_RPN + 1 # ========================================================================= @@ -986,10 +990,10 @@ def card_to_h5group(card, group): elif value is None: next else: - if name not in ['region']: + if name not in ["region"]: group[name] = value - elif name == 'region': + elif name == "region": group[name] = str(value) diff --git a/mcdc/type_.py b/mcdc/type_.py index b17c38d6..43b0e333 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -1375,6 +1375,7 @@ def make_type_RPN_array(input_deck): N_max = max(N_max, N) RPN_array = into_dtype([("values", bool_, (N_max,))]) + def make_type_mesh(card): Nx = len(card["x"]) - 1 Ny = len(card["y"]) - 1