From 6cdb24497ecead3b6a53ea82dca19425cff4d9ed Mon Sep 17 00:00:00 2001
From: Holger Grosshans <holger.grosshans@ptb.de>
Date: Fri, 23 Feb 2024 10:25:16 +0000
Subject: [PATCH] Update 19 files
- /input/input_example.dat
- /src/COPYING
- /src/bc.f90
- /src/electrostatics.f90
- /src/fluid.f90
- /src/main.f90
- /src/makefile
- /src/parallel.f90
- /src/particles.f90
- /src/particlesTransport.f90
- /src/post.f90
- /src/pre.f90
- /src/restart.f90
- /src/timestep.f90
- /src/var.f90
- /src/write_vtk.f90
- /src/writedat_particles.f90
- /clean.sh
- /monitor.plt
---
clean.sh | 4 -
input/input_example.dat | 30 --
monitor.plt | 17 -
src/COPYING | 674 -------------------------------------
src/bc.f90 | 264 ---------------
src/electrostatics.f90 | 274 ---------------
src/fluid.f90 | 387 ---------------------
src/main.f90 | 53 ---
src/makefile | 51 ---
src/parallel.f90 | 340 -------------------
src/particles.f90 | 300 -----------------
src/particlesTransport.f90 | 493 ---------------------------
src/post.f90 | 122 -------
src/pre.f90 | 400 ----------------------
src/restart.f90 | 104 ------
src/timestep.f90 | 52 ---
src/var.f90 | 195 -----------
src/write_vtk.f90 | 450 -------------------------
src/writedat_particles.f90 | 32 --
19 files changed, 4242 deletions(-)
delete mode 100755 clean.sh
delete mode 100755 input/input_example.dat
delete mode 100755 monitor.plt
delete mode 100755 src/COPYING
delete mode 100755 src/bc.f90
delete mode 100755 src/electrostatics.f90
delete mode 100755 src/fluid.f90
delete mode 100755 src/main.f90
delete mode 100755 src/makefile
delete mode 100755 src/parallel.f90
delete mode 100755 src/particles.f90
delete mode 100755 src/particlesTransport.f90
delete mode 100755 src/post.f90
delete mode 100755 src/pre.f90
delete mode 100755 src/restart.f90
delete mode 100755 src/timestep.f90
delete mode 100755 src/var.f90
delete mode 100755 src/write_vtk.f90
delete mode 100755 src/writedat_particles.f90
diff --git a/clean.sh b/clean.sh
deleted file mode 100755
index 7e5544e..0000000
--- a/clean.sh
+++ /dev/null
@@ -1,4 +0,0 @@
-rm output/grid* output/output.dat src/*.o src/*.mod src/pafiX
-rm -r results
-mkdir results
-rm *.e*
diff --git a/input/input_example.dat b/input/input_example.dat
deleted file mode 100755
index 1d49f6e..0000000
--- a/input/input_example.dat
+++ /dev/null
@@ -1,30 +0,0 @@
-pafiX v.1.2.0 input data
-----------------
-dimensions in x,y,z-direction (m,m,m) =
-0.24,0.125,0.04
-boundary conditions in x,y,z-direction [(w)all/(p)eriodic/(i)n-outlet] =
-p,p,w
-number of cells in x,y,z-direction (-,-,-) =
-16,16,16
-grid in x,y,z-direction [(u)niform/(s)tretch] =
-u,u,s
-start time-step (-), time-steps to compute (-) =
-1,10
-interval of writing result files, restart files (-,-) =
-10,10000
-turbulence model [-/Smagorinsky] =
-Smagorinsky
-CFL number (-), max number interations (-), relative tolerance (-) =
-0.4,200,5.E-3
-in/initial flow (m/s), pressure gradient (N/m**3) =
-3.65,-2.9
-fluid density (kg/m**3), fluid kinematic viscosity (m**2/s) =
-1.2,1.46E-5
-particle number density (-/m**3) or rate (-/s), start seed at nt (-) =
-1E7,1
-particle radius (m), material density (kg/m**3), restitution coeff (-) =
-50E-6,1200,0.9
-particle charge: initial, equilibrium (C,C) =
-0.,1E-12
-gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) =
-9.81,0,0
diff --git a/monitor.plt b/monitor.plt
deleted file mode 100755
index 063b67b..0000000
--- a/monitor.plt
+++ /dev/null
@@ -1,17 +0,0 @@
-set multiplot layout 4,4 title 'pafiX -- by H. Grosshans'
-f = 'output/monitor.dat'
-header=system("head -3 ".f." | tail -1")
-
-set key box opaque width 2 t l
-set xlabel offset 0,0.5
-set ylabel offset 2
-set linetype 1 ps 0.5 lc 3
-set xlabel word(header,3) #header columns are shifted by 1 bc of #
-
-do for [i=4:16] {
-set ylabel word(header,i+1)
-p f u 2:i notitle ls 1
-}
-
-pause 20
-reread
diff --git a/src/COPYING b/src/COPYING
deleted file mode 100755
index f288702..0000000
--- a/src/COPYING
+++ /dev/null
@@ -1,674 +0,0 @@
- GNU GENERAL PUBLIC LICENSE
- Version 3, 29 June 2007
-
- Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
- 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.
-
- <one line to give the program's name and a brief idea of what it does.>
- Copyright (C) <year> <name of author>
-
- 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 <https://www.gnu.org/licenses/>.
-
-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:
-
- <program> Copyright (C) <year> <name of author>
- 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
-<https://www.gnu.org/licenses/>.
-
- 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
-<https://www.gnu.org/licenses/why-not-lgpl.html>.
diff --git a/src/bc.f90 b/src/bc.f90
deleted file mode 100755
index 4c12a88..0000000
--- a/src/bc.f90
+++ /dev/null
@@ -1,264 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief Dirichlet (fixed value) boundary condition for velocity field
- subroutine bcUVW(myu,myv,myw)
- use var
- use parallel
- real(kind=pr),dimension(ii,jj,ll) :: myu,myv,myw
- real(kind=pr) :: wval
- integer :: m
-
- wval=0._pr
-
- call sync(myu); call sync(myv); call sync(myw) ! between procs through the sync routines
- if (bcx.eq.'w') then ! x:
- if (myid.eq.0) then
- do 1 m=1,gc
- myu(imin-m,:,:)= wval
- myv(imin-m,:,:)= myv(imin,:,:)-2._pr*(myv(imin,:,:)-wval)
- myw(imin-m,:,:)= myw(imin,:,:)-2._pr*(myw(imin,:,:)-wval)
-1 enddo
- endif
- if (myid.eq.nrprocs-1) then
- do 2 m=1,gc
- myu(imax-1+m,:,:)=wval
- myv(imax+m,:,:)= myv(imax,:,:)-2._pr*(myv(imax,:,:)-wval)
- myw(imax+m,:,:)= myw(imax,:,:)-2._pr*(myw(imax,:,:)-wval)
-2 enddo
- endif
- elseif (bcx.eq.'p') then
- ! between procs through the sync routines
- elseif (bcx.eq.'i') then
- ! see inflow routine
- if (myid.eq.nrprocs-1) then
- do 8 m=1,gc
- myu(imax+m-1,:,:)=myu(imax-1,:,:)
- myv(imax+m,:,:)= myv(imax,:,:)
- myw(imax+m,:,:)= myw(imax,:,:)
-8 enddo
- endif
- endif
-
- if (bcy.eq.'w') then ! y:
- do 3 m=1,gc
- myu(:,jmin-m,:)= myu(:,jmin,:)-2._pr*(myu(:,jmin,:)-wval)
- myu(:,jmax+m,:)= myu(:,jmax,:)-2._pr*(myu(:,jmax,:)-wval)
- myv(:,jmin-m,:)= wval
- myv(:,jmax-1+m,:)=wval
- myw(:,jmin-m,:)= myw(:,jmin,:)-2._pr*(myw(:,jmin,:)-wval)
- myw(:,jmax+m,:)= myw(:,jmax,:)-2._pr*(myw(:,jmax,:)-wval)
-3 enddo
- elseif (bcy.eq.'p') then
- do 4 m=1,gc
- myu(:,jmin-m,:)= myu(:,jmax+1-m,:)
- myu(:,jmax+m,:)= myu(:,jmin-1+m,:)
- myv(:,jmin-m,:)= myv(:,jmax+1-m,:)
- myv(:,jmax+m,:)= myv(:,jmin-1+m,:)
- myw(:,jmin-m,:)= myw(:,jmax+1-m,:)
- myw(:,jmax+m,:)= myw(:,jmin-1+m,:)
-4 enddo
- endif
-
- if (bcz.eq.'w') then ! z:
- do 5 m=1,gc
- myu(:,:,lmin-m)= myu(:,:,lmin)-2._pr*(myu(:,:,lmin)-wval)
- myu(:,:,lmax+m)= myu(:,:,lmax)-2._pr*(myu(:,:,lmax)-wval)
- myv(:,:,lmin-m)= myv(:,:,lmin)-2._pr*(myv(:,:,lmin)-wval)
- myv(:,:,lmax+m)= myv(:,:,lmax)-2._pr*(myv(:,:,lmax)-wval)
- myw(:,:,lmin-m)= wval
- myw(:,:,lmax-1+m)=wval
-5 enddo
- elseif (bcz.eq.'p') then
- do 6 m=1,gc
- myu(:,:,lmin-m)= myu(:,:,lmax+1-m)
- myu(:,:,lmax+m)= myu(:,:,lmin-1+m)
- myv(:,:,lmin-m)= myv(:,:,lmax+1-m)
- myv(:,:,lmax+m)= myv(:,:,lmin-1+m)
- myw(:,:,lmin-m)= myw(:,:,lmax+1-m)
- myw(:,:,lmax+m)= myw(:,:,lmin-1+m)
-6 enddo
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief inflow boundary condition
- subroutine inflow
- use var
- real(kind=pr),dimension(dimj*diml) :: random
- real(kind=pr) :: alpha,dist
- integer :: i,j,l,m,mm
-
- call random_number(random)
- alpha=.1_pr
- mm=0
- do l=lmin,lmax; do j=jmin,jmax
- mm=mm+1
- if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l)))
- if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
- if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
- do m=1,gc
- u(imin-m,j,l)=max(0._pr,ubulk0*(dist/delta)**(1._pr/6._pr) + (random(mm)-.5_pr)*alpha*ubulk0)
- v(imin-m,j,l)=(random(mm)-.5_pr)*alpha*ubulk0
- w(imin-m,j,l)=(random(mm)-.5_pr)*alpha*ubulk0
- enddo
- enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief Neumann (zero-gradient) boundary condition
- subroutine bcNeumann(myvar)
- use var
- use parallel
- real(kind=pr) :: myvar(ii,jj,ll)
- integer :: m
-
- call sync(myvar) ! between procs through the sync routines
- if (bcx.eq.'w') then
- if (myid.eq.0) then
- do 1 m=1,gc
- myvar(imin-m,:,:)= myvar(imin,:,:)
-1 enddo
- endif
- if (myid.eq.nrprocs-1) then
- do 2 m=1,gc
- myvar(imax+m,:,:)= myvar(imax,:,:)
-2 enddo
- endif
- elseif (bcx.eq.'p') then
- ! between procs through the sync routines
- elseif (bcx.eq.'i') then
- if (myid.eq.0) then
- do 7 m=1,gc
- myvar(imin-m,:,:)= myvar(imin,:,:)
-7 enddo
- endif
- if (myid.eq.nrprocs-1) then
- do 8 m=1,gc
- myvar(imax+m,:,:)= myvar(imax,:,:)
-8 enddo
- endif
- endif
-
- if (bcy.eq.'w') then
- do 3 m=1,gc
- myvar(:,jmin-m,:)= myvar(:,jmin,:)
- myvar(:,jmax+m,:)= myvar(:,jmax,:)
-3 enddo
- elseif (bcy.eq.'p') then
- do 4 m=1,gc
- myvar(:,jmin-m,:)= myvar(:,jmax+1-m,:)
- myvar(:,jmax+m,:)= myvar(:,jmin-1+m,:)
-4 enddo
- endif
-
- if (bcz.eq.'w') then
- do 5 m=1,gc
- myvar(:,:,lmin-m)= myvar(:,:,lmin)
- myvar(:,:,lmax+m)= myvar(:,:,lmax)
-5 enddo
- elseif (bcz.eq.'p') then
- do 6 m=1,gc
- myvar(:,:,lmin-m)= myvar(:,:,lmax+1-m)
- myvar(:,:,lmax+m)= myvar(:,:,lmin-1+m)
-6 enddo
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief Dirichlet (fixed value) boundary condition
- subroutine bcDirichlet(myvar,wval)
- use var
- use parallel
- real(kind=pr) :: myvar(ii,jj,ll)
- real(kind=pr) :: wval
- integer :: m
-
- call sync(myvar) ! between procs through the sync routines
- if (bcx.eq.'w') then
- if (myid.eq.0) then
- do 1 m=1,gc
- myvar(imin-m,:,:)= -myvar(imin,:,:)+2._pr*wval
-1 enddo
- endif
- if (myid.eq.nrprocs-1) then
- do 2 m=1,gc
- myvar(imax+m,:,:)= -myvar(imax,:,:)+2._pr*wval
-2 enddo
- endif
- elseif (bcx.eq.'p') then
- ! between procs through the sync routines
- elseif (bcx.eq.'i') then
- if (myid.eq.0) then
- do 7 m=1,gc
- myvar(imin-m,:,:)= myvar(imin,:,:)
-7 enddo
- endif
- if (myid.eq.nrprocs-1) then
- do 8 m=1,gc
- myvar(imax+m,:,:)= myvar(imax,:,:)
-8 enddo
- endif
- endif
-
- if (bcy.eq.'w') then
- do 3 m=1,gc
- myvar(:,jmin-m,:)= -myvar(:,jmin,:)+2._pr*wval
- myvar(:,jmax+m,:)= -myvar(:,jmax,:)+2._pr*wval
-3 enddo
- elseif (bcy.eq.'p') then
- do 4 m=1,gc
- myvar(:,jmin-m,:)= myvar(:,jmax+1-m,:)
- myvar(:,jmax+m,:)= myvar(:,jmin-1+m,:)
-4 enddo
- endif
-
- if (bcz.eq.'w') then
- do 5 m=1,gc
- myvar(:,:,lmin-m)= -myvar(:,:,lmin)+2._pr*wval
- myvar(:,:,lmax+m)= -myvar(:,:,lmax)+2._pr*wval
-5 enddo
- elseif (bcz.eq.'p') then
- do 6 m=1,gc
- myvar(:,:,lmin-m)= myvar(:,:,lmax+1-m)
- myvar(:,:,lmax+m)= myvar(:,:,lmin-1+m)
-6 enddo
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief a particle located on myid=0, i=imax imposes a source on
-!> myid=0, i=imax+1 which needs to be transferred to myid=1, i=imin
-!> required for rho_el,Fsx,Fsy,Fsz
-!> @todo if 2 dimensions periodic, the corners between diagonal
-!> procs not transferred yet
- subroutine bcLEsource(myvar)
- use var
- use parallel
- real(kind=pr) :: myvar(ii,jj,ll)
- integer :: m
-
- call syncLEsource(myvar) ! between procs through the sync routines
-
- if (bcy.eq.'p') then
- do 4 m=1,gc
- myvar(:,jmax+1-m,:)= myvar(:,jmin-m,:) + myvar(:,jmax+1-m,:)
- myvar(:,jmin-1+m,:)= myvar(:,jmax+m,:) + myvar(:,jmin-1+m,:)
-4 enddo
- endif
-
- if (bcz.eq.'p') then
- do 6 m=1,gc
- myvar(:,:,lmax+1-m)= myvar(:,:,lmin-m) + myvar(:,:,lmax+1-m)
- myvar(:,:,lmin-1+m)= myvar(:,:,lmax+m) + myvar(:,:,lmin-1+m)
-6 enddo
- endif
-
- end
diff --git a/src/electrostatics.f90 b/src/electrostatics.f90
deleted file mode 100755
index dfc8028..0000000
--- a/src/electrostatics.f90
+++ /dev/null
@@ -1,274 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief solve E field
- subroutine solveElectrostatics
- use var
- use parallel
- real(kind=pr) :: err,err00
- integer it
-
- if (pnd.eq.0._pr.or.syncMax(maxval(abs(q_el))).eq.0._pr) return
-
- call chargeDensity
-
- do it=1,itmax
- call calcPhi_el(err)
- if (it.eq.1) err00=max(err,1.e-19_pr)
- if (err/err00.lt.tol) exit
- enddo
-
- call calcE_el
-
- if (myid.eq.0) write(*,'(a,i3,a,es8.2e2,a)') 'E-pot. |L2 #',it-1,' = ',err/err00,' |'
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute volumetric charge density
- subroutine chargeDensity
- use var
- real(kind=pr),dimension(2,2,2) :: weight
- real(kind=pr) :: volE
- integer :: n,ibeg,iend,jbeg,jend,lbeg,lend
-
- rho_el=0._pr
-
- do 1 n=1,np
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)
- rho_el(ibeg:iend,jbeg:jend,lbeg:lend)= rho_el(ibeg:iend,jbeg:jend,lbeg:lend) + q_el(n)*partn(n)*weight/volE
-1 enddo
-
- call bcDirichlet(rho_el,0._pr)
- call bcLEsource(rho_el)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief one sweep to relax phi_el poisson equation using Jacobi method
-!> @param err: L2 error norm
-!> @todo Jacobi converges too slow
- subroutine calcPhi_el(err)
- use var
- use parallel
- real(kind=pr),dimension(ii,jj,ll) :: dphi_el
- real(kind=pr) :: diag,dex,dey,dez,err
- integer :: i,j,l
-
- dphi_el=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- diag= ( 2._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
- +2._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
- +2._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
-
- dex= (phi_el(i+1,j,l)*(xc(i)-xc(i-1))+phi_el(i-1,j,l)*(xc(i+1)-xc(i))) &
- /(0.5_pr*(xc(i+1)-xc(i-1))*(xc(i+1)-xc(i))*(xc(i)-xc(i-1)))
- dey= (phi_el(i,j+1,l)*(yc(j)-yc(j-1))+phi_el(i,j-1,l)*(yc(j+1)-yc(j))) &
- /(0.5_pr*(yc(j+1)-yc(j-1))*(yc(j+1)-yc(j))*(yc(j)-yc(j-1)))
- dez= (phi_el(i,j,l+1)*(zc(l)-zc(l-1))+phi_el(i,j,l-1)*(zc(l+1)-zc(l))) &
- /(0.5_pr*(zc(l+1)-zc(l-1))*(zc(l+1)-zc(l))*(zc(l)-zc(l-1)))
-
- dphi_el(i,j,l) = (rho_el(i,j,l)/eps_el + dex + dey + dez)/diag - phi_el(i,j,l)
- enddo; enddo; enddo
-
- phi_el=phi_el + dphi_el
-
- call bcDirichlet(phi_el,0._pr)
-
- err=sqrt(sum((vol*dphi_el)**2))/volTot
- err=syncSum(err)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief calculate E from phi_el, not an iterative process
- subroutine calcE_el
- use var
- integer :: i,j,l
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- Ex_el(i,j,l)= - (phi_el(i+1,j,l)-phi_el(i-1,j,l))/(xc(i+1)-xc(i-1))
- Ey_el(i,j,l)= - (phi_el(i,j+1,l)-phi_el(i,j-1,l))/(yc(j+1)-yc(j-1))
- Ez_el(i,j,l)= - (phi_el(i,j,l+1)-phi_el(i,j,l-1))/(zc(l+1)-zc(l-1))
- enddo; enddo; enddo
-
- call bcNeumann(Ex_el); call bcNeumann(Ey_el); call bcNeumann(Ez_el)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute specific electric forces on particles from E-field
-!> at seeding time-step E is not available, only Coulomb forces
- subroutine forcesGauss
- use var
- real(kind=pr),dimension(2,2,2) :: weight
- real(kind=pr) :: Extot,Eytot,Eztot,partmass,volE
- integer :: n,ibeg,iend,jbeg,jend,lbeg,lend
-
- do 1 n=1,np
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)
-
- Extot= sum(Ex_el(ibeg:iend,jbeg:jend,lbeg:lend)*weight)
- Eytot= sum(Ey_el(ibeg:iend,jbeg:jend,lbeg:lend)*weight)
- Eztot= sum(Ez_el(ibeg:iend,jbeg:jend,lbeg:lend)*weight)
-
- partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3
- fx_el(n)= q_el(n)*Extot/partmass
- fy_el(n)= q_el(n)*Eytot/partmass
- fz_el(n)= q_el(n)*Eztot/partmass
-1 enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute specific electric forces on particles from Coulomb's law
- subroutine forcesCoulomb
- use var
- use parallel
- integer :: n1,n2,m1,m2,i,j,l
-
- rtau_el_max=0._pr
-
- select case (elforceScheme)
-
- case (2) ! only Coulomb law
- fx_el=0._pr; fy_el=0._pr; fz_el=0._pr
- do 1 n1=2,np
- do 2 n2=1,(n1-1)
- if (n1.eq.n2) cycle
- call forcesCoulombN1N2(n1,n2)
-2 enddo
-1 enddo
-
- case (3) ! hybrid scheme, only particles in same cell
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- do 4 m1=2,npic(i,j,l)
- do 5 m2=1,(m1-1)
- n1=nic(i,j,l,m1)
- n2=nic(i,j,l,m2)
- call forcesCoulombN1N2(n1,n2)
-5 enddo
-4 enddo
- enddo; enddo; enddo
-
- end select
-
- rtau_el_max=syncMax(rtau_el_max)
-
- contains
-
-!> @brief compute Coulomb force of n2 acting on n1
- subroutine forcesCoulombN1N2(n1,n2)
- use var
- real(kind=pr) :: dis,dis3,disx,disy,disz,partmass1,partmass2,rtau_el
- integer :: n1,n2
-
- disx=xp(n1)-xp(n2)
- disy=yp(n1)-yp(n2)
- disz=zp(n1)-zp(n2)
- dis=sqrt(disx**2+disy**2+disz**2)
- dis3=dis**3
- partmass1=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3
- partmass2=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3
- rtau_el= sqrt(abs(q_el(n1)*q_el(n2))/4._pr/pi/eps_el/dis**2/min(partmass1,partmass2)/2._pr/dis) ! =sqrt(2*dis/acceleration)
- rtau_el_max=max(rtau_el_max,rtau_el)
- if (rtau_el.gt.(1._pr/dtNext)) dis3=dis3*(dtNext*rtau_el)**2 ! scale dis for stability
-
- fx_el(n1)= fx_el(n1) + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass1
- fy_el(n1)= fy_el(n1) + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass1
- fz_el(n1)= fz_el(n1) + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass1
-
- fx_el(n2)= fx_el(n2) - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass2
- fy_el(n2)= fy_el(n2) - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass2
- fz_el(n2)= fz_el(n2) - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass2
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute particle-wall charging
- subroutine chargeParticleWall(n,direction)
- use var
- real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp,random(2),random_normal
- integer :: n,direction
- character(70) :: filename
-
-
-! pipe:
-! ut2=up(n)**2
-! un2=vp(n)**2+wp(n)**2
-!
- if (direction.eq.1) then
- ut2=vp(n)**2+wp(n)**2
- un2=up(n)**2
- elseif (direction.eq.2) then
- ut2=up(n)**2+wp(n)**2
- un2=vp(n)**2
- elseif (direction.eq.3) then
- ut2=up(n)**2+vp(n)**2
- un2=wp(n)**2
- endif
-
-! John (1979)
-! alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* &
-! ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
-! AoAtot=alpha1/4._pr
-! dqp=Qaccfactor*AoAtot*(qpmax-q_el(n))
-
- call random_number(random)
- random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
- dqp=Qaccfactor*(qpmax-q_elNext(n))*(random_normal+1._pr) ! ...(mu*random_normal+sigma)
- q_elNext(n)=max(min(q_elNext(n)+dqp,qpmax),qp0)
-
-
- write(filename,'(a,i3.3,a)') 'results/particlesWall_p',myid,'.dat'
- open(access='append',unit=10,file=filename)
- write(10,77) nt,t,sqrt(ut2),sqrt(un2),dqp,q_elNext(n)
-77 format(i7,5(x,es11.3e2))
- close(10)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute particle-particle charging
- subroutine chargeParticleParticle(n1,n2)
- use var
- real(kind=pr) urel2,alpha1,dqp,AoAtot,random(2),random_normal
- integer :: n1,n2
- character(70) :: filename
-
-! urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
-!
-! if (urel2.eq.0._pr) then
-! dqp=0.5_pr*(q_el(n2)-q_el(n1))
-! else
-!! Soo (1971)
-! alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 &
-! *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) &
-! *(1._pr-nyp**2)/Ep)**(0.4_pr)
-! dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
-! endif
-! dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1))
-
- call random_number(random)
- random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
- dqp=Qaccfactor*(q_elNext(n2)-q_elNext(n1))*(random_normal+1._pr)
- q_elNext(n1)=max(min(q_elNext(n1)+dqp,qpmax),qp0)
- q_elNext(n2)=max(min(q_elNext(n2)-dqp,qpmax),qp0)
-
-
- write(filename,'(a,i3.3,a)') 'results/particlesPart_p',myid,'.dat'
- open(access='append',unit=10,file=filename)
- write(10,77) nt,t,q_elNext(n1),q_elNext(n2),dqp
-77 format(i7,4(x,es11.3e2))
- close(10)
-
- end
diff --git a/src/fluid.f90 b/src/fluid.f90
deleted file mode 100755
index c30d1e7..0000000
--- a/src/fluid.f90
+++ /dev/null
@@ -1,387 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief solves coupled pressure-velocity and calculates L2 norms
-!> @param err: L2 error norm
- subroutine solveFluid
- use var
- use parallel
- real(kind=pr),dimension(ii,jj,ll) :: du,dv,dw,dp,massRes
- real(kind=pr) :: err(4),err0(4),errL2
- integer itout,it,itmommax,itpCorr1max,itpCorr2max,velCorr1,velCorr2
-
-!> itmax: outer iterations
- itmommax=1; itpCorr1max=1; itpCorr2max=0 ! sequential SIMPLE
- !itmommax>1; itpCorr1max>1; itpCorr2max=0 ! SIMPLE
- !itmommax>1; itpCorr1max>1; itpCorr2max>1 ! PISO
- velCorr1=0; velCorr2=0 ! velocity correction yes/no (1/0)
-
- if ((bcx.eq.'i').and.(myid.eq.0)) call inflow
-
- do 1 itout=1,itmax
-
- if ((mod(it,10).eq.0).and.(npTot.ne.0)) call momentumCoupling
- if (turbModel.eq.'Smagorinsky') call turbSmagorinsky
-
- do 2 it=1,itmommax ! solve momentum eq
- call momx(du); call momy(dv); call momz(dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
- call bcUVW(u,v,w)
-2 enddo
-
- call mass(massRes,u,v,w) ! first pressure correction
- dp=0._pr
- do 3 it=1,itpCorr1max
- call pressureCorrection(dp,massRes)
- call bcNeumann(dp)
-3 enddo
- p= p+dp*urfp
-
- if (velCorr1.eq.1) then ! first velocity correction
- call velocityCorrection(du,dv,dw,dp)
- call bcUVW(du,dv,dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
- endif
-
- if (itpCorr2max.ge.1) then
- call mass(massRes,du,dv,dw) ! second pressure correction
- dp=0._pr
- do 4 it=1,itpCorr2max
- call pressureCorrection(dp,massRes)
- call bcNeumann(dp)
-4 enddo
- p=p+dp*urfp
-
- if (velCorr2.eq.1) then ! second velocity correction
- call velocityCorrection(du,dv,dw,dp)
- call bcUVW(du,dv,dw)
- u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
- endif
- endif
-
- err=sqrt([sum(volu*du),sum(volv*dv),sum(volw*dw),sum(vol*dp)]**2)/volTot
- err(1)=syncSum(err(1)); err(2)=syncSum(err(2)); err(3)=syncSum(err(3)); err(4)=syncSum(err(4))
- !if (myid.eq.0) write(*,'(x,a,4(es9.2e2))') 'res. inner it. u,v,w,p =',err
- if (itout.eq.1) err0=max(err,1.e-20_pr)
- errL2=sum(err/err0)/4._pr
- if (errL2.lt.tol) exit ! converged
-
-1 enddo
-
- if (myid.eq.0) write(*,'(a,i3,3(a,es8.2e2))') 'fluid |L2 #',itout-1,' = ',errL2, &
- ' |mom. = ',sum(err(1:3)/err0(1:3))/3._pr,' |pc. = ',err(4)/err0(4)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief pressure corrrection step
- subroutine pressureCorrection(dp,massRes)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: dp,ddp,massRes
- integer :: i,j,l
-
- ddp=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- ddp(i,j,l)= ( Cb1(i,j,l)*dp(i+1,j,l) + Cb2(i,j,l)*dp(i-1,j,l) &
- + Cc1(i,j,l)*dp(i,j+1,l) + Cc2(i,j,l)*dp(i,j-1,l) &
- + Cd1(i,j,l)*dp(i,j,l+1) + Cd2(i,j,l)*dp(i,j,l-1) &
- + massRes(i,j,l)*rhof/dt ) / Ca(i,j,l) - dp(i,j,l)
- enddo; enddo; enddo
- dp=dp+ddp
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine velocityCorrection(mydu,mydv,mydw,mydp)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: mydu,mydv,mydw,mydp
- real(kind=pr) :: ux,vy,wz
- integer :: i,j,l
-
- mydu=0._pr; mydv=0._pr; mydw=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- if (celltype(i+1,j,l).ne.wall) mydu(i,j,l)= dt/(xc(i+1)-xc(i))/rhof*(mydp(i+1,j,l)-mydp(i,j,l))
- if (celltype(i,j+1,l).ne.wall) mydv(i,j,l)= dt/(yc(j+1)-yc(j))/rhof*(mydp(i,j+1,l)-mydp(i,j,l))
- if (celltype(i,j,l+1).ne.wall) mydw(i,j,l)= dt/(zc(l+1)-zc(l))/rhof*(mydp(i,j,l+1)-mydp(i,j,l))
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute mass residual
- subroutine mass(massRes,myu,myv,myw)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: massRes,myu,myv,myw
- real(kind=pr) :: ux,vy,wz
- integer :: i,j,l
-
- massRes= 0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- if (celltype(i,j,l).ne.active) cycle
- ux= (myu(i,j,l)-myu(i-1,j,l))/(xf(i)-xf(i-1)) ! 2nd order central
- vy= (myv(i,j,l)-myv(i,j-1,l))/(yf(j)-yf(j-1))
- wz= (myw(i,j,l)-myw(i,j,l-1))/(zf(l)-zf(l-1))
- massRes(i,j,l)= -(ux+vy+wz)
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief momentum eq in x-direction, one sweep using Jacobi method
- subroutine momx(du)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: tu,qu,du
- real(kind=pr) :: ua,va,wa,px,uux,vuy,wuz,uxx,uyy,uzz,dudt,Cdudt,Cviscx
- integer :: i,j,l
-
- tu=0._pr; qu=0._pr; du=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- if (celltype(i+1,j,l).eq.wall) cycle
-
-! bilinear interpolation
- ua= u(i,j,l)
- va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i+1,j,l)+v(i+1,j-1,l))
- wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i+1,j,l)+w(i+1,j,l-1))
-
-! convective terms (2nd order central)
- uux= ua*(u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1))
- vuy= va*(u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
- wuz= wa*(u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
-
-! pressure gradient (2nd order central)
- px= (p(i+1,j,l)-p(i,j,l))/(xc(i+1)-xc(i))
-
-! viscous terms (2nd order central)
- uxx= (u(i+1,j,l)*(xf(i)-xf(i-1))+u(i-1,j,l)*(xf(i+1)-xf(i))-u(i,j,l)*(xf(i+1)-xf(i-1)))/ &
- (0.5_pr*(xf(i+1)-xf(i-1))*(xf(i+1)-xf(i))*(xf(i)-xf(i-1)))
- uyy= (u(i,j+1,l)*(yc(j)-yc(j-1))+u(i,j-1,l)*(yc(j+1)-yc(j))-u(i,j,l)*(yc(j+1)-yc(j-1)))/ &
- (0.5_pr*(yc(j+1)-yc(j-1))*(yc(j+1)-yc(j))*(yc(j)-yc(j-1)))
- uzz= (u(i,j,l+1)*(zc(l)-zc(l-1))+u(i,j,l-1)*(zc(l+1)-zc(l))-u(i,j,l)*(zc(l+1)-zc(l-1)))/ &
- (0.5_pr*(zc(l+1)-zc(l-1))*(zc(l+1)-zc(l))*(zc(l)-zc(l-1)))
-
-! time integration (2nd order)
- dudt= ((1._pr+tau01)*u(i,j,l)-(1._pr+tau01+tau02)*u01(i,j,l)+tau02*u02(i,j,l))/ &
- (2._pr*tau01*dt)
-
-! momentum eq -lhs+rhs (m/s**2)
- tu(i,j,l)= - dudt -(uux+vuy+wuz) -(px+dpdx)/rhof + (nuf+nuf_ru(i,j,l))*(uxx+uyy+uzz) + Fsx(i,j,l)
-
-! coefficients of u(i,j,l) of momentum eq
- Cdudt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscx= -2._pr*(nuf+nuf_ru(i,j,l)) * ( 1._pr/((xf(i+1)-xf(i))*(xf(i)-xf(i-1))) &
- + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
- + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
-
-! coefficients lhs-rhs
- qu(i,j,l)= Cdudt - Cviscx
-
- du(i,j,l)=tu(i,j,l)/qu(i,j,l)
-
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief momentum eq in y-direction, one sweep using Jacobi method
- subroutine momy(dv)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: tv,qv,dv
- real(kind=pr) :: ua,va,wa,py,uvx,vvy,wvz,vxx,vyy,vzz,dvdt,Cdvdt,Cviscy
- integer :: i,j,l
-
- tv=0._pr; qv=0._pr; dv=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- if (celltype(i,j+1,l).eq.wall) cycle
-
-! bilinear interpolation
- ua= 0.25_pr*(u(i,j,l)+u(i-1,j,l)+u(i,j+1,l)+u(i-1,j+1,l))
- va= v(i,j,l)
- wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i,j+1,l)+w(i,j+1,l-1))
-
-! convective terms (2nd order central)
- uvx= ua*(v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
- vvy= va*(v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
- wvz= wa*(v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
-
-! pressure gradient (2nd order)
- py=(p(i,j+1,l)-p(i,j,l))/(yc(j+1)-yc(j))
-
-! viscous terms (2nd order)
- vxx= (v(i+1,j,l)*(xc(i)-xc(i-1))+v(i-1,j,l)*(xc(i+1)-xc(i))-v(i,j,l)*(xc(i+1)-xc(i-1)))/ &
- (0.5_pr*(xc(i+1)-xc(i-1))*(xc(i+1)-xc(i))*(xc(i)-xc(i-1)))
- vyy= (v(i,j+1,l)*(yf(j)-yf(j-1))+v(i,j-1,l)*(yf(j+1)-yf(j))-v(i,j,l)*(yf(j+1)-yf(j-1)))/ &
- (0.5_pr*(yf(j+1)-yf(j-1))*(yf(j+1)-yf(j))*(yf(j)-yf(j-1)))
- vzz= (v(i,j,l+1)*(zc(l)-zc(l-1))+v(i,j,l-1)*(zc(l+1)-zc(l))-v(i,j,l)*(zc(l+1)-zc(l-1)))/ &
- (0.5_pr*(zc(l+1)-zc(l-1))*(zc(l+1)-zc(l))*(zc(l)-zc(l-1)))
-
-! time integration (2nd order)
- dvdt = ((1._pr+tau01)*v(i,j,l)-(1._pr+tau01+tau02)*v01(i,j,l)+tau02*v02(i,j,l))/ &
- (2._pr*tau01*dt)
-
-! momentum eq. -lhs+rhs
- tv(i,j,l)= - dvdt - (uvx+vvy+wvz) - py/rhof + (nuf+nuf_rv(i,j,l))*(vxx+vyy+vzz) + Fsy(i,j,l)
-
-! coefficients of v(i,j,l) of momentum eq
- Cdvdt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscy= -2._pr*(nuf+nuf_rv(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
- + 1._pr/((yf(j+1)-yf(j))*(yf(j)-yf(j-1))) &
- + 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
-
-! coefficients lhs-rhs
- qv(i,j,l)= Cdvdt - Cviscy
-
- dv(i,j,l)=tv(i,j,l)/qv(i,j,l)
-
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief momentum eq in z-direction, one sweep using Jacobi method
- subroutine momz(dw)
- use var
- real(kind=pr),dimension(ii,jj,ll) :: tw,qw,dw
- real(kind=pr) :: ua,va,wa,pz,uwx,vwy,wwz,wxx,wyy,wzz,dwdt,Cdwdt,Cviscz
- integer :: i,j,l
-
- tw=0._pr; qw=0._pr; dw=0._pr
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- if (celltype(i,j,l+1).eq.wall) cycle
-
-! bilinear interpolation
- ua= 0.25_pr*(u(i,j,l)+u(i-1,j,l)+u(i,j,l+1)+u(i-1,j,l+1))
- va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i,j,l+1)+v(i,j-1,l+1))
- wa= w(i,j,l)
-
-! convective terms (2nd order central)
- uwx= ua*(w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
- vwy= va*(w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
- wwz= wa*(w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
-
-! pressure gradient (2nd order)
- pz=(p(i,j,l+1)-p(i,j,l))/(zc(l+1)-zc(l))
-
-! viscous terms (2nd order)
- wxx= (w(i+1,j,l)*(xc(i)-xc(i-1))+w(i-1,j,l)*(xc(i+1)-xc(i))-w(i,j,l)*(xc(i+1)-xc(i-1)))/ &
- (0.5_pr*(xc(i+1)-xc(i-1))*(xc(i+1)-xc(i))*(xc(i)-xc(i-1)))
- wyy= (w(i,j+1,l)*(yc(j)-yc(j-1))+w(i,j-1,l)*(yc(j+1)-yc(j))-w(i,j,l)*(yc(j+1)-yc(j-1)))/ &
- (0.5_pr*(yc(j+1)-yc(j-1))*(yc(j+1)-yc(j))*(yc(j)-yc(j-1)))
- wzz= (w(i,j,l+1)*(zf(l)-zf(l-1))+w(i,j,l-1)*(zf(l+1)-zf(l))-w(i,j,l)*(zf(l+1)-zf(l-1)))/ &
- (0.5_pr*(zf(l+1)-zf(l-1))*(zf(l+1)-zf(l))*(zf(l)-zf(l-1)))
-
-! second-order time integration
- dwdt= ((1._pr+tau01)*w(i,j,l)-(1._pr+tau01+tau02)*w01(i,j,l)+tau02*w02(i,j,l))/ &
- (2._pr*tau01*dt)
-
-! momentum eq. -lhs+rhs
- tw(i,j,l)= - dwdt - (uwx+vwy+wwz) - pz/rhof + (nuf+nuf_rw(i,j,l))*(wxx+wyy+wzz) + Fsz(i,j,l)
-
-! coefficients of u(i,j,l) of momentum eq
- Cdwdt= (1._pr+tau01)/(2._pr*tau01*dt)
- Cviscz= -2._pr*(nuf+nuf_rw(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
- + 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
- + 1._pr/((zf(l+1)-zf(l))*(zf(l)-zf(l-1))) )
-
-! coefficients lhs-rhs
- qw(i,j,l)= Cdwdt - Cviscz
-
- dw(i,j,l)=tw(i,j,l)/qw(i,j,l)
-
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief calculate turbulent viscosity with static Smagorinsky model
- subroutine turbSmagorinsky
- use var
- real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,l_s
- integer :: i,j,l
-
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
-
- if (celltype(i+1,j,l).ne.wall) then
- ux= (u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1)) ! (2nd order central derivatives)
- uy= (u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
- uz= (u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
- vx= (v(i+1,j,l)+v(i+1,j-1,l)-v(i,j,l)-v(i,j-1,l))/(2._pr*(xc(i+1)-xc(i)))
- vy= (v(i+1,j,l)+v(i,j,l)-v(i+1,j-1,l)-v(i,j-1,l))/(2._pr*(yf(j)-yf(j-1)))
- vz= (v(i+1,j,l+1)+v(i,j,l+1)+v(i+1,j-1,l+1)+v(i,j-1,l+1)-v(i+1,j,l-1)-v(i,j,l-1)-v(i+1,j-1,l-1)-v(i,j-1,l)-1)/ &
- (4._pr*(zc(l+1)-zc(l-1)))
- wx= (w(i+1,j,l)+w(i+1,j,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(xc(i+1)-xc(i)))
- wy= (w(i+1,j+1,l)+w(i,j+1,l)+w(i+1,j+1,l-1)+w(i,j+1,l-1)-w(i+1,j-1,l)-w(i,j-1,l)-w(i+1,j-1,l-1)-w(i,j-1,l-1))/ &
- (4._pr*(yc(j+1)-yc(j-1)))
- wz= (w(i+1,j,l)+w(i,j,l)-w(i+1,j,l-1)-w(i,j,l-1))/(2._pr*(zf(l)-zf(l-1)))
-
- l_s= C_s*filteru(i,j,l)
- nuf_ru(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
- endif
-
- if (celltype(i,j+1,l).ne.wall) then
- ux= (u(i,j+1,l)+u(i,j,l)-u(i-1,j+1,l)-u(i-1,j,l))/(2._pr*(xf(i)-xf(i-1)))
- uy= (u(i,j+1,l)+u(i-1,j+1,l)-u(i,j,l)-u(i-1,j,l))/(2._pr*(yc(j+1)-yc(j)))
- uz= (u(i-1,j+1,l+1)+u(i,j+1,l+1)+u(i-1,j,l+1)+u(i,j,l+1)-u(i-1,j+1,l-1)-u(i,j+1,l-1)-u(i-1,j,l-1)-u(i,j,l-1))/ &
- (4._pr*(zc(l+1)-zc(l-1)))
- vx= (v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
- vy= (v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
- vz= (v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
- wx= (w(i+1,j,l)+w(i+1,j+1,l)+w(i+1,j,l-1)+w(i+1,j+1,l-1)-w(i-1,j,l)-w(i-1,j+1,l)-w(i-1,j,l-1)-w(i-1,j+1,l-1))/ &
- (4._pr*(xc(i+1)-xc(i-1)))
- wy= (w(i,j+1,l)+w(i,j+1,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(yc(j+1)-yc(j)))
- wz= (w(i,j,l)+w(i,j+1,l)-w(i,j,l-1)-w(i,j+1,l-1))/(2._pr*(zf(l)-zf(l-1)))
-
- l_s= C_s*filterv(i,j,l)
- nuf_rv(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
- endif
-
- if (celltype(i,j,l+1).ne.wall) then
- ux= (u(i,j,l)+u(i,j,l+1)-u(i-1,j,l)-u(i-1,j,l+1))/(2._pr*(xf(i)-xf(i-1)))
- uy= (u(i,j+1,l)+u(i,j+1,l+1)+u(i-1,j+1,l)+u(i-1,j+1,l+1)-u(i,j-1,l)-u(i,j-1,l+1)-u(i-1,j-1,l)-u(i-1,j-1,l+1))/ &
- (4._pr*(yc(j+1)-yc(j-1)))
- uz= (u(i,j,l+1)+u(i-1,j,l+1)-u(i,j,l)-u(i-1,j,l))/(2._pr*(zc(l+1)-zc(l)))
- vx= (v(i+1,j,l)+v(i+1,j,l+1)+v(i+1,j-1,l)+v(i+1,j-1,l+1)-v(i-1,j,l)-v(i-1,j,l+1)-v(i-1,j-1,l)-v(i-1,j-1,l+1))/ &
- (4._pr*(xc(i+1)-xc(i-1)))
- vy= (v(i,j,l)+v(i,j,l+1)-v(i,j-1,l)-v(i,j-1,l+1))/(2._pr*(yf(j)-yf(j-1)))
- vz= (v(i,j,l+1)+v(i,j-1,l+1)-v(i,j,l)-v(i,j-1,l))/(2._pr*(yc(i+1)-yc(i)))
- wx= (w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
- wy= (w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
- wz= (w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
-
- l_s= C_s*filterw(i,j,l)
- nuf_rw(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
- endif
-
- enddo; enddo; enddo
-
- contains
-
- real(kind=pr) function S(ux,uy,uz,vx,vy,vz,wx,wy,wz)
- !> @brief calculate characteristic filtered rate of strain, Pope Eq. (13.74)
- real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,Sij(3,3)
-
- Sij(1,1)= ux !< filtered rate of strain tensor, Pope Eq. (2.70)
- Sij(2,2)= vy
- Sij(3,3)= wz
- Sij(1,2)= 0.5_pr*(uy+vx)
- Sij(2,1)= Sij(1,2)
- Sij(1,3)= 0.5_pr*(uz+wx)
- Sij(3,1)= Sij(1,3)
- Sij(2,3)= 0.5_pr*(vz+wy)
- Sij(3,2)= Sij(2,3)
- S= sqrt(2._pr*sum(Sij**2))
-
- end
-
- end
diff --git a/src/main.f90 b/src/main.f90
deleted file mode 100755
index d042447..0000000
--- a/src/main.f90
+++ /dev/null
@@ -1,53 +0,0 @@
-! Copyright 2015-2021 Holger Grosshans
-!
-! 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 <https://www.gnu.org/licenses/>.
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief pafiX - particle flow simulation in eXplosion protection
- program pafiX
- use var
- use parallel
- use mpi
- implicit none
-
- call initParallel
- if (myid.eq.0) write(*,'(/a/a)') repeat('#',72),version
- if (myid.eq.0) write(*,'(a/a)') repeat('#',72),'Pre-processing'
- call preProcessing
-
- if (myid.eq.0) write(*,'(a/a)') repeat('#',72),'Start computation'
- do 10 nt=ntstart,ntend
- if (myid.eq.0) write(*,'(a/a,i7,x,a,i7,x,a)') repeat('#',72),'nt=',nt,'of',ntend,'time-steps'
-
- call syncCheck
- call prepareTimestep ! nt-1 to nt
- call solveFluid ! nt-1 to nt
- call solveElectrostatics ! nt
- call nextTimestepSize ! nt to nt+1
- call solveParticles ! nt to nt+1
- call postProcessing ! nt
-
- timecom(1)=syncSum(timecom(1))/nrprocs
- timecom(9)=real(mpi_wtime(),kind=pr)-timebeg+timecom(10)
- if (myid.eq.0) write(*,'(3(a,es8.2e2))') &
- 'time (s) |physical= ',t,' |comput. = ',timecom(9),' |MPI com = ',timecom(1)
- call flush
-
-10 enddo
-
- if (myid.eq.0) write(*,'(a/a/)') repeat('#',72),'Calculations completed!'
- call mpi_finalize(mpierr)
-
- end program pafiX
diff --git a/src/makefile b/src/makefile
deleted file mode 100755
index 29d2aed..0000000
--- a/src/makefile
+++ /dev/null
@@ -1,51 +0,0 @@
-#%.o : %.mod
-
-OBJ = var.o parallel.o main.o pre.o bc.o restart.o fluid.o \
- timestep.o post.o particles.o particlesTransport.o \
- electrostatics.o writedat_particles.o write_vtk.o
-INC =
-
-
-# gcc:
-CMP = mpifort
-#O3: optimization
-FLAGS = -O3 -mcmodel=medium
-#debugging:
-#FLAGS = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
-#-Wall
-FFLAGS = -c $(FLAGS)
-
-
-newfu: $(OBJ)
- $(CMP) $(FLAGS) -o pafiX $(OBJ)
-var.o: var.f90 $(INC)
- $(CMP) $(FFLAGS) var.f90
-parallel.o: parallel.f90 $(INC)
- $(CMP) $(FFLAGS) parallel.f90
-main.o: main.f90 $(INC)
- $(CMP) $(FFLAGS) main.f90
-pre.o: pre.f90 $(INC)
- $(CMP) $(FFLAGS) pre.f90
-post.o: post.f90 $(INC)
- $(CMP) $(FFLAGS) post.f90
-bc.o: bc.f90 $(INC)
- $(CMP) $(FFLAGS) bc.f90
-restart.o: restart.f90 $(INC)
- $(CMP) $(FFLAGS) restart.f90
-fluid.o: fluid.f90 $(INC)
- $(CMP) $(FFLAGS) fluid.f90
-timestep.o: timestep.f90 $(INC)
- $(CMP) $(FFLAGS) timestep.f90
-particles.o: particles.f90 $(INC)
- $(CMP) $(FFLAGS) particles.f90
-particlesTransport.o: particlesTransport.f90 $(INC)
- $(CMP) $(FFLAGS) particlesTransport.f90
-electrostatics.o: electrostatics.f90 $(INC)
- $(CMP) $(FFLAGS) electrostatics.f90
-write_vtk.o: write_vtk.f90 $(INC)
- $(CMP) $(FFLAGS) write_vtk.f90
-writedat_particles.o: writedat_particles.f90 $(INC)
- $(CMP) $(FFLAGS) writedat_particles.f90
-
-clean:
- rm *.o *.mod
diff --git a/src/parallel.f90 b/src/parallel.f90
deleted file mode 100755
index cb26786..0000000
--- a/src/parallel.f90
+++ /dev/null
@@ -1,340 +0,0 @@
-module parallel
-implicit none
-
-contains
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief initialize MPI
- subroutine initParallel
- use var
- use mpi
-
- call mpi_init(mpierr)
- call mpi_comm_size(mpi_comm_world,nrprocs,mpierr)
- call mpi_comm_rank(mpi_comm_world,myid,mpierr)
- call mpi_type_create_f90_real(precision,mpi_undefined,mpi_pr,mpierr)
-
- timenow=real(mpi_wtime(),kind=pr)
- timebeg=real(mpi_wtime(),kind=pr)
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief synchronize particles between processors
- subroutine syncPart(np_sendl,np_sendr,np_recvl,np_recvr,myvar)
- use var
- use mpi
- real(kind=pr),allocatable,dimension(:) :: p_sendl,p_sendr
- real(kind=pr),dimension(maxnp) :: myvar,p_recvl,p_recvr
- integer,dimension(maxnp) :: np_sendl,np_sendr,np_recvl,np_recvr
- integer rsleft,rsright,rrleft,rrright,n_sendl,n_sendr,n_recvl,n_recvr,n
-
- timenow=real(mpi_wtime(),kind=pr)
-
- n_sendl=sum(np_sendl(1:np))
- p_sendl=pack(myvar(1:np),np_sendl(1:np).eq.1)
- n_sendr=sum(np_sendr(1:np))
- p_sendr=pack(myvar(1:np),np_sendr(1:np).eq.1)
-
-! send/receive number of particles
- call mpi_isend(n_sendl,1,mpi_integer,prev,2,mpi_comm_world,rsleft,mpierr)
- call mpi_irecv(n_recvl,1,mpi_integer,prev,1,mpi_comm_world,rrleft,mpierr)
- call mpi_isend(n_sendr,1,mpi_integer,next,1,mpi_comm_world,rsright,mpierr)
- call mpi_irecv(n_recvr,1,mpi_integer,next,2,mpi_comm_world,rrright,mpierr)
- call mpi_wait(rsleft,mpistatus,mpierr)
- call mpi_wait(rrleft,mpistatus,mpierr)
- call mpi_wait(rsright,mpistatus,mpierr)
- call mpi_wait(rrright,mpistatus,mpierr)
- if (next.eq.mpi_proc_null) n_recvr=0
- if (prev.eq.mpi_proc_null) n_recvl=0
-
- if (n_sendl.gt.0) call mpi_isend(p_sendl,n_sendl,mpi_pr,prev,2,mpi_comm_world,rsleft,mpierr)
- if (n_recvl.gt.0) call mpi_irecv(p_recvl,n_recvl,mpi_pr,prev,1,mpi_comm_world,rrleft,mpierr)
- if (n_sendr.gt.0) call mpi_isend(p_sendr,n_sendr,mpi_pr,next,1,mpi_comm_world,rsright,mpierr)
- if (n_recvr.gt.0) call mpi_irecv(p_recvr,n_recvr,mpi_pr,next,2,mpi_comm_world,rrright,mpierr)
- if (n_sendl.gt.0) call mpi_wait(rsleft,mpistatus,mpierr)
- if (n_recvl.gt.0) call mpi_wait(rrleft,mpistatus,mpierr)
- if (n_sendr.gt.0) call mpi_wait(rsright,mpistatus,mpierr)
- if (n_recvr.gt.0) call mpi_wait(rrright,mpistatus,mpierr)
-
- np=npp
- np_recvl=0
- np_recvr=0
-
- np_recvl( np+1 : np+n_recvl ) = 1
- myvar( np+1 : np+n_recvl ) = p_recvl( 1 : n_recvl )
-
- np_recvr( np+n_recvl+1 : np+n_recvl+n_recvr ) = 1
- myvar( np+n_recvl+1 : np+n_recvl+n_recvr ) = p_recvr( 1 : n_recvr )
-
- np=np+n_recvl+n_recvr
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> this should be merged with the previous routine using a select type construct
- subroutine syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,myvar)
- use var
- use mpi
- integer, allocatable, dimension(:) :: p_sendl,p_sendr
- integer, dimension(maxnp) :: myvar,p_recvl,p_recvr
- integer, dimension(maxnp) :: np_sendl,np_sendr,np_recvl,np_recvr
- integer rsleft,rsright,rrleft,rrright,n_sendl,n_sendr,n_recvl,n_recvr,n
-
- timenow=real(mpi_wtime(),kind=pr)
-
- n_sendl=sum(np_sendl(1:np))
- p_sendl=pack(myvar(1:np),np_sendl(1:np).eq.1)
- n_sendr=sum(np_sendr(1:np))
- p_sendr=pack(myvar(1:np),np_sendr(1:np).eq.1)
-
-! send/receive number of particles
- call mpi_isend(n_sendl,1,mpi_integer,prev,2,mpi_comm_world,rsleft,mpierr)
- call mpi_irecv(n_recvl,1,mpi_integer,prev,1,mpi_comm_world,rrleft,mpierr)
- call mpi_isend(n_sendr,1,mpi_integer,next,1,mpi_comm_world,rsright,mpierr)
- call mpi_irecv(n_recvr,1,mpi_integer,next,2,mpi_comm_world,rrright,mpierr)
- call mpi_wait(rsleft,mpistatus,mpierr)
- call mpi_wait(rrleft,mpistatus,mpierr)
- call mpi_wait(rsright,mpistatus,mpierr)
- call mpi_wait(rrright,mpistatus,mpierr)
- if (next.eq.mpi_proc_null) n_recvr=0
- if (prev.eq.mpi_proc_null) n_recvl=0
-
- if (n_sendl.gt.0) call mpi_isend(p_sendl,n_sendl,mpi_integer,prev,2,mpi_comm_world,rsleft,mpierr)
- if (n_recvl.gt.0) call mpi_irecv(p_recvl,n_recvl,mpi_integer,prev,1,mpi_comm_world,rrleft,mpierr)
- if (n_sendr.gt.0) call mpi_isend(p_sendr,n_sendr,mpi_integer,next,1,mpi_comm_world,rsright,mpierr)
- if (n_recvr.gt.0) call mpi_irecv(p_recvr,n_recvr,mpi_integer,next,2,mpi_comm_world,rrright,mpierr)
- if (n_sendl.gt.0) call mpi_wait(rsleft,mpistatus,mpierr)
- if (n_recvl.gt.0) call mpi_wait(rrleft,mpistatus,mpierr)
- if (n_sendr.gt.0) call mpi_wait(rsright,mpistatus,mpierr)
- if (n_recvr.gt.0) call mpi_wait(rrright,mpistatus,mpierr)
-
- np=npp
- np_recvl=0
- np_recvr=0
-
- np_recvl( np+1 : np+n_recvl ) = 1
- myvar( np+1 : np+n_recvl ) = p_recvl( 1 : n_recvl )
-
- np_recvr( np+n_recvl+1 : np+n_recvl+n_recvr ) = 1
- myvar( np+n_recvl+1 : np+n_recvl+n_recvr ) = p_recvr( 1 : n_recvr )
-
- np=np+n_recvl+n_recvr
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief synchronize fluid between processors
- subroutine sync(myvar)
- use var
- use mpi
- real(kind=pr), dimension(ii,jj,ll) :: myvar
- real(kind=pr), dimension(gc*jj*ll) :: sendleft,sendright,recvleft,recvright
- integer :: rsleft,rsright,rrleft,rrright
- integer :: n
-
- timenow=real(mpi_wtime(),kind=pr)
-
- n=gc*jj*ll
- sendleft= reshape(myvar(imin:imin+gc-1,1:jj,1:ll),(/n/))
- sendright=reshape(myvar(imax-gc+1:imax,1:jj,1:ll),(/n/))
-
- call mpi_isend(sendleft,n,mpi_pr,prev,1,mpi_comm_world,rsleft,mpierr)
- call mpi_isend(sendright,n,mpi_pr,next,2,mpi_comm_world,rsright,mpierr)
- call mpi_irecv(recvright,n,mpi_pr,next,1,mpi_comm_world,rrright,mpierr)
- call mpi_irecv(recvleft,n,mpi_pr,prev,2,mpi_comm_world,rrleft,mpierr)
- call mpi_wait(rsleft,mpistatus,mpierr)
- call mpi_wait(rsright,mpistatus,mpierr)
- call mpi_wait(rrright,mpistatus,mpierr)
- call mpi_wait(rrleft,mpistatus,mpierr)
-
- if (next.ne.mpi_proc_null) myvar(imax+1:ii,1:jj,1:ll)=reshape(recvright,(/gc,jj,ll/))
- if (prev.ne.mpi_proc_null) myvar(1:gc,1:jj,1:ll)= reshape(recvleft,(/gc,jj,ll/))
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief synchronize LE sources between processors, see 'bcELsource'
- subroutine syncLEsource(myvar)
- use var
- use mpi
- real(kind=pr), dimension(ii,jj,ll) :: myvar
- real(kind=pr), dimension(gc*jj*ll) :: sendleft,sendright,recvleft,recvright
- integer :: rsleft,rsright,rrleft,rrright
- integer :: i,j,l,n
-
- timenow=real(mpi_wtime(),kind=pr)
-
- n=gc*jj*ll
- sendleft= reshape(myvar(1:gc,1:jj,1:ll),(/n/))
- sendright=reshape(myvar(imax+1:ii,1:jj,1:ll),(/n/))
-
- call mpi_isend(sendleft,n,mpi_pr,prev,1,mpi_comm_world,rsleft,mpierr)
- call mpi_isend(sendright,n,mpi_pr,next,2,mpi_comm_world,rsright,mpierr)
- call mpi_irecv(recvright,n,mpi_pr,next,1,mpi_comm_world,rrright,mpierr)
- call mpi_irecv(recvleft,n,mpi_pr,prev,2,mpi_comm_world,rrleft,mpierr)
- call mpi_wait(rsleft,mpistatus,mpierr)
- call mpi_wait(rsright,mpistatus,mpierr)
- call mpi_wait(rrright,mpistatus,mpierr)
- call mpi_wait(rrleft,mpistatus,mpierr)
-
- if (next.ne.mpi_proc_null) myvar(imax-gc+1:imax,1:jj,1:ll)=reshape(recvright,(/gc,jj,ll/))
- if (prev.ne.mpi_proc_null) myvar(imin:imin+gc-1,1:jj,1:ll)=reshape(recvleft,(/gc,jj,ll/))
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute the max of a scalar over all processors
- real(kind=pr) function syncMax(myvar)
- use var
- use mpi
- real(kind=pr) :: myvar
-
- timenow=real(mpi_wtime(),kind=pr)
-
-! mpi_allreduce may not work with intel compiler
- call mpi_allreduce(myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr)
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute the min of a scalar over all processors
- real(kind=pr) function syncMin(myvar)
- use var
- use mpi
- real(kind=pr) :: myvar
-
- timenow=real(mpi_wtime(),kind=pr)
-
-! mpi_allreduce may not work with intel compiler
- call mpi_allreduce(myvar,syncMin,1,mpi_pr,mpi_min,mpi_comm_world,mpierr)
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute the sum of a real scalar over all processors
- real(kind=pr) function syncSum(myvar)
- use var
- use mpi
- real(kind=pr),intent(in) :: myvar
-
- timenow=real(mpi_wtime(),kind=pr)
-
-! mpi_allreduce does not work with intel compiler
- call mpi_allreduce(myvar,syncSum,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr)
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end function
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute the sum of a integer scalar over all processors
- integer function syncSumI(myvar)
- use var
- use mpi
- integer :: myvar
-
- timenow=real(mpi_wtime(),kind=pr)
-
-! mpi_allreduce does not work with intel compiler
- call mpi_allreduce(myvar,syncSumI,1,mpi_integer,mpi_sum,mpi_comm_world,mpierr)
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief checks if processors run synchronous
- subroutine syncCheck
- use var
- use mpi
- integer :: proc,rs,rr,ntproc
-
- timenow=real(mpi_wtime(),kind=pr)
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call mpi_irecv(ntproc,1,mpi_integer,proc,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- if (ntproc.ne.nt) then
- write(*,'(a)') 'Processors not synchronous!'
- stop
- endif
-1 enddo
- endif
-
- if (myid.ne.0) then
- call mpi_isend(nt,1,mpi_integer,0,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
- endif
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief generate m independent random numbers for each processor
-!> @param m number of elements required per processor
-!> @return random random numbers between 0 and 1
- subroutine syncRandom(m,random)
- use var
- use mpi
- real(kind=pr),dimension(m) :: random
- integer :: m,proc,rs,rr
-
- timenow=real(mpi_wtime(),kind=pr)
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call random_number(random)
- call mpi_isend(random,m,mpi_pr,proc,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
-1 enddo
- call random_number(random) !in the end, myid=0 generates numbers for itself
- endif
-
- if (myid.ne.0) then
- call mpi_irecv(random,m,mpi_pr,0,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- endif
-
- timeend=real(mpi_wtime(),kind=pr)
- timecom(1)=timecom(1)+timeend-timenow
-
- end
-
-end module parallel
diff --git a/src/particles.f90 b/src/particles.f90
deleted file mode 100755
index bfcb874..0000000
--- a/src/particles.f90
+++ /dev/null
@@ -1,300 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute weigths for interpolation from Lagrangian to
-!> 8 Eulerian points using a Gauss distribution function
-!> param weight weights of involved cells
-!> param volE volume of the involved cells
-!> param n particle number
-!> param direction 0 if Eulerian variable is stored in the cell
-!> center, 1-3 if on the x,y or z-face
- subroutine weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction)
- use var
- integer,intent(in) :: n,direction
- real(kind=pr),dimension(2,2,2),intent(out) :: weight
- integer,intent(out) :: ibeg,iend,jbeg,jend,lbeg,lend
- real(kind=pr) :: deltaf2,deltaf,dis2,volE
- integer :: i,j,l
-
-! define distribution block and correct for staggered variables
- if ((direction.eq.1).or.(xp(n).lt.xc(ip(n)))) then
- ibeg=ip(n)-1
- else
- ibeg=ip(n)
- endif
- iend=ibeg+1
-
- if ((direction.eq.2).or.(yp(n).lt.yc(jp(n)))) then
- jbeg=jp(n)-1
- else
- jbeg=jp(n)
- endif
- jend=jbeg+1
-
- if ((direction.eq.3).or.(zp(n).lt.zc(lp(n)))) then
- lbeg=lp(n)-1
- else
- lbeg=lp(n)
- endif
- lend=lbeg+1
-
-! volume on the Eulerian grid
- if (direction.eq.0) then
- volE= sum(vol(ibeg:iend,jbeg:jend,lbeg:lend))
- else if (direction.eq.1) then
- volE= sum(volu(ibeg:iend,jbeg:jend,lbeg:lend))
- else if (direction.eq.2) then
- volE= sum(volv(ibeg:iend,jbeg:jend,lbeg:lend))
- else if (direction.eq.3) then
- volE= sum(volw(ibeg:iend,jbeg:jend,lbeg:lend))
- endif
-
- deltaf=xc(iend)-xc(ibeg)
- deltaf2=(deltaf**2)/6._pr
-
-!compute weights
- do l=lbeg,lend; do j=jbeg,jend; do i=ibeg,iend
- if (direction.eq.0) then
- dis2= (xp(n)-xc(i))**2 + (yp(n)-yc(j))**2 + (zp(n)-zc(l))**2
- elseif (direction.eq.1) then
- dis2= (xp(n)-xf(i))**2 + (yp(n)-yc(j))**2 + (zp(n)-zc(l))**2
- elseif (direction.eq.2) then
- dis2= (xp(n)-xc(i))**2 + (yp(n)-yf(j))**2 + (zp(n)-zc(l))**2
- elseif (direction.eq.3) then
- dis2= (xp(n)-xc(i))**2 + (yp(n)-yc(j))**2 + (zp(n)-zf(l))**2
- endif
- weight(i+1-ibeg,j+1-jbeg,l+1-lbeg)=exp(-dis2/deltaf2)
- enddo; enddo; enddo
-
- weight=weight/sum(weight) ! normalized
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief allocate all public particle arrays
- subroutine allocateParticleArrays
- use var
- use parallel
- integer :: nppTot
-
- npTot=syncSumI(np)
- nppTot=syncSumI(npp)
- maxnp=npTot*3 ! one particle can be in domain, gc and sync storage
-
- call allocateParticleArray(nppTot,up)
- call allocateParticleArray(nppTot,vp)
- call allocateParticleArray(nppTot,wp)
- call allocateParticleArray(nppTot,upNext)
- call allocateParticleArray(nppTot,vpNext)
- call allocateParticleArray(nppTot,wpNext)
- call allocateParticleArray(nppTot,uf)
- call allocateParticleArray(nppTot,vf)
- call allocateParticleArray(nppTot,wf)
- call allocateParticleArray(nppTot,dufdy)
- call allocateParticleArray(nppTot,dufdz)
- call allocateParticleArray(nppTot,xp)
- call allocateParticleArray(nppTot,yp)
- call allocateParticleArray(nppTot,zp)
- call allocateParticleArray(nppTot,xpNext)
- call allocateParticleArray(nppTot,ypNext)
- call allocateParticleArray(nppTot,zpNext)
- call allocateParticleArray(nppTot,radp)
- call allocateParticleArray(nppTot,q_el)
- call allocateParticleArray(nppTot,q_elNext)
- call allocateParticleArray(nppTot,fx_el)
- call allocateParticleArray(nppTot,fy_el)
- call allocateParticleArray(nppTot,fz_el)
- call allocateParticleArray(nppTot,fx_d)
- call allocateParticleArray(nppTot,fy_d)
- call allocateParticleArray(nppTot,fz_d)
- call allocateParticleArray(nppTot,fx_l)
- call allocateParticleArray(nppTot,fy_l)
- call allocateParticleArray(nppTot,fz_l)
- call allocateParticleArray(nppTot,partn)
-
- call allocateParticleArrayI(nppTot,ip)
- call allocateParticleArrayI(nppTot,jp)
- call allocateParticleArrayI(nppTot,lp)
- call allocateParticleArrayI(nppTot,wcollnum)
- call allocateParticleArrayI(nppTot,ppcollnum)
- call allocateParticleArrayI(nppTot,nGlob)
-
- contains
-
- subroutine allocateParticleArray(nppTot,myvar)
- use var
- real(kind=pr), allocatable, dimension(:) :: myvar,tmyvar
- integer :: nppTot
-
- if (nppTot.ge.1) then
- allocate(tmyvar(maxnp))
- tmyvar(1:npp)=myvar(1:npp)
- endif
- if (allocated(myvar)) deallocate(myvar)
- allocate(myvar(maxnp))
- myvar=0._pr
- if (nppTot.ge.1) then
- myvar(1:npp)=tmyvar(1:npp)
- endif
-
- end
-
- subroutine allocateParticleArrayI(nppTot,myvar)
- use var
- integer, allocatable, dimension(:) :: myvar,tmyvar
- integer :: nppTot
-
- if (nppTot.ge.1) then
- allocate(tmyvar(maxnp))
- tmyvar(1:npp)=myvar(1:npp)
- endif
- if (allocated(myvar)) deallocate(myvar)
- allocate(myvar(maxnp))
- myvar=0
- if (nppTot.ge.1) then
- myvar(1:npp)=tmyvar(1:npp)
- endif
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief send particles to next processor, all required variables for further transport
- subroutine particlesNextProc
- use var
- use parallel
- use mpi
- integer, dimension (maxnp) :: np_sendr,np_sendl,np_recvl,np_recvr
- integer :: n,m
-
- np_sendl=0
- np_sendr=0
-
- do 1 n=1,np
-
- if (prev.ne.mpi_proc_null) then
- if (xp(n).le.xmin) then
- np_sendl(n)=1
- else
- np_sendl(n)=0
- endif
- endif
- if (next.ne.mpi_proc_null) then
- if (xp(n).ge.xmax) then
- np_sendr(n)=1
- else
- np_sendr(n)=0
- endif
- endif
-
-1 enddo
-
- npp=np
- np_recvl=0
- np_recvr=0
-
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,xp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,yp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,zp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,xpNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,ypNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,zpNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,radp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,up)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wp)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,upNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vpNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wpNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,q_el)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,q_elNext)
- call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,partn)
- call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,wcollnum)
- call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,ppcollnum)
- call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,nGlob)
-
- if ((bcx.eq.'p').and.(myid.eq.nrprocs-1)) then
- do 2 n=npp+1,np
- if (np_recvr(n).eq.1) xp(n)= xp(n) + dimxtot
-2 enddo
- endif
- if ((bcx.eq.'p').and.(myid.eq.0)) then
- do 3 n=npp+1,np
- if (np_recvl(n).eq.1) xp(n)= xp(n) - dimxtot
-3 enddo
- endif
-
- m=0 ! remove sent particles
- do 4 n=1,np
- if ((xp(n).ge.xmin).and.(xp(n).le.xmax)) then
- m=m+1
- call partN2M(n,m)
- endif
-4 enddo
- np=m
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief move particle from storage position n to m, all required variables for further transport
- subroutine partN2M(n,m)
- use var
- integer :: n,m
-
- xp(m)=xp(n)
- yp(m)=yp(n)
- zp(m)=zp(n)
- xpNext(m)=xpNext(n)
- ypNext(m)=ypNext(n)
- zpNext(m)=zpNext(n)
- radp(m)=radp(n)
- up(m)=up(n)
- vp(m)=vp(n)
- wp(m)=wp(n)
- upNext(m)=upNext(n)
- vpNext(m)=vpNext(n)
- wpNext(m)=wpNext(n)
- q_el(m)=q_el(n)
- q_elNext(m)=q_elNext(n)
- partn(m)=partn(n)
-
- wcollnum(m)=wcollnum(n)
- ppcollnum(m)=ppcollnum(n)
- nGlob(m)=nGlob(n)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief locate particles on Eulerian grid, required after moving
-!> particles: nextTimestep, particlesNextProc, particlesSeed
- subroutine particlesToCells
- use var
- integer :: n
-
- npic=0
- do 1 n=1,np
- ip(n)=minloc(xf, dim=1, mask=(xp(n).lt.xf))
- jp(n)=minloc(yf, dim=1, mask=(yp(n).lt.yf))
- lp(n)=minloc(zf, dim=1, mask=(zp(n).lt.zf))
-
- if (celltype(ip(n),jp(n),lp(n)).ne.active) then
- write(*,*) 'particle lost, proc=',myid,', xp=',xp(n),', yp=',yp(n),', zp=',zp(n)
- stop
- endif
-
- npic(ip(n),jp(n),lp(n))=npic(ip(n),jp(n),lp(n))+1
-1 enddo
-
- if (allocated(nic)) deallocate(nic)
- allocate(nic(ii,jj,ll,maxval(npic)))
-
- nic=0
- do 2 n=1,np
- nic(ip(n),jp(n),lp(n),count(nic(ip(n),jp(n),lp(n),:).ne.0)+1)=n
-2 enddo
-
- end
diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90
deleted file mode 100755
index b8428d4..0000000
--- a/src/particlesTransport.f90
+++ /dev/null
@@ -1,493 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief solve particulate phase
- subroutine solveParticles
- use var
-
- if (pnd.ne.0.or.npTot.ne.0) then
- if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
- call particlesDrag
- call particlesLift
- if ((elForceScheme.eq.1).or.(elForceScheme.eq.3)) call forcesGauss
- if ((elForceScheme.eq.2).or.(elForceScheme.eq.3)) call forcesCoulomb
- call particlesVelocityNext
- call particlesCollideNext
- call particlesTransportNext
- if (myid.eq.0) write(*,'(a,i8,2(a,es8.2e2))') &
- 'particles |transp. = ',npTot,' |dt/t_el = ',dtNext*rtau_el_max,' |dt/t_p = ',dtNext/tau_p_min
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief explicit first-order Euler forward integration of particle
-!> velocities to next time-step
- subroutine particlesVelocityNext
- use var
- use parallel
- real(kind=pr),dimension(maxnp) :: dup,dvp,dwp
- real(kind=pr) :: f_g(3)
-
- f_g= (1-rhof/rhop)*g
-
- dup(1:np)= (fx_d(1:np) + fx_l(1:np) + fx_el(1:np) + f_g(1)) * dtNext
- dvp(1:np)= (fy_d(1:np) + fy_l(1:np) + fy_el(1:np) + f_g(2)) * dtNext
- dwp(1:np)= (fz_d(1:np) + fz_l(1:np) + fz_el(1:np) + f_g(3)) * dtNext
-
- upNext(1:np)= up(1:np) + dup(1:np)
- vpNext(1:np)= vp(1:np) + dvp(1:np)
- wpNext(1:np)= wp(1:np) + dwp(1:np)
-
- dup_max=max(maxval(abs(dup(1:np))),maxval(abs(dvp(1:np))),maxval(abs(dwp(1:np))))
- dup_max=syncMax(dup_max)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief specific particle drag force
- subroutine particlesDrag
- use var
- use parallel
- real(kind=pr) :: tau_p
- integer :: n
-
- call fluidVelocity
-
- tau_p_min=1.e10_pr
-
- do 1 n=1,np
- call calcTau_p(tau_p,n)
- fx_d(n)= (uf(n)-up(n))/tau_p
- fy_d(n)= (vf(n)-vp(n))/tau_p
- fz_d(n)= (wf(n)-wp(n))/tau_p
-1 enddo
-
- tau_p_min=syncMin(tau_p_min)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief specific particle lift force (Saffman, 1965, 1968, correction by Mei, 1992)
- subroutine particlesLift
- use var
- real(kind=pr) :: urel,Reyp,Reyf,beta,Cls
- integer :: n
-
- call fluidVelocityGradient
-
- do 1 n=1,np
- urel= sqrt((uf(n)-up(n))**2+(vf(n)-vp(n))**2+(wf(n)-wp(n))**2)
- Reyp= max(1.e-10_pr,2._pr*radp(n)*urel/nuf)
- Reyf= radp(n)**2/nuf*sqrt(dufdy(n)**2+dufdz(n)**2)
- beta= Reyf/Reyp/2._pr
-
- if (Reyp.le.40._pr) then
- Cls= (1._pr-0.3314_pr*sqrt(beta))*exp(-Reyp/10._pr)+0.3314_pr*sqrt(beta)
- else
- Cls= 0.0524_pr*sqrt(beta*Reyp)
- endif
-
- fx_l(n)= 0._pr
- fy_l(n)= 1.54_pr*sqrt(nuf)*rhof/rhop/radp(n)*(uf(n)-up(n))*sign(sqrt(abs(dufdy(n))),dufdy(n))*Cls
- fz_l(n)= 1.54_pr*sqrt(nuf)*rhof/rhop/radp(n)*(uf(n)-up(n))*sign(sqrt(abs(dufdz(n))),dufdz(n))*Cls
-1 enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief particle response time
- subroutine calcTau_p(tau_p,n)
- use var
- real(kind=pr) :: urel,Reyp,Cd,tau_p
- integer :: n
-
- urel= sqrt((uf(n)-up(n))**2+(vf(n)-vp(n))**2+(wf(n)-wp(n))**2)
- urel= max(1.e-10_pr,urel) ! avoid NaN
- Reyp= 2._pr*radp(n)*urel/nuf
-
- if (Reyp.lt.1000._pr) then ! Putnam (1961)
- Cd= 24._pr/Reyp * (1._pr + 1._pr/6._pr*Reyp**(2._pr/3._pr))
- else
- Cd= 0.424_pr
- endif
-
- tau_p= 8._pr*rhop*radp(n)/(3._pr*rhof*urel*Cd)
-
- tau_p_min=min(tau_p_min,tau_p)
- if (tau_p.lt.dtNext) tau_p= dtNext ! stability criteria Euler forward
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief calculate the fluid velocity around each particle
- subroutine fluidVelocity
- use var
- real(kind=pr),dimension(2,2,2) :: weight
- real(kind=pr) volE
- integer :: n,ibeg,iend,jbeg,jend,lbeg,lend
- integer :: direction
-
- do 1 n=1,np
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,1)
- uf(n)=sum(weight*u(ibeg:iend,jbeg:jend,lbeg:lend))
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,2)
- vf(n)=sum(weight*v(ibeg:iend,jbeg:jend,lbeg:lend))
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,3)
- wf(n)=sum(weight*w(ibeg:iend,jbeg:jend,lbeg:lend))
-1 enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief calculate the fluid velocity gradient around each particle
- subroutine fluidVelocityGradient
- use var
- real(kind=pr),dimension(2,2,2) :: weight
- real(kind=pr) volE
- integer :: n,ibeg,iend,jbeg,jend,lbeg,lend,i,j,l,iw,jw,lw
-
- dufdy=0._pr; dufdz=0._pr
-
- do 1 n=1,np
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,1) ! could be extended for 2 and 3
- do l=lbeg,lend; do j=jbeg,jend; do i=ibeg,iend
- iw=i+1-ibeg; jw=j+1-jbeg; lw=l+1-lbeg
- dufdy(n)=dufdy(n)+weight(iw,jw,lw)*(u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
- dufdz(n)=dufdz(n)+weight(iw,jw,lw)*(u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
- enddo; enddo; enddo
-1 enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief second-order Crank-Nicolson integration of particle positions to next time-step
- subroutine particlesTransportNext
- use var
- use parallel
- integer :: n,m
- character(70) :: filename
- logical :: remove(maxnp)
-
- remove=.false.
-
- npTot=syncSumI(np)
-
- xpNext(1:np)= xp(1:np)+(up(1:np)*dt+upNext(1:np)*dtNext)*dtNext/(dt+dtNext)
- ypNext(1:np)= yp(1:np)+(vp(1:np)*dt+vpNext(1:np)*dtNext)*dtNext/(dt+dtNext)
- zpNext(1:np)= zp(1:np)+(wp(1:np)*dt+wpNext(1:np)*dtNext)*dtNext/(dt+dtNext)
-
- do 1 n=1,np ! boundary conditions
-
- if (bcx.eq.'p') then ! x-direction
- ! done in particlesNextProc
- elseif (bcx.eq.'w') then
- if (((myid.eq.0.).and.(xpNext(n).lt.xmin+radp(n))).or. &
- ((myid.eq.nrprocs-1).and.(xpNext(n).gt.xmax-radp(n)))) then
- upNext(n)=-upNext(n)*restRatio
- if (xpNext(n).gt.xmax-radp(n)) then
- xpNext(n)=(xmax-radp(n)) - (xpNext(n)-(xmax-radp(n)))
- elseif (xpNext(n).lt.xmin+radp(n)) then
- xpNext(n)=(xmin+radp(n)) - (xpNext(n)-(xmin+radp(n)))
- endif
- wcollnum(n)=wcollnum(n)+1
- if (qpmax.ne.qp0) call chargeParticleWall(n,1)
- endif
- elseif (bcx.eq.'i') then
- if (myid.eq.nrprocs-1) then
- if (xp(n).gt.xmax) then
- remove(n)= .true.
- else
- remove(n)= .false.
- endif
- endif
- endif
-
- if (bcy.eq.'p') then ! y-direction
- if (ypNext(n).lt.ymin) ypNext(n)=ypNext(n)+dimy
- if (ypNext(n).gt.ymax) ypNext(n)=ypNext(n)-dimy
- elseif (bcy.eq.'w') then
- if ((ypNext(n).lt.ymin+radp(n)).or.(ypNext(n).gt.ymax-radp(n))) then
- vpNext(n)=-vpNext(n)*restRatio
- if (ypNext(n).gt.ymax-radp(n)) then
- ypNext(n)=(ymax-radp(n)) - (ypNext(n)-(ymax-radp(n)))
- elseif (ypNext(n).lt.ymin+radp(n)) then
- ypNext(n)=(ymin+radp(n)) - (ypNext(n)-(ymin+radp(n)))
- endif
- wcollnum(n)=wcollnum(n)+1
- if (qpmax.ne.qp0) call chargeParticleWall(n,2)
- endif
- endif
-
- if (bcz.eq.'p') then ! z-direction
- if (zpNext(n).lt.zmin) zpNext(n)=zpNext(n)+dimz
- if (zpNext(n).gt.zmax) zpNext(n)=zpNext(n)-dimz
- elseif (bcz.eq.'w') then
- if ((zpNext(n).lt.zmin+radp(n)).or.(zpNext(n).gt.zmax-radp(n))) then
- wpNext(n)=-wpNext(n)*restRatio
- if (zpNext(n).gt.zmax-radp(n)) then
- zpNext(n)=(zmax-radp(n)) - (zpNext(n)-(zmax-radp(n)))
- elseif (zpNext(n).lt.zmin+radp(n)) then
- zpNext(n)=(zmin+radp(n)) - (zpNext(n)-(zmin+radp(n)))
- endif
- wcollnum(n)=wcollnum(n)+1
- if (qpmax.ne.qp0) call chargeParticleWall(n,3)
- endif
- endif
-
-1 enddo
-
-
- if ((bcx.eq.'i').and.(myid.eq.nrprocs-1)) then ! remove particles from outlet
- m=0
- do 4 n=1,np
- if (remove(n)) then
- write(filename,'(a,i3.3,a)') 'results/particlesOut_p',myid,'.dat'
- open(access='append',unit=10,file=filename)
- write(10,77) nt,t,q_el(n)
-77 format(i7,3(x,es11.2e2))
- close(10)
- else
- m=m+1
- call partN2M(n,m)
- endif
-4 enddo
- np=m
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief compute momentum source term for fluid phase (habil eq 2.12)
- subroutine momentumCoupling
- use var
- real(kind=pr) :: volE,partmass
- real(kind=pr),dimension(2,2,2) :: weight
- integer :: n,ibeg,iend,jbeg,jend,lbeg,lend
-
- Fsx=0._pr; Fsy=0._pr; Fsz=0._pr
-
- call particlesDrag ! only aerodynamic forces
- call particlesLift
-
- do 1 n=1,np
- partmass= rhop*4._pr/3._pr*pi*partn(n)*radp(n)**3
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,1)
- Fsx(ibeg:iend,jbeg:jend,lbeg:lend)= Fsx(ibeg:iend,jbeg:jend,lbeg:lend) &
- - partmass/rhof/volE*(fx_d(n)+fx_l(n))*weight
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,2)
- Fsy(ibeg:iend,jbeg:jend,lbeg:lend)= Fsy(ibeg:iend,jbeg:jend,lbeg:lend) &
- - partmass/rhof/volE*(fy_d(n)+fy_l(n))*weight
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,3)
- Fsz(ibeg:iend,jbeg:jend,lbeg:lend)= Fsz(ibeg:iend,jbeg:jend,lbeg:lend) &
- - partmass/rhof/volE*(fz_d(n)+fz_l(n))*weight
-1 enddo
-
- call bcLEsource(Fsx); call bcLEsource(Fsy); call bcLEsource(Fsz)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief particles collisions in next time-step using ray casting
- subroutine particlesCollideNext
- use var
- real(kind=pr),dimension(maxnp) :: xpc,ypc,zpc
- real(kind=pr),dimension(3) :: relvelo,midlink,ncontact,vn,ve
- real(kind=pr) :: sumrad,upp1,vpp1,wpp1,upp2,vpp2,wpp2, &
- partdist,r3n1,r3n2,projvtox,srel,projvton, &
- closestrelvelo,fracdt
- integer :: n1,n2,m1,m2,i,j,l,syncSumI
-
- do l=1,ll; do j=1,jj; do i=1,ii
-! condition I only particles in same cell:
- do 1 m1=2,npic(i,j,l)
- do 2 m2=1,(m1-1)
-
- n1=nic(i,j,l,m1)
- n2=nic(i,j,l,m2)
-
-! setting the 2 particles in the rest frame of particle n2
- relvelo= (/ up(n1)-up(n2) , vp(n1)-vp(n2) , wp(n1)-wp(n2) /)
- midlink= (/ xp(n2)-xp(n1) , yp(n2)-yp(n1) , zp(n2)-zp(n1) /)
- projvtox=dot_product(midlink,(relvelo*dtNext))
-
-! condition IIa to check collision (propagation direction):
- if (projvtox.lt.0._pr) cycle
-
- srel= sqrt(dot_product((relvelo*dtNext),(relvelo*dtNext)))
-
-! condition IIb to check collision (no rel velocity, no collision):
- if (abs(srel).lt.1.e-19_pr) cycle
-
- partdist= sqrt(dot_product(midlink,midlink))
- sumrad= radp(n1)+radp(n2)
-
-! condition III to check collision (contact expected):
- if ((partdist**2-(projvtox/srel)**2).gt.sumrad**2) then
- cycle
- elseif (partdist.le.sumrad) then ! if overlapping after seeding
- cycle
- endif
-
- closestrelvelo= sqrt(sumrad**2-(partdist**2-(projvtox/srel)**2))
- fracdt= (projvtox/srel-closestrelvelo)/srel
-
-! condition IV to check collision (contact during next timestep):
- if ((fracdt.gt.1._pr).or.(fracdt.lt.0._pr)) cycle
-
-! all conditions fullfilled -> collision happens
- ppcollnum(n1)= ppcollnum(n1)+1
- ppcollnum(n2)= ppcollnum(n2)+1
- if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2)
-
-! fictitious contact point
- xpc(n1)= xp(n1)+fracdt*up(n1)*dtNext
- ypc(n1)= yp(n1)+fracdt*vp(n1)*dtNext
- zpc(n1)= zp(n1)+fracdt*wp(n1)*dtNext
- xpc(n2)= xp(n2)+fracdt*up(n2)*dtNext
- ypc(n2)= yp(n2)+fracdt*vp(n2)*dtNext
- zpc(n2)= zp(n2)+fracdt*wp(n2)*dtNext
-
- r3n1= radp(n1)**3
- r3n2= radp(n2)**3
-
-! normal vector to the contact point ('Stossnormale')
- ncontact= (/ xpc(n1)-xpc(n2) , ypc(n1)-ypc(n2) , zpc(n1)-zpc(n2) /)
-
-! project the contact velocity on the contact normal
-! these components are zero for n2 bc it is in rest
- projvton= dot_product(ncontact,relvelo)/(sumrad**2)
- vn= projvton*ncontact
-
-! velocity components in the contact plane which do not change during contact
-! these components are zero for n2 bc it is in rest
- ve= relvelo-vn
-
-! new v = rest frame + unchanged components in contact plane + central collision
- upp1= up(n2) + ve(1) + (r3n1-restRatio*r3n2)*vn(1)/(r3n1+r3n2)
- vpp1= vp(n2) + ve(2) + (r3n1-restRatio*r3n2)*vn(2)/(r3n1+r3n2)
- wpp1= wp(n2) + ve(3) + (r3n1-restRatio*r3n2)*vn(3)/(r3n1+r3n2)
-
- upp2= up(n2) + (1._pr+restRatio)*r3n1*vn(1)/(r3n1+r3n2)
- vpp2= vp(n2) + (1._pr+restRatio)*r3n1*vn(2)/(r3n1+r3n2)
- wpp2= wp(n2) + (1._pr+restRatio)*r3n1*vn(3)/(r3n1+r3n2)
-
- if (max(vpp1,wpp1,vpp2,wpp2,.5_pr*upp1,.5_pr*upp2).gt.2._pr*ubulk) then
- write(*,*) 'Warning: post-collision velocity too high'
- write(*,*) up(n1),vp(n1),wp(n1),upp1,vpp1,wpp1,up(n2),vp(n2),wp(n2),upp2,vpp2,wpp2
- endif
-
- if(partn(n1).gt.partn(n2)) then
- upNext(n2)= upp2
- vpNext(n2)= vpp2
- wpNext(n2)= wpp2
- upNext(n1)= (partn(n2)*upp1+(partn(n1)-partn(n2))*up(n1))/partn(n1)
- vpNext(n1)= (partn(n2)*vpp1+(partn(n1)-partn(n2))*vp(n1))/partn(n1)
- wpNext(n1)= (partn(n2)*wpp1+(partn(n1)-partn(n2))*wp(n1))/partn(n1)
- else
- upNext(n1)= upp1
- vpNext(n1)= vpp1
- wpNext(n1)= wpp1
- upNext(n2)= (partn(n1)*upp2+(partn(n2)-partn(n1))*up(n2))/partn(n2)
- vpNext(n2)= (partn(n1)*vpp2+(partn(n2)-partn(n1))*vp(n2))/partn(n2)
- wpNext(n2)= (partn(n1)*wpp2+(partn(n2)-partn(n1))*wp(n2))/partn(n2)
- endif
-
-2 enddo
-1 enddo
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief seed particles
- subroutine particlesSeed
- use var
- use parallel
- real(kind=pr),allocatable,dimension(:) :: randomx,randomy,randomz
- integer :: n,npseed,npseedTot
- character(70) :: filename
-
-!> number of particles to be seeded
- npp=np
- if (bcx.eq.'i') then
- if (myid.eq.0) then
- npseed=int((t-tstart)*pnd)-npstart
- npstart=npstart+npseed
- else
- npseed=0
- endif
- else
- npseed=int(pnd*dimx*dimy*dimz)
-! npseed=2
- endif
- npseedTot=syncSumI(npseed)
- np=npp+npseed
- call allocateParticleArrays
-
- if (npseed.eq.0) return
-
- write(filename,'(a,i3.3,a)') 'results/particlesSeed_p',myid,'.dat'
- open(access='append',unit=10,file=filename)
-
- allocate(randomx(npseed),randomy(npseed),randomz(npseed))
- if (bcx.eq.'i') then
- if (myid.eq.0) then
- call random_number(randomx)
- call random_number(randomy)
- call random_number(randomz)
- endif
- else
- call syncRandom(npseed,randomx)
- call syncRandom(npseed,randomy)
- call syncRandom(npseed,randomz)
- endif
-
- do 1 n=npp+1,np
-!> seed not directly on wall or an active gc:
- if (bcx.eq.'i') then
- xp(n)=randomx(n-npp)*ubulk*dt+xf(imin-1)
- else
- xp(n)=randomx(n-npp)*(dimx-2._pr*prad)+xf(imin-1)+prad
- endif
- yp(n)=randomy(n-npp)*(dimy-2._pr*prad)+yf(jmin-1)+prad
- zp(n)=randomz(n-npp)*(dimz-2._pr*prad)+zf(lmin-1)+prad
-
- nseedLoc=nseedLoc+1
- nGlob(n)=nseedLoc*1000+myid
-
- wcollnum(n)=0
- ppcollnum(n)=0
-
- radp(n)=prad
- partn(n)=1._pr
- q_el(n)=qp0
- q_elNext(n)=q_el(n)
-
-! xp(1)=0.02;yp(1)=0.;zp(1)=0.
-! xp(2)=0.03;yp(2)=0.;zp(2)=0.
-! up(1)=1.;up(2)=-1;vp(1)=0.;vp(2)=0.;wp(1)=0.;wp(2)=0.
-
- write(10,'(3(es10.2e2,x))') radp(n),q_el(n),partn(n)
-
-1 enddo
-
- call particlesToCells
- call fluidVelocity
-
- up(npp+1:np)= uf(npp+1:np)
- vp(npp+1:np)= vf(npp+1:np)
- wp(npp+1:np)= wf(npp+1:np)
-
- close(10)
-
- end
diff --git a/src/post.f90 b/src/post.f90
deleted file mode 100755
index cacd8cd..0000000
--- a/src/post.f90
+++ /dev/null
@@ -1,122 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
- subroutine postProcessing
- use var
-
- if ((mod(nt,int_results).eq.0).or.(nt.eq.1).or.(nt.eq.ntend)) then
- call monitor
- call writevtk_fluid_xy
- call writevtk_fluid_xz
- call writevtk_fluid_yz
- call writevtk_fluid_xyz
- call writevtk_particles
- call writedat_particles
- endif
-
- if ((mod(nt,int_restart).eq.0).or.(nt.eq.ntend)) call saveField
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine monitor
- use var
- use parallel
- real(kind=pr) :: ucl,avu,avv,avw,gradu,graduy,graduz, &
- Rec,Reb,avvp,avwp,avyp,avzp,avvf,avwf,avxp,avqp, &
- rmsv,rmsw,C05,A05
- integer :: m,N05
- logical :: file_ex
- character(70) :: filename2
-
- if (bcy.eq.'w'.and.bcz.eq.'w') then
- ucl= sum(u(imin:imax,int(jj/2),int(ll/2)))/(imax-imin+1)
- elseif (bcy.eq.'w'.and.bcz.eq.'p') then
- ucl= sum(u(imin:imax,int(jj/2),lmin:lmax))/(imax-imin+1)/(lmax-lmin+1)
- elseif (bcy.eq.'p'.and.bcz.eq.'w') then
- ucl= sum(u(imin:imax,jmin:jmax,int(ll/2)))/(imax-imin+1)/(jmax-jmin+1)
- endif
-
- avu= sum(u(imin:imax,jmin:jmax,lmin:lmax)*volu(imin:imax,jmin:jmax,lmin:lmax)) / volTot
- avv= sum(v(imin:imax,jmin:jmax,lmin:lmax)*volv(imin:imax,jmin:jmax,lmin:lmax)) / volTot
- avw= sum(w(imin:imax,jmin:jmax,lmin:lmax)*volw(imin:imax,jmin:jmax,lmin:lmax)) / volTot
-
- rmsv=sum((v(imin:imax,jmin:jmax,lmin:lmax))**2 *volv(imin:imax,jmin:jmax,lmin:lmax)) / volTot
- rmsv=sqrt(rmsv)
- rmsw=sum((w(imin:imax,jmin:jmax,lmin:lmax))**2 *volw(imin:imax,jmin:jmax,lmin:lmax)) / volTot
- rmsw=sqrt(rmsw)
-
- graduy= ( sum(u(imin:imax,jmin,lmin:lmax)*volu(imin:imax,jmin,lmin:lmax)) &
- +sum(u(imin:imax,jmax,lmin:lmax)*volu(imin:imax,jmax,lmin:lmax))) &
- / (sum(volu(imin:imax,jmin,lmin:lmax))+sum(volu(imin:imax,jmax,lmin:lmax))) &
- / (ymax-yc(jmax))
- graduz= ( sum(u(imin:imax,jmin:jmax,lmin)*volu(imin:imax,jmin:jmax,lmin)) &
- +sum(u(imin:imax,jmin:jmax,lmax)*volu(imin:imax,jmin:jmax,lmax))) &
- / (sum(volu(imin:imax,jmin:jmax,lmin))+sum(volu(imin:imax,jmin:jmax,lmax))) &
- / (zmax-zc(lmax))
-
- if (bcy.eq.'w'.and.bcz.eq.'w') gradu=(graduy+graduz)/2._pr
- if (bcy.eq.'w'.and.bcz.eq.'p') gradu=graduy
- if (bcy.eq.'p'.and.bcz.eq.'w') gradu=graduz
-
- if (npTot.gt.0) then
- avyp=syncSum(sum(yp(1:np)))/npTot
- avzp=syncSum(sum(zp(1:np)))/npTot
- avqp=syncSum(sum(q_el(1:np)))/npTot
- if (bcy.eq.'w'.and.bcz.eq.'w') then
- N05=count((min((yp(1:np)-ymin),(ymax-yp(1:np))).le.(dimy*0.05_pr)).or. &
- (min((zp(1:np)-zmin),(zmax-zp(1:np))).le.(dimz*0.05_pr)))
- A05=(dimy*dimz)-(dimy*0.90_pr)*(dimz*0.90_pr)
- elseif (bcy.eq.'w'.and.bcz.eq.'p') then
- N05=count(min((yp(1:np)-ymin),(ymax-yp(1:np))).le.(dimy*0.05_pr))
- A05=(dimy*dimz)-(dimy*0.90_pr)*dimz
- elseif (bcy.eq.'p'.and.bcz.eq.'w') then
- N05=count(min((zp(1:np)-zmin),(zmax-zp(1:np))).le.(dimz*0.05_pr))
- A05=(dimy*dimz)-dimy*(dimz*0.90_pr)
- endif
- C05=real(syncSumI(N05),kind=pr)/(A05*dimxTot)/pnd
- else
- npTot=0
- avyp= 0._pr
- avzp= 0._pr
- avqp= 0._pr
- C05= 0._pr
- endif
-
- ucl= syncSum(ucl)/nrprocs
- ubulk=avu
- avu= syncSum(avu)/nrprocs
- avv= syncSum(avv)/nrprocs
- avw= syncSum(avw)/nrprocs
- rmsv= syncSum(rmsv)/nrprocs
- rmsw= syncSum(rmsw)/nrprocs
- gradu=syncSum(gradu)/nrprocs
-
- if (myid.ne.0) return
-
- Rec= ucl*delta/nuf
- Reb= avu*delta/nuf
- tau_w= muf*abs(gradu)
- u_tau= sqrt(tau_w/rhof)
- Re_tau= u_tau*delta/nuf
- if (Re_tau.gt.0._pr) delta_v= delta/Re_tau
-
- filename2='output/monitor.dat'
- inquire(file=filename2,exist=file_ex)
- if (file_ex.and.nt.ne.1) then
- open(10,file=filename2,access='append')
- else
- open(10,file=filename2)
- write(10,'(a,a)') '# ',version
- write(10,'(a,i10,i14,15(i10))') '#',(m,m=1,16)
- write(10,'(a,a10,a14,14(a10))') &
- '#','nt','t','dt','ucl','av(u)','av(v)','av(w)','rms(v)','rms(w)', &
- 'Rec','Reb','tau_w','u_{tau}','Re_{tau}','C_{05}','av(qp)'
- endif
-
- write(10,'(x,i10,es14.6e2,20(es10.2e2))') &
- nt,t,dt,ucl,avu,avv,avw,rmsv,rmsw,Rec,Reb,tau_w,u_tau,Re_tau,C05,avqp
-
- close(10)
-
- end
diff --git a/src/pre.f90 b/src/pre.f90
deleted file mode 100755
index 44c5afe..0000000
--- a/src/pre.f90
+++ /dev/null
@@ -1,400 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief pre processing
- subroutine preProcessing
- use var
-
- nt=0
- call flowConditions
- call initVariables
- call gridGeneration
- call computeCoefficients
- call writevtk_grid
- if (ntstart.eq.1) then
- call initFlow
- elseif (ntstart.gt.1) then
- call readField
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief read input data, compute flow conditions and write out
- subroutine flowConditions
- use var
- use mpi
- integer :: dimitotTemp
- character :: filename*70,comp_date*8,comp_time*10
-
-! read input file
- open(unit=10,file='input/input.dat',status='old')
-
- read(10,*)
- read(10,*)
- read(10,*)
- read(10,*) dimxtot,dimy,dimz
- read(10,*)
- read(10,*) bcx,bcy,bcz
- read(10,*)
- read(10,*) dimitotTemp,dimj,diml
- read(10,*)
- read(10,*) gridx,gridy,gridz
- read(10,*)
- read(10,*) ntstart,ntime
- read(10,*)
- read(10,*) int_results,int_restart
- read(10,*)
- read(10,*) turbModel
- read(10,*)
- read(10,*) cfl,itmax,tol
- read(10,*)
- read(10,*) ubulk0,dpdx
- read(10,*)
- read(10,*) rhof,nuf
- read(10,*)
- read(10,*) pnd,ntseed
- read(10,*)
- read(10,*) prad,rhop,restRatio
- read(10,*)
- read(10,*) qp0,qpmax
- read(10,*)
- read(10,*) g(1),g(2),g(3)
-
- close(10)
-
-! connect processors
- next=myid+1
- prev=myid-1
-
- select case (bcx)
- case ('p')
- if(prev.eq.-1) prev=nrprocs-1
- if(next.eq.nrprocs) next=0
- case ('w','i')
- if(prev.eq.-1) prev=mpi_proc_null
- if(next.eq.nrprocs) next=mpi_proc_null
- end select
-
-! compute flow conditions
- ntend=ntstart+ntime-1
-
- dimx=dimxtot/nrprocs
- dimi=int(dimitotTemp/nrprocs)
- dimitot=dimi*nrprocs
- if (myid.eq.0.and.dimitot.ne.dimitotTemp) write(*,'(a,i8)') 'dimi set to ',dimitot
-
- npTot=int(pnd*dimxtot*dimy*dimz)
- muf=nuf*rhof
- ubulk=ubulk0
-
- if (bcy.eq.'w'.and.bcz.eq.'w') delta= min(dimy/2._pr,dimz/2._pr)
- if (bcy.eq.'w'.and.bcz.eq.'p') delta= dimy/2._pr
- if (bcy.eq.'p'.and.bcz.eq.'w') delta= dimz/2._pr
-
- Re = ubulk*delta/nuf
- ftt = dimxtot/ubulk
-
-! write output file
- if (myid.ne.0) return
- call date_and_time(date=comp_date,time=comp_time)
-
- write(filename,'(3a,a6,a)') 'output/output_',comp_date,'_',comp_time,'.dat'
- open(unit=20,file=filename,status='replace')
-
- write(20,'(a,x,a8,a,x,a6)') 'Simulation start date:',comp_date,', time',comp_time
- write(20,*)
- write(20,'(a)') 'Conditions:'
- write(20,'(a,i4)') &
- 'Nr. of processors = ',nrprocs
- write(20,'(a,3(es11.2e2))') &
- 'dimensions in x,y,z-direction (m,m,m) = ',dimxtot,dimy,dimz
- write(20,'(a,3x,3(a,x))') &
- 'BCs in x,y,z-direction [(w)all/(p)eriodic] = ',bcx,bcy,bcz
- write(20,'(a,3(i5),a,x,i8)') &
- 'Nr. of cells in x,y,z-direction (-,-,-) = ',dimitot,dimj,diml,', total =',dimitot*dimj*diml
- write(20,'(a,3(i5),a,x,i8)') &
- 'Nr. of cells in x,y,z-direction per proc. (-,-,-) = ',dimi,dimj,diml,', total =',dimi*dimj*diml
- write(20,'(a,3x,3(a,x))') &
- 'Grid in x,y,z-direction [(u)niform/(s)tretch] = ',gridx,gridy,gridz
- write(20,'(a,4(i8))') &
- 'nt: start, compute, results, restart (-,-,-,-) = ',ntstart,ntime,int_results,int_restart
- write(20,'(a,3x,a)') &
- 'turbulence model: = ',turbModel
- write(20,'(a,es11.2e2,i11,es11.2e2)') &
- 'CFL number (-), max number iterations (-), rel. tol. = ',cfl,itmax,tol
- write(20,'(a,2(es11.2e2))') &
- 'in/initial flow (m/s), pressure gradient (N/m**3) = ',ubulk,dpdx
- write(20,'(a,2(es11.2e2))') &
- 'fluid density (kg/m**3), kinematic visc. (m**2/s) = ',rhof,nuf
- write(20,'(a,i11,es11.2,i11)') &
- 'particle number (-), (-/m**3) or (-/s), seed nt (-) = ',npTot,pnd,ntseed
- write(20,'(a,3(es11.2e2))') &
- 'particle rad (m),mat density (kg/m**3),rest coeff (-) = ',prad,rhop,restRatio
- write(20,'(a,3(es11.2e2))') &
- 'particle charge initial, maximum (C,C) = ',qp0,qpmax
- write(20,'(a,3(es11.2e2))') &
- 'gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) = ',g(1),g(2),g(3)
- write(20,'(a,8(es11.2e2))') &
- 'delta (m), Re (-), flow through time (s) = ',delta,Re,ftt
-
- close(20)
-
- end
-
-!####################################################################
-!> @brief allocate and initialize all variables
- subroutine initVariables
- use var
- use mpi
-
- ii=dimi+2*gc; jj=dimj+2*gc; ll=diml+2*gc
- maxnp=0
-
-! parallel
- allocate(mpistatus(mpi_status_size))
-
-! grid
- allocate(celltype(ii,jj,ll))
- allocate(xc(ii),yc(jj),zc(ll),xf(ii),yf(jj),zf(ll), &
- vol(ii,jj,ll),volu(ii,jj,ll),volv(ii,jj,ll),volw(ii,jj,ll), &
- Ca(ii,jj,ll),Cb1(ii,jj,ll),Cb2(ii,jj,ll),Cc1(ii,jj,ll),Cc2(ii,jj,ll), &
- Cd1(ii,jj,ll),Cd2(ii,jj,ll))
-
-! fluid
- allocate(u(ii,jj,ll),v(ii,jj,ll),w(ii,jj,ll),p(ii,jj,ll), &
- rho_el(ii,jj,ll),phi_el(ii,jj,ll), &
- Ex_el(ii,jj,ll),Ey_el(ii,jj,ll),Ez_el(ii,jj,ll), &
- u01(ii,jj,ll),v01(ii,jj,ll),w01(ii,jj,ll), &
- u02(ii,jj,ll),v02(ii,jj,ll),w02(ii,jj,ll), &
- Fsx(ii,jj,ll),Fsy(ii,jj,ll),Fsz(ii,jj,ll), &
- nuf_ru(ii,jj,ll),nuf_rv(ii,jj,ll),nuf_rw(ii,jj,ll))
-
-! particles
- allocate(npic(ii,jj,ll))
- np=0; npp=0; nseedLoc=0; npstart=0
- call allocateParticleArrays
-
-! init
- u=0._pr; v=0._pr; w=0._pr; p=0._pr
- u01=u; v01=v; w01= w
- u02=u; v02=v; w02= w
- Fsx=0._pr; Fsy=0._pr; Fsz=0._pr
- nuf_ru=0._pr; nuf_rv=0._pr; nuf_rw=0._pr
- celltype=active
- t=0._pr; timecom=0._pr
- tau_p_min=0._pr; dup_max=0._pr
- vol=0._pr; volu=0._pr; volv=0._pr; volw=0._pr;
- Ca=0._pr; Cb1=0._pr; Cb2=0._pr; Cc1=0._pr; Cc2=0._pr; Cd1=0._pr; Cd2=0._pr;
-
-! electrostatics
- rho_el=0._pr
- phi_el=0._pr
- Ex_el= 0._pr; Ey_el=0._pr; Ez_el=0._pr
- fx_el= 0._pr; fy_el=0._pr; fz_el=0._pr
- fx_d= 0._pr; fy_d= 0._pr; fz_d= 0._pr
- fx_l= 0._pr; fy_l= 0._pr; fz_l= 0._pr
- q_el= 0._pr; q_elNext=0._pr;
-
- call init_random_seed()
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief initialize random number generator
- subroutine init_random_seed()
- use var
- integer :: i,n,clock
- integer, dimension(:), allocatable :: seed
-
- call random_seed(size=n)
- allocate(seed(n))
- call system_clock(count=clock)
- seed = clock + 37 * (/ (i - 1, i = 1, n) /)
- call random_seed(put=seed)
- deallocate(seed)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief generate grid for gaseous phase
- subroutine gridGeneration
- use var
- real(kind=pr) :: hy_temp,hz_temp
- integer :: i,j,l
-
- imin=1+gc; imax=ii-gc
- jmin=1+gc; jmax=jj-gc
- lmin=1+gc; lmax=ll-gc
- gp=ii*jj*ll; dimgp=dimi*dimj*diml; dimgptot=(dimi*nrprocs)*dimj*diml
-
-!> celltypes, default=active
- do 1 i=1,ii
- select case (bcx)
- case ('p')
- if ((i.lt.imin).or.(i.gt.imax)) celltype(i,:,:)= passive
- case ('w','i')
- if ((i.lt.imin).and.(myid.eq.0)) then
- celltype(i,:,:)= wall
- elseif ((i.lt.imin).and.(myid.ne.0)) then
- celltype(i,:,:)= passive
- elseif ((i.gt.imax).and.(myid.eq.nrprocs-1)) then
- celltype(i,:,:)= wall
- elseif ((i.gt.imax).and.(myid.ne.nrprocs-1)) then
- celltype(i,:,:)= passive
- endif
- end select
-1 enddo
- do 2 j=1,jj
- select case (bcy)
- case ('p')
- if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= passive
- case ('w','i')
- if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= wall
- end select
-2 enddo
- do 3 l=1,ll
- select case (bcz)
- case ('p')
- if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= passive
- case ('w','i')
- if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= wall
- end select
-3 enddo
-
-!> grid
- do 4 i=1,ii
- if (gridx.eq.'u') then
- xf(i)=real(i-gc,kind=pr)*dimx/real(dimi,kind=pr)
- elseif (gridx.eq.'s') then
- !not implemented yet
- endif
-
- if (i.eq.1) then
- xc(i)= xf(1)-(dimx/real(dimi,kind=pr))/2._pr
- else
- xc(i)= (xf(i)+xf(i-1))/2._pr
- endif
-4 enddo
-
- do 5 j=1,jj
- if (gridy.eq.'u') then
- hy_temp=dimy/real(dimj,kind=pr)
- yf(j)=(real(j,kind=pr)-real(jj,kind=pr)/2._pr)*hy_temp
- elseif (gridy.eq.'s') then
- hy_temp= (1._pr-cos(pi/real(dimj,kind=pr)))*dimy/2._pr !width of first active cell
- if (j.lt.jmin) then
- yf(j)= real(j+1-jmin,kind=pr)*hy_temp-dimy/2._pr
- elseif (j.gt.jmax) then
- yf(j)= real(j-jmax,kind=pr)*hy_temp+dimy/2._pr
- else
- yf(j)=-(cos(real(j-gc,kind=pr)*pi/real(dimj,kind=pr)))*dimy/2._pr
- endif
- endif
-
- if (j.eq.1) then
- yc(j)= yf(1)-hy_temp/2._pr
- else
- yc(j)= (yf(j)+yf(j-1))/2._pr
- endif
-5 enddo
-
- do 6 l=1,ll
- if (gridz.eq.'u') then
- hz_temp=dimz/real(diml,kind=pr)
- zf(l)=(real(l,kind=pr)-real(ll,kind=pr)/2._pr)*hz_temp
- elseif (gridz.eq.'s') then
- hz_temp= (1._pr-cos(pi/real(diml,kind=pr)))*dimz/2._pr !width of first active cell
- if (l.lt.lmin) then
- zf(l)= real(l+1-lmin,kind=pr)*hz_temp-dimz/2._pr
- elseif (l.gt.lmax) then
- zf(l)= real(l-lmax,kind=pr)*hz_temp+dimz/2._pr
- else
- zf(l)=-(cos(real(l-gc,kind=pr)*pi/real(diml,kind=pr)))*dimz/2._pr
- endif
- endif
-
- if (l.eq.1) then
- zc(l)= zf(1)-hz_temp/2._pr
- else
- zc(l)= (zf(l)+zf(l-1))/2._pr
- endif
-6 enddo
-
-
- deltax=xf(ii-gc)*real(myid,kind=pr)
- xc=xc+deltax
- xf=xf+deltax
-
- xmin=xf(gc); xmax=xf(ii-gc)
- ymin=yf(gc); ymax=yf(jj-gc)
- zmin=zf(gc); zmax=zf(ll-gc)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief pre-compute coefficients to speed of the solution loops
- subroutine computeCoefficients
- use var
- integer :: i,j,l
-
- do i=imin-1,imax+1; do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
- vol(i,j,l)= (xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
- volu(i,j,l)=(xc(i+1)-xc(i))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
- volv(i,j,l)=(xf(i)-xf(i-1))*(yc(j+1)-yc(j))*(zf(l)-zf(l-1))
- volw(i,j,l)=(xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zc(l+1)-zc(l))
- enddo; enddo; enddo
- volTot= sum(vol(imin:imax,jmin:jmax,lmin:lmax))
-
-! LES filter width
- filteru= volu**(1._pr/3._pr)
- filterv= volv**(1._pr/3._pr)
- filterw= volw**(1._pr/3._pr)
-
-! coefficients for PISO algorithm, not sure which version is better
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- Cb1(i,j,l)=0._pr; Cc1(i,j,l)=0._pr; Cd1(i,j,l)=0._pr
- Cb2(i,j,l)=0._pr; Cc2(i,j,l)=0._pr; Cd2(i,j,l)=0._pr
- if (celltype(i+1,j,l).ne.wall) Cb1(i,j,l)= 1._pr/((xc(i+1)-xc(i))*(xf(i)-xf(i-1)))
- if (celltype(i-1,j,l).ne.wall) Cb2(i,j,l)= 1._pr/((xc(i)-xc(i-1))*(xf(i)-xf(i-1)))
- if (celltype(i,j+1,l).ne.wall) Cc1(i,j,l)= 1._pr/((yc(j+1)-yc(j))*(yf(j)-yf(j-1)))
- if (celltype(i,j-1,l).ne.wall) Cc2(i,j,l)= 1._pr/((yc(j)-yc(j-1))*(yf(j)-yf(j-1)))
- if (celltype(i,j,l+1).ne.wall) Cd1(i,j,l)= 1._pr/((zc(l+1)-zc(l))*(zf(l)-zf(l-1)))
- if (celltype(i,j,l-1).ne.wall) Cd2(i,j,l)= 1._pr/((zc(l)-zc(l-1))*(zf(l)-zf(l-1)))
- Ca(i,j,l)= Cb1(i,j,l)+Cb2(i,j,l)+Cc1(i,j,l)+Cc2(i,j,l)+Cd1(i,j,l)+Cd2(i,j,l)
- enddo; enddo; enddo
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief Initialize the flow to reduce spin-up time (Nelson & Fringer, 2017)
- subroutine initFlow
- use var
- use parallel
- real(kind=pr) yplus,random(gp),uplus,vplus,wplus,dist,alpha,utau
- integer :: i,j,l,m
-
- call syncRandom(gp,random)
-
- alpha=0.25_pr
-
- m=0
- do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
- m=m+1
- if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l)))
- if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
- if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
- u(i,j,l)= (1._pr + (2._pr*random(m)-1._pr)*alpha) * ubulk*dist/delta
- enddo; enddo; enddo
-
- v=0._pr
- w=0._pr
- call bcUVW(u,v,w)
- u01=u; v01=v; w01=w
- u02=u01; v02=v01; w02=w01
-
- end
diff --git a/src/restart.f90 b/src/restart.f90
deleted file mode 100755
index 9130983..0000000
--- a/src/restart.f90
+++ /dev/null
@@ -1,104 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief save fluid and particle field for restart
- subroutine saveField
- use var
- integer :: n
- character(70) :: filename
- integer :: fileformat=-9999003 ! current format identifier
-
- write(filename,'(a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i6.6)') &
- 'restart/fluidField_p',myid,'_',dimi,'_',dimj,'_',diml,'_',nt
- open(unit=11,file=filename,form='unformatted')
- write(11) fileformat
- write(11) u(imin:imax,jmin:jmax,lmin:lmax),u01(imin:imax,jmin:jmax,lmin:lmax), &
- v(imin:imax,jmin:jmax,lmin:lmax),v01(imin:imax,jmin:jmax,lmin:lmax), &
- w(imin:imax,jmin:jmax,lmin:lmax),w01(imin:imax,jmin:jmax,lmin:lmax), &
- p(imin:imax,jmin:jmax,lmin:lmax),t,dtNext,dt,dt01,timecom
- close(11)
-
- if (npTot.gt.0) then
- write(filename,'("restart/particleField_p",i3.3,"_",i6.6)') myid,nt
- open(unit=12,file=filename,form='unformatted')
- write(12) fileformat
- write(12) np
- write(12) (upNext(n),vpNext(n),wpNext(n),xpNext(n),ypNext(n),zpNext(n), &
- radp(n),partn(n),q_elNext(n),wcollnum(n),ppcollnum(n),nGlob(n),n=1,np)
- close(12)
- endif
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief read fluid and particle field for restart, old or new format
- subroutine readField
- use var
- integer :: n,stat
- character(70) :: filename
- integer :: fileformat
-
- if (myid.eq.0) write(*,'(a)') 'read old data'
-
- write(filename,'(a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i6.6)') &
- 'restart/fluidField_p',myid,'_',dimi,'_',dimj,'_',diml,'_',ntstart
- open(unit=11,file=filename,form='unformatted',status='old',iostat=stat)
- if (stat.ne.0) then
- write(*,'(a)') 'old field not existing ',filename
- stop
- endif
- rewind(11)
- read(11) fileformat
- if (fileformat.eq.-9999003) then
- read(11) u(imin:imax,jmin:jmax,lmin:lmax),u01(imin:imax,jmin:jmax,lmin:lmax), &
- v(imin:imax,jmin:jmax,lmin:lmax),v01(imin:imax,jmin:jmax,lmin:lmax), &
- w(imin:imax,jmin:jmax,lmin:lmax),w01(imin:imax,jmin:jmax,lmin:lmax), &
- p(imin:imax,jmin:jmax,lmin:lmax),t,dtNext,dt,dt01,timecom
- call bcUVW(u,v,w)
- call bcUVW(u01,v01,w01)
- call bcNeumann(p)
- timecom(10)=timecom(9)
- elseif (fileformat.eq.-9999002) then ! works with gc=3
- read(11) u,u01,v,v01,w,w01,p,t,dtNext,dt,dt01,timecom
- else ! oldest format, works with gc=3
- rewind(11)
- read(11) u,u01,v,v01,w,w01,p,Fsx,Fsy,Fsz,rho_el,t,dt,dt01
- call nextTimestepSize
- endif
- close(11)
-
- write(filename,'("restart/particleField_p",i3.3,"_",i6.6)') myid,ntstart
- open(unit=12,file=filename,form='unformatted',status='old',iostat=stat)
- if (stat.ne.0) then
- np=0
- call allocateParticleArrays
- else
- rewind(12)
- read(12) fileformat
- if (fileformat.eq.-9999002.or.fileformat.eq.-9999003) then ! newest format
- read(12) np
- call allocateParticleArrays
- read(12) (upNext(n),vpNext(n),wpNext(n),xpNext(n),ypNext(n),zpNext(n), &
- radp(n),partn(n),q_elNext(n),wcollnum(n),ppcollnum(n),nGlob(n),n=1,np)
- else ! oldest format
- rewind(12)
- read(12) np
- call allocateParticleArrays
- read(12) (up(n),vp(n),wp(n),xp(n),yp(n),zp(n),uf(n),vf(n),wf(n), &
- radp(n),partn(n),q_el(n),wcollnum(n),ppcollnum(n),n=1,np)
- xpNext=xp; ypNext=yp; zpNext=zp;
- upNext=up; vpNext=vp; wpNext=wp;
- q_elNext=q_el
- do n=1,np
- nseedLoc=nseedLoc+1
- nGlob(n)=nseedLoc*1000+myid
- enddo
- endif
- endif
- close(12)
-
- ntstart=ntstart+1
- ntend=ntend+1
- tstart=t
-
- end
diff --git a/src/timestep.f90 b/src/timestep.f90
deleted file mode 100755
index 13f2bc3..0000000
--- a/src/timestep.f90
+++ /dev/null
@@ -1,52 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief copy variables between time-steps
- subroutine prepareTimestep
- use var
-
- if (nt.eq.1) then
- dt=cfl*dimx/dimi/10._pr
- dt01=dt
- dt02=dt01
- else
- dt02=dt01
- dt01=dt
- dt=dtNext
- endif
-
-! fixed time-step:
-! dt=cfl*minval(hx)/ubulk
-! dt01=dt
-! dt02=dt
-
- tau01=dt/(dt+dt01)
- tau02=dt01/(dt01+dt02)
- t=t+dt
-
- u02=u01; v02=v01; w02=w01
- u01=u; v01=v; w01=w
- xp=xpNext; yp=ypNext; zp=zpNext;
- up=upNext; vp=vpNext; wp=wpNext;
- q_el=q_elNext
-
- call particlesNextProc
- call particlesToCells
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief next time-step size, needed AFTER flow field and BEFORE particle calculation
- subroutine nextTimestepSize
- use var
- use parallel
- real(kind=pr) umax
-
- umax= syncMax(maxval(abs(u)))
- if (umax.eq.0._pr) then
- dtNext=dt
- else
- dtNext=cfl*dimx/dimi/umax
- endif
-
- end
diff --git a/src/var.f90 b/src/var.f90
deleted file mode 100755
index b54c9a6..0000000
--- a/src/var.f90
+++ /dev/null
@@ -1,195 +0,0 @@
-!> @author Holger Grosshans
-!> @brief definition of public variables, global denotes the complete domain, local only one processor
-!>
-!> how to add a new public variable:
-!> fluid: var,initVariables(allocate+init)
-!> particles: var,allocateParticleArrays,initVariables,particlesNextProc(?),partN2M(?)
-
- module var
- implicit none
-
-! hardcoded parameters which the user cannot change
- character(70) :: version='pafiX v1.1.0 (Copyright 2015-21 by H. Grosshans)'
-
- integer, parameter :: &
- precision=8, & !< number precision
- pr=selected_real_kind(precision), &
- elforceScheme=3, & !< 1=Gauss, 2=Coulomb, 3=Hybrid (Grosshans&Papalexandris, 2017)
- gc= 2 !< number of ghost cells, depends on stencils, but >= 2
-
- real(kind=pr), parameter :: &
- eps_el= 8.85e-12_pr, & !< permittivity vacuum/air (F/m)
- pi= 4._pr*atan(1._pr), &
- urfu= 0.25_pr, & !< under-relaxation factor variable u
- urfv= urfu, &
- urfw= urfv, &
- urfp= 4.0_pr, &
- C_s= 0.17_pr !< Smagorinsky constant
-
- real(kind=pr), parameter :: &
- Ew=1.e11_pr, & !< duct Young's modulus (kg/s**2/m)
- Ep=1.e8_pr, & !< particle Young's modulus (kg/s**2/m)
- nyp=0.4_pr , & !< particle Poisson ratio
- nyw=0.28_pr, & !< duct Poisson ratio
- Qaccfactor= 0.1_pr !< artificially accelerate the charging rate
-
-! input
- real(kind=pr) :: &
- cfl, & !< CFL number (-)
- ubulk, ubulk0, & !< streamwise bulk fluid velocity, initial/inlet (m/s)
- prad, & !< particle radius (m)
- qp0,qpmax, & !< initial and maximum particle charge (C,C)
- dimx,dimy,dimz, & !< dimension of local domain in x-direction (m)
- rhof, & !< fluid density (kg/m**3)
- nuf, & !< fluid kinematic viscosity (m**2/s)
- pnd, & !< particle number density (-/m**3)
- rhop, & !< particle material density (kg/m**3)
- restRatio, & !< particle material restitution ratio
- g(3), & !< gravity vector (m/s**2,m/s**2,m/s**2)
- tol, & !< factor by which L2 shall decrease
- dimxtot, & !< dimension of domain (m)
- delta, & !< duct half width
- tau_w, & !< wall shear stress
- u_tau, & !< friction velocity
- Re, & !< bulk Reynolds number
- Re_tau, & !< friction Reynolds number
- delta_v, & !< viscous length scale (m)
- t_tau, & !< viscous time scale (s)
- ftt !< flow through time (s)
-
- integer :: &
- ntstart, & !< start time step
- ntime, & !< end time step
- dimi,dimj,diml, & !< local number of cells in x-direction
- int_results, & !< interval write out results
- int_restart, & !< interval write out restart files
- itmax, & !< max number of iterations for 1 equation
- dimitot, & !< global number of cells in x-direction
- ntseed !< time step particle seeding
-
- character(11) :: &
- turbModel !< turbulence model
-
-! run
- real(kind=pr) :: &
- dpdx, & !< streamwise pressure gradient (N/m**3)
- dtNext, & !< time-step size from (nt) to (nt+1) (s)
- dt, & !< time-step size from (nt-1) to (nt) (s)
- dt01, & !< time step size from (nt-2) to (nt-1) (s)
- dt02, & !< time step size from (nt-3) to (nt-2) (s)
- tau01,tau02, &
- t, & !< physical time (s)
- tstart !< start time step
-
- integer :: nt, & !< time step
- ntend !< last time step of computation
-
-! grid
- character(1) :: &
- bcx,bcy,bcz, & !< type of boundary condition [(w)all/(p)eriodic]
- gridx,gridy,gridz !< grid [(u)niform/(s)tretch]
-
- real(kind=pr), allocatable, dimension(:,:,:) :: &
- vol,volu,volv,volw !< volume per cell (m**3), centered and staggered variables
-
- real(kind=pr), allocatable, dimension(:) :: &
- xc,yc,zc, & !< center point of a cell (m)
- xf,yf,zf !< face of a cell in pos. direction from xc (m)
-
- real(kind=pr) :: &
- xmin,ymin,zmin, & !< lower local boundary in global coordinates (m)
- xmax,ymax,zmax, & !< upper local boundary in global coordinates (m)
- deltax, & !< shift between processors (m)
- volTot !< total domain volume (m**3)
-
- integer, allocatable, dimension(:,:,:) :: &
- celltype !< celltype
-
- integer, parameter :: active=1, passive=2, wall=3
-
- integer :: &
- gp, & !< local number of grid cells
- dimgp, & !< local number of cells containing fluid
- dimgptot, & !< global number of cells containing fluid
- ii,jj,ll, & !< max index local grid
- imin,jmin,lmin, & !< local min index of cell containing fluid
- imax,jmax,lmax !< local max index of cell containing fluid
-
-! pre-computed coefficients
- real(kind=pr), allocatable, dimension(:,:,:) :: &
- Ca,Cb1,Cb2,Cc1,Cc2,Cd1,Cd2 !< coefficients for SIMPLE algorithm
-
-! fluid
- real(kind=pr) :: muf !< fluid dynamic viscosity (kg/s/m)
-
- real(kind=pr), allocatable, dimension(:,:,:) :: &
- u,v,w, & !< fluid velocity (m/s)
- p, & !< fluid pressure (N/m**2)
- u01,v01,w01, & !< fluid velocity previous time-step
- u02,v02,w02, & !< fluid velocity previous time-step
- Fsx,Fsy,Fsz, & !< source term particles -> fluid (m/s**2)
- nuf_ru,nuf_rv,nuf_rw, & !< turbulent viscosity (m**2/s)
- filteru,filterv,filterw !< LES filter width (m)
-
-! particles
- real(kind=pr) :: &
- tau_p_min, & !< min particle response time in a time-step
- rtau_el_max, & !< max reciprocal Coulomb force time-scale in a time-step
- dup_max !< max particle velocity change in a time-step
-
- integer :: &
- np, & !< local number of particles
- npTot, & !< total number of particles
- npp, & !< local number of particles before some operation
- maxnp, & !< maximum possible local number of particles
- npstart, & !< number of particles seeded since tstart
- nseedLoc !< local counter for seeded particles
-
- real(kind=pr), allocatable, dimension(:) :: &
- up,vp,wp, & !< particle velocity at t (m/s)
- upNext,vpNext,wpNext, & !< particle velocity at t+dt (m/s)
- uf,vf,wf, & !< fluid velocity around a particle (m/s)
- dufdy,dufdz, & !< fluid velocity gradient around a particle (m/s)
- xp,yp,zp, & !< particle position at t (m)
- xpNext,ypNext,zpNext, & !< particle position at t+dt (m)
- radp, & !< particle radius (m)
- q_el, & !< particle charge at t (C)
- q_elNext, & !< particle charge at t+dt (C)
- fx_el,fy_el,fz_el, & !< electrostatic force on particle (m/s**2)
- fx_d,fy_d,fz_d, & !< drag force on particle (m/s**2)
- fx_l,fy_l,fz_l, & !< lift force on particle (m/s**2)
- partn !< number of particles represented by parcel
-
- integer, allocatable, dimension(:) :: &
- ip,jp,lp, & !< Eulerian indice of particle position
- wcollnum, & !< number of particle-wall collisions
- ppcollnum, & !< number of particle-particle collisions
- nGlob !< global particle id
-
- integer, allocatable, dimension(:,:,:) :: &
- npic !< number of particles in a Eulerian cell
-
- integer, allocatable, dimension(:,:,:,:) :: &
- nic !< indices of particles in a Eulerian cell
-
-! parallel
- real(kind=pr) :: timebeg,timenow,timeend, &
- timecom(10) !< 1-8: sync operations, 9: physical time, 10: restart
-
- integer :: &
- mpi_pr, &
- myid, & !< id of current processor
- nrprocs, & !< total number of processors
- mpierr, &
- next,prev !< id of next and previous processor
-
- integer, allocatable, dimension(:) :: &
- mpistatus
-
-! electrostatics
- real(kind=pr), allocatable, dimension(:,:,:) :: &
- rho_el, & !< electric charge density (C/m**3)
- phi_el, & !< electric potential (V)
- Ex_el,Ey_el,Ez_el !< electric field strength (V/m)
-
- end module var
diff --git a/src/write_vtk.f90 b/src/write_vtk.f90
deleted file mode 100755
index bb3af38..0000000
--- a/src/write_vtk.f90
+++ /dev/null
@@ -1,450 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief write out particle field in vtk format
- subroutine writevtk_particles
- use var
- integer :: n
- character(80) :: filename,rowfmt,rowfm2
-
-100 format(es13.4e2)
-101 format(3(es13.4e2))
-102 format(3(es14.5e2))
- rowfmt='(10(1x,es11.4e2))'
- rowfm2='(10(1x,i9))'
-
- if(myid.eq.0) call write_visit_container('particles')
-
- write(filename,'("results/particles_p",i3.3,"_",i6.6,".vtk")') myid,nt
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 3.0'
- write(10,'(a,es14.6e2,a,x,a)') 'particle data: t=',t,' s --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET UNSTRUCTURED_GRID'
- write(10,'(a6,i9,a6)') 'POINTS ',np,' FLOAT'
- write(10,102) (xp(n),yp(n),zp(n),n=1,np)
- write(10,'(a5,i9,i9)') 'CELLS ',np,2*np
- write(10,*) ('1 ',n-1,n=1,np)
- write(10,'(a10,i9)') 'CELL_TYPES ',np
- write(10,*) ('1 ', n=1,np)
- write(10,'(a,i9)') 'POINT_DATA ',np
-
- call writeVar('radp',radp)
- call writeVar('up',up)
- call writeVar('vp',vp)
- call writeVar('wp',wp)
- call writeVar('uf',uf)
- call writeVar('vf',vf)
- call writeVar('wf',wf)
- call writeVar('partn',partn)
- call writeVar('q_el',q_el)
- call writeVar('fx_d',fx_d)
- call writeVar('fy_d',fy_d)
- call writeVar('fz_d',fz_d)
- call writeVar('fx_l',fx_l)
- call writeVar('fy_l',fy_l)
- call writeVar('fz_l',fz_l)
- call writeVar('fx_el',fx_el)
- call writeVar('fy_el',fy_el)
- call writeVar('fz_el',fz_el)
- call writeVarI('wcollnum',wcollnum)
- call writeVarI('ppcollnum',ppcollnum)
-
- write(10,'(a)') 'VECTORS UPVPWP FLOAT'
- write(10,101) (up(n),vp(n),wp(n),n=1,np)
- write(10,'(a)')
-
- close(10)
-
- contains
-
- subroutine writeVar(name,myvar)
- use var
- real(kind=pr) :: myvar(:)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- write(10,fmt=rowfmt) (myvar(n),n=1,np)
- write(10,'(a)')
-
- end
-
- subroutine writeVarI(name,myvar)
- use var
- integer :: myvar(:)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- write(10,fmt=rowfm2) (myvar(n),n=1,np)
- write(10,'(a)')
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine writevtk_fluid_xyz
- use var
- real(kind=pr) :: temp(ii,jj,ll)
- integer :: i,j,l
- character(70) :: filename,rowfmt,vecfmt
-
-100 format(es11.4e2)
- write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
- write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
-
- if(myid.eq.0) call write_visit_container('fluid_xyz')
-
- write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xyz_p',myid,'_',nt,'.vtk'
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 3.0'
- write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a10,3(i4))') 'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),(lmax-lmin+3)
- write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
- write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
- write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
- write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
- write(10,*) 'Z_COORDINATES',(lmax-lmin+3),'float'
- write(10,fmt=rowfmt) (zc(l),l=lmin-1,lmax+1)
- write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3)*(lmax-lmin+3)
-
- temp=0._pr
- temp(2:ii,:,:)=(u(2:ii,:,:)+u(1:ii-1,:,:))/2._pr
- call writeVar('u',temp)
- temp(:,2:jj,:)=(v(:,2:jj,:)+v(:,1:jj-1,:))/2._pr
- call writeVar('v',temp)
- temp(:,:,2:ll)=(w(:,:,2:ll)+w(:,:,1:ll-1))/2._pr
- call writeVar('w',temp)
- call writeVar('p',p)
- call writeVar('rho_el',rho_el)
- call writeVar('phi_el',phi_el)
- call writeVar('Ex_el',Ex_el)
- call writeVar('Ey_el',Ey_el)
- call writeVar('Ez_el',Ez_el)
-
- !temp=sqrt(Ex_el**2+Ey_el**2+Ez_el**2)
- !call writeVar('Emag',temp)
-
-! write(10,'(a)') 'VECTORS uvw float'
-! do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
-! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
-! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
-! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
-! enddo; enddo
-
-! write(10,'(a)') 'VECTORS coordinates float 1'
-! write(10,'(a)') 'LOOKUP_TABLE default '
-! do l=lmin-1,lmax+1
-! do i=imin,imax+1
-! write(10,'(3(i8))') i,j,l
-! enddo
-! enddo
-
- close(10)
-
- contains
-
- subroutine writeVar(name,myvar)
- use var
- real(kind=pr) :: myvar(ii,jj,ll)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do l=lmin-1,lmax+1; do j=jmin-1,jmax+1
- write(10,fmt=rowfmt) (myvar(i,j,l),i=imin,imax+1)
- enddo; enddo
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine writevtk_fluid_xy
- use var
- real(kind=pr) :: temp(ii,jj,ll)
- integer :: i,j,l
- character(70) :: filename,rowfmt,vecfmt
-
-100 format(es11.4e2)
- write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
- write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
-
- l=int(ll/2.)+1
-
- if(myid.eq.0) call write_visit_container('fluid_xy')
-
- write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xy_p',myid,'_',nt,'.vtk'
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 3.0'
- write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a10,i4,i4,a4)') 'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),'1'
- write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
- write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
- write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
- write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
- write(10,*) 'Z_COORDINATES 1 float'
- write(10,100) zc(l)
- write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3)
-
- temp=0._pr
- temp(2:ii,:,:)=(u(2:ii,:,:)+u(1:ii-1,:,:))/2._pr
- call writeVar('u',temp)
- temp(:,2:jj,:)=(v(:,2:jj,:)+v(:,1:jj-1,:))/2._pr
- call writeVar('v',temp)
- temp(:,:,2:ll)=(w(:,:,2:ll)+w(:,:,1:ll-1))/2._pr
- call writeVar('w',temp)
- !call writeVar('uraw',u)
- !call writeVar('vraw',v)
- !call writeVar('wraw',w)
- call writeVar('p',p)
- call writeVar('rho_el',rho_el)
- call writeVar('phi_el',phi_el)
- call writeVar('Ex_el',Ex_el)
- call writeVar('Ey_el',Ey_el)
- call writeVar('Ez_el',Ez_el)
- !temp=sqrt(Ex_el**2+Ey_el**2+Ez_el**2)
- !call writeVar('Emag',temp)
-
-! write(10,'(a)') 'VECTORS uvw float'
-! do j=jmin-1,jmax+1
-! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
-! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
-! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
-! enddo
-
- close(10)
-
- contains
-
- subroutine writeVar(name,myvar)
- use var
- real(kind=pr) :: myvar(ii,jj,ll)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do j=jmin-1,jmax+1
- write(10,fmt=rowfmt) (myvar(i,j,l),i=imin,imax+1)
- enddo
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine writevtk_fluid_xz
- use var
- real(kind=pr) :: temp(ii,jj,ll)
- integer :: i,j,l
- character(70) :: filename,rowfmt,vecfmt
-
-100 format(es11.4e2)
- write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
- write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
-
- j=int(jj/2.)+1
-
- if(myid.eq.0) call write_visit_container('fluid_xz')
-
- write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xz_p',myid,'_',nt,'.vtk'
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 3.0'
- write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a10,i4,a4,i4)') 'DIMENSIONS',(imax-imin+2),'1',(lmax-lmin+3)
- write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
- write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
- write(10,*) 'Y_COORDINATES 1 float'
- write(10,100) yc(j)
- write(10,*) 'Z_COORDINATES',(lmax-lmin+3),'float'
- write(10,fmt=rowfmt) (zc(l),l=lmin-1,lmax+1)
- write(10,*) 'POINT_DATA ',(imax-imin+2)*(lmax-lmin+3)
-
- temp=0._pr
- temp(2:ii,:,:)=(u(2:ii,:,:)+u(1:ii-1,:,:))/2._pr
- call writeVar('u',temp)
- temp(:,2:jj,:)=(v(:,2:jj,:)+v(:,1:jj-1,:))/2._pr
- call writeVar('v',temp)
- temp(:,:,2:ll)=(w(:,:,2:ll)+w(:,:,1:ll-1))/2._pr
- call writeVar('w',temp)
- !call writeVar('uraw',u)
- !call writeVar('vraw',v)
- !call writeVar('wraw',w)
- call writeVar('p',p)
- call writeVar('rho_el',rho_el)
- call writeVar('phi_el',phi_el)
- call writeVar('Ex_el',Ex_el)
- call writeVar('Ey_el',Ey_el)
- call writeVar('Ez_el',Ez_el)
- !temp=sqrt(Ex_el**2+Ey_el**2+Ez_el**2)
- !call writeVar('Emag',temp)
-
-! write(10,'(a)') 'VECTORS uvw float'
-! do l=lmin-1,lmax+1
-! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
-! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
-! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
-! enddo
-
- close(10)
-
- contains
-
- subroutine writeVar(name,myvar)
- use var
- real(kind=pr) :: myvar(ii,jj,ll)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do l=lmin-1,lmax+1
- write(10,fmt=rowfmt) (myvar(i,j,l),i=imin,imax+1)
- enddo
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine writevtk_fluid_yz
- use var
- real(kind=pr) :: temp(ii,jj,ll)
- integer :: i,j,l
- character(70) :: filename,rowfmt
-
-100 format(es11.4e2)
- write(rowfmt,'(a,i4,a)') '(',(jmax-jmin+3),'(1x,es11.4e2))'
-
- i=int(ii/2.)+1
-
- if(myid.eq.0) call write_visit_container('fluid_yz')
-
- write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_yz_p',myid,'_',nt,'.vtk'
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 3.0'
- write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a12,i8,i8)') 'DIMENSIONS 1',(jmax-jmin+3),(lmax-lmin+3)
- write(10,*) 'X_COORDINATES 1 float'
- write(10,100) xc(i)
- write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
- write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
- write(10,*) 'Z_COORDINATES',(lmax-lmin+3),'float'
- write(10,fmt=rowfmt) (zc(l),l=lmin-1,lmax+1)
- write(10,*) 'POINT_DATA ',(jmax-jmin+3)*(lmax-lmin+3)
-
- temp=0._pr
- temp(2:ii,:,:)=(u(2:ii,:,:)+u(1:ii-1,:,:))/2._pr
- call writeVar('u',temp)
- temp(:,2:jj,:)=(v(:,2:jj,:)+v(:,1:jj-1,:))/2._pr
- call writeVar('v',temp)
- temp(:,:,2:ll)=(w(:,:,2:ll)+w(:,:,1:ll-1))/2._pr
- call writeVar('w',temp)
- !call writeVar('uraw',u)
- !call writeVar('vraw',v)
- !call writeVar('wraw',w)
- call writeVar('p',p)
- call writeVar('rho_el',rho_el)
- call writeVar('phi_el',phi_el)
- call writeVar('Ex_el',Ex_el)
- call writeVar('Ey_el',Ey_el)
- call writeVar('Ez_el',Ez_el)
- !temp=sqrt(Ex_el**2+Ey_el**2+Ez_el**2)
- !call writeVar('Emag',temp)
-
- close(10)
-
- contains
-
- subroutine writeVar(name,myvar)
- use var
- real(kind=pr) :: myvar(ii,jj,ll)
- character(*) :: name
-
- write(10,'(3a)') 'SCALARS ',name,' float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do l=lmin-1,lmax+1
- write(10,fmt=rowfmt) (myvar(i,j,l),j=jmin-1,jmax+1)
- enddo
-
- end
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
- subroutine writevtk_grid
- use var
- integer :: i,j,l
- character(70) :: filename,rowfmt
-
-100 format(es11.4e2)
- write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
-
- write(filename,'("output/grid_p",i3.3,".vtk")') myid
- open(10,file=filename)
-
- write(10,'(a)') '# vtk DataFile Version 2.0'
- write(10,'(a,x,a)') 'grid --',version
- write(10,'(a)') 'ASCII'
- write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,*) 'DIMENSIONS',ii,jj,ll
- write(10,*) 'X_COORDINATES',ii,'float'
- write(10,fmt=rowfmt) (xc(i),i=1,ii)
- write(10,*) 'Y_COORDINATES',jj,'float'
- write(10,fmt=rowfmt) (yc(j),j=1,jj)
- write(10,*) 'Z_COORDINATES',ll,'float'
- write(10,fmt=rowfmt) (zc(l),l=1,ll)
- write(10,*) 'POINT_DATA ',ii*jj*ll
-
- write(10,*) 'SCALARS celltype int 1'
- write(10,*) 'LOOKUP_TABLE default '
- do l=1,ll; do j=1,jj; do i=1,ii
- write(10,'(i2)') celltype(i,j,l)
- enddo; enddo; enddo
-
- close(10)
-
- end
-
-!####################################################################
-!> @author Holger Grosshans
-!> @brief write visit container file
- subroutine write_visit_container(filename)
- use var
- integer :: m
- character(*) :: filename
- character(80) :: filename2
- logical :: file_ex
-
- write(filename2,'(3a)') 'results/',filename,'.visit'
- inquire(file=filename2,exist=file_ex)
- if (file_ex.and.nt.ne.1) then
- open(11,file=filename2,access='append')
- else
- open(11,file=filename2)
- write(11,'(a8,x,i3)') '!NBLOCKS',nrprocs
- endif
- do m=1,nrprocs
- write(11,'(2a,i3.3,a,i6.6,a)') filename,'_p',(m-1),'_',nt,'.vtk'
- enddo
- close(11)
-
- end
diff --git a/src/writedat_particles.f90 b/src/writedat_particles.f90
deleted file mode 100755
index 5374c5a..0000000
--- a/src/writedat_particles.f90
+++ /dev/null
@@ -1,32 +0,0 @@
-!####################################################################
-!> @author Holger Grosshans
-!> @brief write out particle field in dat format
- subroutine writedat_particles
- use var
- real(kind=pr) :: partmass
- integer n
- character(80) :: filename
-
-100 format(22(es13.4e2,x),2(i6,x))
-
- write(filename,'("results/particles_p",i3.3,"_",i6.6,".dat")') myid,nt
- open(10,file=filename)
-
- write(10,'(a,es14.6e2,a,x,a)') '# t=',t,' s --',version
-
- do n=1,np
- partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3
- write(10,100) xp(n),yp(n),zp(n), &
- radp(n), &
- up(n),vp(n),wp(n), &
- uf(n),vf(n),wf(n), &
- partn(n),q_el(n),q_el(n)/partmass, &
- fx_el(n),fy_el(n),fz_el(n), &
- fx_d(n),fy_d(n),fz_d(n), &
- fx_l(n),fy_l(n),fz_l(n), &
- wcollnum(n),ppcollnum(n)
- enddo
-
- close(10)
-
- end
--
GitLab