Merge branch 'develop' into feature/openfoam6

Just resolve merge conflicts, does not yet compile with OF6 !

applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C

.circleci/config.yml

@@ -0,0 +1,72 @@
+version: 2
+jobs:
+  build:
+    branches:
+      only:
+        - master
+        - develop
+
+    docker:
+      - image: ubuntu:trusty
+
+    environment:
+      WM_NCOMPPROCS: 2
+
+    working_directory: /root/CFDEM/CFDEMcoupling
+
+    steps:
+      - run:
+          name: Install package dependencies
+          command: sudo apt-get update && sudo apt-get install -y build-essential cmake openmpi-bin libopenmpi-dev python-dev git bc
+          
+      - run:
+          name: Make project and user dir
+          command: mkdir -p /root/CFDEM/CFDEMcoupling && mkdir -p /root/CFDEM/-4.1
+
+      - checkout:
+          path: /root/CFDEM/CFDEMcoupling
+
+      - run:
+          name: Add OpenFOAM package repository
+          command: sudo apt-get install -y software-properties-common wget apt-transport-https && sudo add-apt-repository http://dl.openfoam.org/ubuntu && sudo sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -"
+
+      - run:
+          name: Install OpenFOAM 4.1
+          command: sudo apt-get update && sudo apt-get -y install openfoam4
+
+      - run:
+          name: Clone LIGGGHTS repository
+          command: git clone https://github.com/ParticulateFlow/LIGGGHTS-PFM.git /root/CFDEM/LIGGGHTS
+
+      - run:
+          name: Build LIGGGHTS
+          command: >
+            shopt -s expand_aliases &&
+            source /opt/openfoam4/etc/bashrc &&
+            source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
+            bash /root/CFDEM/CFDEMcoupling/etc/compileLIGGGHTS.sh
+          no_output_timeout: 30m
+
+      - run:
+          name: Build CFDEMcoupling library
+          command: >
+            shopt -s expand_aliases &&
+            source /opt/openfoam4/etc/bashrc &&
+            source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
+            bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_src.sh
+          
+      - run:
+          name: Build CFDEMcoupling solvers
+          command: >
+            shopt -s expand_aliases &&
+            source /opt/openfoam4/etc/bashrc &&
+            source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
+            bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_sol.sh
+
+      - run:
+          name: Build CFDEMcoupling utilities
+          command: >
+            shopt -s expand_aliases &&
+            source /opt/openfoam4/etc/bashrc &&
+            source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
+            bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_uti.sh

.gitignore

@@ -5,6 +5,8 @@
 log_*
 log.*
 *~
+*.swp
+*.swo
 
 **/linux*Gcc*/
 **/.vscode

LICENSE

@@ -0,0 +1,674 @@
+                    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>.

README

@@ -1,80 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 2012-2015 DCS Computing GmbH, Linz
-                                Copyright 2015-     JKU Linz
--------------------------------------------------------------------------------
-License
-    This file is part of CFDEMcoupling.
-
-    CFDEMcoupling 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.
-
-    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code provides models and solvers to realize coupled CFD-DEM simulations
-    using LIGGGHTS and OpenFOAM.
-    Note: this code is not part of OpenFOAM (see DISCLAIMER).
-\*---------------------------------------------------------------------------*/
-
-
-CFDEM(R) coupling provides an open source parallel coupled CFD-DEM framework 
-combining the strengths of the LIGGGHTS(R) DEM code and the Open Source 
-CFD package OpenFOAM(R)(*). The CFDEM(R)coupling toolbox allows to expand 
-standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
-code LIGGGHTS(R). In this toolbox the particle representation within the 
-CFD solver is organized by "cloud" classes. Key functionalities are organised 
-in sub-models (e.g. force models, data exchange models, etc.) which can easily 
-be selected and combined by dictionary settings.
-
-The coupled solvers run fully parallel on distributed-memory clusters. 
-
-Features are:
-
-- its modular approach allows users to easily implement new models
-- its MPI parallelization enables to use it for large scale problems
-- the use of GIT allows to easily update to the latest version
-- basic documentation is provided
-
-The file structure:
-
-- "src" directory including the source files of the coupling toolbox and models
-- "applications" directory including the solver files for coupled CFD-DEM simulations
-- "doc" directory including the documentation of CFDEM(R)coupling
-- "tutorials" directory including basic tutorial cases showing the functionality
-
-
-
-Details on installation are given on the "www.cfdem.com"
-
-The functionality of this CFD-DEM framwork is described via "tutorial cases" showing 
-how to use different solvers and models.
-
-CFDEM(R)coupling stands for Computational Fluid Dynamics (CFD) -
-Discrete Element Method (DEM) coupling.
-
-CFDEM(R)coupling is an open-source code, distributed freely under the terms of the 
-GNU Public License (GPL).
-
-Core development of CFDEM(R)coupling is done by 
-Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
-
-
-/*---------------------------------------------------------------------------*\
-(*) "OpenFOAM(R)" is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.
-This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks.
-\*---------------------------------------------------------------------------*/

README.md

@@ -0,0 +1,33 @@
+# CFDEMcoupling
+
+CFDEM®coupling stands for Computational Fluid Dynamics (CFD) - Discrete Element Method (DEM) coupling. It combines the open source packages OpenFOAM® (CFD) and LIGGGHTS® (DEM) to simulate particle-laden flows. CFDEM®coupling is part of the [CFDEM®project](https://www.cfdem.com).
+
+[![CircleCI](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling.svg?style=shield&circle-token=e4b6af30d3aa7aee109d206116f01600bf9ee9c6)](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling)
+[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)
+
+## Disclaimer
+
+> This is an academic adaptation of the CFDEM®coupling software package, released by the
+[Department of Particulate Flow Modelling at Johannes Kepler University in Linz, Austria.](https://www.jku.at/pfm)
+> LIGGGHTS® and CFDEM® are registered trademarks, and this offering is not approved or
+endorsed by DCS Computing GmbH, the official producer of the LIGGGHTS® and CFDEM®coupling software.
+> This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software via www.openfoam.com, and owner of the OPENFOAM®  and OpenCFD®  trade marks.
+
+## Features
+
+- Documentation and tutorials to get started
+- A modular approach that allows for easy implementation of new models
+- MPI parallelization for large scale problems
+
+## License
+
+[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)
+
+- This software is distributed under the [GNU General Public License](https://opensource.org/licenses/GPL-3.0).
+- Copyright © 2009-     JKU Linz
+- Copyright © 2012-2015 DCS Computing GmbH, Linz
+- Some parts of CFDEM®coupling are based on OpenFOAM® and Copyright on these
+  parts is held by the OpenFOAM® Foundation (www.openfoam.org)
+  and potentially other parties.
+- Some parts of CFDEM®coupling are contributed by other parties, which are
+  holding the Copyright. This is listed in each file of the distribution.

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -6,9 +6,6 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();
-    thCond=particleCloud.thermCondM().thermCond();
-
     addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
       + (
             he.name() == "e"
@@ -35,7 +32,6 @@
       - Qsource
       - fvm::Sp(QCoeff/Cpv, he)
      // thermal conduction of the fluid with effective conductivity
-     // - fvm::laplacian(rhoeps*thDiff,he)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
       // + particle-fluid energy transfer due to work
       // + fluid energy dissipation due to shearing

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -58,10 +58,11 @@ Info<< "Reading thermophysical properties\n" << endl;
            "addSource",
            runTime.timeName(),
            mesh,
-           IOobject::MUST_READ,
+           IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
@@ -94,21 +95,6 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
     );
 
-  /*  Info<< "\nCreating thermal diffusivity field\n" << endl;
-    volScalarField thDiff
-    (
-        IOobject
-        (
-            "thDiff",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
-    );
-  */
     Info<< "\nCreating thermal conductivity field\n" << endl;
     volScalarField thCond
     (
@@ -117,11 +103,12 @@ Info<< "Reading thermophysical properties\n" << endl;
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -6,12 +6,11 @@ volScalarField& he = thermo.he();
 particleCloud.energyContributions(Qsource);
 particleCloud.energyCoefficients(QCoeff);
 
-thCond=particleCloud.thermCondM().thermCond();
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
 // correct source for the thermodynamic reference temperature
-dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-Qsource += QCoeff*Tref;
+// dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+// Qsource += QCoeff*Tref;
 
 fvScalarMatrix EEqn
 (
@@ -49,7 +48,13 @@ fvScalarMatrix EEqn
 
     thermo.correct();
 
+    Info << "Qsource :" << max(Qsource).value() << " " << min(Qsource).value() << endl;
+    Info << "QCoeff :" << max(QCoeff).value() << " " << min(QCoeff).value() << endl;
+    Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
     Info<<  "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-    Info << "he min/max : " << max(he).value() << "  " << min(he).value() << endl;
+    Info << "he max/min : " << max(he).value() << "  " << min(he).value() << endl;
 
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -53,7 +53,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-//    #include "postProcess.H"
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -62,7 +61,6 @@ int main(int argc, char *argv[])
     #include "createRDeltaT.H"
 
     #include "createFields.H"
-    #include "createFieldRefs.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -28,6 +28,8 @@
     }
 
     volScalarField& p = thermo.p();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -51,9 +53,6 @@
     );
 
     //  kinematic fields
-
-
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -82,18 +81,8 @@
         mesh
     );
 
-    volScalarField rhoeps
-    (
-        IOobject
-        (
-            "rhoeps",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        rho*voidfraction
-    );
+    volScalarField rhoeps ("rhoeps", rho*voidfraction);
+
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
@@ -133,11 +122,12 @@
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;
@@ -282,9 +272,18 @@
         mesh,
         dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
     );
+
+    volScalarField dSauter
+    (
+        IOobject
+        (
+            "dSauter",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
+    );
 //===============================
-
-
-
- //   singlePhaseTransportModel laminarTransport(U, phi);
-

applications/solvers/cfdemSolverRhoSimple/Make/files

@@ -1,3 +0,0 @@
-cfdemSolverRhoSimple.C
-
-EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSimple

applications/solvers/rStatAnalysis/Make/files

@@ -0,0 +1,3 @@
+rStatAnalysis.C
+
+EXE=$(CFDEM_APP_DIR)/rStatAnalysis

applications/solvers/rStatAnalysis/Make/options

@@ -0,0 +1,26 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rStatAnalysis/rStatAnalysis.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rStatAnalysis
+
+Description
+    Creates and analyzes a recurrence statistics
+
+\*---------------------------------------------------------------------------*/
+
+// #include "fvCFD.H"
+// #include "singlePhaseTransportModel.H"
+// #include "turbulentTransportModel.H"
+// #include "fvOptions.H"
+
+#include "recBase.H"
+#include "recStatAnalysis.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+
+
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nAnalyzing recurrence statistics\n" << endl;
+
+    recurrenceBase.recStatA().init();
+    recurrenceBase.recStatA().statistics();
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverBase/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverBase.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverBase

applications/solvers/rcfdemSolverBase/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverBase/createFields.H

@@ -0,0 +1,125 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info << "\nCreating fields subject to calculation\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfractionRec
+    );
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        UsRec
+    );
+
+    // write fields for t=t_start
+    voidfraction.write();
+    Us.write();
+//===============================
+
+    Info << "Calculating face flux field phi\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(URec*voidfractionRec) & mesh.Sf()
+    );
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );

applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverRecurrence
+
+Description
+    Solves a transport equation for a passive scalar on a two-phase solution
+    Test-bed for a solver based on recurrence statistics
+
+Rules
+    Solution data to compute the recurrence statistics from, needs to
+        reside in $CASE_ROOT/dataBase
+    Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
+    cfdemCloudRec<cfdemCloud> particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,URec);
+
+        particleCloud.clockM().stop("Coupling");
+
+
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        particleCloud.clockM().start(27,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverBase/readFields.H

@@ -0,0 +1,4 @@
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolVectorField("U",URec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);

applications/solvers/rcfdemSolverCoupledHeattransfer/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverCoupledHeattransfer.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverCoupledHeattransfer

applications/solvers/rcfdemSolverCoupledHeattransfer/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H

@@ -0,0 +1,35 @@
+    volScalarField rhoeps = rhoRec*voidfractionRec;
+
+    particleCloud.energyContributions(Qsource);
+
+    particleCloud.energyCoefficients(QCoeff);
+
+    //K = 0.5*magSqr(URec);
+
+    addSource = fvc::div(phiRec/fvc::interpolate(rhoRec), pRec);
+    // main contribution due to gas expansion, not due to transport of kinetic energy
+    // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K)
+
+    fvScalarMatrix TEqn =
+    (
+        fvm::ddt(rhoeps, T)
+        + fvm::div(phiRec, T)
+        + addSource/Cv
+        - fvm::laplacian(voidfractionRec*thCond/Cv, T)
+        - Qsource/Cv
+        - fvm::Sp(QCoeff/Cv, T)
+        ==
+        fvOptions(rhoeps, T)    // no fvOptions support yet
+    );
+
+    fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.solve();
+
+    particleCloud.clockM().start(31,"postFlow");
+    counter++;
+
+    if((counter - couplingSubStep) % dtDEM2dtCFD == 0)
+        particleCloud.postFlow();
+
+    particleCloud.clockM().stop("postFlow");

applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H

@@ -0,0 +1,230 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure field\n" << endl;
+    volScalarField pRec
+    (
+        IOobject
+        (
+            "pRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-2,0,0,0,0), 0.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+
+    volScalarField rhoRec
+    (
+        IOobject
+        (
+            "rhoRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+    );
+
+    // heat transfer fields
+    Info << "\nCreating heat transfer fields.\n" << endl;
+
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    volScalarField QCoeff
+    (   IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
+    );
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info << "\nCreating fields subject to calculation\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfractionRec
+    );
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        UsRec
+    );
+
+    // write fields for t=t_start
+    voidfraction.write();
+    Us.write();
+//===============================
+
+
+    Info << "Calculating face flux field phiRec\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,0,-1,0,0,0,0), 0.0)
+    );
+    phiRec.write();
+
+
+    Info << "Creating dummy turbulence model\n" << endl;
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+
+
+    const IOdictionary& transportProps = mesh.lookupObject<IOdictionary>("transportProperties");
+    dimensionedScalar molMass(transportProps.lookup("molM"));
+    // need to scale R down with 1e3 because return value of RR in g, not kg
+    dimensionedScalar R("R",dimensionSet(0,2,-2,-1,0,0,0),Foam::constant::thermodynamic::RR / (1e3*molMass.value()));
+    Info << "specific gas constant R = " << R << endl;
+
+    dimensionedScalar Cv(transportProps.lookup("Cv"));
+
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::READ_IF_PRESENT,
+           IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+
+   // place to put weight functions
+    IOdictionary weightDict
+    (
+        IOobject
+        (
+            "weightDict",
+            runTime.constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    );
+    weightDict.add("weights",scalarList(1,1.0));

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rcfdemSolverHeattransfer
+
+Description
+    Solves heat transfer between fluid and particles based on rCFD
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fvOptions.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
+
+#include "cfdemCloudEnergy.H"
+#include "clockModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
+    #include "updateRho.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    // control coupling behavior in case of substepping
+    // assumes constant timestep size
+    label counter = 0;
+    label couplingSubStep = recurrenceBase.couplingSubStep();
+    double dtProp =  particleCloud.dataExchangeM().couplingTime() / runTime.deltaTValue();
+    label dtDEM2dtCFD = int(dtProp + 0.5);
+    Info << "deltaT_DEM / deltaT_CFD = " << dtDEM2dtCFD << endl;
+    if (dtDEM2dtCFD > 1)
+        Info << "coupling at substep " << couplingSubStep << endl;
+
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,URec);
+
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+        #include "updateRho.H"
+        #include "TEqImp.H"
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(32,"ReadFields");
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+        }
+        particleCloud.clockM().stop("ReadFields");
+
+        particleCloud.clockM().start(33,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverCoupledHeattransfer/updateFields.H

@@ -0,0 +1,38 @@
+// get current weights for various databases
+
+// A: triggered over current value of boundary field
+// word boundaryName = "inlet";
+// label myinlet = mesh.boundary().findPatchID(boundaryName);
+// label startIndex = mesh.boundary()[boundaryName].start();
+
+// B: explicitly define weights
+
+scalarList wList(weightDict.lookupOrDefault("weights",scalarList(1,0.0)));
+
+recurrenceBase.recP().updateIntervalWeights(wList);
+
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolScalarField("p",pRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);
+
+Info << "current database weights: = " << wList << endl;
+Info << "current database: " << recurrenceBase.recM().currDataBase() << endl;
+for(int i=0;i<wList.size();i++)
+{
+    scalar w = wList[i];
+    if (recurrenceBase.recM().currDataBase() == i) w -= 1.0;
+    phiRec += w*recurrenceBase.recM().exportSurfaceScalarFieldAve("phi",i)();
+}
+
+{
+    volScalarField& NuField(const_cast<volScalarField&>(mesh.lookupObject<volScalarField> ("NuField")));
+    recurrenceBase.recM().exportVolScalarField("NuField",NuField);
+}
+

applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H

@@ -0,0 +1 @@
+rhoRec = pRec / (T * R);

applications/solvers/rcfdemSolverHeattransfer/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverHeattransfer.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverHeattransfer

applications/solvers/rcfdemSolverHeattransfer/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverHeattransfer/TEqImp.H

@@ -0,0 +1,40 @@
+    volScalarField rhoeps = rhoRec*voidfractionRec;
+
+    particleCloud.energyContributions(Qsource);
+
+    particleCloud.energyCoefficients(QCoeff);
+
+    K = 0.5*magSqr(URec);
+
+    addSource = fvc::ddt(rhoeps, K) + fvc::div(phiRec, K) +
+    fvc::div
+    (
+        fvc::absolute(phiRec/fvc::interpolate(rhoRec), voidfractionRec*URec), pRec
+    );
+
+
+    fvScalarMatrix TEqn =
+    (
+        fvm::ddt(rhoeps, T)
+        + fvm::div(phiRec, T)
+        + addSource/Cv
+        - fvm::laplacian(voidfractionRec*thCond/Cv, T)
+        - Qsource/Cv
+        - fvm::Sp(QCoeff/Cv, T)
+        ==
+        fvOptions(rhoeps, T)    // no fvOptions support yet
+    );
+
+    //TEqn.relax(relaxCoeff);
+
+    fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.solve();
+
+    particleCloud.clockM().start(31,"postFlow");
+    counter++;
+
+    if((counter - couplingSubStep) % dtDEM2dtCFD == 0)
+        particleCloud.postFlow();
+
+    particleCloud.clockM().stop("postFlow");

applications/solvers/rcfdemSolverHeattransfer/createFields.H

@@ -0,0 +1,206 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure field\n" << endl;
+    volScalarField pRec
+    (
+        IOobject
+        (
+            "pRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+
+    volScalarField rhoRec
+    (
+        IOobject
+        (
+            "rhoRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh//,
+        //dimensionedScalar("rhoRec", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // heat transfer fields
+    Info << "\nCreating heat transfer fields.\n" << endl;
+
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
+    );
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info << "\nCreating fields subject to calculation\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfractionRec
+    );
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        UsRec
+    );
+
+    // write fields for t=t_start
+    voidfraction.write();
+    Us.write();
+//===============================
+
+
+    Info << "Calculating face flux field phiRec\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(URec*voidfractionRec*rhoRec) & mesh.Sf()
+    );
+    phiRec.write();
+
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+
+    const IOdictionary& transportProps = mesh.lookupObject<IOdictionary>("transportProperties");
+    dimensionedScalar Cv(transportProps.lookup("Cv"));
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+            "addSource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info << "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(URec));

applications/solvers/rcfdemSolverHeattransfer/rcfdemSolverHeattransfer.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rcfdemSolverHeattransfer
+
+Description
+    Solves heat transfer between fluid and particles based on rCFD
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloudEnergy.H"
+#include "clockModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    // control coupling behavior in case of substepping
+    // assumes constant timestep size
+    label counter = 0;
+    label couplingSubStep = recurrenceBase.couplingSubStep();//5;//3;
+    double dtProp =  particleCloud.dataExchangeM().couplingTime() / runTime.deltaTValue();
+    label dtDEM2dtCFD = int(dtProp + 0.5);
+    Info << "deltaT_DEM / deltaT_CFD = " << dtDEM2dtCFD << endl;
+    if (dtDEM2dtCFD > 1)
+        Info << "coupling at substep " << couplingSubStep << endl;
+
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,URec);
+
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+        #include "TEqImp.H"
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(32,"ReadFields");
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+        }
+        particleCloud.clockM().stop("ReadFields");
+
+        particleCloud.clockM().start(33,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverHeattransfer/updateFields.H

@@ -0,0 +1,14 @@
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolScalarField("rho",rhoRec);
+recurrenceBase.recM().exportVolScalarField("p",pRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportVolVectorField("U",URec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);
+
+{
+    volScalarField& NuField(const_cast<volScalarField&>(mesh.lookupObject<volScalarField> ("NuField")));
+    recurrenceBase.recM().exportVolScalarField("NuField",NuField);
+}
+
+
+#include "updateRho.H"

applications/solvers/rcfdemSolverHeattransfer/updateRho.H

@@ -0,0 +1,32 @@
+// work-around for transient properties
+// needs to be specified for each case
+
+// case 1
+
+forAll(rhoRec,cellI)
+{
+    if (mesh.C()[cellI].z() < 0.00228)
+        rhoRec[cellI] = 1.18+(1.085-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00456)
+        rhoRec[cellI] = 1.18+(1.01-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00684)
+        rhoRec[cellI] = 1.18+(0.98-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
+    else
+        rhoRec[cellI] = 1.18+(0.955-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
+}
+
+
+// case 2
+/*
+forAll(rhoRec,cellI)
+{
+    if (mesh.C()[cellI].z() < 0.00228)
+        rhoRec[cellI] = 1.18+(1.115-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00456)
+        rhoRec[cellI] = 1.18+(1.04-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00684)
+        rhoRec[cellI] = 1.18+(1.005-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
+    else
+        rhoRec[cellI] = 1.18+(0.96-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
+}
+*/

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -6,9 +6,6 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();
-    thCond=particleCloud.thermCondM().thermCond();
-
     addSource =
         (
             he.name() == "e"
@@ -16,7 +13,7 @@
             fvc::div(phi, K) +
             fvc::div
             (
-                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U),
                 p,
                 "div(phiv,p)"
             )
@@ -25,9 +22,6 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
-    // correct source for the thermodynamic reference temperature
-    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (
@@ -35,12 +29,12 @@
       + addSource
       - Qsource
       - fvm::Sp(QCoeff/Cpv, he)
-      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+ //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverRhoSteadyPimple.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverRhoSteadyPimple

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options

@@ -15,9 +15,11 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H

@@ -4,7 +4,7 @@ particleCloud.otherForces(fOther);
 tmp<fvVectorMatrix> tUEqn
 (
     fvm::div(phi, U)
-  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + particleCloud.divVoidfractionTau(U, voidfractionRec)
   + fvm::Sp(Ksl,U)
   - fOther
  ==
@@ -18,13 +18,16 @@ fvOptions.constrain(UEqn);
 
 if (modelType=="B" || modelType=="Bfull")
 {
-    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
 }
 else
 {
-    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
 }
 
+//U.relax();
+#include "limitU.H"
+
 fvOptions.correct(U);
 
 K = 0.5*magSqr(U);

applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H

@@ -51,6 +51,19 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfraction
+    );
+
     volScalarField addSource
     (
         IOobject
@@ -58,10 +71,11 @@ Info<< "Reading thermophysical properties\n" << endl;
            "addSource",
            runTime.timeName(),
            mesh,
-           IOobject::MUST_READ,
+           IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
@@ -102,11 +116,12 @@ Info<< "Reading thermophysical properties\n" << endl;
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;
@@ -158,7 +173,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMax",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             GREAT
         )
@@ -169,12 +184,45 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMin",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             0
         )
     );
 
+    dimensionedScalar pMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMax",
+            pimple.dict(),
+            dimPressure,
+            GREAT
+        )
+    );
+
+    dimensionedScalar pMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMin",
+            pimple.dict(),
+            dimPressure,
+            -GREAT
+        )
+    );
+
+    dimensionedScalar UMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "UMax",
+            pimple.dict(),
+            dimVelocity,
+            -1.0
+        )
+    );
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -189,7 +237,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
 
     mesh.setFluxRequired(p.name());
 
@@ -217,11 +265,11 @@ Info<< "Reading thermophysical properties\n" << endl;
             "Ksl",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh
-        //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
     );
 
 
@@ -239,4 +287,20 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Us
+    );
+
+    
+    dimensionedScalar kf("0", dimensionSet(1, 1, -3, -1, 0, 0, 0), 0.026);
+
 //===============================

applications/solvers/rcfdemSolverRhoSteadyPimple/limitP.H

@@ -0,0 +1,2 @@
+p = max(p, pMin);
+p = min(p, pMax);

applications/solvers/rcfdemSolverRhoSteadyPimple/limitU.H

@@ -0,0 +1,11 @@
+if (UMax.value() > 0)
+{
+    forAll(U,cellI)
+    {
+        scalar mU(mag(U[cellI]));
+        if (mU > UMax.value())
+        {
+            U[cellI] *= UMax.value() / mU;
+        }
+    }
+}

applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H

@@ -7,14 +7,15 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
 if (modelType=="A")
 {
-    rhorAUf *= fvc::interpolate(voidfraction);
+    rhorAUf *= fvc::interpolate(voidfractionRec);
 }
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+//tUEqn.clear();
 
-surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf());
 
 
-if (simple.transonic())
+if (pimple.transonic())
 {
 //  transonic version not implemented yet
 }
@@ -34,7 +35,7 @@ else
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
-    while (simple.correctNonOrthogonal())
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -49,7 +50,7 @@ else
 
         pEqn.solve();
 
-        if (simple.finalNonOrthogonalIter())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
@@ -59,6 +60,8 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
+#include "limitP.H"
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
 rho = max(rho, rhoMin);
@@ -69,13 +72,15 @@ Info<< "rho max/min : " << max(rho).value()
 
 if (modelType=="A")
 {
-    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec);
 }
 else
 {
-    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
 }
 
+#include "limitU.H"
+
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);

applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C

@@ -17,23 +17,31 @@ License
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
 Application
-    cfdemSolverRhoSimple
+    rcfdemSolverRhoSteadyPimple
 
 Description
-    Steady-state solver for turbulent flow of compressible fluids based on
-    rhoSimpleFoam where functionality for CFD-DEM coupling has been added.
-
+    Transient (DEM) + steady-state (CFD) solver for compressible flow using the
+    flexible PIMPLE (PISO-SIMPLE) algorithm. Particle-motion is obtained from
+    a recurrence process.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
-#include "simpleControl.H"
+#include "pimpleControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
+
 #include "cfdemCloudEnergy.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
@@ -60,18 +68,37 @@ int main(int argc, char *argv[])
     #include "createFvOptions.H"
 
     // create cfdemCloud
-    #include "readGravitationalAcceleration.H"
-    cfdemCloudEnergy particleCloud(mesh);
+    //#include "readGravitationalAcceleration.H"
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
     #include "checkModelType.H"
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
 
     turbulence->validate();
+    //#include "compressibleCourantNo.H"
+    //#include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    const IOdictionary& couplingProps = particleCloud.couplingProperties();
+    label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
+    Info << "Solving flow equations every " << nEveryFlow << " steps.\n" << endl;
+    label stepcounter = 0;
+
     Info<< "\nStarting time loop\n" << endl;
 
-    while (simple.loop(runTime))
+    while (runTime.run())
     {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
         particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -80,6 +107,9 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
+        //voidfraction = voidfractionRec;
+        //Us = UsRec;
+
         if(hasEvolved)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
@@ -100,25 +130,52 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
+        volScalarField rhoeps("rhoeps",rho*voidfractionRec);
+        if (stepcounter%nEveryFlow==0)
+        {
+            while (pimple.loop())
+            {
+                // if needed, perform drag update here
+                if (pimple.nCorrPIMPLE() <= 1)
+                {
+                    #include "rhoEqn.H"
+                }
 
-        volScalarField rhoeps("rhoeps",rho*voidfraction);
-        // Pressure-velocity SIMPLE corrector
+                // --- Pressure-velocity PIMPLE corrector loop
 
                 #include "UEqn.H"
-
-
-        // besides this pEqn, OF offers a "simple consistent"-option
-        #include "pEqn.H"
-        rhoeps=rho*voidfraction;
-
                 #include "EEqn.H"
 
-        turbulence->correct();
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    // besides this pEqn, OF offers a "pimple consistent"-option
+                    #include "pEqn.H"
+                    rhoeps=rho*voidfractionRec;
+                }
 
-        particleCloud.clockM().start(32,"postFlow");
-        if(hasEvolved) particleCloud.postFlow();
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+        }
+        stepcounter++;
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
         particleCloud.clockM().stop("postFlow");
 
+        particleCloud.clockM().start(32,"ReadFields");
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+        }
+        particleCloud.clockM().stop("ReadFields");
+
         runTime.write();
 
 
@@ -126,7 +183,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-        particleCloud.clockM().stop("Flow");
+
         particleCloud.clockM().stop("Global");
     }
 

applications/solvers/rcfdemSolverRhoSteadyPimple/updateFields.H

@@ -0,0 +1,8 @@
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);

applications/solvers/recSolverTurbTransport/Make/files

@@ -0,0 +1,3 @@
+recSolverTurbTransport.C
+
+EXE=$(CFDEM_APP_DIR)/recSolverTurbTransport

applications/solvers/recSolverTurbTransport/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/recSolverTurbTransport/TEq.H

@@ -0,0 +1,17 @@
+
+volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + dU2/Sct);
+
+TEqn =
+(
+    fvm::ddt(T)
+    + fvm::div(phiRec, T)
+    - fvm::laplacian(alphaEff, T)
+    ==
+    fvOptions(T)
+);
+
+TEqn.relax(relaxCoeff);
+
+fvOptions.constrain(TEqn);
+
+TEqn.solve();

applications/solvers/recSolverTurbTransport/createFields.H

@@ -0,0 +1,174 @@
+    // dummy fields
+    Info<< "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+    
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+    
+    // recurrence fields
+    Info<< "\nCreating recurrence fields.\n" << endl;
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField U2Rec
+    (
+        IOobject
+        (
+            "U2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info<< "\nCreating fields subject to calculation\n" << endl;
+
+    volScalarField delta
+    (
+        IOobject
+        (
+            "delta",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+       mesh,
+       dimensionedScalar("delta", dimLength, 0.0)
+    );
+
+    delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0);
+    delta.write();
+
+
+    Info<< "\ncreating dU2\n" << endl;
+
+    volScalarField dU2
+    (
+        IOobject
+        (
+            "dU2",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sqrt(0.5*mag(U2Rec - magSqr(URec)))*delta*0.094
+     );
+
+    forAll(dU2, cellI)
+    {
+        if (U2Rec[cellI]-magSqr(URec[cellI]) < 0.0)
+        {
+            dU2[cellI] = 0.0;
+        }
+    }
+
+    dU2.write();
+
+Info<< "Calculating face flux field phiRec\n" << endl;
+surfaceScalarField phiRec
+ (
+     IOobject
+     (
+        "phiRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(URec) & mesh.Sf()
+ );
+
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+  
+    dimensionedScalar Sc("Sc", dimless, laminarTransport);
+    dimensionedScalar Sct("Sct", dimless, laminarTransport);
+  
+ // create concentration field
+ Info<< "Creating scalar transport field\n" << endl;
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
+
+    scalar relaxCoeff(0.0);
+
+ Info<< "reading clockProperties\n" << endl;
+
+    IOdictionary clockProperties
+    (
+            IOobject
+            (
+                "clockProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+    );
+
+    autoPtr<clockModel> myClock
+    (
+            clockModel::New
+            (
+                clockProperties,
+                mesh.time()
+            )
+    );
+

applications/solvers/recSolverTurbTransport/readFields.H

@@ -0,0 +1,14 @@
+
+recurrenceBase.recM().exportVolScalarField("U2Mean",U2Rec);
+recurrenceBase.recM().exportVolVectorField("UMean",URec);
+phiRec=linearInterpolate(URec) & mesh.Sf();
+
+dU2=sqrt(0.5*mag(U2Rec - magSqr(URec)))*delta*0.094;
+
+forAll(dU2, cellI)
+{
+    if (U2Rec[cellI]-magSqr(URec[cellI]) < 0.0)
+    {
+        dU2[cellI] = 0.0;
+    }
+}

applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Turbulent Transport Solver Recurrence
+
+Description
+    Solves a transport equation for a passive scalar on a single-phase solution
+    for a solver based on recurrence statistics
+
+Rules
+	Solution data to compute the recurrence statistics from, needs to 
+		reside in $CASE_ROOT/dataBase
+	Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "clockModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+  
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+    
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+    
+    while (runTime.run())
+    {
+
+        myClock().start(1,"Global");
+
+        runTime++;
+        
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        myClock().start(2,"fieldUpdate");
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	    Info << "Updating fields at run time " << runTime.timeOutputValue()
+	        << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        myClock().stop("fieldUpdate");
+
+        myClock().start(3,"speciesEqn");
+        #include "TEq.H"
+        myClock().stop("speciesEqn");
+
+        runTime.write();
+        
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+        
+        myClock().stop("Global");
+
+    }
+    
+    myClock().evalPar();
+    myClock().normHist();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rtfmSolverSpecies/Make/files

@@ -0,0 +1,3 @@
+rtfmSolverSpecies.C
+
+EXE=$(CFDEM_APP_DIR)/rtfmSolverSpecies

applications/solvers/rtfmSolverSpecies/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rtfmSolverSpecies/TEq.H

@@ -0,0 +1,14 @@
+    TEqn =
+    (
+        fvm::ddt(alpha2Rec, T)
+        + fvm::div(phi2Rec, T)
+        - fvm::laplacian(alpha2Rec*turbulence->nu(), T)
+        ==
+        fvOptions(alpha2Rec, T) // no fvOptions support yet
+    );
+
+    TEqn.relax(relaxCoeff);
+
+    fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.solve();

applications/solvers/rtfmSolverSpecies/createFields.H

@@ -0,0 +1,135 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+    volVectorField U1Rec
+    (
+        IOobject
+        (
+            "U1Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField alpha1Rec
+    (
+        IOobject
+        (
+            "alpha1Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField U2Rec
+    (
+        IOobject
+        (
+            "U2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info << "\nCreating fields subject to calculation\n" << endl;
+    volScalarField alpha2Rec
+    (
+        IOobject
+        (
+            "alpha2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        1-alpha1Rec
+    );
+
+
+    // write fields for t=t_start
+    alpha2Rec.write();
+//===============================
+
+
+    Info << "Calculating face flux field phi\n" << endl;
+    surfaceScalarField phi2Rec
+    (
+        IOobject
+        (
+            "phi2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(U2Rec*alpha2Rec) & mesh.Sf()
+    );
+    phi2Rec.write();
+
+    singlePhaseTransportModel laminarTransport(U2Rec, phi2Rec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U2Rec, phi2Rec, laminarTransport)
+    );
+
+    // transport stuff
+
+    // create concentration field
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
+
+    T.write();
+
+    scalar relaxCoeff(0.0);

applications/solvers/rtfmSolverSpecies/readFields.H

@@ -0,0 +1,6 @@
+recurrenceBase.recM().exportVolScalarField("alpha.air",alpha1Rec);
+alpha2Rec=1-alpha1Rec;
+recurrenceBase.recM().exportVolVectorField("U.air",U1Rec);
+recurrenceBase.recM().exportVolVectorField("U.water",U2Rec);
+recurrenceBase.recM().exportSurfaceScalarField("phi.water",phi2Rec);
+

applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverRecurrence
+
+Description
+    Solves a transport equation for a passive scalar on a two-phase solution
+    Test-bed for a solver based on recurrence statistics
+
+Rules
+    Solution data to compute the recurrence statistics from, needs to
+        reside in $CASE_ROOT/dataBase
+    Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
+    cfdemCloudRec<cfdemCloud> particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Flow");
+        #include "TEq.H"
+        particleCloud.clockM().stop("Flow");
+
+
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+            Info << "Updating fields at run time " << runTime.timeOutputValue()
+                 << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        particleCloud.clockM().start(27,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/Allwclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso recurrenceTurbulence
+
+wclean
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/testTwoFluidRecurrenceTurbulence/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wmake libso recurrenceTurbulence
+
+wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/testTwoFluidRecurrenceTurbulence/CEqn.H

@@ -0,0 +1,29 @@
+// build equation system
+/*
+    Note the use of the effective viscosity, which is provided by the turbulence model
+    The recurrence-based turbulence models are derived from the standard base classes 
+    of OpenFOAM, thus they behave as a normal turbulence model would.
+*/
+
+alphaRhoPhiCarrier = linearInterpolate(alpha2*rhoCarrier)*phi2;
+
+fvScalarMatrix CEqn
+(
+	fvm::ddt(alphaCarrier*rhoCarrier, C)
+  + fvm::div(alphaRhoPhiCarrier, C, "div(alphaRhoPhi,C)")
+  - fvm::Sp(fvc::div(alphaRhoPhiCarrier), C)
+
+  - fvm::laplacian
+	(
+		fvc::interpolate(alpha2)
+	   *fvc::interpolate(carrierPhase.turbulence().muEff()/Sc),
+		C
+	)
+	==
+	fvm::SuSp(alphaCarrier*(1.0 - alphaCarrier)*rhoCarrier*K, C)
+	+ fvOptions(alphaCarrier*rhoCarrier, C)
+);
+
+// solve equations
+fvOptions.constrain(CEqn);
+CEqn.solve();

applications/solvers/testTwoFluidRecurrenceTurbulence/Make/files

@@ -0,0 +1,3 @@
+testTwoFluidRecurrenceTurbulence.C
+
+EXE = $(FOAM_USER_APPBIN)/testTwoFluidRecurrenceTurbulence

applications/solvers/testTwoFluidRecurrenceTurbulence/Make/options

@@ -0,0 +1,31 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/phaseSystems/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialModels/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialCompositionModels/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/twoPhaseCompressibleTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I../../../src/recurrence/lnInclude \
+    -IrecurrenceTurbulence/lnInclude
+
+EXE_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingTwoPhaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -ltwoPhaseReactingTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -L$(FOAM_USER_LIBBIN) \
+    -lrecurrence \
+    -lrecurrenceTwoPhaseTurbulenceModels

applications/solvers/testTwoFluidRecurrenceTurbulence/checkTurbulenceModels.H

@@ -0,0 +1,70 @@
+//===============================
+// recurrence turbulence
+//===============================
+
+
+    // check both phases for turbulence models
+    forAllIter(PtrListDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        
+        Info << "Checking phase " << phase.name() << "'s turbulence model: "
+            << phase.turbulence().type() << endl;
+        
+        /*
+            Check for laminar turbulence. This works with OpenFOAM-4.0 and OpenFOAM-5.0,
+            as the laminar, multi-phase turbulence model is named "laminar" in OF-4.0
+            and "Stokes" in OF-5.0
+        */
+        if (phase.turbulence().type() == "laminar" || phase.turbulence().type() == "Stokes")
+        {
+            // do nothing
+        }
+        else if (isA<Foam::recurrenceTurbulenceModel>(phase.turbulence()))
+        {
+            /*
+                create a reference of the type recurrenceTurbulenceModel
+                register the recurrence model with the recurrenceTurbulenceModel
+            */
+            
+            // get const-reference to the turbulence model
+            const phaseCompressibleTurbulenceModel& turbConstRef = phase.turbulence();
+            
+            // cast away const-ness, the underlying turbulence model is not a const object, so this is bad but fine
+            phaseCompressibleTurbulenceModel& turbRef = const_cast<phaseCompressibleTurbulenceModel&>(turbConstRef);
+            
+            // cast away the wrapper class, to get a reference to the turbulence models' base class
+            PhaseCompressibleTurbulenceModel<phaseModel>& baseTurbRef
+            (
+                static_cast<PhaseCompressibleTurbulenceModel<phaseModel>&>(turbRef)
+            );
+            
+            // casting down the family tree
+            Foam::recurrenceTurbulenceModel& recTurbRef
+            (
+                dynamic_cast<Foam::recurrenceTurbulenceModel&>(baseTurbRef)
+            );
+            
+            
+            // set recurrenceBase pointer
+            recTurbRef.setRecurrenceBasePtr(&recurrenceBase);
+            
+            
+            // check model settings
+            turbRef.validate();
+        }
+        else
+        {
+            /*
+                In a recurrence run, we do not compute any turbulence as we do not solve the fluid flow
+                At this point, the phase is not laminar (i.e. not using turbulence) or 
+                    using recurrenceTurbulence (i.e. taking turbulent quantities from the data base).
+                Hence, abort!
+            */
+            FatalError
+                << "Wrong turbulence model type "
+                << phase.turbulence().type() << " for phase " << phase.name() << nl << nl
+                << "Valid turbulence model types are types derived from recurrenceTurbulenceModel or laminar" << endl
+                << exit(FatalError);
+        }
+    }

applications/solvers/testTwoFluidRecurrenceTurbulence/createTransportFields.H

@@ -0,0 +1,79 @@
+/* --------------------------------------------------------------------------------- */
+/* 								read flotation properties							 */
+/* --------------------------------------------------------------------------------- */
+
+Info<< "Reading scalarTransportProperties\n" << endl;
+IOdictionary scalarTransportProperties
+(
+    IOobject
+    (
+        "scalarTransportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+
+const scalar Sc(scalarTransportProperties.lookupOrDefault<scalar>("Sc",scalar(1.0)));
+
+
+const word carrierPhaseName(scalarTransportProperties.lookup("carrierPhase"));
+
+if (carrierPhaseName != phase1.name() && carrierPhaseName != phase2.name())
+{
+	FatalError << "No valid carrier phase specified" << nl
+		<< "Valid phase names are: " << nl
+		<< phase1.name() << ",  " << phase2.name()
+		<< abort(FatalError);
+}
+
+phaseModel& carrierPhase = (carrierPhaseName == phase1.name()) ? phase1 : phase2;
+
+const word dispersePhaseName = (carrierPhaseName == phase1.name()) ? phase2.name() : phase1.name();
+
+
+volScalarField& rhoCarrier = carrierPhase.thermo().rho();
+volScalarField& alphaCarrier = carrierPhase;
+surfaceScalarField& alphaRhoPhiCarrier = carrierPhase.alphaRhoPhi();
+
+
+volScalarField contErrCarrier
+(
+    "contErrCarrier",
+    fvc::ddt(alphaCarrier, rhoCarrier)
+);
+
+
+
+Info<< "Reading field C\n" << endl;
+volScalarField C
+(
+    IOobject
+    (
+        "C",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField K
+(
+    IOobject
+    (
+        "K",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+
+
+

applications/solvers/testTwoFluidRecurrenceTurbulence/readFields.H

@@ -0,0 +1,9 @@
+// update flow fields
+recurrenceBase.recM().exportVolScalarField("alpha."+carrierPhaseName,alpha2);
+recurrenceBase.recM().exportVolScalarField("alpha."+dispersePhaseName,alpha1);
+
+recurrenceBase.recM().exportVolVectorField("U."+carrierPhaseName,U2);
+
+// update turbulence models
+phase1.correctTurbulence();
+phase2.correctTurbulence();

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Allwclean

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+
+wclean libso
+
+#------------------------------------------------------------------------------

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+
+wmake $targetType
+
+#------------------------------------------------------------------------------

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Make/files

@@ -0,0 +1,10 @@
+recurrenceTurbulenceModel/recurrenceTurbulenceModel.C
+
+
+recurrenceTurbulenceModels.C
+
+recurrenceKEpsilon/recurrenceKEpsilon.C
+recurrenceKOmega/recurrenceKOmega.C
+recurrenceSmagorinsky/recurrenceSmagorinsky.C
+
+LIB = $(FOAM_USER_LIBBIN)/librecurrenceTwoPhaseTurbulenceModels

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Make/options

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/phaseSystems/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialModels/lnInclude\
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialCompositionModels/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I../recurrenceTurbulenceModel/lnInclude \
+    -I../../../../src/recurrence/lnInclude
+
+LIB_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingTwoPhaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -L$(FOAM_USER_LIBBIN) \
+    -lrecurrence \
+    -lreactingTwoPhaseSystem

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/ReadMe.md

@@ -0,0 +1,107 @@
+# Recurrence-based, multi-phase turbulence modelling
+
+This model implements recurrence-based turbulence models, i.e. the fundamental 
+turbulent field quantities are read from the data base and are not solved for.
+All derived field quantities are computed just in the same way as the proper 
+turbulence models do. By deriving the recurrence-based turbulence models from 
+somewhere up the family tree of OpenFOAM's turbulence model class hierarchy, 
+the recurrence-based turbulence models are fully compatible with OpenFOAM's 
+generic treatment of turbulence modelling, i.e. solvers and libraries interact 
+with references to a generic base type of the actual turbulence model. Hence, 
+solvers and libraries may remain blissfully ignorant of the actual turbulence 
+model in use.
+
+For laminar phases no special treatment is necessary, as the *laminar* 
+turbulence model does not compute any fields.
+
+
+## Development notes
+
+The initial development covers only a small number of turbulence models.
+
+
+## Notes on usage
+
+The turbulence model in use for the recurrence run must be the recurrence-based 
+equivalent of the turbulence model used for generating the data base, i.e. if 
+the data base was computed using the *kEpsilon* model, then the recurrence solver 
+is to employ the *recurrenceKEpsilon* turbulence model. This model will read 
+the relevant model coefficients from the *turbulenceProperties* dictionary, and 
+make sure that the turbulent fields `k` and `epsilon` are contained in the data 
+base.
+
+Whenever, the solver or a library calls `turbulence->nut()` to access the 
+turbulent viscosity, the recurrence-based kEpsilon model will compute `nut` 
+according to kEpsilon's relations `nut = Cmu*sqr(k)/epsilon`, with the fields 
+`k` and `epsilon` being from the current snapshot provided by the recurrence model.
+
+Thus, the fundamental turbulent field quantities of the employed turbulence model 
+have to be added to the *volScalarFields* list in the `recProperties` dictionary 
+controlling the recurrence model. This will ensure that the turbulent field 
+quantities are read from the data base.
+
+
+## Notes on the implementation
+
+The base class implements the method `void setRecurrenceBasePtr(recBase*)`, which 
+is used to give the recurrence-based turbulence models a reference (technically 
+a pointer) to the recurrence model. Thus, after construction of the turbulence 
+models and the recurrence model, `setRecurrenceBasePtr()` needs to be called as 
+the pointer to the recurrence model is initialized by the constructor with `NULL`.
+Trying to access the recurrence model from within the recurrence-based turbulence 
+model prior to setting the pointer to the recurrence model with 
+`setRecurrenceBasePtr()` will result in a segmentation fault.
+In order to be able to call `setRecurrenceBasePtr()`, the generic reference to 
+the turbulence model needs to be converted into a reference of the base class' 
+type, i.e. `recurrenceTurbulenceModel`.
+This unfortunate deviation from good standards, i.e. making full use of C++'s 
+polymorphism, should be the only instance of having to use non-pretty hacks.
+However, apart from initialisation, i.e. setting the pointer to the recurrence 
+model, the recurrence-based turbulence models adhere to the generic interface of 
+OpenFOAM's turbulence models, and can be used as any other turbulence model.
+
+
+The concrete implementations, e.g. *recurrenceKEpsilon*, use the method 
+`validate()` to check whether the underlying turbulent quantities are specified 
+for use in the data base in the *volScalarFields* list in the `recProperties` 
+dictionary. This method is part of the signature of the class `Foam::turbulenceModel`, 
+which is the very base class of all turbulence models in OpenFOAM.
+In proper turbulence models, this method is used to check whether the internal 
+fields are properly initialized and to update all derived quantities.
+In the solver, `validate()` must not be called prior to `setRecurrenceBasePtr()`, 
+as validate accesses the recurrence model. The wrong order of function calls will 
+result in a segmentation fault, as the pointer to the recurrence model is 
+initialized by the constructor with `NULL`.
+
+
+The concrete implementations, e.g. *recurrenceKEpsilon*, use the method 
+`correct()` to update the turbulent field quantities from the data base, 
+and in turn update the derived quantities, such as `nut`.
+This method is part of the signature of the class `Foam::turbulenceModel`, 
+which is the very base class of all turbulence models in OpenFOAM.
+In proper turbulence models, this method is used to solve for the next time step.
+
+
+## Compilation
+
+Source OpenFOAM and simply compile with
+
+```bash
+./Allwclean
+./Allwmake
+```
+
+The script `Allwclean` will clear all previous builds. This step is not needed for
+first-time compilation. It is, however, recommended for subsequent compilations, as
+it completely clears the slate. The script `Allwmake` will run the compilation for
+the passive particle model.
+
+
+## Required software
+
+This model has been tested with the following versions of OpenFOAM:
+
+* OpenFOAM-4.0
+* OpenFOAM-5.0
+
+

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKEpsilon/recurrenceKEpsilon.C

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceKEpsilon.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+void Foam::RASModels::recurrenceKEpsilon::correctNut()
+{
+    this->nut_ = Cmu_*sqr(k_)/epsilon_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RASModels::recurrenceKEpsilon::recurrenceKEpsilon
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+
+)
+:
+    eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("k0", dimensionSet(0,2,-2,0,0), 0.0)
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("eps0", dimensionSet(0,2,-3,0,0), 0.0)
+    )
+{
+    if (type == typeName)
+    {
+        printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+Foam::RASModels::recurrenceKEpsilon::~recurrenceKEpsilon()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::RASModels::recurrenceKEpsilon::read()
+{
+    if
+    (
+        eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>::read()
+    )
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::RASModels::recurrenceKEpsilon::correct()
+{
+    // update turbulence fields
+    recurrenceBasePtr_->recM().exportVolScalarField("k."+group_, this->k_);
+    recurrenceBasePtr_->recM().exportVolScalarField("epsilon."+group_, this->epsilon_);
+    
+    // update nut
+    correctNut();
+}
+
+
+void Foam::RASModels::recurrenceKEpsilon::validate()
+{
+    /*
+        Check whether k and epsilon are included in the dataBase.
+        The check only makes sure that these fields are included in the
+            volScalarFields list of recProperties.
+        Whether the fields are actually contained in the dataBase is 
+            done by the recurrenceModel itself.
+    */
+    bool foundK(false);
+    bool foundEpsilon(false);
+    
+    wordList fieldNames(recurrenceBasePtr_->recM().volScalarFieldNames());
+    
+    forAll(fieldNames, i)
+    {
+        word curFieldName = fieldNames[i];
+        
+        if (curFieldName == k_.name())
+        {
+            foundK = true;
+        }
+        
+        if (curFieldName == epsilon_.name())
+        {
+            foundEpsilon = true;
+        }
+        
+        
+    }
+    
+    if (not (foundK and foundEpsilon))
+    {
+        FatalError
+            << "Fields " << k_.name() << " and " << epsilon_.name()
+            << " not specified in the volScalarFields list of recProperties!" << nl
+            << "volScalarFields : " << fieldNames << nl
+            << "Add these fields and make sure they are contained in the dataBase." << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKEpsilon/recurrenceKEpsilon.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceKEpsilon
+
+Description
+    recurrence-based kEpsilon turbulence model
+
+    This model provides kEpslion's turbulence quantities that were computed 
+    elsewhere, i.e. taken from the recurrence dataBase.
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceKEpsilon.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceKEpsilon_H
+#define recurrenceKEpsilon_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceKEpsilon Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceKEpsilon
+:
+    public eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>,
+    public recurrenceTurbulenceModel
+{
+    // Private data
+
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceKEpsilon(const recurrenceKEpsilon&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceKEpsilon&);
+        
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField epsilon_;
+    
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceKEpsilon");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceKEpsilon
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+
+    //- Destructor
+    virtual ~recurrenceKEpsilon();
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+        
+        
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
+        //- Update the turbulent fields
+        virtual void correct();
+        
+        //- Check model settings
+        virtual void validate();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKOmega/recurrenceKOmega.C

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceKOmega.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+void Foam::RASModels::recurrenceKOmega::correctNut()
+{
+    this->nut_ = k_/omega_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RASModels::recurrenceKOmega::recurrenceKOmega
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+
+)
+:
+    eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("k0", dimensionSet(0,2,-2,0,0), 0.0)
+    ),
+    omega_
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("om0", dimensionSet(0,0,-1,0,0), 0.0)
+    )
+{
+    if (type == typeName)
+    {
+        printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+Foam::RASModels::recurrenceKOmega::~recurrenceKOmega()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::RASModels::recurrenceKOmega::read()
+{
+    if
+    (
+        eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>::read()
+    )
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::RASModels::recurrenceKOmega::correct()
+{
+    // update turbulence fields
+    recurrenceBasePtr_->recM().exportVolScalarField("k."+group_, this->k_);
+    recurrenceBasePtr_->recM().exportVolScalarField("omega."+group_, this->omega_);
+    
+    // update nut
+    correctNut();
+}
+
+
+void Foam::RASModels::recurrenceKOmega::validate()
+{
+    /*
+        Check whether k and omega are included in the dataBase.
+        The check only makes sure that these fields are included in the
+            volScalarFields list of recProperties.
+        Whether the fields are actually contained in the dataBase is 
+            done by the recurrenceModel itself.
+    */
+    bool foundK(false);
+    bool foundOmega(false);
+    
+    wordList fieldNames(recurrenceBasePtr_->recM().volScalarFieldNames());
+    
+    forAll(fieldNames, i)
+    {
+        word curFieldName = fieldNames[i];
+        
+        if (curFieldName == k_.name())
+        {
+            Info << "Found " << k_.name()<< endl;
+            foundK = true;
+        }
+        
+        if (curFieldName == omega_.name())
+        {
+            Info << "Found " << omega_.name()<< endl;
+            foundOmega = true;
+        }
+    }
+    
+    if (not (foundK and foundOmega))
+    {
+        FatalError
+            << "Fields " << k_.name() << " and " << omega_.name()
+            << " not specified in the volScalarFields list of recProperties!" << nl
+            << "volScalarFields : " << fieldNames << nl
+            << "Add these fields and make sure they are contained in the dataBase." << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKOmega/recurrenceKOmega.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceKOmega
+
+Description
+    recurrence-based kOmega turbulence model
+
+    This model provides kOmega's turbulence quantities that were computed 
+    elsewhere, i.e. taken from the recurrence dataBase
+    To be used by recurrence solvers
+
+SourceFiles
+    recurrenceKOmega.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceKOmega_H
+#define recurrenceKOmega_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceKOmega Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceKOmega
+:
+    public eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>,
+    public recurrenceTurbulenceModel
+{
+    // Private data
+
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceKOmega(const recurrenceKOmega&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceKOmega&);
+        
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField omega_;
+    
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceKOmega");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceKOmega
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+
+    //- Destructor
+    virtual ~recurrenceKOmega();
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+        
+        
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence specific dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "epsilon",
+                        this->mesh_.time().timeName(),
+                        this->mesh_
+                    ),
+                    Cmu_*k_*omega_,
+                    omega_.boundaryField().types()
+                )
+            );
+        }
+
+        //- Update the turbulent fields
+        virtual void correct();
+        
+        //- Check model settings
+        virtual void validate();
+        
+        
+    // Setters
+        
+        void setRecurrenceBasePtr(recBase* recurrenceBasePtr)
+        {
+            recurrenceBasePtr_ = recurrenceBasePtr;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceSmagorinsky/recurrenceSmagorinsky.C

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceSmagorinsky.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::LESModels::recurrenceSmagorinsky::k
+(
+    const tmp<volTensorField>& gradU
+) const
+{
+    volSymmTensorField D(symm(gradU));
+
+    volScalarField a(this->Ce_/this->delta());
+    volScalarField b((2.0/3.0)*tr(D));
+    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("k", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_
+            ),
+            sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
+        )
+    );
+}
+
+
+void Foam::LESModels::recurrenceSmagorinsky::correctNut()
+{
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    this->nut_ = Ck_*this->delta()*sqrt(k);
+    this->nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LESModels::recurrenceSmagorinsky::recurrenceSmagorinsky
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+)
+:
+    LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    )
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::LESModels::recurrenceSmagorinsky::read()
+{
+    if (LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::LESModels::recurrenceSmagorinsky::epsilon() const
+{
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_
+            ),
+            this->Ce_*k*sqrt(k)/this->delta()
+        )
+    );
+}
+
+
+void Foam::LESModels::recurrenceSmagorinsky::correct()
+{
+    // update nut
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceSmagorinsky/recurrenceSmagorinsky.H

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::recurrenceSmagorinsky
+
+Description
+    The recurrenceSmagorinsky SGS model.
+
+    This model provides Smagorinsky's turbulence quantities.
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceSmagorinsky.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceSmagorinsky_H
+#define recurrenceSmagorinsky_H
+
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceSmagorinsky Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceSmagorinsky
+:
+    public LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>,
+    public recurrenceTurbulenceModel
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        recurrenceSmagorinsky(const recurrenceSmagorinsky&);
+        void operator=(const recurrenceSmagorinsky&);
+
+
+protected:
+
+    // Protected data
+
+        dimensionedScalar Ck_;
+
+
+    // Protected Member Functions
+
+        //- Return SGS kinetic energy
+        //  calculated from the given velocity gradient
+        tmp<volScalarField> k(const tmp<volTensorField>& gradU) const;
+
+        //- Update the SGS eddy viscosity
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceSmagorinsky");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceSmagorinsky
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~recurrenceSmagorinsky()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k(fvc::grad(this->U_));
+        }
+
+        //- Return sub-grid disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Correct Eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModel/recurrenceTurbulenceModel.C

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceTurbulenceModel.H"
+
+
+namespace Foam
+{
+    defineTypeNameAndDebug(recurrenceTurbulenceModel, 0);
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::recurrenceTurbulenceModel::recurrenceTurbulenceModel
+(
+    const word group
+)
+:
+    recurrenceBasePtr_(NULL),
+    group_(group)
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::recurrenceTurbulenceModel::~recurrenceTurbulenceModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModel/recurrenceTurbulenceModel.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceTurbulenceModel
+
+Description
+    recurrence-based turbulence model base class
+
+    This turbulence model class provides a framework for derived turbulence
+    models to provide turbulence quantities that were computed elsewhere, 
+    i.e. taken from the recurrence dataBase.
+    
+    A concrete recurrence-based turbulence model implementation, e.g. recurrenceKEpsilon, 
+    needs to be derived from the base class of the original turbulence model 
+    and this class. This class provides the link to the recurrence model. The original 
+    base class provides the proper turbulence modelling interfaces of OpenFOAM's 
+    turbulence modelling framework.
+    
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceTurbulenceModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceTurbulenceModel_H
+#define recurrenceTurbulenceModel_H
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceTurbulenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceTurbulenceModel
+{
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceTurbulenceModel(const recurrenceTurbulenceModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceTurbulenceModel&);
+        
+
+protected:
+
+    // Protected data
+
+        recBase* recurrenceBasePtr_;
+        const word group_;
+        
+    
+    // Protected Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("recurrenceTurbulenceModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceTurbulenceModel
+        (
+            const word group
+        );
+
+
+    //- Destructor
+    virtual ~recurrenceTurbulenceModel();
+
+
+    // Member Functions
+        
+        void setRecurrenceBasePtr(recBase* recurrenceBasePtr)
+        {
+            recurrenceBasePtr_ = recurrenceBasePtr;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModels.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseCompressibleTurbulenceModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
+
+#include "laminar.H"
+#include "RASModel.H"
+#include "LESModel.H"
+
+#include "recurrenceKEpsilon.H"
+#include "recurrenceKOmega.H"
+#include "recurrenceSmagorinsky.H"
+
+
+
+// Instructions for OpenFOAM-5.0 
+/*makeTurbulenceModelTypes
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
+    phaseModel
+);
+
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, RAS, recurrenceKEpsilon);
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, RAS, recurrenceKOmega);
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, LES, recurrenceSmagorinsky);
+*/
+
+
+
+// Instructions for OpenFOAM-4.0 
+makeBaseTurbulenceModel
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
+    phaseModel
+);
+
+#define makeRASModel(Type)                                            \
+    makeTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
+
+#define makeLESModel(Type)                                            \
+    makeTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
+
+makeRASModel(recurrenceKEpsilon);
+
+makeRASModel(recurrenceKOmega);
+
+makeLESModel(recurrenceSmagorinsky);
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/testTwoFluidRecurrenceTurbulence.C

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    testTwoFluidRecurrenceTurbulence
+
+Description
+    A modified variant of the two-fluid, recurrence model A solver
+    with the extension of recurrence-based, multi-phase turbulence modelling.
+    This application is used to test whether turbulent fields can be provided
+    by the recurrence-based turbulence models.
+    
+    Run this test application in a recurrence case, with turbulence enabled and 
+    the necessary turbulent field quantities present in the data base.
+    
+    Note the initialisation in checkTurbulenceModels.H
+    Updating the turbulence model is done by calling phaseX.correctTurbulence()
+    in the file readFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "twoPhaseSystem.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+/* // uncomment for OpenFOAM-5.0
+namespace Foam
+{
+    tmp<volScalarField> byDt(const volScalarField& vf)
+    {
+        if (fv::localEulerDdt::enabled(vf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
+        }
+        else
+        {
+            return vf/vf.mesh().time().deltaT();
+        }
+    }
+
+    tmp<surfaceScalarField> byDt(const surfaceScalarField& sf)
+    {
+        if (fv::localEulerDdt::enabled(sf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaTf(sf.mesh())*sf;
+        }
+        else
+        {
+            return sf/sf.mesh().time().deltaT();
+        }
+    }
+}
+*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H" // remove for OpenFOAM-5.0
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    
+    #include "createTransportFields.H"
+
+    if (!LTS)
+    {
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+    
+    Switch faceMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
+    
+    recBase recurrenceBase(mesh);
+    
+    #include "checkTurbulenceModels.H"
+
+    #include "pUf/createDDtU.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting recurrence-based time loop\n" << endl;
+    
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        #include "CourantNos.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        
+        #include "CEqn.H"
+        
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	        Info << "Updating fields at run time " << runTime.timeOutputValue()
+	            << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+	        
+            recurrenceBase.updateRecFields();
+	        #include "readFields.H"
+            recTimeIndex++;
+        }
+        
+
+        runTime.write();
+        
+        #include "writeCField.H"
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/writeCField.H

@@ -0,0 +1,12 @@
+/* ----------------------------------------------
+	Write averaged particle volume fraction
+---------------------------------------------- */
+
+// essential information
+Info << "Total mass :";
+Info << tab << sum(C*rhoCarrier*alphaCarrier*mesh.V());
+Info << endl;
+
+Info << "Total Carrier mass :";
+Info << tab << sum(rhoCarrier*alphaCarrier*mesh.V());
+Info << endl;

applications/utilities/mapDataBase/mapFieldData.sh

@@ -0,0 +1,30 @@
+#!/bin/sh
+# Source run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# to be executed from one level above the source directory
+
+if [ $# -eq 0 ]
+  then
+    sourceName="dataBase"
+    targetName="dataBaseCoarse"
+  else
+    sourceName=$1
+    targetName=$2
+fi
+
+cd $sourceName
+
+for time in *
+do
+    if [ $time != "system" ] && [ $time != "constant" ];
+    then
+        cd ../$targetName
+        echo "Found $time."
+        sed -i "/^startTime/c\startTime \t$time;" ./system/controlDict
+        grep 'startTime' ./system/controlDict
+        mapFields ../$sourceName -sourceTime $time -consistent
+        cd ../$sourceName
+    fi
+done
+

applications/utilities/rBaseMirror/mirrorProperties

@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          mirrorProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//===========================================================================//
+// sub-models & settings
+
+
+refPoint (0 0 0);
+
+refDirection (1 0 0);
+
+fieldName U;
+
+dataBaseName dataBase;
+
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/files

@@ -0,0 +1,3 @@
+rBaseMirrorScalar.C
+
+EXE=$(CFDEM_APP_DIR)/rBaseMirrorScalar

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/options

@@ -0,0 +1,16 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude 
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorScalar/createFields.H

@@ -0,0 +1,17 @@
+    IOdictionary mirrorProperties
+    (
+        IOobject
+        (
+            "mirrorProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    vector refPoint(mirrorProperties.lookup("refPoint"));
+    vector refDirection(mirrorProperties.lookup("refDirection"));
+
+    word fieldName(mirrorProperties.lookup("fieldName"));
+    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorScalar/rBaseMirrorScalar.C

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rBaseMirror
+
+
+Description
+    Read time series and extend it by mirrored fields if geometry possesses
+    the same symmetry
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+   // read in start and end time from controlDict
+
+    scalar startTime=runTime.startTime().value();
+    scalar endTime=runTime.endTime().value();
+    scalar origTimeRange = endTime - startTime;
+
+    Info << "start time = " << runTime.startTime() << endl;
+    Info << "end time = " << runTime.endTime() << endl;
+
+    // check which time directories are present
+    // instantList timeDirs = timeSelector::select0(runTime, args);
+    // runTime.setTime(timeDirs[0], 0);
+
+    #include "createMesh.H"
+
+    #include "createFields.H"
+
+    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
+    instantList timeDirs(recTime.times());
+    recTime.setTime(timeDirs[0],0);
+
+    #include "readFields.H"
+
+    Info << fieldName << endl;
+
+    volScalarField transformedField = origField;
+
+    scalar t;
+
+    label shiftedTimeI = 0;
+
+    // check number of time directories
+    label shift = 0;
+    forAll(timeDirs, timeI)
+    {
+
+       if (recTime.timeName() == "constant") continue;
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       shift++;
+    }
+
+    scalar dt = origTimeRange / (shift - 1.0);
+    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
+
+    label cellI_transformed = -1;
+    forAll(timeDirs, timeI)
+    {
+
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       Info << "time = " << t << ", time index = " << timeI << endl;
+
+       #include "readFields.H"
+
+       forAll(transformedField, cellI)
+       {
+           vector position = mesh.C()[cellI];
+           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+           cellI_transformed = mesh.findCell(transformedPosition);
+           if(cellI_transformed < 0)
+           {
+               Info << "Couldn't find transformed cell. Stopping." << endl;
+               return 0;
+           }
+
+           scalar value = origField[cellI_transformed];
+           scalar transformedValue = value;
+
+           transformedField[cellI] = transformedValue;
+       }
+
+       shiftedTimeI = timeI + shift;
+       t = recTime.value() + origTimeRange + dt;
+       runTime.setTime(t, shiftedTimeI);
+       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
+       transformedField.write();
+    }
+
+    Info << "\nEnd" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorScalar/readFields.H

@@ -0,0 +1,13 @@
+
+    volScalarField origField
+    (
+        IOobject
+        (
+            fieldName,
+            recTime.timePath(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    );

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/files

@@ -0,0 +1,3 @@
+rBaseMirrorVec.C
+
+EXE=$(CFDEM_APP_DIR)/rBaseMirrorVec

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/options

@@ -0,0 +1,16 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude 
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorVec/createFields.H

@@ -0,0 +1,17 @@
+    IOdictionary mirrorProperties
+    (
+        IOobject
+        (
+            "mirrorProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    vector refPoint(mirrorProperties.lookup("refPoint"));
+    vector refDirection(mirrorProperties.lookup("refDirection"));
+
+    word fieldName(mirrorProperties.lookup("fieldName"));
+    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorVec/rBaseMirrorVec.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rBaseMirror
+
+Description
+    Read time series and extend it by mirrored fields if geometry possesses
+    the same symmetry
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+   // read in start and end time from controlDict
+
+    scalar startTime=runTime.startTime().value();
+    scalar endTime=runTime.endTime().value();
+    scalar origTimeRange = endTime - startTime;
+
+    Info << "start time = " << runTime.startTime() << endl;
+    Info << "end time = " << runTime.endTime() << endl;
+
+    // check which time directories are present
+    //instantList timeDirs = timeSelector::select0(runTime, args);
+    //runTime.setTime(timeDirs[0], 0);
+
+    #include "createMesh.H"
+
+    #include "createFields.H"
+
+    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
+    instantList timeDirs(recTime.times());
+    recTime.setTime(timeDirs[0],0);
+
+    #include "readFields.H"
+
+    Info << fieldName << endl;
+
+    volVectorField transformedField = origField;
+
+    scalar t;
+
+    label shiftedTimeI = 0;
+
+    // check number of time directories
+    label shift = 0;
+    forAll(timeDirs, timeI)
+    {
+
+       if (recTime.timeName() == "constant") continue;
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       shift++;
+    }
+
+    scalar dt = origTimeRange / (shift - 1.0);
+    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
+
+    label cellI_transformed = -1;
+    forAll(timeDirs, timeI)
+    {
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       Info << "time = " << t << ", time index = " << timeI << endl;
+
+       #include "readFields.H"
+
+       forAll(transformedField, cellI)
+       {
+           vector position = mesh.C()[cellI];
+           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+           cellI_transformed = mesh.findCell(transformedPosition);
+           if(cellI_transformed < 0)
+           {
+               Info << "Couldn't find transformed cell. Stopping." << endl;
+               return 0;
+           }
+
+           vector value = origField[cellI_transformed];
+           vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
+
+           transformedField[cellI] = transformedValue;
+       }
+
+       shiftedTimeI = timeI + shift;
+       t = recTime.value() + origTimeRange + dt;
+       runTime.setTime(t, shiftedTimeI);
+       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
+       transformedField.write();
+    }
+
+    Info << "\nEnd" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorVec/readFields.H

@@ -0,0 +1,14 @@
+
+volVectorField origField
+    (
+        IOobject
+        (
+            fieldName,
+            recTime.timePath(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    );
+

applications/utilities/rSmoother/Make/files

@@ -0,0 +1,3 @@
+rSmoother.C
+
+EXE=$(CFDEM_APP_DIR)/rSmoother

applications/utilities/rSmoother/Make/options

@@ -0,0 +1,18 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lfiniteVolume \
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+    -lsampling \
+    -lfvOptions \
+$(CFDEM_ADD_LIBS)

applications/utilities/rSmoother/createFields.H

@@ -0,0 +1,28 @@
+    IOdictionary recProperties
+    (
+        IOobject
+        (
+            "recProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    scalar threshold(readScalar(recProperties.lookup("threshold")));
+
+    volVectorField U_smooth
+    (
+        IOobject
+        (
+            "U_smooth",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector(0.0, 0.0, 0.0))
+    );

applications/utilities/rSmoother/rSmoother.C

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rSmoother
+
+Description
+    Loops over all recurrence times and averages fields over given similarity range
+
+\*---------------------------------------------------------------------------*/
+
+ #include "fvCFD.H"
+ #include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+
+
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    // set time step to that of recurrence database
+    runTime.setDeltaT(recurrenceBase.recM().recTimeStep());
+
+    // check start and end time
+    if (abs(runTime.startTime().value() - recurrenceBase.recM().recStartTime()) > 1e-5)
+    {
+        Info << "Stopping. Start time and database start time are different." << endl;
+        Info << "Start time = " << runTime.startTime().value() << endl;
+        Info << "Database start time = " << recurrenceBase.recM().recStartTime() << endl;
+        return 0;
+    }
+
+    if (runTime.endTime().value() > recurrenceBase.recM().recEndTime())
+    {
+        runTime.setEndTime(recurrenceBase.recM().recEndTime());
+        Info << "End time set to database end time." << endl;
+    }
+
+    label index = -1;
+
+    Info << "\nSmoothing recurrence statistics\n" << endl;
+
+    while (runTime.run())
+    {
+        // runtime can't be larger than recurrence database size
+        index = runTime.timeIndex();
+        #include "updateFields.H"
+        runTime++;
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/rSmoother/updateFields.H

@@ -0,0 +1,4 @@
+Info << "averaging fields for time = " << runTime.value() << ", time index = " << index << endl;
+
+recurrenceBase.recM().exportAveragedVolVectorField(U_smooth, "UMean", threshold, index);
+U_smooth.write();

doc/_themes/lammps_theme/breadcrumbs.html

@@ -0,0 +1,34 @@
+<div role="navigation" aria-label="breadcrumbs navigation">
+  <ul class="wy-breadcrumbs">
+    <li><a href="{{ pathto(master_doc) }}">Docs</a> &raquo;</li>
+      {% for doc in parents %}
+          <li><a href="{{ doc.link|e }}">{{ doc.title }}</a> &raquo;</li>
+      {% endfor %}
+    <li>{{ title }}</li>
+      <li class="wy-breadcrumbs-aside">
+        {% if pagename != "search" %}
+          {% if display_github %}
+            <a href="https://{{ github_host|default("github.com") }}/{{ github_user }}/{{ github_repo }}/blob/{{ github_version }}{{ conf_py_path }}{{ pagename }}{{ source_suffix }}" class="fa fa-github"> Edit on GitHub</a>
+          {% elif display_bitbucket %}
+            <a href="https://bitbucket.org/{{ bitbucket_user }}/{{ bitbucket_repo }}/src/{{ bitbucket_version}}{{ conf_py_path }}{{ pagename }}{{ source_suffix }}" class="fa fa-bitbucket"> Edit on Bitbucket</a>
+          {% elif show_source and source_url_prefix %}
+            <a href="{{ source_url_prefix }}{{ pagename }}{{ source_suffix }}">View page source</a>
+          {% elif show_source and has_source and sourcename %}
+            <a href="{{ pathto('_sources/' + sourcename, true)|e }}" rel="nofollow"> View page source</a>
+          {% endif %}
+            <a href="https://www.cfdem.com">Website</a>
+        {% endif %}
+      </li>
+  </ul>
+  <hr/>
+  {% if next or prev %}
+    <div class="rst-footer-buttons" style="margin-bottom: 1em" role="navigation" aria-label="footer navigation">
+      {% if next %}
+        <a href="{{ next.link|e }}" class="btn btn-neutral float-right" title="{{ next.title|striptags|e }}" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>
+      {% endif %}
+      {% if prev %}
+        <a href="{{ prev.link|e }}" class="btn btn-neutral" title="{{ prev.title|striptags|e }}" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+      {% endif %}
+    </div>
+  {% endif %}
+</div>

doc/_themes/lammps_theme/footer.html

@@ -0,0 +1,36 @@
+<footer>
+  {% if next or prev %}
+    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
+      {% if next %}
+        <a href="{{ next.link|e }}" class="btn btn-neutral float-right" title="{{ next.title|striptags|e }}" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>
+      {% endif %}
+      {% if prev %}
+        <a href="{{ prev.link|e }}" class="btn btn-neutral" title="{{ prev.title|striptags|e }}" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>
+      {% endif %}
+    </div>
+  {% endif %}
+
+  <hr/>
+
+  <div role="contentinfo">
+    <p>
+    {%- if show_copyright %}
+      {%- if hasdoc('copyright') %}
+        {% trans path=pathto('copyright'), copyright=copyright|e %}&copy; <a href="{{ path }}">Copyright</a> {{ copyright }}.{% endtrans %}
+      {%- else %}
+        {% trans copyright=copyright|e %}&copy; Copyright {{ copyright }}.{% endtrans %}
+      {%- endif %}
+    {%- endif %}
+
+    {%- if last_updated %}
+      {% trans last_updated=last_updated|e %}Last updated on {{ last_updated }}.{% endtrans %}
+    {%- endif %}
+    </p>
+  </div>
+
+  {%- if show_sphinx %}
+  {% trans %}Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>{% endtrans %}.
+  {%- endif %}
+
+</footer>
+