From bd7e8b040b85c4f3c068b6c4817c716dcd8beb1f Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Fri, 9 Oct 2020 08:45:07 +0200 Subject: [PATCH] First commit of code. New functionality: - DSMC module: - 2D cyl geometry - GMSH file format - Elastic cross-section for Argon-Argon collisions. - Basic boundary conditions: reflection, absorption and axis symmetry. Bugs fixed: Other comments: - Still searching for name. --- .gitignore | 6 + COPYING | 674 ++++++++++++++++++++++++++ README.md | 10 + data/collisions/Ar-Ar_Elastic.dat | 302 ++++++++++++ doc/PPartiC_UserManual.pdf | Bin 0 -> 82316 bytes doc/logos/PPartiC.eps | 177 +++++++ doc/logos/PPartiC.svg | 159 ++++++ makefile | 34 ++ src/DSMC_Neutrals.f95 | 110 +++++ src/makefile | 17 + src/modules/makefile | 49 ++ src/modules/moduleBoundary.f95 | 36 ++ src/modules/moduleCaseParam.f95 | 8 + src/modules/moduleCollisions.f95 | 152 ++++++ src/modules/moduleCompTime.f95 | 10 + src/modules/moduleConstParam.f95 | 14 + src/modules/moduleErrors.f95 | 43 ++ src/modules/moduleInject.f95 | 155 ++++++ src/modules/moduleInput.f95 | 366 ++++++++++++++ src/modules/moduleList.f95 | 79 +++ src/modules/moduleMesh.f95 | 198 ++++++++ src/modules/moduleMeshCyl.f95 | 534 ++++++++++++++++++++ src/modules/moduleMeshCylBoundary.f95 | 80 +++ src/modules/moduleMeshCylRead.f95 | 428 ++++++++++++++++ src/modules/moduleOutput.f95 | 118 +++++ src/modules/moduleRefParam.f95 | 10 + src/modules/moduleSolver.f95 | 113 +++++ src/modules/moduleSpecies.f95 | 66 +++ src/modules/moduleTable.f95 | 121 +++++ 29 files changed, 4069 insertions(+) create mode 100644 .gitignore create mode 100644 COPYING create mode 100644 README.md create mode 100644 data/collisions/Ar-Ar_Elastic.dat create mode 100644 doc/PPartiC_UserManual.pdf create mode 100644 doc/logos/PPartiC.eps create mode 100644 doc/logos/PPartiC.svg create mode 100644 makefile create mode 100644 src/DSMC_Neutrals.f95 create mode 100644 src/makefile create mode 100644 src/modules/makefile create mode 100644 src/modules/moduleBoundary.f95 create mode 100644 src/modules/moduleCaseParam.f95 create mode 100644 src/modules/moduleCollisions.f95 create mode 100644 src/modules/moduleCompTime.f95 create mode 100644 src/modules/moduleConstParam.f95 create mode 100644 src/modules/moduleErrors.f95 create mode 100644 src/modules/moduleInject.f95 create mode 100644 src/modules/moduleInput.f95 create mode 100644 src/modules/moduleList.f95 create mode 100644 src/modules/moduleMesh.f95 create mode 100644 src/modules/moduleMeshCyl.f95 create mode 100644 src/modules/moduleMeshCylBoundary.f95 create mode 100644 src/modules/moduleMeshCylRead.f95 create mode 100644 src/modules/moduleOutput.f95 create mode 100644 src/modules/moduleRefParam.f95 create mode 100644 src/modules/moduleSolver.f95 create mode 100644 src/modules/moduleSpecies.f95 create mode 100644 src/modules/moduleTable.f95 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0a9a292 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +DSMC_Neutrals +mod/ +obj/ +doc/user_manual/ +doc/coding_style/ +json-fortran-8.2.0/ diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..f288702 --- /dev/null +++ b/COPYING @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/README.md b/README.md new file mode 100644 index 0000000..61e1a8e --- /dev/null +++ b/README.md @@ -0,0 +1,10 @@ +# Introduction +Welcome to PPartiC (Plasma Particle Code), a modern object oriented code for kinetic simulations of plasma. This code works by simulating charged and neutral particles, following their trajectories, collisions and boundary conditions imposed by the usser. + +This code is currenlty in very early steps of development. + +The code aims to be easy to maintain and easy to use, allowing its application from complex problems to easy examples that can be used, for example, as teaching exercies. + +Parallelization techniches as OpenMP, MPI will be used, as well as making the code run in GPU architectures. + +PPartiC uses Finite Element meshes to adap to complex geometries. diff --git a/data/collisions/Ar-Ar_Elastic.dat b/data/collisions/Ar-Ar_Elastic.dat new file mode 100644 index 0000000..34a8118 --- /dev/null +++ b/data/collisions/Ar-Ar_Elastic.dat @@ -0,0 +1,302 @@ +# A. V. Phelps et. al, J. Phys. B, 33(2000) 2965-2981 +# Relative energy (eV) cross section (m^2) +0.010 1.325E-17 +0.020 1.004E-17 +0.030 8.538E-18 +0.040 7.610E-18 +0.050 6.960E-18 +0.060 6.471E-18 +0.070 6.084E-18 +0.080 5.767E-18 +0.090 5.502E-18 +0.100 5.275E-18 +0.110 5.078E-18 +0.120 4.904E-18 +0.130 4.749E-18 +0.140 4.611E-18 +0.150 4.485E-18 +0.160 4.371E-18 +0.170 4.266E-18 +0.180 4.170E-18 +0.190 4.081E-18 +0.200 3.998E-18 +0.210 3.921E-18 +0.220 3.848E-18 +0.230 3.780E-18 +0.240 3.717E-18 +0.250 3.656E-18 +0.260 3.600E-18 +0.270 3.546E-18 +0.280 3.494E-18 +0.290 3.446E-18 +0.300 3.399E-18 +0.310 3.355E-18 +0.320 3.313E-18 +0.330 3.272E-18 +0.340 3.233E-18 +0.350 3.196E-18 +0.360 3.160E-18 +0.370 3.126E-18 +0.380 3.093E-18 +0.390 3.061E-18 +0.400 3.030E-18 +0.410 3.000E-18 +0.420 2.971E-18 +0.430 2.944E-18 +0.440 2.917E-18 +0.450 2.891E-18 +0.460 2.865E-18 +0.470 2.841E-18 +0.480 2.817E-18 +0.490 2.794E-18 +0.500 2.771E-18 +0.510 2.750E-18 +0.520 2.728E-18 +0.530 2.708E-18 +0.540 2.688E-18 +0.550 2.668E-18 +0.560 2.649E-18 +0.570 2.630E-18 +0.580 2.612E-18 +0.590 2.594E-18 +0.600 2.577E-18 +0.610 2.560E-18 +0.620 2.543E-18 +0.630 2.527E-18 +0.640 2.511E-18 +0.650 2.496E-18 +0.660 2.480E-18 +0.670 2.466E-18 +0.680 2.451E-18 +0.690 2.437E-18 +0.700 2.423E-18 +0.710 2.409E-18 +0.720 2.396E-18 +0.730 2.383E-18 +0.740 2.370E-18 +0.750 2.357E-18 +0.760 2.345E-18 +0.770 2.332E-18 +0.780 2.320E-18 +0.790 2.309E-18 +0.800 2.297E-18 +0.810 2.286E-18 +0.820 2.275E-18 +0.830 2.264E-18 +0.840 2.253E-18 +0.850 2.242E-18 +0.860 2.232E-18 +0.870 2.221E-18 +0.880 2.211E-18 +0.890 2.201E-18 +0.900 2.192E-18 +0.910 2.182E-18 +0.920 2.173E-18 +0.930 2.163E-18 +0.940 2.154E-18 +0.950 2.145E-18 +0.960 2.136E-18 +0.970 2.127E-18 +0.980 2.118E-18 +0.990 2.110E-18 +1.000 2.101E-18 +1.010 2.093E-18 +1.020 2.085E-18 +1.030 2.077E-18 +1.040 2.069E-18 +1.050 2.061E-18 +1.060 2.053E-18 +1.070 2.046E-18 +1.080 2.038E-18 +1.090 2.031E-18 +1.100 2.023E-18 +1.110 2.016E-18 +1.120 2.009E-18 +1.130 2.002E-18 +1.140 1.995E-18 +1.150 1.988E-18 +1.160 1.981E-18 +1.170 1.974E-18 +1.180 1.967E-18 +1.190 1.961E-18 +1.200 1.954E-18 +1.210 1.948E-18 +1.220 1.941E-18 +1.230 1.935E-18 +1.240 1.929E-18 +1.250 1.923E-18 +1.260 1.917E-18 +1.270 1.911E-18 +1.280 1.905E-18 +1.290 1.899E-18 +1.300 1.893E-18 +1.310 1.887E-18 +1.320 1.882E-18 +1.330 1.876E-18 +1.340 1.870E-18 +1.350 1.865E-18 +1.360 1.859E-18 +1.370 1.854E-18 +1.380 1.849E-18 +1.390 1.843E-18 +1.400 1.838E-18 +1.410 1.833E-18 +1.420 1.828E-18 +1.430 1.823E-18 +1.440 1.818E-18 +1.450 1.813E-18 +1.460 1.808E-18 +1.470 1.803E-18 +1.480 1.798E-18 +1.490 1.793E-18 +1.500 1.788E-18 +1.510 1.784E-18 +1.520 1.779E-18 +1.530 1.774E-18 +1.540 1.770E-18 +1.550 1.765E-18 +1.560 1.761E-18 +1.570 1.756E-18 +1.580 1.752E-18 +1.590 1.748E-18 +1.600 1.743E-18 +1.610 1.739E-18 +1.620 1.735E-18 +1.630 1.730E-18 +1.640 1.726E-18 +1.650 1.722E-18 +1.660 1.718E-18 +1.670 1.714E-18 +1.680 1.710E-18 +1.690 1.706E-18 +1.700 1.702E-18 +1.710 1.698E-18 +1.720 1.694E-18 +1.730 1.690E-18 +1.740 1.686E-18 +1.750 1.682E-18 +1.760 1.679E-18 +1.770 1.675E-18 +1.780 1.671E-18 +1.790 1.667E-18 +1.800 1.664E-18 +1.810 1.660E-18 +1.820 1.657E-18 +1.830 1.653E-18 +1.840 1.649E-18 +1.850 1.646E-18 +1.860 1.642E-18 +1.870 1.639E-18 +1.880 1.635E-18 +1.890 1.632E-18 +1.900 1.629E-18 +1.910 1.625E-18 +1.920 1.622E-18 +1.930 1.619E-18 +1.940 1.615E-18 +1.950 1.612E-18 +1.960 1.609E-18 +1.970 1.606E-18 +1.980 1.602E-18 +1.990 1.599E-18 +2.000 1.596E-18 +2.010 1.593E-18 +2.020 1.590E-18 +2.030 1.587E-18 +2.040 1.584E-18 +2.050 1.581E-18 +2.060 1.578E-18 +2.070 1.574E-18 +2.080 1.572E-18 +2.090 1.569E-18 +2.100 1.566E-18 +2.110 1.563E-18 +2.120 1.560E-18 +2.130 1.557E-18 +2.140 1.554E-18 +2.150 1.551E-18 +2.160 1.548E-18 +2.170 1.546E-18 +2.180 1.543E-18 +2.190 1.540E-18 +2.200 1.537E-18 +2.210 1.534E-18 +2.220 1.532E-18 +2.230 1.529E-18 +2.240 1.526E-18 +2.250 1.524E-18 +2.260 1.521E-18 +2.270 1.518E-18 +2.280 1.516E-18 +2.290 1.513E-18 +2.300 1.511E-18 +2.310 1.508E-18 +2.320 1.505E-18 +2.330 1.503E-18 +2.340 1.500E-18 +2.350 1.498E-18 +2.360 1.495E-18 +2.370 1.493E-18 +2.380 1.490E-18 +2.390 1.488E-18 +2.400 1.486E-18 +2.410 1.483E-18 +2.420 1.481E-18 +2.430 1.478E-18 +2.440 1.476E-18 +2.450 1.474E-18 +2.460 1.471E-18 +2.470 1.469E-18 +2.480 1.467E-18 +2.490 1.464E-18 +2.500 1.462E-18 +2.510 1.460E-18 +2.520 1.457E-18 +2.530 1.455E-18 +2.540 1.453E-18 +2.550 1.451E-18 +2.560 1.448E-18 +2.570 1.446E-18 +2.580 1.444E-18 +2.590 1.442E-18 +2.600 1.440E-18 +2.610 1.438E-18 +2.620 1.435E-18 +2.630 1.433E-18 +2.640 1.431E-18 +2.650 1.429E-18 +2.660 1.427E-18 +2.670 1.425E-18 +2.680 1.423E-18 +2.690 1.421E-18 +2.700 1.419E-18 +2.710 1.417E-18 +2.720 1.415E-18 +2.730 1.413E-18 +2.740 1.411E-18 +2.750 1.409E-18 +2.760 1.407E-18 +2.770 1.405E-18 +2.780 1.403E-18 +2.790 1.401E-18 +2.800 1.399E-18 +2.810 1.397E-18 +2.820 1.395E-18 +2.830 1.393E-18 +2.840 1.391E-18 +2.850 1.389E-18 +2.860 1.387E-18 +2.870 1.385E-18 +2.880 1.383E-18 +2.890 1.382E-18 +2.900 1.380E-18 +2.910 1.378E-18 +2.920 1.376E-18 +2.930 1.374E-18 +2.940 1.372E-18 +2.950 1.371E-18 +2.960 1.369E-18 +2.970 1.367E-18 +2.980 1.365E-18 +2.990 1.364E-18 +3.000 1.362E-18 diff --git a/doc/PPartiC_UserManual.pdf b/doc/PPartiC_UserManual.pdf new file mode 100644 index 0000000000000000000000000000000000000000..a7b4465235aa3f8fd87c1cd77468ad64192dacd8 GIT binary patch literal 82316 zcma%?L$oN&lBJJr+qP}nwyk??+qP}nwr$(@*r)&Zs;YbN23@mUk)vFhEB1~rB1siQ z#Aq4mSfEJ9uHy=!n3)I|2<(llpm=!bWlZhNoh=BMng0nX(TiEyIGZ{W(2Lm^I-81^ z8rz$g^6^1AIXjvf+CsT!muO7Il5!&MoT(qjc85#Zgym@?U|AL;!sv#B8Uc!`vX@sX zzTU^bdCoQr#c997VJ}DZ{+WauJK`Y9)^~N<+_}&D`F(2r{%SjXjpMh`XVXsPJ;IAl zV7_Zb8QS|?JNth9Je#@3XFrqW#7&i8GKc!y3Uv1~YZG=YGa0ALX1tMe2EFeOfO}ry?IB?DgZ+9=v1xv!Ea{nl-Ar>$DbB zsUtelLF7!r4vj$JE!CclVRgAgV?Py*4*0koJ2Z7-JY8B(Q2|y|5_PvTO{Dr*mSB@W zH8d{KBvu{9Pj8p2j3-`dh=)>3xJpunk^xEpdAdc*0Di;?C<-o-8Az^X`lQc~hybkX#rh;pqd9z0NOo;2I>;%Voev!8 zUW~hbej&!q74x@xeqqn>9ZuyJ)gwQYg*Qh#xff}!$C{C(BJqm{TBAf?dkfbxw2C3H zrQwydK;D;Z=mg5W&xOULood64t+NP+(eiKq;H1h>;!{ z;FITU(Ci3wK;bsIJ(oCJBIj`xE_d^vVLa<8}Xl_>E%*qhd+GL736a-z( z%}ovQ0f9fu-1xu;(-2(WW8lmPbHgjl_#3DB%otoY2;%7 z+7$#d8CSkTWV}|!jVYK)qI)5fN^-3Bi%oRs_3ICf+1QM^Vpd~_QQDvHDR~70+=1$C zuJKg5&U1Y@%}KF3dX8;kR(97JB%&_bvv4=kusan;__*DBh4v>EvKz?i1Od+FhKmF) zdQgzA>y&48!BEfMQWXvv#YCq&0Ti*h*A@Sn54Ou}l)Dk7V{4i)E1@UGg?}095ui|H z`DLfbupsq-vk5X%yr=5&XiLP*q;UA^hcVt#xw94+nn@EwN~0_R6ch2sd5p+?u`F6p znoWv*Cxs@L7larvOmCuAl_1M&otZ3vODK7B2hqhxdt>vr?AoOeE4+=1P^P#BlZ#gj z=iVv%QaPvJy)yGex@LXI1~;Nh-;kthU#Zij7B5_UVGecTn@dJSvkW3u@`I)_sJWg% zI>$5QJ0~;%7$?13eTmV-GgVdR>>6#u(i|6Jm-hOw4mxn4$sz92L?Wx{Mx|A@V%Me1 z9s3>x!!?P826cxyCC#0paoCvj*UPiprrjm5qmcKe@36zJUPJ>~vDk2MOGL)%C|t0f z!f{YcdF+%cpekC<;%T4?A~aS9Kcv!N+;ez|4Uw;Zp}q@j>YaF&S>+USCcGhEdk0m= z>#vCX6??RS&A_yXJni7@rNNt2VG(eD_Xy>r=d^k(KFj9T;GyDRcqI&VwJoe;HA9OS zEGX3~aJ!CaZ;-+9dW`iT``Y}TiQPHc6Y9=j?ZG6d=LE!-&8wfi?=Rej=z)aTi+%d& zi$!yd<)Hp79o>0SDEub86Yx<`uZbfyYUe;}MP?kHvbgFLJ5?3v-ggC|Etjzw1|kZb z_8h=;ybg7Z@%)S+J^xO;{V~HJeAn;I4$?b*t40ecoc5QG_RA}mz3c%ezLtn`C@EMf zg8*sfV7L~2X5a9eY7tPDPyG$*4zU*da(Jp)T;10`nE6&R!lhS~zHYUz5CWe8^?s(D zy#}k*s*sJ;TfKLg%8t9+d{21l3-9$N+UxHY((QU9-jm^HuP^3T+xOjt?e7S2 zA8lyP9cYyJ@zK;O6>aG2ov63D?d(OhyORB81Y+#BC6;Ho?j=AfzFYRuXhhm z|HWYv&aVZ35$u@$dy_=t$=k_E0KQ+O2WhxH-~HE2g|uAmw5-<`0XEmX0}pi39y+iJEp8(#aO3C8}@_AKR{o9 z?T#mX0MJJ}njXiK_5(w)2#bbiS7`R5!UEQo-<%M#%ot+<#5QG6D0tA>5!o!5uWLW_lCc@zezqxBa>2zK^8>EVA%!RRdMU3iC8(-gzc77sbn5!NM%udFnirQw&z_hhP|veJ<;@H%}~&5bDI!5^l^iHqXcHpWrTD?^kM7L}DJgW4O7pwx`ewbzZEe1} zXUViz-lN@&_YT{WBq_zN+eTbf`)_?`v0Zm9!=y1BQ_Fof=OO{Y1IO4r^W-NC@?%?f z0YuUaVr*(MSIN{CGQ1&~Jz9p`>TvHVv4gY$bg~4;jzo@`&P{@|Zv9z7*YN?fgQmpr z0w##r?y1GxF)_a2e#(Hq)B?LX?`n?^(1xWFeBN6n+NUM{ z1SKki=>i`?7Tt%?6Xp+7p-r)%rWQc@&M%@=FH8p=6z)m3| zTtXZ4*whHojIA?8hJz`r!Iv)|4@Bo3zZCUjSCARvBJ@|If23|t#))%FAbPnN-;Ix* zjG|*);JO(Vsf6~5gOX7(L1A-h)O z6e%A{m8ntV?oc?^Gglekoc$L4JOD-DEOE zb{4@J+PTRy`QIi|I{cGCjEfYL7h1!2nNvd(ouclPkpsG;O?b(2;2!$aIAlhs18GDV zK=B6!x@bxA@GV-70UqQoV4CPlL1a+DMxZ|dx&fy-)7yyj!>IurB@OKma;s74$A2t& z$1&syrO4wzHO-~tdu?vc;{oAZR?o0IVlqQrUV zys76=k-7kMlOG!MrrWnP!)YJN-1GlAa`5J!KtM>~!%ks;>!(x}IRmtaXk#CU`mjmN z^fEf>%ANGBk@NINbqo|i7CqX@J6cGpY4m*F$!?3k{+|5)$%~8nIq+$%e*=e&$kja0 zxT?mQ0z*cIG`eAX^KSB(>FBSmUf8RJ-W=)K{<53MlmRy-yfnY^6D|G6t_E7%$cDbm zN2W)WB@uTpv;GZ64MOgIsjOsN)DFuy8hS%_L)m#K{z%zSu1$hO)Ty8)GKxzljQiH} ztfTbfyIm|s$2$FH5|^-;HE&qM2Tz0(ubwfxaCVxpS6oPTQSTv!C)ifTLp^0he_qNM ztpEpd2GY!k(4%dDtU+sDgE&Q39VX*USMy^uZSeGbC#(Y18FE2Dqj-sP)3uk?{x;C! znlzkO&T%`}y&{hrdf`p_J05{?4dj$P8mYuPTw@*Wp_c`8S2gWRAF7NOF};Ink&%8T z!Xd`lO@;EK{}C(OhZ&E;uLMVkKLvglXr^ds%AwifqZR^0$2&?s)S_hW?q4kJQJ2_q z1PJ)Xd1)Ki3AE869u7=%pKs;8#$_`Kj!(d?1m<-ZsD%f>fnJrrKRKd=wHIbRW6Ht@ zXau|>#0kpR%#GA6j&+iH@B53B zSJ0wXp6m(naa0%x5o5Q9$N1GoV5lX{14wN+W(b8 zs1pa{|3;l!w56SJ+7Nr0e-L?QOUEI1;rEbK4;?4vI3bgwQh-U7mN}%QQf0V5UN7tb zg22Fr!jRD6vo}WoVE3+nC1D$-9X`uo0qpKm7aeoAL9X2y{Jcse~CMHJi1KlZJ) zmHD;%BaY@!F>;q}UyNiWJyr6gOimvxTe4m`c>eKv_c~t<)RIpF zXN)q)S#tEuWuy_VInN?f^z_l-CVOz1ExKG7;y-s#e!xiHc+{H@Ko#@0MgnD3z%dTy1 z_`8Q&oL7^8gvQn81i^?~7gWT+w5YWzMqUQ$kRP4B^W><-Z^DfG2r?j#|Mym=-akxx z-7(x9F72%|>%jl9a8l+urMuXPjP?Q~d+-746$Sw~BX;e+m(^hj4w2G~Uc@2s7T*=C zyW}%Y3y25M`-v!f{OJLJI6F`~^zkD+sb6R@k%w*CVOUcK9-poa-1J~OkQm}*V|{he z0|2S?War}mx7J|akp-C#`cyUja-5j1(m-!ctal_uqy*`c~4?&3Gxf&)tIEuO+5PIv-3km&6FC{2KE?N1B zQRWQ7*R%+NW5FXA(izy5%p|k^wju#$Tm9gjvP6gzuxmMN2QvR2VCYScuXp;>}p0KN)j?HLBs4y&CQ zi6lQJ=Xy|Yh!^N3=-!6g{yu)7dx(%+f+J4Lzh6<&KY;?mB+AY%R;{L78bVh54<`zB zM2#RW8rxliVxTu6$N|k{hl(VI8FiGXVx-}f?FePi0r@f^{HiEBRZJTs^0zktH7)8k zq!t&jazLs+T!UXqb+|w0!RCMgq%Bcb)YXAzLy-jv&?yo#U$2xJ5Y2^45d$0^Y$0qf zsH}CbnKUZ0;S=5Eq;?yJ3b5b0dQJ1dcLfRoih7NZy28#{f(&)7r;4+0Sxdl$nu|>X z6kJyEZbGC1Hc!xvDr+iC70%^Vm}Jn>I%Ypj__GYU%*5U*qUw6^eDh`G{2Hq32;#tz zYLGPSD_m_xH<49o+*ivdS86~wFOvy50-N)}8>MU|=my!mW596zKQ$b17pV6Pgh5iw z1ZLIW1X#0OyML~A=+jf)rK<)PDAHWup}}0%QTpPoc#I1$tM8v}Tz8Q{Eh29$iMw9) zf+EN1FJ^fr5-IyIB0mYce(;d)((xxgZ1w%kg8b%&a)Y?!``+zjIPY1OIGCUwRqWu< z7WVvRoIe2g-=mSO^@|3MLe|{BLBX7!BF1FQwYdH9T^)YX^I^(vO<>+Esjc2~Bmydu z;LehcDyFP~)Sm2}ZSM?`Sk1l^1`4bOZ)oU33Ki3Dtm_-5xDimt0VwGVKB5EA?QrUa z5D}wMp_BRJuBY&I2f}wC7zsxoF1YS%@W}}<_#D8bIlHhV8Ym*IkcA8wJMTQO=yU_f z-mc{ZJ%rL-4NQ~V?Sdk}!JZ7RQ%8a3KPi^fVHX{6kql}n4i)qVEnxsurdpxKAnpQd zTvD5qA;?@tHrS({ULwY(BzQmW0FNrdfi2GUTlZPSQu_;RJcLjS1A#G&!jbQUwk&X8 zI)8YZEzdBa=Ms{_s5|BWK|s{40S1UcM|;!dtFtB*fqQnvV57^62EbUXlx64%?XiQg z`eIm_?3bWesPg<hClBfWc zDYFY6WwA&g_<;For(hW1LE^dkEx@Bag>lu;1rNs=V1~xEsrMJ@3`iPIRHo#vJ8Xss zX5_R&Sy(j>XdAQ(LJBckgx@wiVLoJ4g|^g@uCWk%zYwo1UErys*5BC*xnT*ELpD+@yk+a^AcN9SqR;(1%ks+A~9Ov_f2@ zl?URn6xY?x)t-voxkVk(F$gePcx3g>`c6g;G6qwo>v3kV!&@~J=)%y-X7^CR zG^Hpzv+2DvO_XxV>Ic18x& zfu((^g^po1u}c_ePUZa zs!nSRwQ>pwJ4|0G_}4w`*&hweF7~mMMd(msE#nkoTvhu>LkHDXJ)Mkv%D`LI^@0oi zmEl8+yP1qwiD8Q4HxY8NAp=E8)do9#zTA^)0*HpU|RpTxgDpgeoCng`d$dTZ%{f}2%t zR9Wjz8gqNj5l@6aRH`<^zI{iFb1i^-)n5DpU)ojrvfIPMtS2ZX!e@n_;`=aYFyB){ zlPk)EATpX(oK{X+R_RM0*F=Ja3RH>|@F*pcOpLa;VHBjS1;!6~$sq<}1gZjoN2myp z4~AMK*1NosWd2HIr@YdZ6sCp=-dLh3sHi1S73MQQx=T=piKe!7+DmSQ@$*%Q@7vtZ+0=-ga_P zbP`v81*Vp^?Qnwf;pvPA z8gh*I%Moc`0CTwyi;$8ZcVib&oz&;%yywEQ8KP@*DyQMrkgY7U&c#SysGR@hE+=r^ z2rcvKu|i3xAHt#6h;hl-cB6F?BRtxnj*cu5NX%!x~9| zzYJ$7&D+6?z<%LRE{&>Rqg^{Gss|UEM8jF1$uGWkG^qW1$kU*;K}4T`lr`a-CBBT2-Gx06c z+$;hEj62uCnMw!D%I9!9*zU2{8cynLT=vKAsI-OT9me&Uz`wzsTVG5N-rgN^)mf1S z-<$@E2ThY-fR9bhdPk6gN52?`!8_cRsHC93y5%vpZ10$f`?|5P&dL2`OYn~qH11^< zKvg&_WvKES`&&g$17tldaVjL+MyssTa(VtbnQRo0@dv;TOZ)&oF3>IaCfg=+B89p( z?IxD;Q-hKcI*s^-9K;40d)I3BSSQ1*YQmpJ$&M7Y-V6x)r45HL+hJGupTEQpX;C!DfzDC+o?bNV{}Lw9WQyJ+}!B| zcY`zQbn1ZN{d#q=ZdTH=8S((+CDNa+g-n!9Q^v1RoQ^QsgmUkhos%@JVEhO-sf)WE zjO*Aq_)=b0@Gqlrnkd9;zs7u}AQ4)CxgLU6yx==#&xxX7QZ;Y%UxsMc$B(_Ke z%B5}A0O`WAN2Jg4LoVWZg+d>lU@=7q@jxe)yDb+NI|!9CbFQsw_cZ3WLEdsO>+hb0 zA9iCEtplK;#qjNW$*cBNg@YyKPu3em)>_z&W(P6RNtk5ejIXxsuWiNi?`GiC$HKbTuOYoL-Hv1lO z+ISC?0ie${Rf|`5^muOpFfWH(qBrzYV>tKwb>v-4-N~?PW`ZiH)MH#vFM4<47mY0q z)*O>Qt@0P_g*YqeJzxHT~4I5 z@N9NR;02Bh@h5}qgU$tNO2z(O|KQWh{XqHsg8l<1{r>}7{2NF73tO--au!++E*4!V4e*5WbU5BjCESWY zZj>Iq~N0^$_m4Yc7;~Ec3myW83R)t?l}%U)VW0{SGZ{)BWDu zao)3UApWGQ!jW;t@>7poDP9r*Y;F-Re_)%o3d5}|vl z3uUCh?|%7Kv+-F-)rT))RWr5s#1gAhn%sD?!Z2QbOZdRwnDxA|`JF;CbVO^B3)mlU zw;{OnGX>vp2l+8{h~!wli<)sVB)$;^I53EK9(y4NIE(=q-Yhl(Alq{CHz@07PB`mP z=5*)GgaEdk#MAs<$wwqHM|nMCNQ4?tKIaF3I7y7^wg^()RPtLT=rEtY&wXefu$C+9 z6Ue0W#Q*4OeJ%iP&%slMY2`l?mkzzYW*hlAGDv=%IQ8{wO9rf35G%^oUyRhqc>&AW z#({nr^s5;m0%U{Qhk|Rd?b)v`a`f-^=+xTSI!}37zRRRu$0op+O$6l$waA(BpaqEg z9bNNGnF*Br)(;{S4Q56n#nUpp)7s7o05A-v{?sNbU(H1X0m+wRVZDiDdmJ@U8lB13 zf{Gt*x4Kvvf;!I+YxPWrVgT%eLUl37BPOgWu4qTg=KA9GuVx#3-4Kw(tg-!tq$n6r zyhq6lglh1M60M>Lbc{xSi(})Q=>^UYpqH9`p6N6A@EE$n1oziBAm)q!kPA8f3=-(J z)%b(-ZRRd6JE|!;=j93tPRTRFSvxXwgBqxv^N@CkpO8D?cus0ZfQ@4HskEJXrGC@0 zf886x>rjJNiIO}1qRUG{-yNy*T85=DK3NtLMu&m)5IN52S@jcoJN{S+e zsBu?B9QGY9z$K0*N+dzReLF zy5m*sLvcH?OfI!8A*Oi)ql_|}v?`{t-GnV6W%DW3zIX&CCCiz!=>EiGg6MR6HuG}D zp3Y<;f5-g%IR)BXRnY{9YjG|Fa%PU2Is97tq>^2!dyJnwwlDm+D#n2aB5KD%6fyEk z#qOj=jp>S51bIVS`dAfE_`oJA=p@)Un$0^}?v&5=24dF%B2?JFBWcRBLvfa^rWCP? zJR z${f4;dn1EE`ZGwET+Ns4JBfH>ZLc2eOEh?^bZFjGmC0cJmw|ux1!D$0laZwyqsE{DZG; z%zjHBVNGZHT_=-^8;=xkQ@}z8;kEXt47{mOEnS(u-k}jVH#{G~~5eZ+oxU ziCs9ij&a8tfKdDz>nDjr2lrX_bZ22f7$<`BeGk7;T~T}A!dXK%pU{8Wd)^L&O54I7 z_M*Cm@M{?c+k=pboF;K=^>mNrPXwu$LxYyAmJb|U-2U#|j~$iZgfO$g)y1V*sYqGz zetlAiAyzwS!QD+Bm52V~oLHzh|8bA#D@2=x=f|UOc`br6RABip82nu_ZXkx=O@b*i zyr{0Puk-!VbzuKyp|WU|R|-#@+L-=?)b%7Z#!`b(qjYO1T;IY|J1u zv!#GA=1DxmfVW~%(yxQGD0L!k9$Ix&f|&w6XtrZj1XBUWWYx_+^jH*%B}x!-)WQ(? zBAb`twf>w*5UV}`VGCQ%RIhA^p{Lrfk9+_(!w>ITM50p}D>No2 zh1PxAtCO1wx(77vH61;-^4x`$Cf@`@)+DJ{&J%Li^5vsc?kPcz(sOZZ>$%tBrLx+C zUM2%e_%5xW?z~{YB$VFaV+L{G0HqdF8DI&`>m#+5D!64Gz=SkJc*}&3O2j7o{`3;; zUa1uy_nPzyikqTw(&QBZG!)84~DoYSo0?*-DAq|+OIs-j-G z1hjkThUMBnBbxn5FJ{#82UQu0x&OA(WwpZY^z+ClpLkasfN=IQ5@7)vF~56lXXFwF zG6JBFy$ z76V#eHvVsItpnWQJ7K(D8QRPe;r^UAbP}a+d8f0H|PZHB4%hX{2e7M})hB267s=kGKtl0nlD%pECzMB>R?!OB=G)k z^1bg{#h`pA$vPWZZCY_O3RfJ6zVp5ZOf_U7vCMC^UZ*n=rq^l|@}sE)n|9ybQVzPd zJ_B%KP@Jl_`jBO2WcdEd+LsM->Ll-zr_1s&VY6n&08jHTYDcVr7<8Df+rZaMmgy>d z4L=o4*ga_=w|^jif-3NH;letnAO_>>2p~1tq07pT8UxO41n1yZ;J~m4QnmPC@cSw< z8rb1YWq^g54@_ZOx`<3^)It3d%#laSuz>fI#{&OSfHgE$t`SN-JuL}y#@dR(#Xc%r zQ7-nzf?piLVK^WV`)ou9fSE6>wPX|}{_3qjLefJQNOGP6=Z7f6H;TCYs}j45A*QxGR;QC zSO$%&)`%Xkdg5`bbtnB=6EesWqDB0Zf7x88F*zDI!}U+0hI6X*z_xNJ_NK$U*M{gK z-QsO#B@jdeRG{2n&A%x{{J`_bbuDwjPQEig#Tftuiy!2J6sVHKhI?cx2>vP)HIeV_ zmRyYao2FkUnS#IR@Pe*Eezc$RTS2IEHSyfB=&hzta&0wW$gp{R)Pls z8{2_Kn*W#v4tf)7zU7c~HSZaGM%(S>vTtQ>WPnpqOnwUm`5}kE4RAf#SqU`2MP)X@ zIDu@8m2o(MP1UFO6Ge~yp7B}ah1#m~`omrNSa;j};$kC75B8}#=@bL7Op6TU?e%`7 z6O$>z>5!iCe%{)oa@OCFlHatmZ;p*5ZwlMLCrl<5& zigpw-2_g1XlF~EV9JJaUf02&alhcN9lsTq?WUVPq%fyu#no6>zdaZ`mDo)LU|4iT6 z-n*f4CCGAmn|+@#s#T5FU%^VVS-lqLgt9ZrXJbp}YbzbW*LZ^sCON)n#-JwO93el6 z?ip6IPdGRG@3Gr9x3w<5=lhVu=ZR-)eQ1;EPA4dq*5p)SVbUT8y{P~P-|sU}wQhLn$rFZ%?%xR06v)~#7&G;21H$(j zZ_J-k=qNYY>?e0TkyN?L8y+)&+5<=o!s;t-Y=BNr!N-#dmd1?%K$9^ARM}OB9gOp4 zOSeL>SEf|M{RF&4ouIA4G)*si-pL6i|L74x%t^zV;>t`s=}ELe8-^(}eQfF%++CUL zp>}ZE(79^&`VFOrbhw?qcz*OT!8!=w&mq)2N{(t|HOB^Rku{tVQ@%j2_WDP^ye5p5s^y zl()w~r1yuR$G;`YGm(dMPXAJ(7s)2xY2wzOYLkH`ex_?7@x#1jsocqcIBE*FTsXP{ z=`M9wNVGv$f(aDaF6}kin`AlhKG(F>FdG zW6-JI4`*i2MEb6AU#*=~t@1Wx!vfo{HZAuV!{lypB&4_Lw_u$ZZDllz-4{U|NNOdN z!_zIp<)9&baxRpBP`)zytH6`=)oAY3?l8g9lr`h)F18+&cVaNgR@Bl~Ed4DToTq`d zSC5!bN}C3W_Qcf`*mS?#5Lq2bztMEfXUD$lBvq<18vmqCSVq^iBfr`+liHiCO+7=wDjthr{0L~;#02RsSh*paMz{_W-+(ES%S>e_R=!$&e# z(lC|SCOSY~!7L~2LFcqGj$`ok{5$xI=Q{#PUjR%>ln-rT(QQk?ceYUwciPleQ~koX z>o)K^bagra@Wvd~$Lz;Mqcz18{GpIf>B#%U$XZ7BWYtqvx^jx8Lv8dB^VBpFm&DtX zjznZ3IkEI&ncz5){==kd7i9(?R4{w1q`&z3G_pn6wR1`6C~J=}_3_B%MS4t<$eG69yDOD&H91_gcyie;H;CVj{)x;18zmDSvXMN*?0dmPY_y z!>kQg2YD4<>AZR4Vn^|}foD;kV8vQx4@Zt(%94_9%O{#LQoj%R&IPT_g}C{_x&p2% zaH9*yo(_lH&KVA)QHmcPFy$uCDzw5On;0I&^+$2?Hsl$WrBb39J8Z(H^{W+n$su96 ztltZGHEaxN#iy->Zl`V-T5h7*UBWp!w!cXYFtzj%rP64Q{bO9I&WGr!&CFEXqH7Va zsf5GST4s3~f;w~M9Mv!t>>N2omN~Cn?ZrSQMq%@`{KA~VtKUH}3o75vZOizUrO`O~ zIUX2P3}2@HE?@VuQhCW96sLDb@o+_RHOy+tjVb*N^tGG#%4{Vp^MvpfT>`CCZw{uW)T{0TYmpE76?L_Tz7bN1#+)9Cl?85*9L zx22Vly}x&U7d9rV|50xLR~7xQdYhSp_5T7)Mn?AkAefAd9RC}HIm0uNXwvqFr_FhB zBaBOJLpg38of8F4u~7J^-bT?*GOPe@$kMU&Iy2*xM(62yRt5k*;?-Yk{U%-*(P+k+ zCSsd?e=8>aSZx3I!B4-t_x&ku(?+ki4>2o^eQ#IdTNpp1AbdgRw(nv?`9aKMMbIX=-9&<=fLN9{s(V#lfZla`byva7goUx zRXk(W@Byyg)eGFB(YcYr6wi8FOp1efI-uc`$Gwh+hCn5=2my^U*a2-Guzg-{;Lr1= zlfr`Si5-0c4;34-AM_ZERiM5D3n%2w%!_#XNsu@X?p0^l$Kwo~nLF*xORO)TmU<)GX188=BTq}b>U;oNfl z_G5;eLv#<+A0$0{{nGD1`g~5@R{OKKx&T077OPMk^$&pbCbskOHlhB&^x(y+yTYwV%GjOZTMLv0msR8dN)ud4X--V`3*{u`baQxTYh0G&Z}sX z_UUoQS%grPU{?PH486{^17)_fH^DY|OgrN&TT@XU(?BBrkS2X8&O%mR% zA5=;}FRCu0d`E#n7>as5unB2Yj_+{VY%ofnEJj2lGlF=IH;4y9;rK0Qwi4hYP7nm= z;9yiXf*XHO4zK|iea=r0RJc)(2wY+SFM~iwlR~w)AvFWzjO6YbV7}EDLF;u%2-VVgb=NH^ZaJ5F4DJl$s3ILgZpZ=tpuOE_`M+N#jOP z4A6BrdDn;_$mE`B*VqQ^LIKcuVla4UGsIik|o##_Mf%Y6SxAaKcc9^r@>VRQnGT8;r?$s*8b zG&wQ-P-;FLT^5YE$QeKO=q!FXc~WqJ1nE*45eyivBKw zt?gm&B>~J?KFnWSJh?B7PzewFeq22)KPICuIDR?2mSJNCvD_oZ(BUwbSnHtYB4_!W zzW{^9KFgI;P;xS(({c)Vkjtzx`(%J5iEK4*z*yb;GH#CyGu-;(BT2D2))Nm5M0Cf~ z!aNj`t1#{(pcpDyQsoT!l2}{QoYF~?!N8g4))qiK$Y_ZXFFsJC!?pOB1W!UM_eteW z3AL5&Z=h7JW!(zIxnw}P0sMD_FWf==POv`=v!LYzJIo`QG#U95(?0a_2Qa{Q4u!CUrxC7~SHYgZ+2AL;KW$bx;-`}T0 zU2*z*-XCc(-S1CsZ+>~-2kiHGWqmu|j*st4Gw4Gw(}|>xP+)9A+nNx>F;blAAyV-~ zt&vOzf&st?x-=c25eMm5<6Y#Y;I)?J>0NTyX^ZthQVWwaaO44rWMuHca~So7J#sxp zrzPcUd4M%J%|w%U2#7iiX0yH_2Rt8lu^v9tHJ%I5keB1^+ztoE0BLtJD)uNh+O0Vk zJ>TZ*{Q1VJ5*kU`*G*i=X#yl zl?B>U0>js&@curW^AehjgFy8pPa62~2S4FNmI&&MgJ~3%aNIawung)WKdx7rH4+Ll(8>9@G6k-P5FrJumYS0k3O4$`T zxwt1J=HzDoo!%06)sU%ydLkKmnR9B1|NfS+e5Yiil#%HM0;^g4XXsx!1B^Y>DI@>h zz~_Ek&V-6N>$yaU#dnjQ;uAV*xb`Ij6`0xvSL;W=c-R(rE#QO~b3Z+rMM=3L#*`cL z`K}YC!r8`JJ)idg-!VCu1<4!S6uk4n-6Z1gkoGChor|*r0=JnqQE`N~sn}(wjsEcZ z8hM8G>7$Zo+f*CsV4Um)v7dL#K0I5_rMNyZA3-1r8FET*W8i7lM@m7jT6}$X5Ve6} z)1;SDIV0oRbj_LL#T=*CG_am)75`Z4EyyhxM6?fuxtWl$?p9bPXF)fO$usda+XP}gN`5Y~_=s6-K?Fv?aN0o`^{Y39m&VQe3`CwKQYnTZrhf9WbXd?n zp8M6fdaEw_3sCC|kE*GjUFBs(NNf!?zNsm>PSM9Nr6Dfw%#yFqWl*@vg*YOKOQ^e zNY>ecra$L(AWo1#wmsdL&+BP@IQzs2wqo(+ zP^0~%2DT+R@k&rMHPHy$UxK4;mpO$ibaKLZvRc^7LSoPakWHd3hW#4Ix5BDdzSWzw z?Ebj%YT?USNXJDoh2V1{_W%U>vQrW0dG^w*ba*2yMss-u5|;JF!+e)%Ui*W|B=cA= zfmgLcZb!GEf7C!ea{0)O_b$)^U+xJ{QKW24mwXuYqv1S=i+Ev-Y;DrN>S!7ASr~D7 z=#A^3<>P&>f|*_k48X*!aVWf{Y=hZnev|pd{He<3k*@}RyPg-Xu@t#GEnO0~1UnvW zRcF-QvViT&U|_&Hn8^rWIpR9f^`Zn$1>m?najqEg$^stqX39$@K2=-oh#c0aQu*yRk(p1l@W@vs2L&wr~d#YCo@NbC@J)Y_=To6dzW zt~AE%hJ&RYOD$y=9+s}PS6hlZF%C}&Wbx+1dQ(HQ&wY~7qvgJ;D0|-m9Yg!NXN((W z=%pKCHl~YhWFu6xpC74<89rrga#HNEJfvp&eqmf-43j#7&23kLd7?{9!+p8BD6bkL zYa#gm@b)F}P=4S4rIKuA&sIYcVfGn2*>_p9W#4yMqLLEXOIfp&r6`5SzGjy!p%5h^ zN=m6jME^T8WQ5P>^ZR_i|L;rF%(LA4oO{muyzjZ^-sgGpJ{V-)-(MieBXsUqM$`EI zy?PH-t)2y8UfH;66#ndsyzDr3-oSNA%=;t!Q6dAxa#-_g63S`tM776wQkfpwjpAFf z5oGTkZnJ8?XD=u}a`CEELXX08pL%+a-@YACs2lAPEj7S(LcC<+Kp@tGyQ`x*`t9x4 z%@`LeaxdEK@Gz!Kld-=H*{=3&yQp znblonK6s%2SaxNNJ}TM}y}$osXrV-6_Y&#i{0?tbZ&Y@MJeld_VBNO+uk!{RThg$7 zaEMM#Qh?g_*!x>}S{5a_w>9sNRyLG<+O}Y$zCSvGNjPI~T@OQH)Ad_Foj)j;n)Zob zQ}BIq-u8ViwOwThmPgt9DDAz|nOK7+o5v+`8Vzytk*rVrMBH=-159H>zX-DGDi|jp zj`|WKXfjryd|5MSmm}ZK@~n=;`W@a5OIg|TkFsVU*>x*#FeXn6Eq*6H#zV$r0 zb zsQsKVs_Z6BEO9&S{VXN>hlYyfixLYNSmV+!vhuVP2qtz9<`m4P?G^Oek(1(?GZrrO zI_t&f7v;M{Q5cig$7i2DiGjx^RZs;$Ap2g7dD{P42%NG#GdAvWKXYinTqM__Kl*y_ zm|vrxiP#~EEp?~Q3~sps9WsgEN;+9yI3sUivemUcfDGusP#d={#NX@*`@rRB7SZY(T98Ko-S(Ey9!1ozoVTTx}PIke(DA9)%c|0 zlgFzM+n46*I!_sN>)I#fE}Y@#X~ZVWK9TO<*^_*Vg4IOMJ0PS~X+DPEWN1ncHRF4- z=k;zm^zN_GJJq)Be$xG_C@xIV&){|12aM6--Wvsvoo{7(3p3S@)i!@G`}wi+$h_3a zK*X@WK}JO->1~MtpVS-V^o1i&G(WP|YrAsul434w^@yq4I&K2zd;M+KxjR-*_c2^} z)5rJqU4|LFBN{m^dqnDS_wEFlGyGzS=3A~Q8$a;&@ZCe^!qksWY+E=W`u?3;9}mLX zED}EDm43q{q&nBq1b&{bnyyo)2(FaPzx`@caj3?HT1pqkC-6sU@zP)Ysi3D`1V9Q} zcFULrvNbPR4@sO`(m2zjvD^RtO|=irmnZ`l_NgiBZc()|S~z4k2>Ef=Y_TK0OffY6 z7?k$d*71?0XjQIz`#7}otuzMiT!JoyhNiRMjTa3=?QR>Js(5>cETC4UF`{*7Jf=`< zbYI|0R*~Rup?BayyM_V$Z<|E*X?zqZtaD@wmDq>Y+Ar00`iES8p&#lsT0HcC8R;w?_S#t)amew|k=?mtXIQ@- zZbKw++kZ%Y-|l#wnrqUCX7W`*EY_(%GCd{nc#2 zj_%`3eygr-+7W$;GWp4ij)VOk9Q+CnC6!%;2ivw{&WS`+X@9@X-e0qKv~&JRrQVo* zPZfXtQmEgO_Mo8Atzo=VKEE^Kfoh>$CdT9?*7wu=q*^aRS+HwV zJZFMK$kLM8jp5ce+cf2GyQYe#s?GdMB0bOE|6uC9wbeVF23gqvLE9;`dij2h=vE$` zFz=sN%4%+%d~y1LAG==}o_jFC-$5Fe&px|H!*Ej8naM3QaI83G-z-vD zrMg$-RCzY+jn3Es-Ov`=Lhqx7wztpMJQ#EQRQZ`PZp`@M>;jvOdJyd{)WH-DSN{i` zN8CpzoXfi*Tkaucx>F0HrS~$A{V;$H%=t31$?60?r+TKfh3-~N2D8B4)(Ve266v!h z2_stb-=Ado=e+3s0zJ65?hB*xbc_LCO_@nN-;6*txejgRyE*OLZ43@C8oEmFlJS$5 zy3WVvDmHti9G>*84V{nl41lS;8SqcxK^@8uH_ZqZ5U-PiSX~1%{Kx2*?n6cQod&zK zFHUw;aO_M7(OVLKGEOh7Sh6*Y$7=}rLE-Hs5|RQ(bMuEEzl)frCWL0^{~B-^NENxX)#8zq68qb9ASEqNNWM32(^<1E@3JeB>_4O)ZN4nr5Xi)OYxG zIYMtpBTocdvS?j3JjxPO>$w(aiJ`}@0X3U&tkdhsyw(tY{VlzFL~)w923AvpQcma#e0nqy3}&sZ)} zlANr+uNx~SbmkwEKj^kN9C(3ZW{&51j#bLYv?u}p3Yf9BDz#yYA&5VXV1CjBFkNUvK}W}i+k%Qf0dqDoP9XdKeOwuXAxue z%y9HjL`m~~?(eKT0x^3(Oa~ZNSp8U(60709Ze0?Su6EWgFP64Jzmexj>v>n}i;!{U z>-21Ib0C(lq_eNbF@1=a`IXL(IFlz#QII}Ycq%1U#>?OO-Hx4&eWPDTxz1d;qM6zB zd|WCSD-gPG;ZqCyP}@Nk%a(Yk$Z^N_+b-<#=jRvmcx|Pu-*{?UOY8X4=^r~ButshZ zLcw>LmDT;ONxggr9+>7RoNtFHnIoalOwI_Cu9`dc**w6;MFXVb;jNa6i4#R`_&hYu9!q)Z?gQU}os6sxRavVOus6a!1**}j>ZetFqN z6sme*+b@f1TZAUXv?nQh>N+HNN}1X>wG~y2D(JMR{-|{Gd>`oBr_*8H{s09NhtZ$3 zFKrIu)_pB=n)+&;zV+KvEvghnf)(ALqPSxAGAW!b|2a6Me50*PYUHqT;V+8QS9NYS z45)GrsqOpTb+~;nPW4^YEy_XuJZrua0@+6$5rr1e6HyPLgZq6CkQVIKD=E^OzAGEQaBW7Jv!~2sh8w0q` zmic=5Z`s-YBQNR8?(`-uW+n9eG(%25rL^V4^cVL_<0YLq_P(8^SbAM;U^F?N_cmbr z1t-_LyoO50U(N_IZ54VD>TT!jx)XBl$^Ho{E7|ABGaPT7qXLV|TSiY&9gF<<>TFl> zh43C(X|O3W{l98S^;pZ~W`t>A`E4&;6L6 zDtXR-Gs8U?g7Kr+*&o)D`}~y0-1w8=A$F6gy!MWRhQ;MWl@$e2T7g@-6de!8UG7e^ zI!E@ZmHP4H?1~a^=8F5rsn1nW1e;P^G(Fptt@r$@>&Q{ee6d#S$(|mBLOFKP)?nhP zrO6=k2WIEO*oDPdw{LmQDfHL5D)pa`4Y7T_%g&K4Vjqy?dFDb2+t>HfD!Sc9ZHxOE zY{+OHY``~cfvQgyvA;Iurg#?5jmr4+4n}3cDKFttwP0BX?`TjO+azq8HC=x>|m#HEbJw#fk zC>3@3%>}O9Ev2RH+ps6LM0%@#iJaG1=#SsR(fo~?@2u)OGCmhZ>VWfBTZPiPUC&c~ zYK=@0eds)Y<5zQq-!tx$KU-S9zvmQxyyszB@?&nVuZJyU=E@nu-bu-rW$mEm58*#i z6QkuWop`{vwx#D5+zcwk%PJk_oEdZJNE2^An8v3sM06>5a-+_+uEs0$SBoD!;@!oiY+y<{SC~+^=bf9S1!qjF(vgh4e4dCy z6lagv=F>fI>Ac_`mL>fC_Mve;i9K%e$o#@$t<&Ol3yT->pO5Xg+Rf(zal3QkUQ~f& z<^eDLxLj#u&H>-&+VApm{h8ShjIny9YKJ_?H|(>6+n}K=Htt#U9b-p$%VItXu;%_U zMLvv*L+SiDYN$S){=j;PC65sqoiMzM!z5$r!5x)qMbGa+G*4)lJVhSs{b0C$;wob| zlIz3vqvY~*v@B(GP2prmh9~O7uEa@I#?cjqqY)(>k0kqzYqMLokKb{;_7V}$cllI# z-9sk_nTFhBy{PW&JMIlz{4CVB9Uy1ln?hq6f6Xf{#F-<<>Ue%^+-#T@8p7XySGv*t+Es(#EDuDbOr=NrbK6Ht#amALBnLBd zR@`iws3S2N<7DjT6mWZx$7$=SdO?A#z5G;hl=Ln$x@~WRML19_&yQ<;)+tnKLGw0# zqic%L@5*hXp-J-1(0J@GZ^0B3+Ds`HIZSrKB7%$sqZqi);zzpUt;hR3IKy;5y`iU| zZFHSD5=L`8p6u;#BJbhjy8ZbxHW$^tzA0BND->sS-z9jGy-Qb*Q}QGf6@QM)exxGJTp2Tg{L9Q#FWlU_z1vHR?aj zVX!Nxtu!W9vz0&UxQ>a}d9s8g{<*VF(IKBg%&Tdh*EQ2&x$So-xOqnG8gOQMk1;o* zy$9v3=E!>yrcE)Ko{64tOQqlPR8JP^WGA9{h&*X}a!Zrt1#JTn_wF#0ck_{n3M95@ z?uKX?BcyqdYXVF4Q(t35Z=R@AKko_3Iw{4GyzQyCxE}G23rdaYYEl&lVr!dQ&(;Cs5v>;|d~C$Ro!KdC zc5vnHu$gn;+G<0-v?WN(!K^2Oqw5|WnfAN84%%(@wIN?wb*;3ZjDfd;{ZpS?@fE&# zXV1qf-=J)DD1`kM;?Ch<9@*h^i2?In{b>4(y)7RFUdzL8M34&a>!I0q_wi#Zi&;fR zIW&NF}b!tGFxRC^ywDH0{x5t03yw z)?3hPsdwq3n8V-1otFupP55@8Alj^4&qhb9&dLZaE%kK|KSpzbj6NpS(fyZU@;IfB zqX#OjU#HgL(5?cVo?@E>QVLklOpJ0w6J>PJK&3$`dS>T~vM61?(k}u(xeA!;G@s5v zXY7N-T7!fn3k=jCymv@FEvw(eaR{hutG8aeTx&q_N|`1mZ-O&PXl~(1GPgw(O>ZkE zuz-cv^6`$kVgJ6%u_|rX-Zb{8RfbJ};LRn6lq*U_Vw9~}7}KIGpEb~?!P8=B4|qTV zHKMXZrm`FsA`ExA#3hHklYaQU~RedK~_LKP8tccPDT9Nc=4)SiE$J@I`*mg6xg3htUtwf>&_ zaKZh_UGTy4>M0yGHC%fzT;hECM(M(2#%&1Aux%}2TP_(souLwjzcZZ(kP>fzXS(fE z+IdThpZ$uo5t&K~R=;+5SK@r7|Fa6S(f)6ze!dpzShmaU`bADE2sK))){**0VsB#}I8LaWao4an-LaN&iq=TVgB^!oeJk0@kLBVU zXOZDd{3>vlicQ5Q=5f1#Kc_kcXI{_2r@HO*@6Pbavxa&dvyvsD>kl-kkI0BZ`uHC) z)V4;S6djq+BFWU~O1*zrOXqQkJ>=rY2bA^d;gpP(2Pw9{rN};m zsys7EcKf@i3cM~YPt5Hgdu^}YaLD_znwJeH43s!u zxok~_N2BORsb($nbe|obD*CYNWudiUQP}AS(&zRajY(deH)OfExv&qm8|=|`P}!4t z_?xPE#G`g*)mCqdPHFx9w?nWA)TvHky3e*%^}gN(8=tr~MMBED+aFat{eXm(%l{sB z{LR*u`ihQiD)tRYPvf+*Mlq$K?HGCwLxgHbgIRbLC#ivjj_2h&JWll2HHU*;M{aE; zlNdekt2&S$uIv(H{buWDGNPu9=RmHO*Xj9Sb=^WMFwG5fVFR29% zBsP;#Qq!^il#i4+hFOvh&%7hImA{i$+>8GCRk2r8f#@`;r&=kHj`vX_68oY_GQBEZ zjTTZ-(Ck6nr5L~E#@(!ItYd!I|9UqYhXaLDJeh{&iT6+Z4!pJNCMpE!sm+9Lg_3uoE_ZKtDZX8u@BKwQIC{$dhGYW351IZaZ>a%!*nQ=eg{ znVBK$q@Uj2`z+>ybGxrxM3kS)z%e7G_ax(d|a0+jAakXE4@Fea_|jI3iZ1m&x}zJB4&$2V~pN~%ESOBAIx9D0wf`3Drege@kLjN58hs64uL({$R> zP)s*8pOy(+pv63%ysK$%Rn2^)IxpjL!T|RWT08i*75k{sUT^a1yWi(o5=HKW$6dB4 zNWO7X^Vd6va<+Eq!Z#OtUUg{SkK-G9qFov|DmwV+bD6Xgmw)K#+ujR@z`D&IEfk~~S-HhjcpjG61>@yj!>Nl<-{gK|a~8ag|NOU-ZODo!6OsrH95 z8O^so|202hx9+bO3rTX4-WwXuKQOQ0a!+39$qFM4n3xWtWi#nAaJY361}&c?dLG)EcglzcJx9Hnkd;nb3oV0W<8D#x};cDE1xw7Ib7SKo(DI0Tey7WRp>iWGC*5f^uGW4?Z#V-jk}~2|N%wqnrLV%z*gti_ zGt%y)uz*#3TItS9-m(sz_Zk^L@-f$^9m*{|Cm#EK@Q1D#&sclUi;O1=D#uEVGP+Y3 z6JSDyxzAsJpK{mzaigy=`Ro3&SC1}!PG{-pA<<&~(^}p)#gZ#ha|zv_knrsIFN>OfaT|$q z9+N$KzASaaWgl{244w?%gsyZvI8&E%UzFc_*k>nt;T}_WoDxfIZR}T3w~rm4Q}%A_ zoiO&*ETf5q9b@=Px&NHW^uq1~pRpd4DV5GYuP&VH#4?^rFu(oa;I7HtqK%@JtnWr| zyLOoR4oQgwkG4-(A=DG{L z+w{$u1+k+Q5niXZ(pi|hy0rd!-Msai406dPd32ft>H{m0OTk2C*jAN>bZTN?-pOOuHyRXzOOSP!kkur>%1gP z)5sOqY_d65Qyud~6_`2#dAKbnKS@*NoC3UHXkSeDtvoyVMp>90Sr4fMK+6UiP-a1#j ze@~RM%e0E*o;@5MmG9|pCF(MdmNdTPS$ankm{{y%S<$4V^7Ga1r(dWDmw3%V}_a+s>@bM?*!Pbv`W1X29M2 z>(Eg@^v_`Jwi9<{ABB~RUWCaPey)l&RXcB!yesPYlkBAr)hCpg(YAVht@Rf2B=>Kg zmJtewo!yDlJ26cQPjfecpNMLBYNR4~s4wQBSmv?f@p7(XhOx3Rr6c#Yo1T?L-ow0Q zo-<$|ZH|d?HgOR=G1(vc?Cv9n3M8u$J1Z4a-43HOq41QC{TE&sid~~ALd2fxC{gp+ z`gV%@)}@aI2bgq!X-gVR952<}?(wsXWd9SFUnI^l79|NYDKdo~93LNKdE1%|Fr`<+M}1@?V*`5@6h%Sy3-Zg@Sq@R-3* zCL!RzZ|cWxa8p0_AKlb{R;$&80YNz!ITyk&u1Ghcz;=d$M7B|;=EGseJ1-e;nJ7E7Ll&yFghwe}7CZ{E9w zJSwvb7#eNJu|{s2&~MB6692>LDWuGxK37WY&fAjupWoPC$!f|wQc4CWQEG~XoYpTq z{IkE{_T~EtyI$*Bb&-eGMM9NxV#xc>`0Tws*F}EjjK=vWgP8IQ0}qta&d|_`KjY@P zqVSj@T(HJ&YnC;GPah^|Th}bdtE6<+mtH}g>BsdR>MK@8%U#WXVL4nF@X$n-BH>lQ zf%jcFJ`P@I3j2F$@aQ# z#+!aibL)Qnna{!ka|e$vGFCW@>fAIi3R-aAC(-m{XEtZX>GfE^n;3@$G*b9C76doI zg5W>Gf|`b?iX4LS)dC-~xGq$!JG15PI|x*{R*LP-UA-&eCi{O4efJ@4xf10Et3F`G z^~j^RN;L4P5pUiva^~-^_u4I9bpA-~*1I#E$4y@M*gM$kNaZ8PJC9t{q@3q?+5EL| zkBvajsH5TD1SKws%!vNGcQt-K@93BorC~fdGbd&OH8A?A8eSOW#5nnk)?Rntd{11O zS%ExCX4-u4zKmp9WPA&=(%{Gwx5M%}zPh#^{VoRg)Mp|jt~_9J@n#S>i5_Bhz~xN)-e6%heOSk_Z9^XSrLu zctY?~x?s1phpU&nrL_lml97hHtEINJrx6Ui?Fa(XvG(x<<>Y)k6}0h1QNp;}gUDg> zuAn#?SGyKDc)1Zw25f^GJs|kkBH{M|`SJG0J<1vP162OMkdT0q459M*uumTBH*hcE?Zj|uYDq&;vd>YEHMK9KNcB*SoS)+3n38rIBTsB5gvrI zMoI(@TmS#ONXZ${hV%yhjzAIcOK^C^viq;`c*1RjYox^!!~Z{N*Rc{}*{heOL=#Yf zlX4jc2*T@5h=4#?&q>sQe-2YO_Wu{+|l0d#|r6dq6&efKs#5p2CMu2?Kc|!2Lrw0V4VsGUEF(N0p zCVH8y302}5PK%sC8{ix7nBaT4FGS!$1h$tu3W6fE$c-U@*S$cW2_oZF0ap2Z|LxqAOz82>gLX1 zLrjbuCTr~h`fAW46IzPzxQs#l;2-R$o3)FSrKi2C3&aSo0h~yjzXgVg$tqup+aauGy(~6@p5tkk}ns=In-~9+Iu=#^GI2^dU--L zK)~Z^FT+b{Igf`mzI{*x6ajWE2Sdw5rw&3MU?ot|8n;{J27ncO&1&ZME)WGTdn;?a zlGoRbTvxY-yY)%PveOf;L%_xKBGxr* zjkDnusAFdh!P~{s$r=I>VdrYKBEgD*ki-T8*2WnHM$p{veshI#a6A5Nn*+hn5LX+B zi?uhz4Tl17HIMbzAu`s=b@0|&>8%MDB}8~c=2}JTin)-o5Sdc|5mxT@mgY{sFXUyh z+9uSRR;{aZwX0vLQ^U){&fdlLcV)})MM80QU$zk9w}sZqPten{(w2=!pcn)`P+c(? zA!;RR8ZkK$5D&%))}A$ z!rIf@+S&!;;bv`VzmW+sn_YDcWC`v^bk+5op>1jI3Glz6#<81SgqS{%f6@nj$69d8 zIRPuWyMn;W#TtMTVqNwxwfU9wexS@Sn+Qh=9keU z_!~XT%EE(xO`8ZFz3l4*w^+u``sS{eMaIe11t=8a37p>D+RfF&-qY269dy=s26Fj| zc-KP``m|NM;diWUuY-rHi-3)*yQjOk3uqdyllx8eC@7JKuhfihB$Aj`kblw&e#hF& z$hbPY*#i$>=Z`47KQ04*?KPHrWI|K^^mycYkJs{Y`K>P83xwACW0Iw$o*jTA~1r$V-g5D4mD;`Z~3xPJ4Q9)4U3M$q$b3LP|S$o*62(ZGBXkvCm zF4s@E*q@AnT+bK^)~?Rhp6ii?USt2|S_$?irU&GoZiU~m*8a-cN7aG5f=&xy&Dq>@ zeFKQPCduH@nd0v7NN>EAYAu?)`Qn{~&8>WLkO5Hx~32tJ%d;v2bQ znLiteT7?K~eb6NgM$v*uEL2bsfdZrDwGkI;wc916cYf!JKqRT520{Fe&G2JKLR(>4 zC^R{~q}HlK6Av1-P=e$ryvX2G3^@u9C8AST`2{a8zTFVmKaCQ9&w_ve7}wlCivMMM z@gvRUR^W3m_z@;KOcAC;4wH6uw*m{<@mV(nh?r#{M!34+xYQe9d$5H>L509$ZZTLi z6a(@*NH`KIjD%nWF;EN~1x~|r){4P#!Kt>Fg(ohrr)uxw_*<3u?A5;)0wqKc0)vEN z!O#tf6vhjIMM1I1-$JaoBOylKSd2ey|9{iPiWq1(K@1eojSvKbL2N7rj=wh+W0gJ$ z5d3#yAcdfULO>g63=}8@gA#@cqBqb6&bK!fW7U`cjTkG)fMcOTpwmI35KxdG#2^HL zRKLYoLB{5~SmoG%BgC>v(86ef5Ex(*2nLQJmSF`LgxGvzRjl&ze4w!tr7tgu$Rw2!a+uY%Ip6z`!SV{=K(iaX>&r&=?_t zHo!L^Av6{^Es-|Xk1sZ~%)0E+zw#9{7=OcY%z=eN36jCr69*EC>##PKf{4HW1&OYp z6ony30W5+u3>pjFPz>C-U}G^>$E5!WK7i*UaKr%@go4-zjRFD8dg1_7yQvue6&gh& z08`_{z@Tybg+>DZ+(-=EfMR29{AV<}%wH%WVJH&76(J0G49pP0fl%w{Yjaom&p31$ z7$`wOC<;d$0PQFtAt;K-PnL(ygn`D!dibwUC<==q$bbN>g@vF5u^WrADS6>1j{ZGs ztf&GoI9?S1oInZ?c5JAO%?bW0Y&5Z0>pSI5Otu}vju8%w!o*H|Gg$wfQb+W zF%}M(fV0pL6bgJ{*CBIL%37Du{#V3_!kGc*9YUZ>1@=G)K{2cDvC^dxMnD^a6F+76 z?{%?^3ZyWIad1Hf=-GupdKtU16q{lbKV$iyOCf{{IDr%>JaHmH(6OO1HU%d>q5tp2 zSOx}Q2}0lp2hM>dL>sHbxe{@0+KK$0lm~-}jfOS&iHT()wJ=C>3_jfe(?ScA6XvPF zPb@haKO$Wn;Rq5%P+*9HpZV2V9pT`ID4Qk01&Mriv$C5F=Y)uQ=1t1t0(2okL6C-D zI}#Kk>PTYWdwRLqlZSG`!`iPgCGt|=685f-3&H|pW zag>6e@n2Ev23Z8cB7l_>gfZk`@e7zW_&dEL;A<7XJioW9`M>ly1bi*_H!i)I&mrJM zQ4g@%hNUsUb%asEXfPqX)@%s;RLgRMH!F<-+$IFFD1R3H2aQGumF z0AX#ziYbVxZw2=-Ed*Q`%%$R7377~07X%ylw9hINBJk4`MC~KYxUO2}?|cFQT$k9V zmZxl2EBMQN514h_yfhju3`GmWL3awlVnDoyLLmXY|Dj}}>@Low*Ib^s0+0p)PJ>=o z{10#tH=npT5*KfR$2b4!)(Cv+Zxj6ySAYNk5gepFuzyN`&m{hR1sGv4#t{}qu478X zAL+%-G=KqNI~L^T{?H$OMcF^NJQidL;DSi7pytmC*7gvaULN$PSddx!v-Dr+A2%h9 zLk|iU!5}c$b=gq_K8d-B@GGT}V8kd0;?Z?IDq-=>a_J2{Y)wve6SCVN&H%24LXm^d zvfm+v$PxHlI?f4^NU(t;2K)vbY~W|J{*kZZQ>jGFC9I10!__v6TM!^DU*o77mfn;_ z@k@M`uT7j(AWFh+n513qzrgT782CCC4bm(4pcEJd1NtuvF7b`)zwoJ98+%(XcWdZh zMmNItc?jX;NUjv8#ErfyT*;{$qc^H9FIv`1Mu(lnnAH#h%{R8I+;+0fLrk8Ft{l^ITtV!52pQr5H3)L*|V+$h?L=8hCR~rTr3OEDcK0s`|=s58=`CE4!WPWkkUvNG!DjqN(K3@jT z3T*!F8xOPP8gMe;YQn<=sE-&>SfCNyBdu48J%21e2`2jg2)}LkMm(F#-!` zVKz1}8+(`yZu#bl%MiyVzk#Ksp(`z~2!JK6zX7svI@^ppgo%(v>@uJ;22<{9K0wh# zdkCh${f|MT)*wp|zwQV)5uk_FvxKlKXJPRR=5a^BSPTcgHGepJ;(KnKf8tKXX&-k6 z?l%fQP{o350PdGCh!5~*2x7o}Cb*hF#0X{q!3PA0U=NYS{wX67ho-;Pt9ndLLt2%f zUgV#Si;qP&Q7-~Xq+TTAx8veqk8kxF#jVl9-yTN+W1w&=jc^(R%(dC)k6OXmLU;@b zfwDrNoFLE~7!^WjTMWP?qpN2A_|@ZG-Yew)&?&NEmhh85YPZsVgYp z^s8;Kp5qXNrBWMrpYT-$!&h7o*H00N24F-O>2C@nqDj`iyu3zDGBxp;e$muaMM-t z#^TlY=b3mY;RkWJBWQdQ7x!yLBEmNYw+5d8B5<=I3b;KAe1pAf?SspFXl>r%5+z3d zt$lz3Pk6cFmNnbjc@kNPkofu2Oj5dPT2k`BOmZ9axd_54$&K-dSkoGl>cpI-`3%>1~+Z7(H-KLmPufxI_!VkG*~xpdJH%2AMvHQ(!ZE$L+gD?B9^Z# zzSGh0f^^tpwT8sI4y=}j^6J0DYfbjHA4O6ek#0>q&03IpH|<%Z*%B{`os6w-Z%YP= zNewd6>)W8f+!Z5Vx-;8+jE?c+;^AefS8i`1x`hDtyIa;*Gt`h$!gm`R0G1HN>$vK= z*yD5pO!XJzOQL=REwln)av^*GLTvs&e*ecHLm`FNQX!6a|AG~8hnLO%Z{f7uRe`_* z=H_niY)uR+{9M$Zyr8e7p{<9He%5Q5aF5l-T!<#>v@qxurGhXK1R(hI5l)%l7m)(M zN&hIhK$s_6C+IOHeFI!9ueL#U0RaZBn~s$>8--)>Qw(^c<4q0L9fC^!{KJK@tLOZA zlxTYsxGWEM@gwy=13jYB|Isym*F#PS8PWK>HEu5&3=V)DKu8KE^FT7?4}K%ePHxI? zM6-+k7LY*s=1y*Q=H#Fs!4n-h68tdtBF6~Bz@)e{G4bK)me>K-22uE!3U{OD1}Qgb z@JI--HeCRW{DnbCBZ$Bn<2cjS12=mq2nSR29+tSrLVyEc?Yg;}qBSm7AZ+5oG~6*z zm>i}A?jEtXlyb2J_dtMR+Md?V$G~t(kQ}CF@8JPbWcWTGH#`FMfGi$(pac(03MK=S zg~{Po%EQ!P>ae3Q4H)oqm<~)A=IUY%12-$#yIR3Kyb=d&U7zmk_jsQ$c!ijo<*kikNwDnm0rmr?fb3!qRs?H{!qzr;J8Heq1M($8@ z_MR0xCa1b%d+n9;N31EhK8>B*<)0Vw;G*0?p5bH7msQv!MRHVBM4gKUCLTSzaqQmB zoYO;-YEKfv^K-%zSwnlD}9c}t%R_f-Q5LLN=dYr%DZ-Plz*Oh zdi(NIHA;%4+@$Mw1-TA3xOO!*vC_?CP0N2dAXkDG;81dH>%Da2svSC0Au!0s+FC8- zX#2G*HqXCVJQ6ZK7KP9(kxfYxw0bZ2emuh&@t(2q_Ioc0H#*<_n)X=}+&5pPhO51} zBj}XI<1qZ;PT51#J?s@cF4IvNt#L=mQWP~0+!@lWeWrLgkVQrhM^{n6-q+4iAp@%P%kNSJU(Z~$PMeO0lnA2l@ZOp%%&b1xO+qpnJVFaC)a;Nd> z+%K->8IW##LqSE?y)f*4N$hx?qghd+$F%ULcn59Yk=`8;`T_Dl>0(VWC&}$pZc5KB zMZs&OA7HJ!S~C=mn+V4&waun{E@mz?PdjAPlYy)IX>aH581}(S~Kl)u%ii~7_hA?CbC?-`B&XwrO8<) z+2RNKg`aJBXbT&&$_m~#HRc#b70^C+WAh=E>)R6clZ|S4$#*)%)^%$*$rB6mr#7s@ zPeLz3dA5m`rar4BZOItE%CPUE*cjs_={E+4Rg|edTr|dBN*X$NeplE<&ZJXOeL2N8 zpMPYNQ;`}I_y17G7cs4ATuEwv?*Xs;3b$kBE4pguQ zYn)up-d!K|ayzr4MQWy=x9n&^kr}-ew;S?nlsLsA(42ZW*5AP2_P)vGdZkY?H$;9O z(z*Y^99lAV|Ivu)s7EQ|IL}>emTMm^pCS)!6Acv(GA`|Pf_wJ`jRv)EQ$HGh6u>7BstflnN%)SX{M zCej}Fn;6a@@?4rGG`YXP{P$s{>7Cs7bf+k$T`kiwNdDoYk-X2AZ&qjel8VQzKDPlHgHe!ihC$*ADEzv}T3?>8E;r*5mE zqa--1KUoTTP&(x3`d*t28E}B#rwB6bl!??C?~SqvGR-d@g%={8RX2NLk0T!$%!8#mWyk zSIM-tYlV=vH!o^PvIWt7Mt;GL#A|HXcb@xZksX8Kr3XbOIfD9_!y_pA@`r4fKFQP1 zb}X{8yKY^WJ#)fGU6yoTxll|vGhLbazP##K@y@VU-+68Y|4Q$Y_?Am9hkPe!_gUcC z+vxtj&mL?`<8Kugd)`l76?s?^#oAW%bploLRp}yYS(`ltZT+K$;H$cr#FyJUl3YB) zYVuzOohgiotkxJjAz8V-i(l2<%wnLhP2oRQ!P{B~QDO$vJOi|rL! z>2}6Edj2=c-!3q__w(vuy=Mzw4!tgnOQLXY{85Q@wbs(2&6}L7-;W7b>ujs#u)~5E zlqxwra?8KzeJAh2wcd_y`TPfII@;#0hw|(FK9j*O_0XuCkOVI=kdx+giX=HPc&9Qc zZfTn9kg}7s)YX$8quRdfQSW(Pn&PtEugNO%)A8n0Ce8w7H%HG_7o>_fs9& zN>!({h=9)2*yz1fKbs-ans8jx>1eJ}|L#Nme1m$Y_nh_Z{mf>BdE-D)d;FtY*va5H zSymseml;P2+2mC+8enGyjy<0}!mf0q`nvGV?u(RueEB2uCVcX3v`$*Z@cTDA9|zbW zKLm2~6ky2JZ!{G4Ghr_|wAx6x3gw-y(+pqaDLg-3>_a(KZ*C$Z)Qy6@DSycsGjLB4 zycnU3>fmdUNm7CPak+Q=p6{*7Yl_=0ACiB2skx~#;*-@d&9%#H6p4YIr*Gv3u_&;6 z8<3wG*t+!4wfS_cFZPX;Kk4pjQL219G7F6@S{R1_4OTH7-yb3J_mlOCSRIhj-P+oh zikgb6ZuBk=Ktc}%OnzZ=)_5tM(y5jAGyAcLltq8-6wd|Uk1wyeDW0;}pW{Z-QLO&G z_nrW+^eZ1TZC{mVGU$_Jr~O7Mv)bB7;$J9Pcrf)Djgp`KJiL#&<4UmIqM9KxRhY~A z;z;(jkHW|9_1NT3=iN98;lspFT=e_3c)1RG-{3N@>)ZgXl6cu5- zKX&A*@o}8f%pN$)^lQM2`uXlzf$Nfogh_Tc`Ln4vIg(QdoXRsk@B<~qT|%!LUaS*a z#$1^;Fhuux8a1%xVHTDD*+B+>fBG5MZX>=&>j!)IdBA@XhUrf+0sFBgs-9^ba7v4oC7$IxMaoEG`T?ryDth1>duhCH3|kjQ!2Hw-lj^v>7= z4iPNv&l33ie?CxnaD8N~{%+}=>2MawK@Z7f55Kpw=MeB-p@>>DRYg^GQ+CN;(_acY zJ>~71N;?%)J1@UvH$d&$!ffg1_9=3NDYR`wskBAx;(6FOZ=K|sGS&FV?=tY3?cq+Jn65WJt7 z3b}l!BRu-CM_gmaYeq5#k||QZeWtl@Y(4ap4wAc&A8mVkxwi4rS=tN>gX#;gx7;cP zKiqc(nO~(nb%Qe{KB1lI1B0dh;bIzdN7z~W%(AIf(p?81cq-|<4M41|?Nuz9KwAdqY@8J}^4z{=~ z@wDs9?Ul_BR1%V?WnU7S2ajrOKE@S{ilN8zV|R zkfiu2Av$D-%-oZSrH{KrwRBnP;1H)FDq)ROjc{Fq#2K?DFY3Hr-DC4PfeUwSAzS98 zuGsZ(>t0MeuAYL3*OrH$DEw~W>U^`^@<6)uF2nuh?iH_|kTt%t>yC z=iaw31bNuK|8c$Vww$E&`QrRDo>d~x1H-LUY{zdK9L);VKK71;Vyf#IFGF{_F|Ukp zU*`vE`kvWgR&N_Bn#IUMk>*#6=FzTm>a8b?i=+4GzOzoPFY1v~yR_4>VgY$HJM)l@ zU;f1qzwNEi*uuC`&3Ar6Cdpfrzacx{gsbW>CAL)W?|QfWs5A4OmJFI^%{D(K{hG(q z30IF%)@#>La=ZVJ6$ zsLr>hmerC{6(2($_MEeiYNJS>gP`7!so5Ic<;&5yC^|E|l`;pdqj+*Q0PpieQXLBXJ zj4yQE;Z(~+t<-Ttxck*B*-}BT_%p*obBUANrY;Yp<#isqLmz9Z9^<5wLq}4*+wl0W zn*;TCo>tlX(g~!x=^F4f-8JnYH1GP0NWRGGFTC}4eh3G126QUYbX5%=)(o%{X}@)- z%j2Cv_|6PfdSR#h;4Fv7_r7tn%O!;DJ%r2f{u|Dscx5w>R!ywLT zW?8evUMJZmwp8}zOEpAyMIUr*(fToUE9ZRVkJHLuW%r^yzhmsW?bJ>j>8;z=@zvz= z%p9pmi=MZM>7Z<`tHtDOCcoOwp9jyJ6gpCtL0u$6au}}3Fz&`KRvcpGSBOrIrO~^R z9}xOYvE3!ORd?~mi{N3}=bx*^5B2&FhL*r)C(blO?KswidfQuTbr)`k=;Ea@kadF;%HAvm-Ry?vW7FEJoK6S%)JqQ59bmq zwdL#mF7420>-Tj!dM>ymqwKPs*!y|j;}em)BF4 z`$zhBWaWwL&D9%Sll!rxk^cWEd+VS`+AZ6=ad+3o9SVm=8))3!-QC@xac$h)-QC^Y z-QC^!IPaM`_nmuYZp`;zW>#foL`78XXRr0!>p|I9#ePF#p8P>bZ|cIx5K_N{^J7{Y z1uI1N69}LI%l7Ye^)C$Ne^XZl{%@A$ud+(}1*Fn4ei=kxR>qHia2x+wTK`W~{kr(~ zs`}UIzpCojt^ZS5W%(+&|6W%ARY~~({Qsez2myruQc%PI691;6C;^lKx&U2=ud=K6 zm6~I2_!WR<^mX>nin0;12D98)U^P75gGr4GEAKezsL-L z(N``g!1$|xo4DGT7+L{L0cHU6|ET6)`6Pw_>wii&fUiU#L#w}@_SexrpYThr8^z|wBrY>LW00$F0L&JZ_Jpf056TlhZ0&oSm0o)Aj{;U4;&y2~xEfD^? z@c&!?Vf*so|F%G=cXd^fUtIVI`I9COQPwBVw~A8MFtM<_!eO3noF@1OE4grpD;PAU zjAZN<5txsjAS)O|ak`{_-mO&2d|s)%e0EXIvfK%8ZSD2lDfrGU^KxVLGU_sD5`^r> z;aU5F+hXg6vv7~Y<$wqY^75Kt>0m_t$EBui)-R9R)Te~GSWK;=75THQL{vA&&DR4& z0EDfCHuFQbVKemqm)Ufx!{PWjh8M8p{ATg`R{+(`EQ}Zdz`G@XL?N zhsZ1V&F%AjE?yG$8P{g_=tld)g8_$+3>;+ClJGltXD#%deLVaz`Z2%dZhQWNzzLXF z3J)A}{{fJrg!&iPf%w?)u0I1P+Cuba5XLN9vHr7P7` z;Roz__t7H;{lSP#dRsZ^K-W~&ZdfBFMu8;7G0riw1I>W|tdTO^KnC;V?J*7`(Ao^F zVGErDdYt6uF%u&crYNjl3nL!A4gCWq{baB)b|b>t;N&hEBQi$uB+)TBBgWeB?t7BSfZ%<9<+v_GBF+aQ(?U-4{mq@U-K>t1+Ay{4stt%mU8b zSlDP(@sAB!sTK8-xX^X~mdyrr)5lH5*fXF~eeVM!+aRse&lPR=-;ZfiwgGCy)r@WH_3L@n+IOnH@a`TY(gEf^So2Gzx7LFdrT*N zut10)aB-KMb-GfSHXJ-_56Ts-(V*mi*ZXy5DLE*Y$)q_wI9cJKzC0MnT5e_lqpz7D zVPKIup4aB@te%TM*YQg_;r98e#V7pa;`DL|U*5v%R2l`ViEr_L5vux`S|H_Kwx|Qs z%CM$kIkkU$7fQUe3D4)`7ru~xsu+ENx(A%iPq1pz1p2x~fIHn_Qr6g4MNL4+*_=U4 z#t3$)n9tPY+%dC@Wdq3yyL(uFX1nO?a%f`#=oe zu#<1u+f{zroFw$?0~|9ohD6imj!&q%%fWOS`TIUmVA#6LT+RK|YJuLoTSx*J9=@;wom6wHD zq%*+%jCTIyo`w}EXKr9hyJ8Qvqe%g`3(|NwRT%TbMz&?5a`rHfYcZifN0a)U*^{Ey zGtQuYqC0+nT9R}GwtGKg?E^y*SG=*>ltOYNdqMK}3oRj#m%^X(%&1Mn25n{4xT)W> zs(@{5!I;Nejg(t;~taif7?_l&>qBS>#Yab5Qf8{pCnQu5{kH**b1 zNr2Bwl~;wqI%SBFbkUOIJnYy8ykS`M-n;V(46j%3Y3620e{)2#A=N38*^+8Nt%;qT zB{WT}vTcyw4b&&@7X64WBZwUWt%(3%5?WBTi4%1&gJR?4>I^1rd>Q@PxuhDy+m}?J ztiCC&+)jR|^Ze?R>O@O3wFBF!*H*H(^jp%vK~1tPt$8Ony!*!k(+$P<>m&!|q3@Cq z&W#kL#svLbLRIyHQ49DG7bC_f!Fb2-({alg8v`KOp}(Q0eE$>>gV3BKhk*ug=DdZ} zWqXo+SLt_}uPXIhH(-YY!wsaiWx)#G9g+Y$eq=N5s14rX!C}rA;t2{H*}a9Y7)uA8 z^%cogezFRM)@(n~OmV6{G29~_lNAbTZ>;*Q_pW|@Kj~7l#Vbsep;WV>f4&!GhRi+T z$sTQE9YtUr-4V4R72{{PvcojF zDKf_$Tkug8X5lh^E7jRH7E3q<&WN(ZHT^R5mRc3&U^QU0XEIdVR)}Xayh9BR%GD1xmd6t(7s2zDZ*L5uO{5bThA@$`~KpU zz+2>?6^*P)>H4zsP4B14<9myeCmgc?HFIhCkEQO3g*`SZ8(7QOR&C=`zEX_2l_C|#U$G*RjLh>%mHWI=Ui^;7BOmb zdX2&e`)?Nl`iFiWXA~3`GFcDsuku>TUGh#ZQK)3P3!^^n!#1Evnqt8}K^iau3}V@v z*RnToKNtp*hz&^;n#~;=kuG!E6om30ZjM z3#kDZcWzUO@!#!u%gVECjH2H5E6M#HNYE(WPTTC(>?!1%_W(3pOj+44U3g~+Kj>rt zM{;(22pWV23DjQ47Ol*1B&5(Mcj=VJW|!Cq!G-g+k6K4N%E&e9#nGGPFr~8ZEWQVW zWp#O&3~NLkv%|o#moF@LiGyK}6W0MFMw1>}c7%m0zV=J+SjpOu=nKNFLMWga^jIEW z->R?U8Mi=9lpYkdHR*Z?7pVcm@|@2gCBUML_x1JN*5JrHFTkgzo;=7GMOBI{_PQ^< zY4;Z+kTUh&wenW#NqgpumaPQqt_*RY%9UsCUSk-$PmV(a`|rVSbVkg?z|lH+|Fob& z`;ki#TJaYy39~j0`>#E5IXvLRNy8@jo3yExgw_tx*UXRm_*ONGm*_r5*iLRMp%l@R zGu>G8g{+gnywYXcE!<*6R_On{$CiHZ!J+bMPxCb*RPYOsdl4;RrTfqhu*$s|E6egr zuu`g<_`flH=JflrJ=xO04>#6UeLo-J5nBcMLn`oTg29-X-sIYyMtNLr3GDTF6MLG4 zUr9i`bWD#f40UrPZ&Ot8U_XKX`4F?Y(MBzp&*CNU0hhgkz08cAh#}#QrOM*od&S7; zc*9IM)^dpBF!;JXFv=ZNJv*1~ULjV22%$oiVil7Hu&X@-a;AqmS#{ zc`=W%yI0^+AR+y>DvDCFe#Isx6^|f{L}vV&!4=MCDEG0V$vuTv$Dh$npICGnWT242 z0L=nLt^zBXAS>A%A z*@e-6ZqcS)vS zwonXVs%y#{mC=cmZsy3JS+S>RL#$2DYVD%&l)Y_*z@0?RiTfK_B4{I6twBWfOhP&Q zfw9SJZi^T1e|*mVGUrwst!_-%A4ErK(^5!f`c#5^1my@$@|MCy-!eR3Bb5}@u@O)k z)lhH?pk@N69m;3PSCBQbheh+mn5~TUQp}be29G=Q_NEOCfAHDbR#e}@pA$*qo_$l+ zETYv+{8?zJd7X3>nJYRuUaDl~!lP|o2TbrHAm3*yS_JVR=Fq#^**(|~qyxd0WnhE#p0*y4D6H$Sl@ zzm|m)KSed2uF%Apnw%F~b*<*Q)K_!I;7Ok-7ncpMp?RQR##XK;7XoHKYZS2TeG;2S z#SR;xtGF?bA2`1S#FK*7t?mr!_F5TGzNL;f{p&syXoSLjRaWViIn*y=zN0phb?mU$ zaQKo$#+N-xX^OSFpZiRU;zhQoamZs{|CqSFi0qk-ge+mW47h(?!AV_VW}t%>jeH8c zh{st&8QA;I5w__hzTrisjxX5S|Y zlja=kag#rS%w%e_w>a`pNI1-Nm~<;(=K92@y9=BkTc)Eme8Btn%7J?Qrg_HrBA z8*?jmYs$<~nXx+Z@>;%s083Uz0#}!qA4MP6(0{WugM~U~KNVs6W>lRiI-@r*X<82v ze0kOUtaP#k?aUYcvtlS_n1r*I=5+Cf>`Z9iUqa}F!>y3SAABYaxyJOzRYDo6vW#XT zzuTD&F)mz5P?`AyyQAMrK2ZdQ_UCWyzMUy=iOk#`>q`DK(3DTN{_k=E&aAi@x-HooSLmaLbL7H6t(^$G(m#kds*jBbMXI1^kMtvr=9txG7}o z<&^K@k668qk>;yeB}jz}A`Bg#T$yrN^MbPhMxYBLq1{Ea3rZ!9YQYB_BcFN5`4&bg zQmNUaW0EtbwZ7DTXLbQDaqDp9(vgAxmY%~{ZC z`gU1TH2eH6U@h1rrwb(I4@N4|AMk$`x05z%xxcedSg!)FNFiDJ2$ zkGMm-bIpo2uwOexiu$bsxSq5^?faX#?9LFrdwUGYz!vxO%ehGq{kH!@D%)xv7NID~ zbAQT{(qalSfjKT{sV}s zIA2DmUg+m9gyQPBL9x{I7V2IodXre{DWMFV`+>e(WaZ3e%@CPmnOSs4v_lO`6MN2C z?GJaxz+De&t_?*y39NyVpO7}4&JxAVs}itqyq9wP-${Cf8Vz3F*dL@!!{z@LeRxFv;b)K<^& z&reQ-WMC_5z4jj-FJ_smjlxutbsgjqOGdU?5192PY`4eK^Y8lrdgRPol`ei#)FfEf z7d3idXN1N#E=f-j z`eV~j%=5=mC#4j4EDrgm zvTg@n0dI-0w%@few2+XCVBug`0#wb$&>S6z>M3c@_04gm6ob+S0h5)vP|bdwp^;q2 z#;QtZj>r@2njFu|j&D=HHa1g8j})>-yiMn9#0ycx=FMW|z&hvo>hU24F(55k3Z-!~?WImEdYWhb8lMK}uGfam~8BZS;DkKD2-J5FfZhmrdv(UcY zQ{r(YxMXW9ICyKTBzR?|q(tqc*~4U^P-*LoOf(nE{z5NOb%m$ElWNnq#v~_Cz66pS zZE;tGBDZR6Rv}M-@s}Mc$|~WuWjB+IQ`Z$WZuVzJkFS)3OO3U_wAQIZ9ws%_=un{# z9mI=|-X{@b@vPTsA8BB{NQQPa&u6ae38nLV55CC)tNF(?h4k`mrUuU34lsk%boKG{ zT?_Vl;R|#Ll6if7!>Pwls15R9Tnu)zds9^fXG0C1)I|FWm#~hW5^n$b5JcdPDv3-T z9NhhtHUp90MZfbjiwzK~ol%KDn`7lD^j=m$)m!_imUdij9)b~sn`7{eG-!i8Kj*9; zjD`vJAU#2A9`SR^6zDMtF;afC^(^%2o^svJTP;{X@Gg}-dKzXpvhuh7W(>LWL*+Ho z9_yr17GLI#apboKXTr;stgH%kJVZP9FCwzoq_Ip4W_!*lSg+Ssi$yHDK8rgmz0+-& zU#UaWRiHbUD`lpt$(kMCORA7{Jj+-C z?22O!eFz*i7S_oR0%V8rxAmtY6`28 z{B7hGaAl!v@n^&9-V5iukD|Pc(PF*y0dHUoL~NIYV%~26w@i;W>0eYsoajE2fkorE z4WG#Ng$BcRiDLr>~bQ%ZJToluoVs?{g=my4`g#*a~(yDfA%tW#;QN$qc zA;25kPCX06!do^OHfaIUu$g;lV?{JA_DA;WQ%{MrDq{LP=&q2Z2!i;v_;S)njaTXq zd6xll7qtL2I_)2wUKeJS1{J#1R#aP=fVcOZ@gXFpph`bRy!i6cW;tD?0kb+Re2$H}Xkg-|FBfaSZ; zavfje)u%qhQ2?7fO#fN}3|20q);gwcRkO_C8QJDtBJ@_b5`4Zhu}{gl1`=0JZtd#8 zfM$f5`sv%CO7WWu&R*!hLgP>%ELsEB)exhrQ5$S7t-mwyRox8R!@ANn-)e-K|8oEE zh6S>F*M|)Z&&@x#zsBik!VMoV5s`6tBM}b{EGn|2PzwI~xp@3Jk983o z^d8=^%tV!y-4lUD189WD3>xRTq~_hkCEi!(=K(lCx2U?Nee55$anJlcVal`UZI+M0 z=O=7of<5x9d>zDzd5}PZ5oEwBp3m`qY-g0prdh zgOi>JZOaesNx|ZCGMu4LBqV>_n$+gIYXj7>rQ-pyIoN1Jb}XcDnIEaX2Qe6ib{lmy zSC>zIW~1K1hX%H`hUP32IwO6lQT+O~f3d7ArX-7SxS|S@z-^f9IP%#vO+T*QR;q9t z=6z9e+x?@98Pv!|WmoUGOCJT&c}o6{NqD_lUQUHxEJs-)8lS5I>MwH0yo|Y8C{0qK z7qNtI(?Y~?7*K#t>^e~VUk_Hf_opBlFI`kBk=wMkklrMm(*k@{pSN$w`82j7A+s~z zZDKZO%$TF;qwjulF7I-VXI`4FVVRp4LD`V}lQn&ZLi4-K z9-krFBkP|gSfw}X+@IPo(r>0__+FlSrEBZaa>RLggXh2LgNH)^itTDo8Fg1EEAlQk z35U6zQP|rd*28Dl@D=RD!-I^J3p>ysu$1ckPIdti6jE=}8C5f>#p}}~%ib{dR-*lu zE&)N-Gac8jW~0*b;%mBKV^Fu4DD=FR93GIhj01Q{u2=rlvPy0*Ejv?Dju=dpy>O9F zG1P4sWz2(baUSUSU4=IIbbDzugBi;9<1z>sdLu#m&0P{IVx;qad!*O7N)VRkX(Y_^ zk_iHJhXS23dgkuD9ng1wdTK~f2);G_ToCy_uQcw6I!lseimU!=Y{vw3( z5^^{QO%CU{%_&QYZrX@+9Y^PK4}9r}b?!w4m{gm>Ghh;}!&?uF%BrAwGE`hvy_1&6 zn6p8iCBP_DvL1GR$K0vmdF|X-!>s3{5jT>1!=n4Bb)mTlz429yFGiq7ovM0{LPvD)Yp zv~!L&mM*gE0{Z!K==}v{gz{%cb-iOj2sY|n^VOAmanrzyCAxXGM&Iy-*P7dtcE$+v z%ekZK@j-e987Nz`K_D<&sLd(VZKtK_-SezJuGrw*&}wy*RfC6Z9T~QRUqZzEFV*)5 zF^mq*0FHE2ksZj<6>7up+*efn%0-Fhb~jDthm&#oPXzYD zbIAY3ss6`Y!vDsp{`cYde-U6`1OH$1%U_`cUvqc=!TB?N;aUF)`Tqsy{|oYG{qaTb z|I4ZVNsm;OTh|6xn(|IKVR`paJahqr9nuoNEgDBevmjCKhJABfE&xfeyL@`xl>w#_T z!SE8E!5%w(^xI&}1b3YFhUs#AlieuB*=Qow{$TF7;-uoU(*4r8^L1SvdcX)VSB7HA zkeVCVs96<|zCR7^I!54Hh8(tzISXmo9mm#74%cGYGl%4sL}#L4g*U)F&{i zrzgKd+~epLAwGmUUgmmxAh9^+WS?57ij9W;9SIq?TotQqF+dRP^f)R##gTy!N%sXc zViY)C>pQ3S!Pws39#*$36qDXgFNS1tg=pl|DRYM|lvhl1aOfqK!h*LA*PY zR(eAJuYR(V-W?Wdeg58^T&zH(9DTjrcI3N|@1Q}y1Irf#%fi;M_IZsj+WKpaI?q4g zd&`U#=*^QI!Ed zr{RzBTJ%jy(`(l&gF<@gMv@;q;R&{D#kbKR6sCzG6(dii+z4;?QW6u&8Qu$uH?GU`+a;y-OfkzVx#YHyi&H+|UoJTb5mXnvTK8qTLL-Vi9% zC0ZNS4Vi`I;BWSEdB{_;=`%;=8u1X@4{0=t_JzQHj$m!}f+1vX;{g8<^Ki)lH^6u5 zc}9WLJ6>+M!8h@5n1j}I1uY92Pvb-kECR!Zy^gg>6I9ZGHnnvT?H%f7!<&C(K=&RM zC>~ajo0)l$wF0wYMrfXxr}zRP*LrKZYX2x#vmq z)Ml@jqtxvO1LN3rNu%6W5nX1hu?E`wY)4Q2*hCU4H5+jQ$D0Bv@`(fjF<@58(Em*O zUZ)UO#aF84MsHG5g+wir+vA*NH_gD96sz!XGXNXL)v}RYO#8&NCftc@iry=`(|G8L z;n$~m8Zg5tIelsmI*UKMZMRiIP$>82M;q`KM-0A%M%9OQYA4_mM#SX99vIi3aE`9- zAc^JtXm=u<5YB08E<~FQ6|%cy)p&^lpNJdS-bMN|P>TfG2?Lis**H*t-w%o^(lEJ| zh}VyC>%<0pYwuxCNHe0S14j)fxwR^(HgVCgZtGL5sY&8G&B#4C&x0;%cbHP?mj1eS zk`u6rPG05kl_=IWJ%NS8+Rk=(uSZoZ1y=kCObvv$)JL=Ck=Bba zETg#rx%CE-IqXC#`NpMx6QwCK9b@aTnxgh4$7l0xXF^NKp7&bX&uuZt2Bzgg7Vrr| z^8AN6nu@gHzSvbkw8S!5o*)`uUPm2WG7??!Jt|miw86*|Li6qkUQL@+Qc(5YHjq@T zXTETHen%$CA7dNRp^+nn1fhBK@^=?6%z8-;QBjG6KeKH-RA_cgugJsnYbTk!N67o1 z#zxi)4n<|{4JWt75Q$ouSaoYE;0r1JD?J83Me$4Hm^KiDnIcnmu6=pF^gA$-ax}{| ztF(S;5c77YBW;2{&Dcm}*sUheYB? zsA-!XV*KYtB3b6RZ(7ga$k3JF;!vLZH>qQGE+$LdQ`lY(FEAy^qP$bXiHSQiBI$D6 zCKZn=N=gr|P;g@U3-T@0JCjoCA#j~)ap=63Xciho!5nR{)QG5aOdCol3Oy#a-jtyo zuhjF~BA^unfznDwo7%6_n$~xmXTo^8*FA}__=LFALoc_^9umo=F8q31;ezCu`trs* z{%kn(vlDZqV@xh3rAk`v8}nG0QkH<>!W!+E%uEw~*t?l$mMGfdY3}DAm7JM(SGul= zl5M*nOUng|Xbd?a`(ZNz#oDq#vqSj@c;VQ^YlGmktnzHyo^KgneXB5j@7V2BtK2B+ z;E&&%O-evWQ8EDd?QgjnI8`m8BTV+eb-^9IQU5}#G_A;^*=w&W@fO|%SqMA=dwV!^ ztr}DZq1{Jn`LQOz#t4cf>l@KR&C+dc`JKpx7SiVdGlYbi9;t!VbB|I_N5?B2kdI9p~?NO|nPQkw9??6f(k=dK(F?soA+ampTnC#eRvs z!hrl{oboYI+Q#g(w`9sGST;IRD#fc3EGI+5!YV{fG_8hRmk18 ziS(yw*^#SYD+B91#OSL+OpaN2gx2$}VDeQnxx_(B8KeCVT!_dXy*aFgpO|I8uYmdu zOS5szW`|=w9En2;+xj1(VyI|l(fkNd(S59R#z7h!aj^un-fmxAjsT%+b7$}KTMBT? zV9kTm0uMML`0rn+W+LN`l3R0H7%}Q%{-)e~9Ou4PM;7AP)X*P-c>G}diH);73;m94 zpV<-zl1U{wG!HF}G)ua%b_DVxnHrh;4R`EnRUsh|TicsAOuH!YADu`oDn}D16R|he zaoVO9c=U60o7PT8URB%h8BMj@zV+VL2`Eoma8Vaf!i)Y1}Pc8lqBd8!)nwEJNAQo|zWc!h`HXt3uDc$K0@M@AfZsb*`SphO9^t%m7We$3o0X}l7 z#f{p&;t1Iz9%4sIHWp(IajSf~yDW>d46eXBws>hfJ+bg@oH4X> zodW4JuC-(|E#!?I-agj1= zYHid3)VL-F%#rvYvCAvB?0YT*J4!2PX4+`CU%!UBpnfN!1C4Ifr?2AWhC4=({7UYt zbJVm}oireBac?)_#*ZHjYD_<$un3z9wmMnfh^>iqaspnF#i>`q=>hguc{kj2@ZY=4 zztX+_H(lnxXN3I|i2eeqX&Jv#kH0oG{{o`_PoMcUB=q06nP2buTbucH>;KedenoQq zKikZIOsN0gljZ-1jsCCm<^O0u_5Py){Ttr=|6rT{g*E?kCjB3SOaBwP>WxVrHjj{~NEIt^%Q;u(;5%+qFhbw17i8&57tzs7x-*8bSo`F{R~L5UkPSz>_BL90r4MB#W+b_a*iRs_Me$b);Cv5zWmy&v0R zU&q|hV<@Zt{d&;LU|qoG<*uiOl!PZ;lk6yO&k0L%oaC%rm-ooIKgtoSP}klte^&IUEEAdeH-~>kPYi zL>>a3xC<@;E~Y*0TzF}Izyf$m#UO3CEfs9lr-r&RT@zI@fCzv4=)y#@Xx${rTM&Qr z3ynIJUDO&92h*jReI&C-5d`B&MXgIq=j#X6-BzbT4&Kjii&#vSG{zY)q8cRcjP9WzLc85lt^1*neD z=oc`ICLfPMA6Sox4J51fkLiIk_$PlGh&=ASGX|#z!OC!rNEQ{xpxk&4#jk0SMRe;`{TYLzgoL))CqUIKU>DaEl zGT!LEW^9<5*1`1pa$kF_yy30#c?Sdn5sp=5@pGjM9|QsREiwwF-swE5%P*F8=Y}Y` zsos+{QzdH^0~mXXlw1f|U~sfML^ca%#;lGK5nLlr%y#ODP=Oz8?YG?ere%=o68)85 zR@odZJ@#zMb$z!=S(5~ry0i>sNkq8}%={HfoO4Q~@h}KOsBq+*dTcv%tCM+EFsl*P zh-{#`EExuZF3e~9E7pvDcI_mWHeVDss)3YWeQ0|ee1@3ax@dm_wsjX_ixj+WlVN^!@7y27)tUDKAsqSQi z3_WjO%we;PHt?m}3(6;V-H4jJi(kHAJU%Ju&wPR#ThLZ(Gz7TYV66`am@qBq+c$a+ zjYD+Vxh5vjjf3^lhvI}R5=o6V=g0Ymx7`^e(IEz)%TjB$uvq>a4*clnQ)IluRpfN= zdiNtPj2TxI@;W~xs8!sT3l}cAvcB=|ZE!hmkDMDP7jCoaSl{c`3}I#hEqBHHRZEEa zi!7+y{8q(PAa>{5HObl7mSv1|?r{vgsS+9w>e`-4)ZP_ACF3EnZn6lzx|PZJYPz*| zR;oXv`0SBSpY;pBR=6O}_5A^`WW?xJD%o=JX*Y$!JcX|2^+Pu7@Q6WMh_05ez*_SQ z(3Sbw^4)P(8oAw&p9!TXDKQM<+REn$@Lp131NJQW2>_TuSsS-D{6ocO3`k2=5|3b z#(IaL)4m<-H~&7wM)k!sUcUyor4L4yIGG2eHiT1By@ zhG;o04&3qD5nIOoQY!EYj`~*3|92iyyDpo}b{3o>a>6=TNb|AK5H&J9ddAfc2rwOQ z{8z)=_hTLVU%oT%!Ij7FV`!?EhTWzYCjt~g&8F1@6qq@W;#t&zb8`9X)RBH0<3taa z{BH`l5`%9eLR>#CW~Av4TQe1qT>CU@w+zT<$SL%U zLv7yax5cA@pAVmvIp*32tc%o9Po9V|&>E};2bAZ1hIMnQFcY*f@4rs7o&9m!`;oZd z8>6?Il**@~oN&0+M$ex@6{;hh1+3T1+HAy4lv-;>fq~yQ{mt|J)~SC;o_GRsgv!d$-2Lm|BG z{78n9D3dJjBf(2#6+7O6+u%+Y3ySs!`Ei*TS5V?-QYtFTO!V-%Si^~}R3vwSOzWo83Z zO@jb>V-T{4jcOGB96K_&1Q>}ko);>p*y}ba*2+mxsVxmi96QTU>K>XXjX#?UWn%wm zLj^T6WsrVs@TAb*@@)FR!h3jJH}|Ct_Y+zFsdo=;f4+O0y9?_Ql?^Xoa0oiOFt5fk zCQw6oiDTbS-TS6nz3>78Wk^?r!s-Lh>ZDd*RSfl-GrJ=n%P2sA<*@Rr-HA?>JBeO} zVmas}LB^fe#^1-vJ#A_lz|BZo4SH5fvGIUWuiL%F8k2}b)?N?rRgiZBGA37x89UPQ zogJ?sW~RQ7OO)GI{kZlM@3WiH_v=?l^H==2J`9}6X}$^IEdXOlE3fe z$%cWgOuc?58_`T6y>2mnwfku7d)Eg^)lSk*Ql6{(@UXxmb;~}-kL=1u!T*y@m+Ji7 zO#Rq=Ei9c$Ymznj0{m4ysLP<1EAE_EaStKHbzYC#N#_QI4Qf@}%>KlSY>)U8{Ir=n@F7kf)}# zt?8kwxb-(nyL1}*Eckc~;ihztEPFn5vxC}_aVxd%9OZR@!;mg+SM0BA6D~<1_IwO& zwfaR=#DYbL!e1vucYnG9N>K9$037bSiR*quaq&dRG1y53Y(JsMGuWzBHm>co%P8r- zpWSxHH3MJ}8Xqw&qA7mL>UsbPAGracX$H_WG^hpxaNw^IzfZigPRBfw6<=hE($=@O zG}h<_tWu>HJz7Vm3E>V^6*5hIycP4`+U8S-Np&u?3?Zuy5^J@a2{50-5oyuY1| zb}5vgcIsvK>98j$>CiBu%J3wH+Q18nBDR5Kr&CZGQg=Zbu8g4i)hyuw^-?g;$G{WJYvy0 zu2Dw!D!at)w0eGdB2G?taAM^s?^?4kMIaxJC;NOa9!#KUuHNffMvf3jPFPMwr1CHn zCr;0gsFNXoQkLe3`wh4-GL&YF?cYeZij701oic=24v%YU78l3{ zE}<;~X#Iltn~Gbjo9o_@V^6OJ4K?~6%EEGQ)t6hRdel*|Iu9|1^sy3h$fT_DXfyl` z3cGWIjin4P!DT0Ofpy6jBiAe%6PfTG-~wa(Zr+^qf4WVX zbfW0PH*PU7EvURWALOFrLBK5$sU?1FXMMwcYOkyfye7RIWp{%jm+O3tl6yb07xA&pR|Sk*Ny>Eq zekN)r)h39^v_C?Wao304f;1tKeqJNRo8DGl*ISAT(3_27*Wk~^6>iZ3`O*-LQhe6E34M5p1$7A6B( z8`@48qb)B-V2$PL=Qm?Nf@uvPZmEiimYh4;-2s~?le`|-)u%j#Pvs04N?5R%ST5x{ zb)z^IJgTj-&=h_0JsKfc>MHKpBvhIWzM+qlK)eHs#=fpiay&+^1s4!1h>nl9&1!qw zNobv(>0)}`ZD{c3;3cZDJGnwBZr8lHbKzd$%HuJX%{+Z8;pTdKSrkLkDNYRNOD#_50zpXfO(Y*yD3U2hN4^FE`_>Q^`zVyXuq7`!RQ0U z4wVb7v7Q*Yf)S;5n-P84%tb{anP}nr<6HCoXv=P+?dTyeP>kvvK55GI1o;8v+yc75 zex1rm!?MK7(m!2GdO}V7Lq@xLH5u~HCI-YCcl*~#< zN-0<iFS1i(D6ziAD!Tqk@IOE-s?VwiW|Sz&;UI{Q?N+B);C)L z{V_GgvLnm~+K2s6P>?rpV3Yn0Sv|F73gwDZaK>@TI6cVo!CbX1UeBfRtO$L_NB!~F8-pjgkro1o|i1X zOfmkc)5Fbo+Cf^57spZIQAV})8-EMPJKWj4KOL!-=*F>%Ktk}x!lpW^S4+qi*av;X z>KVSeG$v!6KcUpzR_4qYeV|s`SL4#w8!Ro$ZoheJhP_VJUVSUAfHV}!Dp0^-?T z8#R9^Ldq=;@7>dpDd^^Bg?T%q38UeMzwJq8HNnlOVVC1eSsYpBdWqN}X?7}kPB(G0J8;VI0T z33RdY)J}+5K5PaTof(<&x-ThW>hOt_Mez3FkvGn^gc@(COMSGKe((D~AZBI*KwK5f z*?^^n_)13`o24oSl#Q)(EUTHTa>)gkm~8Pf`2%_Nr?D*=PixEwcB;0%W8eXIuRVlD z6CViV*B>@&;fkUfv;mH3!{6$aKi~tOpfPIq$=ai0);?+vcYaa*blc7utg*02*j3c% z<0Tu9Re`#WR@>)T$5=nmc7t!rA2*hnOE~N&f*OH9fzesb^4hhliF0r^WY57Ip@_k_!VL7bbMK zI}dH?k~5mo^t|Me;rq$|{5gWokplV(e%j{1 z@B-->N+${v65Vso$|tMLUd$fT5!*MJZ`F=jcvYTIz2hO+%p-|X?9a^rYGHEKnL=1(ZMBGje|T1v4?MO?EqG*<}88U4IKw*>_$VtK-h-)%r_#}64O< zH=O?MH&Ea5TDNAYmCn_Y!rco2$E9UTx#3-JiAmbXLFn~ zv4#ex($dzJ(wBAT$>zg!0r?9nRLR+zsiC8;ANtlzeTOE~tCAiWKsmUQ*qXaQ8L%!b zpVO})9h6PGlw&EsFC_{vQ!349v<0{Pz;54A?SX;ODK^IBOlM-_ z%mnk~QVT}p%<`io`d61lkEtxHd&t=(OUG4@F11xu*wvHc**Mi>TZR^;Ws@o-<$sml zD5L-@*65iwsFeoG85N`Mj*sJe2O1gg@0lxN3mQ1ph~*Yggen`G_tb%HkKBurmc$n) z-@Y3`h{zPh&rzLNjs%5f>lLdiK%$A!W^d#s$N|kUQpKbQN#>EsO3TU@6&4$lKvOu6 z{9Mmho^4$7SSz#24UcTT>{#R6MaCy0zF>J*(_LeHXZ=U|vJh zhu5zWfheQF|1_W0Y=!~keCzu62I;-IK5tpTwFU5aXC50Xu_a_QsTCAFtHA8o>vPNI z&DOWiMvgr`73;mH(X+9!uv{O!vo4s^edug@&}_nqDGda|$GPn)rw4#Btn5N~=w6;I zuj?z5Wq%A@x7RyjfMTSYYe1`30crW-#@30fk8k{cy?uE+Rp0kFDTF3Ul;M_yGTeEP zp;975hRpN5X71%0D`ZH+)ojd^h)^_;CQ~KK)I_07rNOMsl4qZr%Dt!0_xXIE=l4AR z@VaN6J*>U5z{Nb?|(1KyJPYD3K@B&Dh?#lv$11X{z~ z!>;Zl^_XK`eJM6xJ8Btu+s;q_)nfGo&A4+2LGa$d;s`B{DBH;v?@p;ZK@X3JE-A-t z4bHRLu6Snp3jZg|FGnIidpaY{k}KZ(u1-9;V)8~$ypi)Fk=qZb+60}1Ja7=7y2|*S zA=e*k@WxklmqxUDhR5I5pbDkYpF8cXIT|t{-tIze4NZ@tUw6w`icaUiZ~so*=%037 zk^kw1w`S(%Q-ZFkyWe{2Q*>}!-xZ!Lv5RwgCgQk`2VKs)_f&jxWL9o};>3bUid#y( z@g4D7+e%z5ldu@sYpKDqad$uK4|HB}k0zRb4~nj3yf@jtO~GXHMD5_?>Vq-W_2$K! zvg}MA=E$e#-KY^zwVLL+_oC>|-5BYJGlbGR!QAg6Q5ELbAKmDk65Sl!!Iv*7b_3Ie zap{)!9Frc)>fhs{z0hmh4&B==7wLt?MMdM;4>a#Mw%gMsX-D?tE$G&kKSq#Gy*2a#MYd+vTFU0Wp$dLOgr<^?R5cQiwH;+UL~YM6=>68RyG1-s_;qFRQl$<0wITL>weLoEE=z7XoOaN>; zzAa~*?D#0#TZ&2u>kg?EN1v$@4Hzq znz-pH@2>RfGW6Q6?Uf;uUifR3uDUJnv3DdjCd()2%VP)JhKD_>=cArZ=ndyAuTdxw zI&qHnpv2L3pyl9xO7W0al|=&o-fN=A6{fO$_}^T(T6#6-@V)F$5gnxu1!g?=+^KAA zBy$aHtbHc@aO$jB&>3YNTaM&&2=&e6j+fjb_5?HQMRpfMv55hg{P+eRm@Yf6y^9 zF7uI{$Ewi$3j}lPn6xCvCxb1`_pj#Lyt!M%?RE$gmhZE3;OV2lpxE>lUK1;dHO2b< zGZS)EZaHtTV`Rlq&bFK)*H?!v^qTV5A9?B3D@H%Jp!f65eRo57*Je&wMA!}?ttS0? zo=tTc91yz_CEXJ<*}P&(3$e2Ht?-MRnE}Jy_cLf*C&m0ydM00ta3yzDPk-1{#`&Ob zr&fGQ*t%mgy|kiH-_}a~yW={d0lVE3?@6r*UndY((RJhU`Yz$UG59zSWx3`lE1pc= zhSFNQS=ojy7qSehYs{0ruBW**ty}Vb2hVtQZPuFI5s|t@wXSQzmR;&qBvu`XdBQDv zYdl(SaejwVz_W)tm%Xisw%3dIh#T2CW>> zE?qek8`Ef&(C{{3S@}N8s`_^+JLCF%q>R@tthx052zS`@*k{36)Aw!BHtS>-XO>Ip zXPOr*ND{d#bTes<2l}FilePZRJ|rbsYVlizfBEd$+@T8fH*(QWM2LXGyn0-FDMp;ze)Z3sWuZ%hxO4i zbWAY9~MEHAIEPOxq1O13pHk4a43d-1}VD>-&9uCDnRc4Rxd zoQ{MfidC^N`DmY?OuJU z)wM@9#yO5GyYP~5ws7pFEZN)LpJHVvq(5eD+WIz&Cq4NaoqxBEA!FsiZ$=<>A_D%x z0!6umq@)-hPc1D?qnDX_I^7?CzvY|$=%C3(`bYnUG}~rt=~-tfOA|^!(B#_z($|BP zjL;S(^0b2LiD1lYCohts@b$v|67_XAg8~mF#D{TmnAlp6z_o{oO!RM*2L^BSYnuY$ zOM;{O0YLVzck#h_tVqJ27bU=*`6|S`{G~z)(1HAaKXe}jYQdH0_eO=~JXQ^GVw?du zZQmYjP+XZ?lRt<+e?6R)2tn|$C>%uiuN4l3BhJ^j;Fw?op@2120ov7Ghad0nx47y@ zAnXcOIhAru9tjs>$fWCe|78RkLPMh)7}{ConH9>{Hme0exO#ZELg35M$PuT;VuKLp8Hv& zm+u7!-lRM8Vb}ITzKe|y_1(OFEa-Jg_ulBFwmrw?ZsrxmW@ESF zm0n8MdW5?bnkuC?Zj^sHvw!mI4r{~f&kYltWuISKb7j$EhX|7(*%zL;?p159`oFu_ zxM-jIm)#|$_wTDT8CVIs-|(=#8JL)V|Cz-F9YUO_!j==0cUKx6Gbw&z;PL8Rn!?x* zT(hy2+s!o!ca4oopBrx98*^(zv!PRwsNo4qJ`6+H@Rif1DDc6kLE4Dh#K^GzhZ|ld2}=vQly?S17_H3f zo$QN@SMpH1&;9Xi(GLB%vsPzM?Z2@wl%e1LAgCG}kzLPu!tY`S;_l&w<$Q0lwH}yw zJoL8;>fRE4|9S$ib-3}8{Ec>&TYFX(PP||Fr5{CB%%LrAqP3TfcUP&Uz4`pUVjoqp zH@DxVf3#oycE&wgzvMYSw?!*lm)OtB%z7Vo;eBjJJ1%!;%}Kd5#K}Y7J+`RQMBJC= z4TazLIiHhzp6c7wh*u5}aK013sbd`E>gpX3uoTBDIWXWi_NF!WfakRoTnwe>Jwq|X z{?pplihho!b;r9G9JBB9dVa5=e{Ar5iJtgO+sA>ibs@oDLjM4eqkBOR^hG8s9+keR^=Y(7i8%8r1>0w4D1#HLOp6 zMmFYqSFU_?m3BQgROX@XqUE?7pQ6-nKTazGCzS8D#l(E@3$5H7wO8cRDHWtZ>gm@$ zFTZ-&3EpUumTsE;f;~kL`8uT9d9oaF`it=x#g|J2a5=h&_4W0dYKO~yO^D1OR z$)}qQKdH)fsEb}r?&h_lT$t_|BXHli>oo(jGM5}m)ua0E$aU8=BL{-R5y^+-03PFQ!g90CVbe{=+Z0Mth4R* zCyUA{YFZK$|8QBn@tKQH8a{$u(8)c8ejRkq{XcUS1=VcelS8()(xL9K%rKHEsK>Mu&-3Tqabh8x#seU?f7eadO_C5$iAA+IB26ne3F6AVJX*c)@?VvRU%| z+rgm86;2MdTQad`*4)=*W4W4gqRJ}r;^g9Jud0JjXFps^6dWBKO89ObxKwr+HT7-l z4sDy~5!wV7UEa$`HL*%^r_-c9A$vbW>7Qv0o6+Z37@C)Qb; zl;}|hmSBflo*ufNze?}N*rBKoSX6~*x@<^Pi`VfG*=J(n;%FWX}{*NaLpT8;pvyavA*GF!on_vJ9@e8b@FjY4W~M2Yr1(k zEIxfMV|j8~a-xM{&IS8(dkqckE*03?*yP72SYC;@jq%FP$uLMYKbyWh#v<9mIPUz# zj<___iyiR^Nx2-}DH4%B^$iNL3fc;WLv{Mv3Qsrh*5yo0iAk(4FE2MMFMm^6{;0aT z>w8ATpy`rlJp5#{mQR=;*ErETP5s74bW_G1?+NH6KFSn*SMn}c?B)6FSFei*xWXqh ziNk}B25Ti1%DPuMnb;~fY_{+h&`-1rZ7QU7ef`Hr(7C{ad_H zf77&iDOe`H(?8)Y!N3E5D;mI-0+8GD64ZY-4uCj6>SJF1QM7+jGXn8HYuQf_DW(0a zO2!0#o3CX55kc>tfRykg+&m={Kjo<<H3J{HTH+i$b-wusw$ftusX!)V>t3~f4)vj; z+Ye5Lv>%)}=&hR?-bpOIX;4shVY>a~@guZh+VQ|0+&136H$8Y{$rnzyx2nX2^^9e- zlz=m?Hdc8v8ef}KUsJ5NYA-*RsJSF)y}E6-G@b4Pa;4bH#1q{7IfS7lD?K}kHHFNd ziqx$fT#w;j6PT)+n9@%p-Vi9nF* zEwb8RmcEo&n%JpVy!*mAmt$LPGmDi(jd#XZ*{E@^Ur~5{kFE8VJTK>x*IPVKyay{ond35%w z46hI^!tYL@vwp96{g>Nj9yyC7JMVDKY)33#V<7!0>pggCJ!so@P;7dWiedH^2ZwY-s37}7i=|D7VBLn#8D>&Ejl1oS7T9?tm@{<1XTsZ{w1)swRaVwir! z#@p!L{~tf%pCpms%jQR9pVsP8?Ck6ukal)vkAWYlt**WpI|3X@eEIkd)#!H75DVn& zl}Ve8$_Kp1yQHPxnS2ekxn{CDVD)N+o>ZxZ@tJ|RMe2N?Fatl@8omd&o)UZhguZ`c zfPR-*Y#C9eNGpG=vYI43hY0zN~lGQB^MJtKIl{p-u-aaw<8 z8z$`t&(R=J5suV_!B768fJ_+4f2DwISd_oK42WZJC@erV{{0vffcoKqqnJxkfgAI+ zka_t_3&A%W{}9j+{<9iFxP8pi+}S|;m}uH8iG%+mR2LizW*?`2EQ??W0? zS`hPOdW5W_96&=Q5#$kQ=qUefp5WoEqrDeG&mO=pdfNL!fWZhOs-r8#HwYoC;^ON| zQ&v(62nbNL_oORQeVofOiwKA_Md+Uf)syB2fN2qW z0DaZRi-UuM3y(r`aIEGK(y5)z=@=c9LCLvG@8#b|NRUvUY}&YoSc-j%a!(JveO%;Wu`@@-KXQn zm#Tu}2A7V09$!5!wDHUI;w`rOZ2Tm~#ikopnKQ_j{fiIgM3gY9c>5}wTJPirIb}*@ zP@?S_dl`lddm5SU;vQ(h;04Su3MMiOiqkm9&%Wgve?Pf=Ro1a+s)UX%^0)(o*O@1S zQO9_Y&HMc<4~WxcBq+|@DbeH0PF`c=GRDD^iQ^>#QUw$q8AietNvN@o@xIZ|1)NVu z$0PcTcuXXHOVCyyHvOHBPUbUQWO& znca$t&A6OgkP?}1Dv|r5A^k3ccd$%V8LOx1?kqWdcse{RoNSc z7z)SC5)+U(xFZ0>`6kG`oM(dn{4xTt-2bqtU_yr~{e1Jl86X`)P{=fp&jSFWbAh1$ zPjCemPXUllS#x_#!hh{4zyTxkzf(GK?trf4Yw@bV6Jg04q7|-Q*=XTYka5|1{fhIx zhI?R3 zUG&R+l$qsrRxjz6R^y$|W2fFN>_0VAtbVedG@IN$K0O(-OgqyPp>_6cf~mpjP`lZh z_mxhjYtT)NeYqN^>t7!rk3CQ z1wOWQFZboci*zh6Na@Ph*An=QPW$LR_>_2$3ng;2aJ$ipNaYMVk$}~$>d}lPaM(u{ zcNGs&Q&;DU$!*HXH-34wxJJmtfG=+My~bQOwMB}rS{Ah_o;t#BI85bMOX6HP=w*ad zyp<7V<8&#V>{D3ayx&27B@cE>rb@7df=A&dm&OCqJ)5>p*5z$e3Jl(>%usLGIAp5+ zh$uKallb9DyrF5vxiEin|E=*I(+{}TJ?+)+6#r2=WZWKM_d;@QiQ1#X8jezF4|qok z1Gm1kp)ORtu~WXQwj%M%`C)#Un&U;nxt^u_zjn>8Yx&@i)WFy--?#Z-g4KY<1FA~= z&~6EBYVU!o+!>{%q>kuo#t+0^YT(X>%9sbF9HQ&{so(Z;pRUZttsm`7Rgpl&;H9F(b8_@VEq| zBXIYV*Nz2kDe>*ej_q7|;fL#VVoQ{eNqF?FRkGcbj)mO1*E>1l5ABV)vW(y+Fm?V# zL!GIqlf=Sc<3hJApT*A@_so>vQh4Iuz3ZgCbZho@TykM|EBQFT&zZdSc~0Bv>r@jh zc{JsB#cepqRc6vN$gLDk4$oLM*{NNS1zdhc1RQ>bW2InfH*`3(|W}lB+H*LRtv~AMjTX(GSLEjiRwRS@* z?>z;EuVS|g9qetYk=fPhHlssAyY8Cq6wq zcrI`{cVw*V%~0*az=3z?6DM0S1DGMH=f$m^KCP{CddHqPEtD!8!9G@qSi^a)m(C1-Z=QxGB93?zL~MWf`W| z#X_*hB)vbp<=b;rYeKM1+WNY~lwO#yPo!(bT`B;EJsQyIp~4W|Z6@fy`+W^pdtJym ztJ~izRW2Ovm_A(TxjOoQ%YhRy3b>NIY9k5LmGNALNsm+O7Omd(W8LHP?PoiV?G&{A zM9>~ikA5SN8KrBNemF+}W4yh2r(e;ng~^`<-(=VqsEcI^mzHc^ur*Kgz~L%!g*`&P z9?lwguT997_Zc^uY;0$Z(o`P59@Gg`b-FaMwf+IhKKq*98->27=E`3O3nZ6@j5=@f z_sZQD6Xf2thg$tb$ouw=^-W!FXE<{YmkCHeve@zEQy_hT)DU9u60(&~(R~YLW^3?= zE?3{K%Ok}{mfc)+Xpig1_|=+7R8U0ekI(+Ih1H^wgdYZrobNC3ZL+YXgzGOkou?9a zrQ3ml73AeX@IAl9v36pTd!VgJ%qN0>xqI==u|qdow|z_Qm^8Q2JLxa;rub&c*in6t zpy;O$-alwws5T{Z>>FvZ;x{c*0n<&5P7(VE7t)SD{+?nT)6hE{$PuIxfBvMLeoiB{ z=6IY{^}3G|lLMCbk3G%(G;v>RLln*`3USIK`OSA4e-TL`06~9ALt@S$esC$2?nA)6 zIOfl}SqLR#R{(${tLf+K=R-l5QhnUzm{~kuSw|Op8t{^eI6yFfkc$AwHRKJ!l@G-q zfo7)T`2M6y{yg>KQ_^+yboHepEcs-2P#qPFeC>UF<@jXjehzesBXnSuB3conh?3(2 znC?(C$W{`8fi{|ke6rdAGaD2GAdT~`(cqf-$ZhUE23#{YbbP;9-UQXbGG7z?O@l}T zCXC`A5@-p9WGe<*L*bZTC;usj1vETcIUJx7=4cQs)Vxt4Q0O^+S#ChRq30N9xdC-K zCrw$F8&JRKpE75?0d>xtL?9!U8hGf&obXs~5TF}#!(zEXgl^1n<>cxNlmN^PkQEC1 zbADxM_RbWhJm+dNufbS;LXPwOybp#zntvbS966r=aSWdBGG~?90>mkHj#DN9;*be8 z55~kKK%C*{w?P5~!W`k39-$!Ch%h3gBXg7UrzFJWoHTzCe-0Hn$Jie#G*fdKv4*^rztqOdVm#+j-b#TRo816>T0fs z3Ir+#deS@r8;C$p4==iMpsEbBiZb{Ok(Fc+%p$(-sxoGJ1_)m&)x+J@7lBqJDWVim z7`&9yaXOa>Ut%OD^fD#Jy zATMA0Km{+l6zGr1h^yJ}YF>0jCQpiBhoauB}kR23ID2jsu68tHkNk|72fq=v~U>tsRN1MkcOMl=7=nw-s zQ~^Jg4Sc9%KSzp>8rk2KMn`N$=(*|`Q_K*%40I8@T&Q&4pX?w=ifEjQ5=#v}W&=8u zARsI zsDTt8ABxk@4w<#)3jIG9mpx5S?>`I7(-X1|sGKtJ=9{ACRL0rcR@IfJ_zhm`tFnK*zyNUnuBNW18_c z58t0ZJX|5ORQ9mXA`yy9pg)Ay9F0IALBSvk4Nb&DAtfsfPauHUgM~)G zp>e>8v(SL=fMQ5i8j%3SBCIqL27(!7qamR{kd4MXYL<Z`?Bf@DQq-P%oc#8pwo>}XH;1mSCY&0~QU48@{bLYZV4hsd}EHrR|VwV*X zjRwz8Sj%BR+|E7@EZEVo*Toazx`9OkDqDLXSqfngv+x4?fOci`%*f>Pcr<*-QR9*c#BBjKT)Bnu6L0?5YfG!g+eUhu3F1?P_l zZD3g21I5_&6r6ws*Hf^9u=9*Tqu??{gVlt+4>T6Mg~8S@8V9p|3=lMIY~W;Q_BqF( zNiaRYV308T#bD4F_PNAhFnIRyVlX(kyhW2V~bxCkJV@a^}6zDFRy${IE!NvjQ&l zW>8q#12mZafzL0XR1aG@u!ynIK<)~*UZF_L3<-N(*t&~F}w*1fQ9K5puzM7tU~blB>REO1)n1V9?IIWi~|EyoL%;49Ecv-%As*Y z=v5K+x-eVCfS0Ra^FbiN*B9`O4a^4+k+8KBME6ksm8D;v4v&KlpJn9* zYzJX-#sitM>mLD$hV2aqNFo4DVQmk1Jec1npm8w&O8_pOofixV57Tih5=Ua^1qk73AA461iVvTP3PQ=q zbw35#$Rm^tsZ?JCuyTYULP^`piHZQNK+n*$cOWcPH1XR|SmZXsHUd#|I}%GEf-^cW mni%wUB8ddJ#jE}AOE8n~bm##aGkriHf(H4drFZIU^8GJI@rRfI literal 0 HcmV?d00001 diff --git a/doc/logos/PPartiC.eps b/doc/logos/PPartiC.eps new file mode 100644 index 0000000..5517180 --- /dev/null +++ b/doc/logos/PPartiC.eps @@ -0,0 +1,177 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%Creator: cairo 1.16.0 (https://cairographics.org) +%%CreationDate: Wed Oct 07 12:01:13 2020 +%%Pages: 1 +%%DocumentData: Clean7Bit +%%LanguageLevel: 2 +%%BoundingBox: 6 7 205 57 +%%EndComments +%%BeginProlog +50 dict begin +/q { gsave } bind def +/Q { grestore } bind def +/cm { 6 array astore concat } bind def +/w { setlinewidth } bind def +/J { setlinecap } bind def +/j { setlinejoin } bind def +/M { setmiterlimit } bind def +/d { setdash } bind def +/m { moveto } bind def +/l { lineto } bind def +/c { curveto } bind def +/h { closepath } bind def +/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto + 0 exch rlineto 0 rlineto closepath } bind def +/S { stroke } bind def +/f { fill } bind def +/f* { eofill } bind def +/n { newpath } bind def +/W { clip } bind def +/W* { eoclip } bind def +/BT { } bind def +/ET { } bind def +/BDC { mark 3 1 roll /BDC pdfmark } bind def +/EMC { mark /EMC pdfmark } bind def +/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def +/Tj { show currentpoint cairo_store_point } bind def +/TJ { + { + dup + type /stringtype eq + { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse + } forall + currentpoint cairo_store_point +} bind def +/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore + cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def +/Tf { pop /cairo_font exch def /cairo_font_matrix where + { pop cairo_selectfont } if } bind def +/Td { matrix translate cairo_font_matrix matrix concatmatrix dup + /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point + /cairo_font where { pop cairo_selectfont } if } bind def +/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def + cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def +/g { setgray } bind def +/rg { setrgbcolor } bind def +/d1 { setcachedevice } bind def +/cairo_data_source { + CairoDataIndex CairoData length lt + { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def } + { () } ifelse +} def +/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def +/cairo_image { image cairo_flush_ascii85_file } def +/cairo_imagemask { imagemask cairo_flush_ascii85_file } def +%%EndProlog +%%BeginSetup +%%EndSetup +%%Page: 1 1 +%%BeginPageSetup +%%PageBoundingBox: 6 7 205 57 +%%EndPageSetup +q 6 7 199 50 rectclip +1 0 0 -1 0 71 cm q +0 g +49.766 27.305 m 49.766 29.898 49.078 31.813 47.703 33.039 c 46.328 34.27 + 44.297 34.883 41.609 34.883 c 36.625 34.883 l 36.625 47.945 l 32.438 47.945 + l 32.438 14.352 l 41.516 14.352 l 44.234 14.352 46.285 14.953 47.672 16.148 + c 49.066 17.348 49.766 19.242 49.766 21.836 c h +45.594 21.398 m 45.594 20.18 45.313 19.301 44.75 18.758 c 44.195 18.219 + 43.313 17.945 42.094 17.945 c 36.625 17.945 l 36.625 31.289 l 42.094 31.289 + l 43.313 31.289 44.195 31.008 44.75 30.445 c 45.313 29.883 45.594 29 45.594 + 27.789 c h +45.594 21.398 m f +75.641 27.305 m 75.641 29.898 74.953 31.813 73.578 33.039 c 72.203 34.27 + 70.172 34.883 67.484 34.883 c 62.5 34.883 l 62.5 47.945 l 58.313 47.945 + l 58.313 14.352 l 67.391 14.352 l 70.109 14.352 72.16 14.953 73.547 16.148 + c 74.941 17.348 75.641 19.242 75.641 21.836 c h +71.469 21.398 m 71.469 20.18 71.188 19.301 70.625 18.758 c 70.07 18.219 + 69.188 17.945 67.969 17.945 c 62.5 17.945 l 62.5 31.289 l 67.969 31.289 + l 69.188 31.289 70.07 31.008 70.625 30.445 c 71.188 29.883 71.469 29 71.469 + 27.789 c h +71.469 21.398 m f +99.891 47.945 m 99.473 47.945 99.016 47.82 98.516 47.57 c 98.023 47.313 + 97.602 46.957 97.25 46.508 c 96.508 47.469 95.535 47.945 94.328 47.945 +c 90.047 47.945 l 87.617 47.945 85.852 47.461 84.75 46.492 c 83.645 45.516 + 83.094 43.836 83.094 41.461 c 83.094 40.508 l 83.094 35.551 85.426 33.07 + 90.094 33.07 c 96.047 33.07 l 96.047 30.664 l 96.047 29.676 95.742 28.91 + 95.141 28.367 c 94.535 27.816 93.641 27.539 92.453 27.539 c 85.016 27.539 + l 85.016 23.945 l 91.969 23.945 l 94.75 23.945 96.801 24.555 98.125 25.773 + c 99.457 26.984 100.125 28.938 100.125 31.633 c 100.125 43.148 l 100.125 + 43.629 100.305 43.988 100.672 44.227 c 101.047 44.469 101.664 44.586 102.531 + 44.586 c 102.531 47.945 l h +93.844 44.352 m 94.801 44.352 95.406 44.129 95.656 43.68 c 95.914 43.234 + 96.047 42.656 96.047 41.945 c 96.047 36.664 l 90.047 36.664 l 89.211 36.664 + 88.523 36.941 87.984 37.492 c 87.441 38.035 87.172 38.719 87.172 39.539 + c 87.172 41.945 l 87.172 42.781 87.363 43.391 87.75 43.773 c 88.133 44.16 + 88.742 44.352 89.578 44.352 c h +93.844 44.352 m f +114.047 44.352 m 114.047 27.539 l 110.547 27.539 l 110.547 23.945 l 117.891 + 23.945 l 117.891 25.242 l 118.953 24.379 120.234 23.945 121.734 23.945 +c 126.672 23.945 l 126.672 27.789 l 121.25 27.789 l 120.227 27.789 119.461 + 28.117 118.953 28.773 c 118.441 29.43 118.172 30.094 118.141 30.758 c 118.141 + 44.352 l 123.359 44.352 l 123.359 47.945 l 110.547 47.945 l 110.547 44.352 + l h +114.047 44.352 m f +146.359 47.945 m 144.723 47.945 143.473 47.422 142.609 46.367 c 141.754 + 45.305 141.328 43.91 141.328 42.18 c 141.328 27.539 l 136.766 27.539 l +136.766 23.945 l 141.328 23.945 l 141.328 18.18 l 145.406 18.18 l 145.406 + 23.945 l 152.406 23.945 l 152.406 27.539 l 145.406 27.539 l 145.406 42.43 + l 145.406 43.129 145.578 43.625 145.922 43.914 c 146.273 44.207 146.82 +44.352 147.563 44.352 c 152.406 44.352 l 152.406 47.945 l h +146.359 47.945 m f +172.719 14.352 m 172.719 18.664 l 168.156 18.664 l 168.156 14.352 l h +161.344 27.539 m 161.344 23.945 l 172.469 23.945 l 172.469 44.352 l 179.578 + 44.352 l 179.578 47.945 l 168.391 47.945 l 168.391 27.539 l h +161.344 27.539 m f +197.25 47.945 m 194.469 47.945 192.332 47.27 190.844 45.914 c 189.352 44.551 + 188.609 42.508 188.609 39.789 c 188.609 22.508 l 188.609 19.789 189.352 + 17.754 190.844 16.398 c 192.332 15.035 194.469 14.352 197.25 14.352 c 204.063 + 14.352 l 204.063 17.992 l 196.766 17.992 l 195.578 17.992 194.613 18.348 + 193.875 19.055 c 193.145 19.754 192.781 20.598 192.781 21.586 c 192.781 + 40.695 l 192.781 41.688 193.145 42.535 193.875 43.242 c 194.613 43.953 +195.578 44.305 196.766 44.305 c 204.063 44.305 l 204.063 47.945 l h +197.25 47.945 m f +0.831373 0 0 rg +17.277 53.707 m 17.277 56.629 14.91 58.996 11.988 58.996 c 9.066 58.996 + 6.695 56.629 6.695 53.707 c 6.695 50.785 9.066 48.414 11.988 48.414 c 14.91 + 48.414 17.277 50.785 17.277 53.707 c f +0 0 0.501961 rg +23.906 61.84 m 23.906 62.969 22.992 63.883 21.863 63.883 c 20.734 63.883 + 19.82 62.969 19.82 61.84 c 19.82 60.715 20.734 59.801 21.863 59.801 c 22.992 + 59.801 23.906 60.715 23.906 61.84 c f +0.596078 g +0.751181 w +0 J +0 j +[] 0.0 d +4 M q 1 0 0 1 0 0 cm +16.941 50.758 m 29.465 47.813 l S Q +27.848 47.258 m 30.453 47.578 l 28.262 49.023 l 28.523 48.41 28.355 47.699 + 27.848 47.258 c h +27.848 47.258 m f* +0.137102 w +1 j +q -1 0.235294 -0.235294 -1 0 0 cm +-15.851 -50.987 m -18.248 -51.872 l -15.849 -52.753 l -16.234 -52.23 -16.233 + -51.519 -15.851 -50.987 c h +-15.851 -50.987 m S Q +0.6 g +0.751178 w +0 j +q 1 0 0 1 0 0 cm +23.234 59.613 m 30.398 49.301 l S Q +28.824 49.969 m 30.973 48.461 l 30.313 51 l 30.098 50.375 29.496 49.957 + 28.824 49.969 c h +28.824 49.969 m f* +0.115667 w +1 j +q -0.694805 1 -1 -0.694805 0 0 cm +20.193 -42.855 m 18.17 -43.597 l 20.191 -44.342 l 19.87 -43.904 19.87 -43.302 + 20.193 -42.855 c h +20.193 -42.855 m S Q +Q Q +showpage +%%Trailer +end +%%EOF diff --git a/doc/logos/PPartiC.svg b/doc/logos/PPartiC.svg new file mode 100644 index 0000000..9966134 --- /dev/null +++ b/doc/logos/PPartiC.svg @@ -0,0 +1,159 @@ + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + PPartiC + + + + + + diff --git a/makefile b/makefile new file mode 100644 index 0000000..749dad2 --- /dev/null +++ b/makefile @@ -0,0 +1,34 @@ +# compiler +FC := gfortran + +TOPDIR = $(PWD)# top directory +MODDIR := $(TOPDIR)/mod# module folder +OBJDIR := $(TOPDIR)/obj# object folder +SRCDIR := $(TOPDIR)/src# source folder +JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build +JSONLIB := $(JSONDIR)/lib/libjsonfortran.a +JSONINC := $(JSONDIR)/include + +INCDIR := +INCDIR += $(JSONINC) + +# compile flags +FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall + +#Output file +OUTPUT = DSMC_Neutrals + +export + +$(OUTPUT): src.o + +src.o: + @mkdir -p $(MODDIR) + @mkdir -p $(OBJDIR) + $(MAKE) -C src all + +clean: + rm -f $(OUTPUT) + rm -f $(MODDIR)/*.mod + rm -f $(OBJDIR)/*.o + diff --git a/src/DSMC_Neutrals.f95 b/src/DSMC_Neutrals.f95 new file mode 100644 index 0000000..ab537cd --- /dev/null +++ b/src/DSMC_Neutrals.f95 @@ -0,0 +1,110 @@ +! PPartiC neutral DSMC main program. +PROGRAM DSMC_Neutrals + USE moduleMesh + USE moduleCompTime + USE moduleSolver + USE moduleSpecies + USE moduleInject + USE moduleInput + USE moduleErrors + USE moduleList + USE moduleOutput + USE moduleCaseParam + IMPLICIT NONE + + ! t = time step; n = number of particle; e = number of element; i = number of inject + INTEGER:: t, n, e, i + CHARACTER(200):: arg1 + CHARACTER(:), ALLOCATABLE:: inputFile + REAL(8), EXTERNAL::omp_get_wtime, omp_get_wtick + + !Starts threads for OpenMP parallelization + !TODO: make this an input parameter + CALL OMP_SET_NUM_THREADS(8) + + !Gets the input file + CALL getarg(1, arg1) + inputFile = TRIM(arg1) + !If no input file is declared, a critical error is called + IF (inputFile == "") CALL criticalError("No input file provided", "DSMC_Neutrals") + + !Reads the json configuration file + CALL readConfig(inputFile) + + t = 0 + !Output initial state + CALL mesh%printOutput(t) + CALL printTime(t, .TRUE.) + CALL mesh%printColl(t) + PRINT *, "t/tmax: ", t, "/", tmax + PRINT *, "Particles: ", n_part_old + + !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) + DO t = 1, tmax + tStep = omp_get_wtime() + !Insert new particles + !$OMP SINGLE + tPush = omp_get_wtime() + DO i=1, nInject + CALL inject(i)%addParticles() + END DO + !$OMP END SINGLE NOWAIT + + !$OMP DO + DO n=1, n_part_old + CALL push(part_old(n)) + + END DO + !$OMP END DO + + !$OMP SINGLE + tPush = omp_get_wtime() - tPush + tReset = omp_get_wtime() + + !Reset particles + CALL resetParticles() + + tReset = omp_get_wtime() - tReset + tColl = omp_get_wtime() + !$OMP END SINGLE + + !Collisions + !$OMP DO + DO e=1, mesh%numVols + mesh%vols(e)%obj%nColl = 0 !Reset number of collisions + CALL mesh%vols(e)%obj%collision() + END DO + !$OMP END DO + + !$OMP SINGLE + tColl = omp_get_wtime() - tColl + + !Weight + tWeight = omp_get_wtime() + !$OMP END SINGLE + + CALL scatterGrid() + + !$OMP SINGLE + tWeight = omp_get_wtime() - tWeight + tStep = omp_get_wtime() - tStep + !Output data + counterOutput=counterOutput+1 + IF (counterOutput>=triggerOutput .OR. t == tmax) THEN + counterOutput=0 + CALL mesh%printOutput(t) + CALL printTime(t) + CALL mesh%printColl(t) + PRINT *, "t/tmax: ", t, "/", tmax + PRINT *, "Particles: ", n_part_old + WRITE (*, "(2(5X,A26))") "total time (ms)", "avg t/particle (ns)" + WRITE (*, "(2(F31.3))") 1.D3*tStep, 1.D9*(tPush+tReset+tColl+tWeight)/DBLE(n_part_old) + + END IF + !$OMP END SINGLE + + END DO + !$OMP END PARALLEL + +END PROGRAM DSMC_Neutrals + diff --git a/src/makefile b/src/makefile new file mode 100644 index 0000000..ecc2529 --- /dev/null +++ b/src/makefile @@ -0,0 +1,17 @@ +OBJECTS = $(OBJDIR)/$(OUTPUT).o \ + $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \ + $(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \ + $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ + $(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \ + $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ + $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o + + +all: $(OUTPUT) + $(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(JSONLIB) + +$(OUTPUT): modules.o $(OUTPUT).f95 + $(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95 + +modules.o: + $(MAKE) -C modules all diff --git a/src/modules/makefile b/src/modules/makefile new file mode 100644 index 0000000..5b59144 --- /dev/null +++ b/src/modules/makefile @@ -0,0 +1,49 @@ + +OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\ + moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\ + moduleMeshCylRead.o moduleOutput.o moduleInput.o \ + moduleSolver.o moduleCollisions.o moduleTable.o + +all: $(OBJS) + +moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleMeshCyl.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleInput.o: moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleMesh.o moduleMeshCylRead.o moduleErrors.o moduleSpecies.o moduleInput.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleInject.o: moduleSpecies.o moduleSolver.o moduleMesh.o moduleMeshCyl.o moduleInject.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleList.o: moduleErrors.o moduleList.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleMesh.o: moduleOutput.o moduleList.o moduleSpecies.o moduleMesh.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleMeshCyl.o: moduleRefParam.o moduleCollisions.o moduleOutput.o moduleMesh.o moduleMeshCyl.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleMeshCylRead.o: moduleBoundary.o moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleRefParam.o: moduleConstParam.o moduleRefParam.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleSpecies.o: moduleCaseParam.o moduleList.o moduleSpecies.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +moduleTable.o: moduleErrors.o moduleTable.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ + +%.o: %.f95 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/moduleBoundary.f95 b/src/modules/moduleBoundary.f95 new file mode 100644 index 0000000..a2e579c --- /dev/null +++ b/src/modules/moduleBoundary.f95 @@ -0,0 +1,36 @@ +MODULE moduleBoundary + + TYPE, PUBLIC:: boundaryGeneric + INTEGER:: id = 0 + CHARACTER(:), ALLOCATABLE:: name + INTEGER:: physicalSurface = 0 + CHARACTER(:), ALLOCATABLE:: boundaryType !TODO: substitute for extended types + + END TYPE boundaryGeneric + + TYPE:: boundaryCont + CLASS(boundaryGeneric), ALLOCATABLE:: obj + + END TYPE boundaryCont + + + INTEGER:: nBoundary = 0 + TYPE(boundaryCont), ALLOCATABLE:: boundary(:) + + CONTAINS + FUNCTION getBoundaryId(physicalSurface) RESULT(id) + IMPLICIT NONE + + INTEGER:: physicalSurface + INTEGER:: id + INTEGER:: i + + id = 0 + DO i = 1, nBoundary + IF (physicalSurface == boundary(i)%obj%physicalSurface) id = boundary(i)%obj%id + + END DO + + END FUNCTION getBoundaryId + +END MODULE moduleBoundary diff --git a/src/modules/moduleCaseParam.f95 b/src/modules/moduleCaseParam.f95 new file mode 100644 index 0000000..a70ce6b --- /dev/null +++ b/src/modules/moduleCaseParam.f95 @@ -0,0 +1,8 @@ +!Problems of the case +MODULE moduleCaseParam + !Maximum number of iterations and number of species + INTEGER:: tmax + REAL(8):: tau + +END MODULE moduleCaseParam + diff --git a/src/modules/moduleCollisions.f95 b/src/modules/moduleCollisions.f95 new file mode 100644 index 0000000..f59084b --- /dev/null +++ b/src/modules/moduleCollisions.f95 @@ -0,0 +1,152 @@ +MODULE moduleCollisions + USE moduleTable + + TYPE, ABSTRACT:: collisionBinary + TYPE(table1D):: crossSec + CONTAINS + PROCEDURE(initBinary_interface), PASS, DEFERRED:: init + PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide + + END TYPE collisionBinary + + ABSTRACT INTERFACE + SUBROUTINE initBinary_interface(self, crossSectionFilename) + IMPORT:: collisionBinary + CLASS(collisionBinary), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename + + END SUBROUTINE + + SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrel, part_i, part_j) + USE moduleSpecies + IMPORT:: collisionBinary + + CLASS(collisionBinary), INTENT(in):: self + REAL(8), INTENT(in):: sigmaVrelMax + REAL(8), INTENT(out):: sigmaVrel + TYPE(particle), INTENT(inout):: part_i, part_j + + END SUBROUTINE + + END INTERFACE + + !Container for collisions + TYPE:: collisionCont + CLASS(collisionBinary), ALLOCATABLE:: obj + + END TYPE collisionCont + + TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic + CONTAINS + PROCEDURE, PASS:: init => initBinaryElastic + PROCEDURE, PASS:: collide => collideBinaryElastic + + END TYPE collisionBinaryElastic + + TYPE:: interactionsBinary + INTEGER:: amount + TYPE(collisionCont), ALLOCATABLE:: collisions(:) + CONTAINS + PROCEDURE, PASS:: init => initInteractionBinary + + END TYPE interactionsBinary + + !Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory + TYPE(interactionsBinary), ALLOCATABLE:: interactionMatrix(:) + !Folder for collision cross section tables + CHARACTER(:), ALLOCATABLE:: pathCollisions + + CONTAINS + SUBROUTINE initInteractionMatrix(interactionMatrix) + USE moduleSpecies + IMPLICIT NONE + + TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:) + INTEGER:: nInteractions + + nInteractions = (nSpecies*(nSpecies+1))/2 + ALLOCATE(interactionMatrix(1:nInteractions)) + + END SUBROUTINE initInteractionMatrix + + FUNCTION interactionIndex(i,j) RESULT(k) + + INTEGER:: i, j + INTEGER:: p + INTEGER:: k + + k = i + j + p = (k + ABS(i - j))/2 + k = k + (p*(p-3))/2 + + END FUNCTION interactionIndex + + SUBROUTINE initInteractionBinary(self, amount) + IMPLICIT NONE + + CLASS(interactionsBinary), INTENT(inout):: self + INTEGER, INTENT(in):: amount + INTEGER:: k + + self%amount = amount + + ALLOCATE(self%collisions(1:self%amount)) + DO k= 1, self%amount + ALLOCATE(collisionBinaryElastic:: self%collisions(k)%obj) + + END DO + + END SUBROUTINE initInteractionBinary + + SUBROUTINE initBinaryElastic(self, crossSectionFilename) + USE moduleTable + USE moduleRefParam + USE moduleConstParam + IMPLICIT NONE + + CLASS(collisionBinaryElastic), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename + + !Reads data from file + CALL self%crossSec%init(crossSectionFilename) + + !Convert to no-dimensional units + CALL self%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) + + END SUBROUTINE + + SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrel, part_i, part_j) + USE moduleConstParam + USE moduleSpecies + USE moduleTable + IMPLICIT NONE + + CLASS(collisionBinaryElastic), INTENT(in):: self + REAL(8), INTENT(in):: sigmaVrelMax + REAL(8), INTENT(out):: sigmaVrel + TYPE(particle), INTENT(inout):: part_i, part_j + REAL(8):: vRel !relative velocity + REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities + REAL(8):: alpha + + !v relative (in units of [m][L]^2[s]^-2 + vRel = species(1)%obj%m*NORM2(part_i%v - part_j%v)**2 + sigmaVrel = self%crossSec%get(vRel)*vRel + IF (sigmaVrel/sigmaVrelMax > RAND()) THEN + !Applies the collision + v_i = NORM2(part_i%v) + v_j = NORM2(part_j%v) + vp_j = (v_i + v_j)*(1.D0+DSQRT(3.D0))/2.D0 + vp_i = (v_i + v_j)*(DSQRT(3.D0)-1.D0)/2.D0 + alpha = PI*RAND() + part_i%v(1) = v_i*DCOS(alpha) + part_i%v(2) = v_i*DSIN(alpha) + alpha = PI*RAND() + part_j%v(1) = v_j*DCOS(alpha) + part_j%v(2) = v_j*DSIN(alpha) + + END IF + + END SUBROUTINE + +END MODULE moduleCollisions diff --git a/src/modules/moduleCompTime.f95 b/src/modules/moduleCompTime.f95 new file mode 100644 index 0000000..a45a07d --- /dev/null +++ b/src/modules/moduleCompTime.f95 @@ -0,0 +1,10 @@ +!Information to calculate computation time +MODULE moduleCompTime + REAL(8):: tStep=0.D0 + REAL(8):: tPush=0.D0 + REAL(8):: tReset=0.D0 + REAL(8):: tColl=0.D0 + REAL(8):: tWeight=0.D0 + +END MODULE moduleCompTime + diff --git a/src/modules/moduleConstParam.f95 b/src/modules/moduleConstParam.f95 new file mode 100644 index 0000000..a335a53 --- /dev/null +++ b/src/modules/moduleConstParam.f95 @@ -0,0 +1,14 @@ +!Physical and mathematical constants +MODULE moduleConstParam + IMPLICIT NONE + + PUBLIC + + REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi + REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1 + REAL(8), PARAMETER:: eV2J = 1.60217662D-19 !Electron volt to Joule (elementary charge) + REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI + ! REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 (Not used in Neutrals) + +END MODULE moduleConstParam + diff --git a/src/modules/moduleErrors.f95 b/src/modules/moduleErrors.f95 new file mode 100644 index 0000000..0a95db8 --- /dev/null +++ b/src/modules/moduleErrors.f95 @@ -0,0 +1,43 @@ +!moduleErrors: Manages the different type of errors the program can produce. +! By errors we understand critical errors (that stop the program), +! warnings (that only output a message with a WARNING tag), +! o verbose outputs that can be used for debugging porpouse. +MODULE moduleErrors + CONTAINS + SUBROUTINE criticalError(msg, pgr) + IMPLICIT NONE + + CHARACTER(*), INTENT(in):: msg, pgr + CHARACTER(:), ALLOCATABLE:: errorMsg + + errorMsg = "CRITICAL error in " // pgr // " with message:" // NEW_LINE('A') // msg + + ERROR STOP errorMsg + + END SUBROUTINE criticalError + + SUBROUTINE warningError(msg) + IMPLICIT NONE + + CHARACTER(*), INTENT(in):: msg + CHARACTER(:), ALLOCATABLE:: errorMsg + + errorMsg = "WARNING: " // msg + + WRITE (*, '(A)') errorMsg + + END SUBROUTINE warningError + + SUBROUTINE verboseError(msg) + IMPLICIT NONE + + CHARACTER(*), INTENT(in):: msg + CHARACTER(:), ALLOCATABLE:: errorMsg + + errorMsg = msg + + WRITE (*, '(A)') errorMsg + + END SUBROUTINE verboseError + +END MODULE moduleErrors diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 new file mode 100644 index 0000000..9caf61b --- /dev/null +++ b/src/modules/moduleInject.f95 @@ -0,0 +1,155 @@ +!injection of particles +MODULE moduleInject + + TYPE:: injectGeneric + INTEGER:: id + CHARACTER(:), ALLOCATABLE:: name + REAL(8):: v !Velocity (module) + REAL(8):: T(1:3) !Temperature + REAL(8):: n(1:3) !Direction of injection + INTEGER:: nParticles !Number of particles to introduce each time step + INTEGER:: pt !Species of injection + INTEGER, ALLOCATABLE:: edges(:) !Array with edges + CONTAINS + PROCEDURE, PASS:: init => initInject + PROCEDURE, PASS:: addParticles => addParticlesMaxwellian + + END TYPE injectGeneric + + INTEGER:: nInject + TYPE(injectGeneric), ALLOCATABLE:: inject(:) + + CONTAINS + SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface) + USE moduleMesh + USE moduleRefParam + USE moduleConstParam + USE moduleSpecies + IMPLICIT NONE + + CLASS(injectGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: i + REAL(8), INTENT(in):: v, n(1:3), T(1:3) + INTEGER, INTENT(in):: pt, physicalSurface + REAL(8), INTENT(in):: flow + INTEGER:: nEdges, e, et + INTEGER:: phSurface(1:mesh%numEdges) + + self%id = i + self%v = v/v_ref + self%n = n + self%T = T/T_ref + self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight) + self%pt = pt + + !Gets the edge elements from which particles are injected + !TODO: Improve this A LOT + DO e = 1, mesh%numEdges + phSurface(e) = mesh%edges(e)%obj%physicalSurface + + END DO + + nEdges = COUNT(phSurface == physicalSurface) + ALLOCATE(inject(i)%edges(1:nEdges)) + et = 0 + DO e=1, mesh%numEdges + IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN + et = et + 1 + self%edges(et) = mesh%edges(e)%obj%n + + END IF + + END DO + nPartInj = nPartInj + self%nParticles + + END SUBROUTINE + + !Random velocity from maxwellian distribution + REAL(8) FUNCTION vBC(u, vth) + USE moduleConstParam, ONLY: PI + REAL(8), INTENT (in):: u, vth + REAL(8):: x, y + vBC=0.D0 + + x=RAND() + y=RAND() + + vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + + END FUNCTION vBC + + SUBROUTINE addParticlesMaxwellian(self) + USE moduleSpecies + USE moduleSolver + USE moduleMesh + USE moduleMeshCyl + IMPLICIT NONE + + CLASS(injectGeneric), INTENT(in):: self + INTEGER:: randomEdge + REAL(8):: randomPos + REAL(8):: vVec(1:3), vTh(1:3) + !Edge nodes + INTEGER:: n1 = 0, n2 = 0 + !Edge nodes coordinates + REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0 + INTEGER:: n + + vVec = self%v * self%n + vTh = DSQRT(self%T/species(self%pt)%obj%m) + + !Insert particles + DO n = 1, self%nParticles + !Select edge randomly from which inject particle + randomEdge = self%edges(INT(RAND()*(SIZE(self%edges)-1)+1)) + + !Get coordinates of edge nodes + SELECT TYPE(edge => mesh%edges(randomEdge)%obj) + CLASS IS (meshEdgeCyl) + n1 = edge%n1%n + n2 = edge%n2%n + + END SELECT + SELECT TYPE(node => mesh%nodes(n1)%obj) + TYPE IS (meshNodeCyl) + p1(1) = node%z + p1(2) = node%r + + END SELECT + SELECT TYPE(node => mesh%nodes(n2)%obj) + TYPE IS (meshNodeCyl) + p2(1) = node%z + p2(2) = node%r + + END SELECT + + !Volume associated to the edge: + IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN + part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n + + ELSE + part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n + + END IF + + part_inj(n)%pt = self%pt + part_inj(n)%n_in = .TRUE. + part_inj(n)%v = (/ vBC(vVec(1), vTh(1)), & + vBC(vVec(2), vTh(2)), & + vBC(vVec(3), vTh(3)) /) + + !Random position in edge + !TODO: Use edge coordinates and transformations for this process + randomPos = RAND() + part_inj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1)) + part_inj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2)) + part_inj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3)) + + !Push new particle + CALL push(part_inj(n)) + + END DO + + END SUBROUTINE addParticlesMaxwellian + +END MODULE moduleInject diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95 new file mode 100644 index 0000000..28ea16f --- /dev/null +++ b/src/modules/moduleInput.f95 @@ -0,0 +1,366 @@ +! moduleInput: Reads JSON configuration file +MODULE moduleInput + USE json_module + IMPLICIT NONE + + CONTAINS + !Main routine to read the input JSON file + SUBROUTINE readConfig(inputFile) + USE json_module + USE moduleErrors + USE moduleBoundary + USE moduleInject + IMPLICIT NONE + + CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile + TYPE(json_file):: config + + !Initialize the json file variable + CALL config%initialize() + + !Loads the config file + CALL config%load(filename = inputFile) + + !Reads reference parameters + CALL readReference(config) + + !Reads case parameters + CALL readCase(config) + + !Reads output parameters + CALL readOutput(config) + + !Read species + CALL readSpecies(config) + + !Read interactions between species + CALL readInteractions(config) + + !Read boundaries + CALL readBoundary(config) + + !Read Geometry + CALL readGeometry(config) + + !Read injection of particles + CALL readInject(config) + + END SUBROUTINE readConfig + + !Reads the reference parameters + SUBROUTINE readReference(config) + USE moduleRefParam + USE moduleConstParam + USE moduleErrors + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + CHARACTER(:), ALLOCATABLE:: object + + object = 'reference' + !Mandatory parameters that define the case and computes derived parameters + CALL config%get(object // '.density', n_ref, found) + IF (.NOT. found) CALL criticalError('Reference density not found','readReference') + + CALL config%get(object // '.mass', m_ref, found) + IF (.NOT. found) CALL criticalError('Reference mass not found','readReference') + + CALL config%get(object // '.temperature', T_ref, found) + IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference') + + CALL config%get(object // '.radius', r_ref, found) + IF (.NOT. found) CALL criticalError('Reference radius not found','readReference') + + !Derived parameters + sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section + L_ref = 1.D0/(sigma_ref*n_ref) !mean free path + ti_ref = L_ref/v_ref !reference time + Vol_ref = L_ref**3 !reference volume + v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity + + END SUBROUTINE readReference + + !Reads the specific case parameters + SUBROUTINE readCase(config) + USE moduleRefParam + USE moduleErrors + USE moduleCaseParam + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + CHARACTER(:), ALLOCATABLE:: object + REAL(8):: time !simulation time in [t] + + object = 'case' + !Time parameters + CALL config%get(object // '.tau', tau, found) + IF (.NOT. found) CALL criticalError('Required parameter tau not found','readCase') + + CALL config%get(object // '.time', time, found) + IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase') + + !Convert simulation time to number of iterations + tmax = INT(time/(ti_ref*tau)) + + END SUBROUTINE readCase + + !Reads configuration for the output files + SUBROUTINE readOutput(config) + USE moduleErrors + USE moduleOutput + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + CHARACTER(:), ALLOCATABLE:: object + CHARACTER(8) :: date_now='' + CHARACTER(10) :: time_now='' + + object = 'output' + CALL config%get(object // '.path', path, found) + CALL config%get(object // '.trigger', triggerOutput, found) + IF (.NOT. found) THEN + triggerOutput = 100 + CALL warningError('Using default trigger for output file of 100 iterations') + + END IF + + !Creates output folder + !TODO: Add option for custon name output_folder + CALL DATE_AND_TIME(date_now, time_now) + folder = date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' & + // time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6) + CALL SYSTEM('mkdir ' // path // folder ) + + CALL config%get(object // '.cpuTime', timeOutput, found) + CALL config%get(object // '.numColl', collOutput, found) + + END SUBROUTINE readOutput + + !Reads information about the case species + SUBROUTINE readSpecies(config) + USE moduleSpecies + USE moduleErrors + USE moduleRefParam + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + CHARACTER(2):: iString + CHARACTER(:), ALLOCATABLE:: object + CHARACTER(:), ALLOCATABLE:: speciesType + REAL(8):: mass + LOGICAL:: found + INTEGER:: i + + !Gets the number of species + CALL config%info('species', n_children = nSpecies) + !Zero species means critical error + IF (nSpecies == 0) CALL criticalError("No species found", "configRead") + + ALLOCATE(species(1:nSpecies)) + + !Reads information of individual species + DO i = 1, nSpecies + WRITE(iString, '(I2)') i + object = 'species(' // TRIM(iString) // ')' + CALL config%get(object // '.type', speciesType, found) + + SELECT CASE(speciesType) + CASE ("neutral") + ALLOCATE(species(i)%obj, source=speciesNeutral()) + !TODO: move to subroutine + CALL config%get(object // '.name', species(i)%obj%name, found) + CALL config%get(object // '.mass', mass, found) + CALL config%get(object // '.weight', species(i)%obj%weight, found) + species(i)%obj%pt = i + species(i)%obj%m = mass/m_ref + + CASE DEFAULT + CALL warningError("Species " // speciesType // " not supported yet") + + END SELECT + + END DO + + !Set number of particles to 0 for init state + !TODO: In a future, this should include the particles from init states + n_part_old = 0 + + END SUBROUTINE readSpecies + + !Reads information about interactions between species + SUBROUTINE readInteractions(config) + USE moduleSpecies + USE moduleCollisions + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + CHARACTER(2):: iString, kString + CHARACTER(:), ALLOCATABLE:: object + CHARACTER(:), ALLOCATABLE:: species_i, species_j + CHARACTER(:), ALLOCATABLE:: crossSecFile + CHARACTER(:), ALLOCATABLE:: crossSecFilePath + INTEGER:: nInteractions, nCollisions + INTEGER:: i, k, ij + INTEGER:: pt_i, pt_j + + CALL initInteractionMatrix(interactionMatrix) + + CALL config%get('interactions.folderCollisions', pathCollisions) + + CALL config%info('interactions.collisions', n_children = nInteractions) + DO i = 1, nInteractions + WRITE(iString, '(I2)') i + object = 'interactions.collisions(' // TRIM(iString) // ')' + CALL config%get(object // '.species_i', species_i) + pt_i = speciesName2Index(species_i) + CALL config%get(object // '.species_j', species_j) + pt_j = speciesName2Index(species_j) + CALL config%info(object // '.crossSections', n_children = nCollisions) + ij = interactionIndex(pt_i,pt_j) + CALL interactionMatrix(ij)%init(nCollisions) + + DO k = 1, nCollisions + WRITE (kString, '(I2)') k + CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile) + crossSecFilePath = pathCollisions // crossSecFile + CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath) + + END DO + + END DO + + + END SUBROUTINE readInteractions + + !Reads boundary conditions for the mesh + SUBROUTINE readBoundary(config) + USE moduleBoundary + USE moduleErrors + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + CHARACTER(2):: istring + CHARACTER(:), ALLOCATABLE:: object + LOGICAL:: found + INTEGER:: i + + CALL config%info('boundary', n_children = nBoundary) + ALLOCATE(boundary(1:nBoundary)) + DO i = 1, nBoundary + WRITE(istring, '(i2)') i + object = 'boundary(' // trim(istring) // ')' + + ALLOCATE(boundaryGeneric:: boundary(i)%obj) + + CALL config%get(object // '.type', boundary(i)%obj%boundaryType, found) + CALL config%get(object // '.physicalSurface', boundary(i)%obj%physicalSurface, found) + boundary(i)%obj%id = i + + END DO + + END SUBROUTINE readBoundary + + !Read the geometry (mesh) for the case + SUBROUTINE readGeometry(config) + USE moduleMesh + USE moduleMeshCylRead + USE moduleErrors + USE moduleOutput + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + CHARACTER(:), ALLOCATABLE:: geometryType, meshType, meshFile + CHARACTER(:), ALLOCATABLE:: fullPath + + !Selects the type of geometry. + CALL config%get('geometry.type', geometryType, found) + SELECT CASE(geometryType) + CASE ("2DCyl") + !Creates a 2D cylindrical mesh + ALLOCATE(meshCylGeneric:: mesh) + + !Gets the type of mesh + CALL config%get('geometry.meshType', meshType, found) + SELECT CASE(meshType) + CASE ("gmsh") + !Gets the gmsh file + CALL config%get('geometry.meshFile', meshFile, found) + + CASE DEFAULT + CALL criticalError("Mesh type " // meshType // " not supported.", "readGeometry") + + END SELECT + !Reads the mesh + fullpath = path // meshFile + CALL mesh%readMesh(fullPath) + + CASE DEFAULT + CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry") + + END SELECT + + END SUBROUTINE readGeometry + + !Reads the injection of particles from the boundaries + SUBROUTINE readInject(config) + USE moduleSpecies + USE moduleErrors + USE moduleInject + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + INTEGER:: i + character(2):: istring + character(:), allocatable:: object + logical:: found + CHARACTER(:), ALLOCATABLE:: speciesName + CHARACTER(:), ALLOCATABLE:: name + REAL(8):: v + REAL(8), ALLOCATABLE:: T(:), normal(:) + REAL(8):: flow + INTEGER:: physicalSurface + INTEGER:: pt + + CALL config%info('inject', n_children = nInject) + ALLOCATE(inject(1:nInject)) + nPartInj = 0 + DO i = 1, nInject + WRITE(istring, '(i2)') i + object = 'inject(' // trim(istring) // ')' + + !Find species + CALL config%get(object // '.species', speciesName, found) + pt = speciesName2Index(speciesName) + + CALL config%get(object // '.name', name, found) + CALL config%get(object // '.v', v, found) + CALL config%get(object // '.T', T, found) + CALL config%get(object // '.n', normal, found) + CALL config%get(object // '.flow', flow, found) + CALL config%get(object // '.physicalSurface', physicalSurface, found) + + CALL inject(i)%init(i, v, normal, T, flow, pt, physicalSurface) + + END DO + + !Allocate array for injected particles + IF (nPartInj > 0) THEN + ALLOCATE(part_inj(1:nPartInj)) + + END IF + + END SUBROUTINE readInject + +END MODULE moduleInput diff --git a/src/modules/moduleList.f95 b/src/modules/moduleList.f95 new file mode 100644 index 0000000..830688b --- /dev/null +++ b/src/modules/moduleList.f95 @@ -0,0 +1,79 @@ +!Linked list of particles +MODULE moduleList + IMPLICIT NONE + TYPE lNode + INTEGER:: n = 0 + TYPE(lNode), POINTER:: next => NULL() + END TYPE lNode + + TYPE listNode + INTEGER:: amount = 0 + TYPE(lNode),POINTER:: head => NULL() + TYPE(lNode),POINTER:: tail => NULL() + CONTAINS + PROCEDURE,PASS:: add => addToList + PROCEDURE,PASS:: get => getFromList + PROCEDURE,PASS:: erase => eraseList + END TYPE listNode + + CONTAINS + !Adds element to list + SUBROUTINE addToList(this,n) + INTEGER,INTENT(in):: n + CLASS(listNode):: this + TYPE(lNode),POINTER:: temp + + ALLOCATE(temp) + temp%n = n + NULLIFY(temp%next) + IF (.NOT. ASSOCIATED(this%head)) THEN + !First element + this%head => temp + this%tail => temp + this%amount = 1 + ELSE + !Append element + this%tail%next => temp + this%tail => temp + this%amount = this%amount + 1 + END IF + + END SUBROUTINE addToList + + FUNCTION getFromList(self, iObj) RESULT(n) + use moduleErrors + IMPLICIT NONE + + CLASS(listNode):: self + INTEGER:: iObj + INTEGER:: n + INTEGER:: i + TYPE(lNode), POINTER:: tempNode + + IF (iObj > self%amount) CALL criticalError('Accessing to element list outisde range','getFromList') + tempNode => self%head + DO i = 1, iObj - 1 + tempNode => tempNode%next + + END DO + n = tempNode%n + + END FUNCTION getFromList + + !Erase list + SUBROUTINE eraseList(this) + CLASS(listNode):: this + TYPE(lNode),POINTER:: current, next + + current => this%head + DO WHILE (ASSOCIATED(current)) + next => current%next + DEALLOCATE(current) + current => next + END DO + IF (ASSOCIATED(this%head)) NULLIFY(this%head) + IF (ASSOCIATED(this%tail)) NULLIFY(this%tail) + this%amount = 0 + END SUBROUTINE eraseList + +END MODULE moduleList diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95 new file mode 100644 index 0000000..3e54b02 --- /dev/null +++ b/src/modules/moduleMesh.f95 @@ -0,0 +1,198 @@ +!moduleMesh: General module for Finite Element mesh +MODULE moduleMesh + USE moduleList + USE moduleOutput + IMPLICIT NONE + + !Parent of Node element + TYPE, PUBLIC, ABSTRACT:: meshNode + !Node index + INTEGER:: n = 0 + !Node volume + REAL(8):: v = 0.D0 + !Output values + TYPE(outputNode), ALLOCATABLE:: output(:) + CONTAINS + PROCEDURE(initNode_interface), DEFERRED, PASS:: init + PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates + + END TYPE meshNode + + ABSTRACT INTERFACE + !Interface of init a node (3D generic coordinates) + SUBROUTINE initNode_interface(self, n, r) + IMPORT:: meshNode + CLASS(meshNode), INTENT(out):: self + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: r(1:3) + + END SUBROUTINE initNode_interface + + !Interface to get coordinates from node + FUNCTION getCoord_interface(self) RESULT(r) + IMPORT:: meshNode + CLASS(meshNode):: self + REAL(8):: r(1:3) + + END FUNCTION getCoord_interface + + END INTERFACE + + !Containers for nodes in the mesh + TYPE:: meshNodeCont + CLASS(meshNode), ALLOCATABLE:: obj + + END TYPE meshNodeCont + + !Parent of Edge element + TYPE, PUBLIC, ABSTRACT:: meshEdge + !Element index + INTEGER:: n = 0 + !Connectivity to vols + CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL() + !Normal vector + REAL(8):: normal(1:3) + !Physical surface in mesh + INTEGER:: physicalSurface + !id for boundary condition + INTEGER:: bt = 0 + CONTAINS + PROCEDURE(initEdge_interface), DEFERRED, PASS:: init + + END TYPE meshEdge + + ABSTRACT INTERFACE + SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface) + IMPORT:: meshEdge + CLASS(meshEdge), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + INTEGER, INTENT(in):: bt + INTEGER, INTENT(in):: physicalSurface + + END SUBROUTINE initEdge_interface + + END INTERFACE + + !Containers for edges in the mesh + TYPE:: meshEdgeCont + CLASS(meshEdge), ALLOCATABLE:: obj + + END TYPE meshEdgeCont + + !Parent of Volume element + TYPE, PUBLIC, ABSTRACT:: meshVol + !Volume index + INTEGER:: n = 0 + !Maximum collision rate + REAL(8):: sigmaVrelMax = 1.D-15 + !Volume + REAL(8):: volume = 0.D0 + !List of particles inside the volume + TYPE(listNode):: listPart_in + !Number of collisions per volume + INTEGER:: nColl = 0 + CONTAINS + PROCEDURE(initVol_interface), DEFERRED, PASS:: init + PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter + PROCEDURE(collision_interface), DEFERRED, PASS:: collision + PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell + + END TYPE meshVol + + ABSTRACT INTERFACE + SUBROUTINE initVol_interface(self, n, p) + IMPORT:: meshVol + CLASS(meshVol), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + + END SUBROUTINE initVol_interface + + SUBROUTINE scatter_interface(self, part) + USE moduleSpecies + + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + CLASS(particle), INTENT(in):: part + + END SUBROUTINE scatter_interface + + SUBROUTINE collision_interface(self) + IMPORT:: meshVol + CLASS(meshVol), INTENT(inout):: self + + END SUBROUTINE collision_interface + + SUBROUTINE findCell_interface(self, part) + USE moduleSpecies + + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + CLASS(particle), INTENT(inout):: part + + END SUBROUTINE findCell_interface + + END INTERFACE + + !Containers for volumes in the mesh + TYPE:: meshVolCont + CLASS(meshVol), ALLOCATABLE:: obj + + END TYPE meshVolCont + + !Abstract type of mesh + TYPE, PUBLIC, ABSTRACT:: meshGeneric + INTEGER:: numEdges, numNodes, numVols + !Array of nodes + TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) + !Array of boundary elements + TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) + !Array of volume elements + TYPE(meshVolCont), ALLOCATABLE:: vols(:) + !Global stiffness matrix + REAL(8), ALLOCATABLE, DIMENSION(:,:):: K + !Global load vector + REAL(8), ALLOCATABLE, DIMENSION(:):: F + + CONTAINS + PROCEDURE(readMesh_interface), PASS, DEFERRED:: readMesh + PROCEDURE(printOutput_interface), PASS, DEFERRED:: printOutput + PROCEDURE(printColl_interface), PASS, DEFERRED:: printColl + + END TYPE meshGeneric + + ABSTRACT INTERFACE + !Reads the mesh from a file + SUBROUTINE readMesh_interface(self, filename) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(out):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + + END SUBROUTINE readMesh_interface + + !Prints output variables + SUBROUTINE printOutput_interface(self, t) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + + END SUBROUTINE printOutput_interface + + !Prints unmber of collisions + SUBROUTINE printColl_interface(self, t) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + + END SUBROUTINE printColl_interface + + END INTERFACE + + !Generic mesh + CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh + +END MODULE moduleMesh diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 new file mode 100644 index 0000000..42dae52 --- /dev/null +++ b/src/modules/moduleMeshCyl.f95 @@ -0,0 +1,534 @@ +!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format. +! x == z +! y == r +! z == theta (unused) +MODULE moduleMeshCyl + USE moduleMesh + IMPLICIT NONE + + !TODO: Move this, this is horrible + REAL(8), PARAMETER:: w(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) + REAL(8), PARAMETER:: cor(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) + + TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl + !Element coordinates + REAL(8):: r = 0.D0, z = 0.D0 + CONTAINS + PROCEDURE, PASS:: init => initNodeCyl + PROCEDURE, PASS:: getCoordinates => getCoordCyl + + END TYPE meshNodeCyl + + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdgeCyl + !Element coordinates + REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() + CONTAINS + PROCEDURE, PASS:: init => initEdgeCyl + PROCEDURE(boundary_interface), PASS, DEFERRED:: fBoundary + + END TYPE meshEdgeCyl + + ABSTRACT INTERFACE + SUBROUTINE boundary_interface(self, part) + USE moduleSpecies + + IMPORT:: meshEdgeCyl + CLASS (meshEdgeCyl), INTENT(inout):: self + CLASS (particle), INTENT(inout):: part + + END SUBROUTINE + + END INTERFACE + + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl + CONTAINS + PROCEDURE, PASS:: collision => collision2DCyl + + END TYPE meshVolCyl + + TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad + !Element coordinates + REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() + !Connectivity to adjacent elements + CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() + !Maximum collision rate + !TODO: Move to other more appropiate section + REAL(8):: phi(1:4) = 0.D0 + REAL(8):: Ke(1:4,1:4) = 0.D0 + REAL(8):: arNodes(1:4) = 0.D0 + + CONTAINS + PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: locKe => localKeQuad + PROCEDURE, PASS:: detJac => detJQuad + PROCEDURE, PASS:: invJ => invJQuad + PROCEDURE, PASS:: area => areaQuad + PROCEDURE, NOPASS:: weight => weightQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: dVal => dValQuad + PROCEDURE, PASS:: scatter => scatterQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: findCell => findCellCylQuad + + END TYPE meshVolCylQuad + + CONTAINS + + !Inits node element + SUBROUTINE initNodeCyl(self, n, r) + USE moduleSpecies + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshNodeCyl), INTENT(out):: self + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: r(1:3) + + self%n = n + self%r = r(1)/L_ref + self%z = r(2)/L_ref + !Node volume, to be determined in mesh + self%v = 0.D0 + + !Allocates output: + ALLOCATE(self%output(1:nSpecies)) + + END SUBROUTINE initNodeCyl + + FUNCTION getCoordCyl(self) RESULT(r) + IMPLICIT NONE + + CLASS(meshNodeCyl):: self + REAL(8):: r(1:3) + + r = (/self%z, self%r, 0.D0/) + + END FUNCTION getCoordCyl + + !Edge functions + SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface) + IMPLICIT NONE + + CLASS(meshEdgeCyl), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + INTEGER, INTENT(in):: bt + INTEGER, INTENT(in):: physicalSurface + REAL(8):: r1(1:3), r2(1:3) + + self%n = n + self%n1 => mesh%nodes(p(1))%obj + self%n2 => mesh%nodes(p(2))%obj + !Get element coordinates + r1 = self%n1%getCoordinates() + r2 = self%n2%getCoordinates() + self%z = (/r1(1), r2(1)/) + self%r = (/r1(2), r2(2)/) + !Normal vector + self%normal = (/ self%r(2)-self%r(1), & + self%z(2)-self%z(1), & + 0.D0 /) + !Boundary index + self%bt = bt + !Phyiscal Surface + self%physicalSurface = physicalSurface + + END SUBROUTINE initEdgeCyl + + !Vol functions + !Quad Volume + SUBROUTINE initVolQuadCyl(self, n, p) + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + REAL(8):: r1(1:3), r2(1:3), r3(1:3), r4(1:3) + + self%n = n + self%n1 => mesh%nodes(p(1))%obj + self%n2 => mesh%nodes(p(2))%obj + self%n3 => mesh%nodes(p(3))%obj + self%n4 => mesh%nodes(p(4))%obj + !Get element coordinates + r1 = self%n1%getCoordinates() + r2 = self%n2%getCoordinates() + r3 = self%n3%getCoordinates() + r4 = self%n4%getCoordinates() + self%z = (/r1(1), r2(1), r3(1), r4(1)/) + self%r = (/r1(2), r2(2), r3(2), r4(2)/) + + !Assign node volume + CALL self%area() + self%n1%v = self%n1%v + self%arNodes(1) + self%n2%v = self%n2%v + self%arNodes(2) + self%n3%v = self%n3%v + self%arNodes(3) + self%n4%v = self%n4%v + self%arNodes(4) + + CALL self%locKe() + + self%sigmaVrelMax = sigma_ref/L_ref**2 + + END SUBROUTINE initVolQuadCyl + + FUNCTION fpsi(xi1,xi2) RESULT(psi) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi1,xi2 + REAL(8):: psi(1:4) + psi(1) = (1.D0-xi1)*(1.D0-xi2) + psi(2) = (1.D0+xi1)*(1.D0-xi2) + psi(3) = (1.D0+xi1)*(1.D0+xi2) + psi(4) = (1.D0-xi1)*(1.D0+xi2) + psi = psi*0.25D0 + + END FUNCTION fpsi + + !Derivative element function (xi1) + FUNCTION dpsiXi1(xi2) RESULT(psiXi1) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi2 + REAL(8):: psiXi1(1:4) + psiXi1(1) = -(1.D0-xi2) + psiXi1(2) = (1.D0-xi2) + psiXi1(3) = (1.D0+xi2) + psiXi1(4) = -(1.D0+xi2) + psiXi1 = psiXi1*0.25D0 + + END FUNCTION dpsiXi1 + + !Derivative element function (xi2) + FUNCTION dpsiXi2(xi1) RESULT(psiXi2) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi1 + REAL(8):: psiXi2(1:4) + psiXi2(1) = -(1.D0-xi1) + psiXi2(2) = -(1.D0+xi1) + psiXi2(3) = (1.D0+xi1) + psiXi2(4) = (1.D0-xi1) + psiXi2 = psiXi2*0.25D0 + + END FUNCTION dpsiXi2 + + !Transforms physical coordinates to element coordinates + FUNCTION phy2logQuad(this,r) RESULT(xN) + IMPLICIT NONE + + CLASS(meshVolCylQuad):: this + REAL(8):: r(1:2) + REAL(8):: xN(1:2), xO(1:2), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2) + REAL(8):: conv + + !Iterative newton method to transform coordinates + conv=1.D0 + xO=0.D0 + + DO WHILE(conv>1.D-4) + dPsi(1,:) = dpsiXi1(xO(2)) + dPsi(2,:) = dpsiXi2(xO(1)) + j = this%invJ(xO(1),xO(2),dPsi) + psi = fpsi(xO(1), xO(2)) + f(1) = DOT_PRODUCT(psi,this%z)-r(1) + f(2) = DOT_PRODUCT(psi,this%r)-r(2) + detJ = this%detJac(xO(1),xO(2),dPsi) + xN=xO - MATMUL(j, f)/detJ + conv=MAXVAL(DABS(xN-xO),1) + xO=xN + + END DO + + END FUNCTION phy2logQuad + + !Computes element local stiffness matrix + SUBROUTINE localKeQuad(this) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylQuad):: this + REAL(8):: r, xi1, xi2, dPsi(1:2,1:4) + REAL(8):: j(1:2,1:2), detJ + INTEGER:: l, m + + this%Ke=0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 3 + DO m = 1, 3 + xi1 = cor(l) + xi2 = cor(m) + dPsi(1,:) = dpsiXi1(xi2) + dPsi(2,:) = dpsiXi2(xi1) + detJ = this%detJac(xi1,xi2,dPsi) + j = this%invJ(xi1,xi2,dPsi) + r = DOT_PRODUCT(fpsi(xi1,xi2),this%r) + this%Ke = this%Ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ + + END DO + END DO + this%Ke = this%Ke*2.D0*PI + + END SUBROUTINE localKeQuad + + !Computes element Jacobian determinant + FUNCTION detJQuad(this,xi1,xi2,dPsi_in) RESULT(dJ) + IMPLICIT NONE + + REAL(8),OPTIONAL:: dPsi_in(1:2,1:4) + REAL(8):: dPsi(1:2,1:4) + REAL(8):: dz(1:2), dr(1:2) + REAL(8):: xi1,xi2, dJ + CLASS(meshVolCylQuad):: this + + IF(PRESENT(dPsi_in)) THEN + dPsi=dPsi_in + + ELSE + dPsi(1,:) = dpsiXi1(xi2) + dPsi(2,:) = dpsiXi2(xi1) + + END IF + + dz(1) = DOT_PRODUCT(dPsi(1,:),this%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),this%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),this%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),this%r) + dJ = dz(1)*dr(2)-dz(2)*dr(1) + + END FUNCTION detJQuad + + FUNCTION invJQuad(this,xi1,xi2,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant) + IMPLICIT NONE + + REAL(8),OPTIONAL:: dpsi_in(1:2,1:4) + REAL(8):: dPsi(1:2,1:4) + REAL(8):: dz(1:2), dr(1:2) + REAL(8):: xi1,xi2, j(1:2,1:2) + CLASS(meshVolCylQuad):: this + + IF(PRESENT(dPsi_in)) THEN + dPsi=dPsi_in + + ELSE + dPsi(1,:) = dpsiXi1(xi2) + dPsi(2,:) = dpsiXi2(xi1) + + END IF + + dz(1) = DOT_PRODUCT(dPsi(1,:),this%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),this%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),this%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),this%r) + j(1,:) = (/ dr(2), -dz(2) /) + j(2,:) = (/ -dr(1), dz(1) /) + + END FUNCTION invJQuad + + SUBROUTINE areaQuad(this)!Computes element area + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylQuad):: this + REAL(8):: r, xi1, xi2 + REAL(8):: detJ, fpsi0(1:4) + + this%volume = 0.D0 + this%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + xi1 = 0.D0 + xi2 = 0.D0 + detJ = this%detJac(xi1,xi2)*8.D0*PI !4*2*pi + fpsi0 = fpsi(xi1,xi2) + r = DOT_PRODUCT(fpsi0,this%r) + this%volume = r*detJ + this%arNodes = fpsi0*r*detJ + + END SUBROUTINE areaQuad + + FUNCTION weightQuad(xi1_p, xi2_p) RESULT(w) !Computes weights in the four vertices. '_p' references particle position in the logical space + IMPLICIT NONE + + REAL(8):: xi1_p, xi2_p + REAL(8):: w(1:4) + + w = fpsi(xi1_p, xi2_p) + + END FUNCTION weightQuad + + FUNCTION insideQuad(xi1_p, xi2_p) RESULT(ins) !Checks if a particle is inside a quad element + IMPLICIT NONE + + REAL(8):: xi1_p, xi2_p + LOGICAL:: ins + + ins = (xi1_p >= -1.D0 .AND. xi1_p <= 1.D0) .AND. & + (xi2_p >= -1.D0 .AND. xi2_p <= 1.D0) + + END FUNCTION insideQuad + + FUNCTION dValQuad(this,xi1,xi2) RESULT(EF)!Computes electric field in xi1,xi2 + IMPLICIT NONE + + CLASS(meshVolCylQuad):: this + REAL(8):: xi1, xi2, dPsi(1:2,1:4) + REAL(8):: j(1:2,1:2), detJ + REAL(8):: EF(1:3) + + dPsi(1,:) = dpsiXi1(xi2) + dPsi(2,:) = dpsiXi2(xi1) + detJ = this%detJac(xi1,xi2,dPsi) + j = this%invJ(xi1,xi2,dPsi) + EF(1) = -DOT_PRODUCT(dPsi(1,:), this%phi)*j(2,2)/detJ + EF(2) = -DOT_PRODUCT(dPsi(2,:), this%phi)*j(1,1)/detJ + EF(3) = 0.D0 + + END FUNCTION dValQuad + + SUBROUTINE scatterQuad(self, part) + USE moduleOutput + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(particle), INTENT(in):: part + TYPE(outputNode), POINTER:: vertex + REAL(8):: w_p(1:4) + REAL(8):: tensorS(1:3,1:3) + + w_p = self%weight(part%xLog(1), part%xLog(2)) + tensorS = outerProduct(part%v, part%v) + + vertex => self%n1%output(part%pt) + vertex%den = vertex%den + w_p(1) + vertex%mom(:) = vertex%mom(:) + w_p(1)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(1)*tensorS + + vertex => self%n2%output(part%pt) + vertex%den = vertex%den + w_p(2) + vertex%mom(:) = vertex%mom(:) + w_p(2)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(2)*tensorS + + vertex => self%n3%output(part%pt) + vertex%den = vertex%den + w_p(3) + vertex%mom(:) = vertex%mom(:) + w_p(3)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(3)*tensorS + + vertex => self%n4%output(part%pt) + vertex%den = vertex%den + w_p(4) + vertex%mom(:) = vertex%mom(:) + w_p(4)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(4)*tensorS + + + END SUBROUTINE scatterQuad + + RECURSIVE SUBROUTINE findCellCylQuad(self, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(particle), INTENT(inout):: part + REAL(8):: xLog(1:2) + REAL(8):: xLogArray(1:4) + CLASS(*), POINTER:: nextElement + INTEGER:: nextInt + + + xLog = self%phy2log(part%r(1:2)) + IF (self%inside(xLog(1), xLog(2))) THEN + !Checks if particle is inside of current cell + part%e_p = self%n + part%xLog = xLog + ELSE + !If not, searches for a neighbour and repeats the process. + xLogArray = (/ -xLog(2), xLog(1), xLog(2), -xLog(1) /) + nextInt = MAXLOC(xLogArray,1) + !Selects the higher value of directions and searches in that direction + NULLIFY(nextElement) + SELECT CASE (nextInt) + CASE (1) + nextElement => self%e1 + CASE (2) + nextElement => self%e2 + CASE (3) + nextElement => self%e3 + CASE (4) + nextElement => self%e4 + END SELECT + + !Defines the next step + SELECT TYPE(nextElement) + CLASS IS(meshVolCyl) + !Particle moved to new cell, repeat find procedure + CALL nextElement%findCell(part) + + CLASS IS (meshEdgeCyl) + !Particle encountered an edge, execute boundary + CALL nextElement%fBoundary(part) + + CLASS DEFAULT + WRITE(*,*) "ERROR, CHECK findCellCylQuad" + + END SELECT + END IF + + END SUBROUTINE findCellCylQuad + + SUBROUTINE collision2DCyl(self) + USE moduleRefParam + USE moduleConstParam + USE moduleCollisions + USE moduleSpecies + USE moduleList + IMPLICIT NONE + + CLASS(meshVolCyl), INTENT(inout):: self + INTEGER:: Npart !Number of particles inside the cell + INTEGER:: Ncoll !Maximum number of collisions + REAL(8):: Fn !Specific weight + REAL(8):: Pmax !Maximum probability of collision + INTEGER:: i,j !random particles index + INTEGER:: rnd !random index + TYPE(particle), POINTER:: part_i, part_j + INTEGER:: n !collision + INTEGER:: ij, k + REAL(8):: sigmaVrel + REAL(8):: sigmaVrelMaxNew + + Fn = species(1)%obj%weight!Check how to do this for multiple species + Npart = self%listPart_in%amount + Pmax = Fn*self%sigmaVrelMax*tau/self%volume + Ncoll = INT(REAL(Npart*(Npart-1))*Pmax*0.5D0) + self%nColl = Ncoll + + DO n = 1, NColl + !Select random numbers + rnd = 1 + FLOOR(Npart*RAND()) + i = self%listPart_in%get(rnd) + part_i => part_old(i) + rnd = 1 + FLOOR(Npart*RAND()) + j = self%listPart_in%get(rnd) + part_j => part_old(j) + ij = interactionIndex(part_i%pt, part_j%pt) + sigmaVrelMaxNew = 0.D0 + DO k = 1, interactionMatrix(ij)%amount + CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrel, part_i, part_j) + sigmaVrelMaxNew = sigmaVrel + SigmaVrelMaxNew + + END DO + + !Update maximum cross section times vrelative per each collision + IF (sigmaVrelMaxNew > self%sigmaVrelMax) self%sigmaVrelMax = sigmaVrelMaxNew + + + END DO + + CALL self%listPart_in%erase() + + END SUBROUTINE collision2DCyl + +END MODULE moduleMeshCyl diff --git a/src/modules/moduleMeshCylBoundary.f95 b/src/modules/moduleMeshCylBoundary.f95 new file mode 100644 index 0000000..fafe56d --- /dev/null +++ b/src/modules/moduleMeshCylBoundary.f95 @@ -0,0 +1,80 @@ +!moduleMeshCylBoundary: Edge elements for Cylindrical mesh. +MODULE moduleMeshCylBoundary + USE moduleMeshCyl + + TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylRef + CONTAINS + PROCEDURE, PASS:: fBoundary => reflection + + END TYPE meshEdgeCylRef + + TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAbs + CONTAINS + PROCEDURE, PASS:: fBoundary => absorption + + END TYPE meshEdgeCylAbs + + TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAxis + CONTAINS + PROCEDURE, PASS:: fBoundary => symmetryAxis + + END TYPE meshEdgeCylAxis + + CONTAINS + SUBROUTINE reflection(self, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdgeCylRef), INTENT(inout):: self + CLASS(particle), INTENT(inout):: part + REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2) + + edgeNorm = DSQRT((self%r(2)-self%r(1))**2 + (self%z(2)-self%z(1))**2) + cosT = (self%z(2)-self%z(1))/edgeNorm + sinT = DSQRT(1-cosT**2) + + rp(1) = part%r(1) - self%z(1); + rp(2) = part%r(2) - self%r(1); + + rpp(1) = cosT*rp(1) - sinT*rp(2) + rpp(2) = sinT*rp(1) + cosT*rp(2) + rpp(2) = -rpp(2) + + vpp(1) = cosT*part%v(1) - sinT*part%v(2) + vpp(2) = sinT*part%v(1) + cosT*part%v(2) + vpp(2) = -vpp(2) + + part%r(1) = cosT*rpp(1) + sinT*rpp(2) + self%z(1); + part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + self%r(1); + part%v(1) = cosT*vpp(1) + sinT*vpp(2) + part%v(2) = -sinT*vpp(1) + cosT*vpp(2) + + END SUBROUTINE reflection + + !Absoption in a surface + SUBROUTINE absorption(self, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdgeCylAbs), INTENT(inout):: self + CLASS(particle), INTENT(inout):: part + + + !TODO: Add scatter to mesh nodes + part%n_in = .FALSE. + + END SUBROUTINE absorption + + SUBROUTINE symmetryAxis(self, part) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshEdgeCylAxis), INTENT(inout):: self + CLASS(particle), INTENT(inout):: part + + !Do to the 2D Cyl pusher, this function should never be called. + + END SUBROUTINE symmetryAxis + + +END MODULE moduleMeshCylBoundary diff --git a/src/modules/moduleMeshCylRead.f95 b/src/modules/moduleMeshCylRead.f95 new file mode 100644 index 0000000..c3e3b69 --- /dev/null +++ b/src/modules/moduleMeshCylRead.f95 @@ -0,0 +1,428 @@ +MODULE moduleMeshCylRead + USE moduleMesh + USE moduleMeshCyl + USE moduleMeshCylBoundary + + TYPE, EXTENDS(meshGeneric):: meshCylGeneric + CONTAINS + PROCEDURE, PASS:: readMesh => readMeshCyl + PROCEDURE, PASS:: printOutput => printOutputCyl + PROCEDURE, PASS:: printColl => printCollisionsCyl + + END TYPE + + INTERFACE connected + MODULE PROCEDURE connectedVolVol, connectedVolEdge + + END INTERFACE connected + + CONTAINS + SUBROUTINE readMeshCyl(self, filename) + USE moduleRefParam + USE moduleBoundary + IMPLICIT NONE + + CLASS(meshCylGeneric), INTENT(out):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + REAL(8):: r, z + INTEGER:: p(1:4) + INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0 + INTEGER:: totalNumElem + INTEGER:: boundaryType + + !Read msh + OPEN(10, FILE=TRIM(filename)) + !Skip header + READ(10, *) + READ(10, *) + READ(10, *) + READ(10, *) + !Read number of nodes + READ(10, *) self%numNodes + !Allocate required matrices and vectors + ALLOCATE(self%nodes(1:self%numNodes)) + ! ALLOCATE(self%K(1:numNodes,1:numNodes)) + ! ALLOCATE(self%F(1:numNodes)) + !Read nodes cartesian coordinates (x=z, y=r, z=null) + DO e=1, self%numNodes + READ(10, *) n, z, r + ALLOCATE(meshNodeCyl:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /)) + + END DO + !Skips comments + READ(10, *) + READ(10, *) + !Reads Totalnumber of elements + READ(10, *) TotalnumElem + !counts edges and volume elements + self%numEdges = 0 + DO e=1, TotalnumElem + READ(10,*) eTemp, elemType + IF (elemType==1) THEN + self%numEdges=e + END IF + END DO + !Substract the number of edges to the total number of elements + !to obtain the number of volume elements + self%numVols = TotalnumElem - self%numEdges + !Allocates arrays + ALLOCATE(self%edges(1:self%numEdges)) + ALLOCATE(self%vols(1:self%numVols)) + + !Go back to the beggining to read elements + DO e=1, totalNumElem + BACKSPACE(10) + END DO + + !TODO: symplify to read all elements independently if they are edges or vols + !Reads edges + DO e=1, self%numEdges + READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) + !Associate boundary condition procedure. + !TODO: move to subroutine + bt = getBoundaryId(boundaryType) + SELECT CASE(boundary(bt)%obj%boundaryType) + CASE ('reflection') + ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj) + + CASE ('absorption') + ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj) + + CASE ('axis') + ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj) + + END SELECT + + CALL self%edges(e)%obj%init(n, p, bt, boundaryType) + + END DO + + !Read and initialize volumes + !TODO: Extend to triangular elements + DO e=1, self%numVols + READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) + ALLOCATE(meshVolCylQuad:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p) + + END DO + + CLOSE(10) + + !Build connectivity between elements + DO e = 1, self%numVols + !Connectivity between volumes + DO et = 1, self%numVols + CALL connected(self%vols(e)%obj, self%vols(et)%obj) + + END DO + + !Connectivity between vols and edges + DO et = 1, self%numEdges + CALL connected(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + ! !Compute global stiffness matrix + ! GlobalK=0.D0 + ! DO e=1, numElem + ! DO i=1, 4 + ! DO j=1, 4 + ! GlobalK(elems(e)%p(i),elems(e)%p(j)) = GlobalK(elems(e)%p(i),elems(e)%p(j)) + elems(e)%Ke(i,j) + ! END DO + ! END DO + ! END DO + ! ! Apply Dirichlet boundary conditions to GlobalK + ! DO n=1, numNodes + ! IF (nodes(n)%bound == 1) THEN + ! GlobalK(n,:)=0.D0 + ! GlobalK(n,n)=1.D0 + ! END IF + ! END DO + + END SUBROUTINE + + SUBROUTINE connectedVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + + SELECT TYPE(elemA) + TYPE IS (meshVolCylQuad) + SELECT TYPE(elemB) + TYPE IS(meshVolCylQuad) + CALL connectedQuadQuad(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectedVolVol + + SUBROUTINE connectedVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS (meshVolCylQuad) + SELECT TYPE(elemB) + CLASS IS(meshEdgeCyl) + CALL connectedQuadEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectedVolEdge + + SUBROUTINE connectedQuadQuad(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n1%n == elemB%n4%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n4%n) THEN + elemA%e2 => elemB + elemB%e4 => elemA + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3) .AND. & + elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4) .AND. & + elemA%n4%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + END IF + + END SUBROUTINE connectedQuadQuad + + SUBROUTINE connectedQuadEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n4%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (elemA%n4%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e4 => elemB + elemB%e1 => elemA + + END IF + + END IF + + END SUBROUTINE connectedQuadEdge + + SUBROUTINE printOutputCyl(self, t) + USE moduleRefParam + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshCylGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n, i + TYPE(outputFormat):: output(1:self%numNodes) + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations + + time = DBLE(t)*tau*ti_ref + + DO i = 1, nSpecies + WRITE(tstring, '(I6.6)') t + fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' + PRINT *, "Creado archivo de datos:", fileName + OPEN (60, file = path // folder // '/' // fileName) + WRITE(60, "(A)") '$MeshFormat' + WRITE(60, "(A)") '2.2 0 8' + WRITE(60, "(A)") '$EndMeshFormat' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Density (m^-3)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) + WRITE(60, "(I6,ES20.6E3)") n, output(n)%density + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Velocity (m/s)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 3 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Pressure (Pa)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Temperature (K)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature + END DO + WRITE(60, "(A)") '$EndNodeData' + CLOSE (60) + + END DO + + END SUBROUTINE printOutputCyl + + SUBROUTINE printCollisionsCyl(self, t) + USE moduleRefParam + USE moduleCaseParam + USE moduleOutput + IMPLICIT NONE + + CLASS(meshCylGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations + + + IF (collOutput) THEN + time = DBLE(t)*tau*ti_ref + WRITE(tstring, '(I6.6)') t + + fileName='OUTPUT_' // tstring// '_Collisions.msh' + PRINT *, "Creado archivo de datos:", fileName + OPEN (60, file = path // folder // '/' // fileName) + WRITE(60, "(A)") '$MeshFormat' + WRITE(60, "(A)") '2.2 0 8' + WRITE(60, "(A)") '$EndMeshFormat' + WRITE(60, "(A)") '$ElementData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Collisions"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numVols + DO n=1, self%numVols + WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl + END DO + WRITE(60, "(A)") '$EndElementData' + + CLOSE(60) + + END IF + + END SUBROUTINE printCollisionsCyl + +END MODULE moduleMeshCylRead diff --git a/src/modules/moduleOutput.f95 b/src/modules/moduleOutput.f95 new file mode 100644 index 0000000..7ccff50 --- /dev/null +++ b/src/modules/moduleOutput.f95 @@ -0,0 +1,118 @@ +!Contains information about output +MODULE moduleOutput + IMPLICIT NONE + !Output for each node + TYPE outputNode + REAL(8):: den, mom(1:3), tensorS(1:3,1:3) + + END TYPE + + !Output in dimensional units to print + TYPE outputFormat + REAL(8):: density, velocity(1:3), pressure, temperature + + END TYPE + + CHARACTER(:), ALLOCATABLE:: path + CHARACTER(:), ALLOCATABLE:: folder + INTEGER:: triggerOutput, counterOutput + LOGICAL:: timeOutput = .FALSE. + LOGICAL:: collOutput = .FALSE. + + CONTAINS + + FUNCTION outerProduct(a,b) RESULT(s) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3):: a, b + REAL(8):: s(1:3,1:3) + + s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) + + END FUNCTION outerProduct + + FUNCTION tensorTrace(a) RESULT(t) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3,1:3):: a + REAL(8):: t + + t = 0.D0 + t = a(1,1)+a(2,2)+a(3,3) + + END FUNCTION tensorTrace + + SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn) + USE moduleRefParam + USE moduleSpecies + IMPLICIT NONE + + TYPE(outputNode), INTENT(in):: rawValues + TYPE(outputFormat), INTENT(out):: formatValues + REAL(8), INTENT(in):: nodeVol + CLASS(speciesGeneric), INTENT(in):: speciesIn + REAL(8), DIMENSION(1:3,1:3):: tensorTemp + REAL(8), DIMENSION(1:3):: tempVel + REAL(8):: tempVol + + !Resets the node outputs + formatValues%density = 0.D0 + formatValues%velocity = 0.D0 + formatValues%pressure = 0.D0 + formatValues%temperature = 0.D0 + tempVol = speciesIn%weight/(nodeVol*Vol_ref) + IF (rawValues%den > 0.D0) THEN + tempVel = rawValues%mom(:)/rawValues%den + tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel)) + formatValues%density = rawValues%den*tempVol + formatValues%velocity(:) = tempVel + IF (tensorTrace(tensorTemp) > 0.D0) THEN + formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0 + formatValues%temperature = formatValues%pressure/(formatValues%density*kb) + + END IF + END IF + + formatValues%velocity = formatValues%velocity*v_ref + formatValues%pressure = formatValues%pressure*m_ref*v_ref**2 + formatValues%temperature = formatValues%temperature*m_ref*v_ref**2 + + END SUBROUTINE calculateOutput + + SUBROUTINE printTime(t, first) + USE moduleSpecies + USE moduleCompTime + IMPLICIT NONE + + INTEGER, INTENT(in):: t + LOGICAL, INTENT(in), OPTIONAL:: first + CHARACTER(:), ALLOCATABLE:: filename + + filename = 'cpuTime.dat' + + IF (timeOutput) THEN + IF (PRESENT(first)) THEN + IF (first) THEN + OPEN(20, file = path // folder // '/' // filename, action = 'write') + WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting" + + ELSE + + END IF + OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write') + + ELSE + OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write') + + END IF + + WRITE (20, "(I10, I10, 5(ES20.6E3))") t, n_part_old, tStep, tPush, tReset, tColl, tWeight + + CLOSE(20) + + END IF + + END SUBROUTINE + +END MODULE moduleOutput + diff --git a/src/modules/moduleRefParam.f95 b/src/modules/moduleRefParam.f95 new file mode 100644 index 0000000..152f103 --- /dev/null +++ b/src/modules/moduleRefParam.f95 @@ -0,0 +1,10 @@ +!Reference parameters +MODULE moduleRefParam + USE moduleConstParam + !Parameters that define the problem (inputs) + REAL(8):: n_ref, m_ref, T_ref, r_ref, sigma_ref!, q_ref=1.6022D-19 + !Reference parameters for non-dimensional problem + REAL(8):: L_ref, v_ref, ti_ref, Vol_ref + +END MODULE moduleRefParam + diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 new file mode 100644 index 0000000..48f3e7e --- /dev/null +++ b/src/modules/moduleSolver.f95 @@ -0,0 +1,113 @@ +MODULE moduleSolver + + CONTAINS + + SUBROUTINE scatterGrid() + USE moduleSpecies + USE moduleRefParam + USE moduleMesh + USE moduleOutput + INTEGER:: n, e, k + + !Cleans previous output + !$OMP DO PRIVATE(k) + DO e = 1, mesh%numNodes + DO k= 1, nSpecies + mesh%nodes(e)%obj%output(k)%den = 0.D0 + mesh%nodes(e)%obj%output(k)%mom = 0.D0 + mesh%nodes(e)%obj%output(k)%tensorS = 0.D0 + + END DO + + END DO + !$OMP END DO + + !Loops over the particles to scatter them + !$OMP DO + DO n=1, n_part_old + CALL mesh%vols(part_old(n)%e_p)%obj%scatter(part_old(n)) + + END DO + !$OMP END DO + + + END SUBROUTINE scatterGrid + + SUBROUTINE push(part) + USE moduleSpecies + USE moduleMesh + IMPLICIT NONE + TYPE(particle), INTENT(inout):: part + REAL(8):: v_p_oh_star(2:3) + TYPE(particle):: part_temp + REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha, alpha + + part_temp = part + !z + part_temp%v(1) = part%v(1) + part_temp%r(1) = part%r(1) + part_temp%v(1)*tau + !r,theta + v_p_oh_star(2) = part%v(2) + x_new = part%r(2) + v_p_oh_star(2)*tau + v_p_oh_star(3) = part%v(3) + y_new = v_p_oh_star(3)*tau + r = DSQRT(x_new**2+y_new**2) + part_temp%r(2) = r + alpha = 0.D0!ATAN2(y_new,x_new) !0 in 2D problem + part_temp%r(3) = part%r(3) + alpha + IF (r > 0.D0) THEN + sin_alpha = y_new/r + cos_alpha = x_new/r + ELSE + sin_alpha = 0.D0 + cos_alpha = 1.D0 + END IF + part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3) + part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3) + part_temp%pt = part%pt + part_temp%e_p = part%e_p + !Assign cell to particle + part=part_temp + CALL mesh%vols(part%e_p)%obj%findCell(part) + + END SUBROUTINE push + + SUBROUTINE resetParticles() + USE moduleSpecies + USE moduleList + USE moduleMesh + IMPLICIT NONE + + INTEGER:: nn, n + TYPE(particle), ALLOCATABLE:: partTemp(:) + + IF (n_part_old > 0) THEN + n_part_new = COUNT(part_old%n_in) + COUNT(part_inj%n_in) + ELSE + n_part_new = COUNT(part_inj%n_in) + END IF + + CALL MOVE_ALLOC(part_old, partTemp) + ALLOCATE(part_old(1:n_part_new)) + + nn = 0 + DO n = 1, nPartInj + IF (part_inj(n)%n_in) THEN + nn = nn + 1 + part_old(nn) = part_inj(n) + CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn) + END IF + END DO + DO n = 1, n_part_old + IF (partTemp(n)%n_in) THEN + nn = nn + 1 + part_old(nn) = partTemp(n) + CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn) + END IF + END DO + n_part_old = n_part_new + + END SUBROUTINE resetParticles + +END MODULE moduleSolver + diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 new file mode 100644 index 0000000..b5860e4 --- /dev/null +++ b/src/modules/moduleSpecies.f95 @@ -0,0 +1,66 @@ +!Contains the information about species (particles) +MODULE moduleSpecies + USE moduleCaseParam + USE moduleList + IMPLICIT NONE + + TYPE, ABSTRACT:: speciesGeneric + CHARACTER(:), ALLOCATABLE:: name + REAL(8):: m=0.D0, weight=0.D0 + INTEGER:: pt=0 + END TYPE speciesGeneric + + TYPE, EXTENDS(speciesGeneric):: speciesNeutral + + END TYPE speciesNeutral + + TYPE, EXTENDS(speciesGeneric):: speciesCharged + REAL(8):: q=0.D0, qm=0.D0 + + END TYPE speciesCharged + + TYPE:: speciesCont + CLASS(speciesGeneric), ALLOCATABLE:: obj + + END TYPE + + INTEGER:: nSpecies + TYPE(speciesCont), ALLOCATABLE:: species(:) + + TYPE particle + REAL(8):: r(1:3) !Position + REAL(8):: v(1:3) !Velocity + INTEGER:: pt !Particle species id + INTEGER:: e_p !Index of element in which the particle is located + REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p. + LOGICAL:: n_in !Flag that indicates if a particle is in the domain + + END TYPE particle + + INTEGER:: n_part_old, n_part_new + INTEGER:: nPartInj + !Arrays that define the particles + TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: part_old + TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_new + TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_inj !array of inject particles + + CONTAINS + FUNCTION speciesName2Index(speciesName) RESULT(pt) + IMPLICIT NONE + + CHARACTER(:), ALLOCATABLE:: speciesName + INTEGER:: pt + INTEGER:: n + + pt = 0 + DO n = 1, nSpecies + IF (speciesName == species(n)%obj%name) THEN + pt = species(n)%obj%pt + EXIT + END IF + END DO + !TODO: add error handling when species not found + + END FUNCTION speciesName2Index + +END MODULE moduleSpecies diff --git a/src/modules/moduleTable.f95 b/src/modules/moduleTable.f95 new file mode 100644 index 0000000..1ea89e5 --- /dev/null +++ b/src/modules/moduleTable.f95 @@ -0,0 +1,121 @@ +MODULE moduleTable + + TYPE:: table1D + REAL(8):: xMin, xMax + REAL(8):: fMin, fMax + REAL(8), ALLOCATABLE, DIMENSION(:):: x, f, k + CONTAINS + PROCEDURE, PASS:: init => initTable1D + PROCEDURE, PASS:: get => getValueTable1D + PROCEDURE, PASS:: convert => convertUnits + + END TYPE table1D + + CONTAINS + SUBROUTINE initTable1D(self, tableFile) + USE moduleErrors + IMPLICIT NONE + + CLASS(table1D), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(IN):: tableFile + CHARACTER(100):: dummy + INTEGER:: amount + INTEGER:: i + INTEGER:: stat + INTEGER:: id = 20 + + OPEN(id, file = tableFile) + amount = 0 + DO + READ(id, '(A)', iostat = stat) dummy + !Skip comment + IF (INDEX(dummy,'#') /= 0) CYCLE + !If EOF or error, exit file + IF (stat /= 0) EXIT + !Add row + amount = amount + 1 + + END DO + + IF (amount == 0) CALL criticalError('Table ' // tableFile // ' is empty', 'initTable1D') + IF (amount == 1) CALL criticalError('Table ' // tableFile // ' has only one row', 'initTable1D') + + !Go bback to initial point + REWIND(id) + + !Allocate table arrays + ALLOCATE(self%x(1:amount), self%f(1:amount), self%k(1:amount)) + self%x = 0.D0 + self%f = 0.D0 + self%k = 0.D0 + + i = 0 + DO + READ(id, '(A)', iostat = stat) dummy + IF (INDEX(dummy,'#') /= 0) CYCLE + IF (stat /= 0) EXIT + !Add data + !TODO: substitute with extracting information from dummy + BACKSPACE(id) + i = i + 1 + READ(id, *) self%x(i), self%f(i) + + END DO + + CLOSE(10) + + self%xMin = self%x(1) + self%xMax = self%x(amount) + self%fMin = self%f(1) + self%fMax = self%f(amount) + + DO i = 1, amount - 1 + self%k(i) = (self%f(i+1) - self%f(i))/(self%x(i+1) - self%x(i)) + + END DO + + END SUBROUTINE initTable1D + + FUNCTION getValueTable1D(self, x) RESULT(f) + IMPLICIT NONE + + CLASS(table1D), INTENT(in):: self + REAL(8):: x + REAL(8):: f + REAL(8):: deltaX + INTEGER:: i + + IF (x <= self%xMin) THEN + f = self%fMin + + ELSEIF (x >= self%xMax) THEN + f = self%fMax + + ELSE + i = MINLOC(x - self%x, 1) + deltaX = x - self%x(i) + IF (deltaX < 0 ) THEN + i = i - 1 + deltaX = x - self%x(i) + + END IF + + f = self%f(i) + self%k(i)*deltaX + + END IF + + END FUNCTION getValueTable1D + + SUBROUTINE convertUnits(self, data_x, data_f) + IMPLICIT NONE + + CLASS(table1D), INTENT(inout):: self + REAL(8):: data_x, data_f + + self%x = self%x * data_x + self%f = self%f * data_f + self%k = self%k * data_f / data_x + + END SUBROUTINE convertUnits + +END MODULE moduleTable