commit bd7e8b040b85c4f3c068b6c4817c716dcd8beb1f Author: Jorge Gonzalez Date: Fri Oct 9 08:45:07 2020 +0200 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. 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 0000000..a7b4465 Binary files /dev/null and b/doc/PPartiC_UserManual.pdf differ 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