diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c2608e2..d69c110 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -43,9 +43,9 @@ pages: lint: stage: static_analyze before_script: - - apk add clang-extra-tools + - apk add clang-extra-tools openmp-dev script: - mkdir lint && cd lint - - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*, modernize-*" -DTUG_ENABLE_TESTING=OFF .. + - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*" -DTUG_ENABLE_TESTING=OFF .. - make tug when: manual diff --git a/CMakeLists.txt b/CMakeLists.txt index ec9da4c..afed992 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,15 @@ endif() option(TUG_ENABLE_TESTING "Run tests after succesfull compilation" - ON) + OFF) + +option(TUG_HANNESPHILIPP_EXAMPLES + "Compile example applications" + OFF) + +option(TUG_NAAICE_EXAMPLE + "Enables NAAICE examples with export of CSV files" + OFF) add_subdirectory(src) @@ -41,9 +49,9 @@ if(TUG_ENABLE_TESTING) add_subdirectory(test) endif() -add_subdirectory(examples) - -option(TUG_NAAICE_EXAMPLE "Enables NAAICE examples with export of CSV files" OFF) +if(TUG_HANNESPHILIPP_EXAMPLES) + add_subdirectory(examples) +endif() if (TUG_NAAICE_EXAMPLE) add_subdirectory(naaice) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..29ba8d7 --- /dev/null +++ b/LICENSE @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Lesser General Public License instead.) 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 +this service 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 make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. 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. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute 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 and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +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 +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the 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 a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, 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. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE 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. + + 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 +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + tug + Copyright (C) 2023 naaice + + 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 2 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, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision 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, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This 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. diff --git a/README.org b/README.org index 2909fb8..c84711a 100644 --- a/README.org +++ b/README.org @@ -60,7 +60,7 @@ need to manually link your application/library to =tug=. into according =CMakeLists.txt= file: #+BEGIN_SRC - target_link_libraries(your_libapp tug) + target_link_libraries( tug) #+END_SRC 6. Build your application/library with =CMake=. @@ -121,6 +121,3 @@ If you can't get access to this =gitlab= instance: Also provide us the SHA of the commit you've downloaded. Thank you for your contributions in advance! - -* License -TODO? diff --git a/cmake-build-debug/_deps/doctest-src b/cmake-build-debug/_deps/doctest-src deleted file mode 160000 index b7c21ec..0000000 --- a/cmake-build-debug/_deps/doctest-src +++ /dev/null @@ -1 +0,0 @@ -Subproject commit b7c21ec5ceeadb4951b00396fc1e4642dd347e5f diff --git a/docs_sphinx/_build/html/_images/activity_diagram_run.svg b/docs_sphinx/_build/html/_images/activity_diagram_run.svg deleted file mode 100644 index 4f67009..0000000 --- a/docs_sphinx/_build/html/_images/activity_diagram_run.svg +++ /dev/null @@ -1,4 +0,0 @@ - - - -
execute run
execute run
Call FTCS routine
Call FTCS routine
Call BTCS routine
Call BTCS routine
[FTCS method]
[FTCS method]
[BTCS method]
[BTCS method]
[i<iterations]
[i<iterations]
[maximum
iterations reached]
[maximum...
[i<iterations]
[i<iterations]
[maximum
iterations reached]
[maximum...Text is not SVG - cannot display \ No newline at end of file diff --git a/docs_sphinx/_build/html/_images/class_diagram.svg b/docs_sphinx/_build/html/_images/class_diagram.svg deleted file mode 100644 index 3392a01..0000000 --- a/docs_sphinx/_build/html/_images/class_diagram.svg +++ /dev/null @@ -1,4 +0,0 @@ - - - -
Grid
Grid
- dim:int
- col:int
- row:int
- domain_col: int
- domain_row: int
- delta_col: double
- delta_row: double
- dim:int...
+ Grid(col: int)
+ Grid(row: int, col: int)
+ setConcentrations(concentrations: MatrixXd)
+ getConcentrations(): MatrixXd
+ setAlpha(alpha: MatrixXd)
+ setAlpha(alpha_x: MatrixXd, alpha_y: MatrixXd)
+ getAlphaX(): MatrixXd
+ getAlphaY(): MatrixXd
+ getDim(): int
+ getRow(): int
+ getCol(): int
+ setDomain(domain_col: int)
+ setDomain(domain_row: int, domain_col: int)
+ getDeltaCol(): double
+ getDeltaRow(): double
+ Grid(col: int)...
Boundary
Boundary
- grid: Grid
- type: BC_TYPE
- left: VectorXd
- right: VectorXd
- top: VectorXd
- bottom: VectorXd
- grid: Grid...
+ Boundary(grid: Grid, type: BC_TYPE)
+ getBoundaryConditionType(): BC_TYPE
+ setBoundaryConditionValue(side: BC_SIDE, value: VectorXd)
+ getBoundaryConditionValue(): VectorXd
+ Boundary(grid: Grid, type: BC_TYPE)...
Simulation
Simulation
- timestep: double
- iterations: int
- csv_output: CSV_OUTPUT
- grid: GRID
- bc: Boundary
- approach: APPROACH
- timestep: double...
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH
+ setOutputCSV(csv_output: CSV_OUTPUT)
+ setTimestep(timestep: double)
+ getTimestep(): double
+ setIterations(iterations: int)
+ getIterations(): int
+ printConcentrationsConsole()
+ printConcentrationsCSV(ident: string, appendMode: bool)
+ run()
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH...
FTCS
FTCS
BTCS
BTCS
functionalities
functionalities
TUG API
TUG API
BoundaryElement
BoundaryElement
- type: BC_TYPE
- value: double
- type: BC_TYPE...
+ BoundaryElement()
+ BoundaryElement(value: double)
+ setType(type: BC_TYPE)
+ setValue(value: double)
+ getType(): BC_TYPE
+ getValue(): double
+ BoundaryElement()...Text is not SVG - cannot display \ No newline at end of file diff --git a/docs_sphinx/_build/vis.html b/docs_sphinx/_build/vis.html deleted file mode 100644 index 8a11674..0000000 --- a/docs_sphinx/_build/vis.html +++ /dev/null @@ -1,224 +0,0 @@ - -
- - -
- - -
- - - \ No newline at end of file diff --git a/docs_sphinx/conf.py b/docs_sphinx/conf.py index 34477fc..f193aa9 100644 --- a/docs_sphinx/conf.py +++ b/docs_sphinx/conf.py @@ -22,7 +22,7 @@ subprocess.call('doxygen Doxyfile.in', shell=True) # -- Project information ----------------------------------------------------- project = 'TUG' -copyright = 'MIT' +copyright = 'GPL2' author = 'Philipp Ungrund, Hannes Signer' @@ -94,4 +94,4 @@ breathe_projects = { "Tug": "_build/xml/" } breathe_default_project = "Tug" -breathe_default_members = ('members', 'undoc-members') \ No newline at end of file +breathe_default_members = ('members', 'undoc-members') diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index fb04966..1a53738 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -1,48 +1,51 @@ +#include #include +using namespace Eigen; + int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** + // ************** + // **** GRID **** + // ************** - // create a linear grid with 20 cells - int cells = 20; - Grid grid = Grid(cells); + // create a linear grid with 20 cells + int cells = 20; + Grid grid = Grid(cells); - MatrixXd concentrations = MatrixXd::Constant(1,20,0); - concentrations(0,0) = 2000; - // TODO add option to set concentrations with a vector in 1D case - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); + concentrations(0, 0) = 2000; + // TODO add option to set concentrations with a vector in 1D case + grid.setConcentrations(concentrations); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set the number of iterations + simulation.setIterations(100); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); -} \ No newline at end of file + // run the simulation + simulation.run(); +} diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index c7b1a17..9a103cb 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -1,84 +1,85 @@ +#include #include +using namespace Eigen; + int main(int argc, char *argv[]) { - // EASY_PROFILER_ENABLE; - // profiler::startListen(); - // ************** - // **** GRID **** - // ************** - // profiler::startListen(); - // create a grid with a 20 x 20 field - int row = 40; - int col = 50; - Grid grid = Grid(row,col); + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // create a grid with a 20 x 20 field + int row = 40; + int col = 50; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - // grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(10,10) = 2000; - grid.setConcentrations(concentrations); - + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(10, 10) = 2000; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value - // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value - // grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value + // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value + // grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + // bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + // bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + // bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideClosed(BC_SIDE_LEFT); - bc.setBoundarySideClosed(BC_SIDE_RIGHT); - bc.setBoundarySideClosed(BC_SIDE_TOP); - bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - // bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - // bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - // bc.setBoundarySideConstant(BC_SIDE_TOP, 0); - // bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + // (optional) set boundary condition values for one side, e.g.: + // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values + // VectorXd bc_zero_values = VectorXd::Constant(20,0); + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + // VectorXd bc_front_values = VectorXd::Constant(20,2000); + // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // (optional) set boundary condition values for one side, e.g.: - // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - // VectorXd bc_zero_values = VectorXd::Constant(20,0); - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - // VectorXd bc_front_values = VectorXd::Constant(20,2000); - // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the number of iterations + simulation.setIterations(300); - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_XTREME); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // **** RUN SIMULATION **** - // set the number of iterations - simulation.setIterations(300); + // run the simulation - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_XTREME); - - // **** RUN SIMULATION **** - - // run the simulation - - // EASY_BLOCK("SIMULATION") - simulation.run(); - // EASY_END_BLOCK; - // profiler::dumpBlocksToFile("test_profile.prof"); - // profiler::stopListen(); -} \ No newline at end of file + // EASY_BLOCK("SIMULATION") + simulation.run(); + // EASY_END_BLOCK; + // profiler::dumpBlocksToFile("test_profile.prof"); + // profiler::stopListen(); +} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a4ebaaf..0a3173d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -2,6 +2,7 @@ add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp) add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp) add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp) add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) +add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp) add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) add_executable(profiling_openmp profiling_openmp.cpp) @@ -9,6 +10,7 @@ target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(FTCS_2D_proto_example tug) target_link_libraries(BTCS_1D_proto_example tug) target_link_libraries(BTCS_2D_proto_example tug) +target_link_libraries(CRNI_2D_proto_example tug) target_link_libraries(reference-FTCS_2D_closed tug) target_link_libraries(profiling_openmp tug) diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp new file mode 100644 index 0000000..39262e7 --- /dev/null +++ b/examples/CRNI_2D_proto_example.cpp @@ -0,0 +1,27 @@ +#include +#include + +using namespace Eigen; + +int main(int argc, char *argv[]) { + int row = 20; + int col = 20; + Grid grid(row, col); + + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(10, 10) = 2000; + grid.setConcentrations(concentrations); + + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + + Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); + simulation.setTimestep(0.1); + simulation.setIterations(50); + simulation.setOutputCSV(CSV_OUTPUT_XTREME); + + simulation.run(); +} diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index 3ac0df6..e97f0be 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -1,47 +1,50 @@ #include "tug/Boundary.hpp" +#include #include +using namespace Eigen; + int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** + // ************** + // **** GRID **** + // ************** - // create a linear grid with 20 cells - int cells = 20; - Grid grid = Grid(cells); + // create a linear grid with 20 cells + int cells = 20; + Grid grid = Grid(cells); - MatrixXd concentrations = MatrixXd::Constant(1,20,20); - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); + grid.setConcentrations(concentrations); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 1); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 1); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // (optional) set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // (optional) set the number of iterations + simulation.setIterations(100); - // (optional) set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); - // (optional) set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); -} \ No newline at end of file + // run the simulation + simulation.run(); +} diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 4ed9599..845d36f 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -1,85 +1,90 @@ /** * @file FTCS_2D_proto_closed_mdl.cpp * @author Hannes Signer, Philipp Ungrund, MDL - * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary condition; optional command line argument: number of cols and rows - * + * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary + * condition; optional command line argument: number of cols and rows + * */ +#include #include #include #include +using namespace Eigen; + int main(int argc, char *argv[]) { - int row = 64; + int row = 64; - if (argc == 2) { - // no cmd line argument, take col=row=64 - row = atoi(argv[1]); - } - int col=row; + if (argc == 2) { + // no cmd line argument, take col=row=64 + row = atoi(argv[1]); + } + int col = row; - std::cout << "Nrow =" << row << std::endl; - // ************** - // **** GRID **** - // ************** + std::cout << "Nrow =" << row << std::endl; + // ************** + // **** GRID **** + // ************** - // create a grid with a 20 x 20 field - int n2 = row/2-1; - Grid grid = Grid(row,col); + // create a grid with a 20 x 20 field + int n2 = row / 2 - 1; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(n2,n2) = 1; - concentrations(n2,n2+1) = 1; - concentrations(n2+1,n2) = 1; - concentrations(n2+1,n2+1) = 1; - grid.setConcentrations(concentrations); + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(n2, n2) = 1; + concentrations(n2, n2 + 1) = 1; + concentrations(n2 + 1, n2) = 1; + concentrations(n2 + 1, n2 + 1) = 1; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value - MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value - grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); - // create a boundary with constant values - Boundary bc = Boundary(grid); + // (optional) set boundary condition values for one side, e.g.: + bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - // (optional) set boundary condition values for one side, e.g.: - bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values - bc.setBoundarySideClosed(BC_SIDE_RIGHT); - bc.setBoundarySideClosed(BC_SIDE_TOP); - bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the timestep of the simulation + simulation.setTimestep(10000); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set the number of iterations + simulation.setIterations(100); - // set the timestep of the simulation - simulation.setTimestep(10000); // timestep + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** + // run the simulation + simulation.run(); - // run the simulation - simulation.run(); - - return 0; + return 0; } diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index e65d425..5f5e520 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -1,92 +1,91 @@ /** * @file FTCS_2D_proto_example.cpp * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach - * and constant boundary condition - * + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS + * approach and constant boundary condition + * */ +#include #include + +using namespace Eigen; // #include // #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); - int main(int argc, char *argv[]) { - // EASY_PROFILER_ENABLE; - // profiler::startListen(); - // ************** - // **** GRID **** - // ************** - // profiler::startListen(); - // create a grid with a 20 x 20 field - int row = 20; - int col = 20; - Grid grid = Grid(row,col); + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // create a grid with a 20 x 20 field + int row = 20; + int col = 20; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - // grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(0,0) = 1999; - grid.setConcentrations(concentrations); - + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(0, 0) = 1999; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value - // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value - // grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value + // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value + // grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - bc.setBoundarySideConstant(BC_SIDE_TOP, 0); - bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + // (optional) set boundary condition values for one side, e.g.: + // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values + // VectorXd bc_zero_values = VectorXd::Constant(20,0); + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + // VectorXd bc_front_values = VectorXd::Constant(20,2000); + // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // (optional) set boundary condition values for one side, e.g.: - // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - // VectorXd bc_zero_values = VectorXd::Constant(20,0); - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - // VectorXd bc_front_values = VectorXd::Constant(20,2000); - // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the number of iterations + simulation.setIterations(10000); - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // **** RUN SIMULATION **** - // set the number of iterations - simulation.setIterations(10000); + // run the simulation - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - - // **** RUN SIMULATION **** - - // run the simulation - - // EASY_BLOCK("SIMULATION") - simulation.run(); - // EASY_END_BLOCK; - // profiler::dumpBlocksToFile("test_profile.prof"); - // profiler::stopListen(); -} \ No newline at end of file + // EASY_BLOCK("SIMULATION") + simulation.run(); + // EASY_END_BLOCK; + // profiler::dumpBlocksToFile("test_profile.prof"); + // profiler::stopListen(); +} diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 80d4112..2c63c0c 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -1,77 +1,81 @@ /** * @file FTCS_2D_proto_example.cpp * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach - * and constant boundary condition - * + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS + * approach and constant boundary condition + * */ +#include #include +using namespace Eigen; + int main(int argc, char *argv[]) { - - // ************** - // **** GRID **** - // ************** - // create a grid with a 20 x 20 field - int row = 64; - int col = 64; - int n2 = row/2-1; - Grid grid = Grid(row,col); + // ************** + // **** GRID **** + // ************** - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // create a grid with a 20 x 20 field + int row = 64; + int col = 64; + int n2 = row / 2 - 1; + Grid grid = Grid(row, col); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(n2,n2) = 1; - concentrations(n2,n2+1) = 1; - concentrations(n2+1,n2) = 1; - concentrations(n2+1,n2+1) = 1; - grid.setConcentrations(concentrations); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set alphas of the grid, e.g.: - MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value - MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value - grid.setAlpha(alphax, alphay); + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(n2, n2) = 1; + concentrations(n2, n2 + 1) = 1; + concentrations(n2 + 1, n2) = 1; + concentrations(n2 + 1, n2 + 1) = 1; + grid.setConcentrations(concentrations); + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); - // ****************** - // **** BOUNDARY **** - // ****************** + // ****************** + // **** BOUNDARY **** + // ****************** - // create a boundary with constant values - Boundary bc = Boundary(grid); + // create a boundary with constant values + Boundary bc = Boundary(grid); - // (optional) set boundary condition values for one side, e.g.: - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - bc.setBoundarySideConstant(BC_SIDE_TOP, 0); - bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + // (optional) set boundary condition values for one side, e.g.: + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // (optional) set the timestep of the simulation + simulation.setTimestep(1000); // timestep - // (optional) set the timestep of the simulation - simulation.setTimestep(1000); // timestep + // (optional) set the number of iterations + simulation.setIterations(5); - // (optional) set the number of iterations - simulation.setIterations(5); + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); - - // **** RUN SIMULATION **** + // **** RUN SIMULATION **** - // run the simulation - simulation.run(); + // run the simulation + simulation.run(); - return 0; + return 0; } diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index 1c83542..5dc13ef 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -1,61 +1,70 @@ -#include -#include -#include +#include #include +#include +#include +#include +#include + +using namespace Eigen; +using namespace std; int main(int argc, char *argv[]) { - int n[4] = {100, 500, 1000, 2000}; - int threads[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - int iterations[1] = {5}; - int repetition = 1; + int n[] = {2000}; + int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + int iterations[1] = {1}; + int repetition = 10; + for (int l = 0; l < size(threads); l++) { + // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv"; ofstream myfile; - myfile.open("testLarge.csv"); + myfile.open("speedup_1000.csv", std::ios::app); + myfile << "Number threads: " << threads[l] << endl; - for (int i = 0; i < size(n); i++){ - cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - for(int j = 0; j < size(iterations); j++){ - cout << "Iterations: " << iterations[j] << endl; - myfile << "Iterations: " << iterations[j] << endl; - for (int k = 0; k < repetition; k++){ - cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); - grid.setDomain(1, 1); + for (int i = 0; i < size(n); i++) { + cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; + // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl; + for (int j = 0; j < size(iterations); j++) { + cout << "Iterations: " << iterations[j] << endl; + // myfile << "Iterations: " << iterations[j] << endl; + for (int k = 0; k < repetition; k++) { + cout << "Wiederholung: " << k << endl; + Grid grid = Grid(n[i], n[i]); + grid.setDomain(1, 1); - MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); - concentrations(n[i]/2,n[i]/2) = 1; - grid.setConcentrations(concentrations); - MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); + MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); + concentrations(n[i] / 2, n[i] / 2) = 1; + grid.setConcentrations(concentrations); + MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); - Boundary bc = Boundary(grid); + Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setSolver(THOMAS_ALGORITHM_SOLVER); + Simulation sim = Simulation(grid, bc, BTCS_APPROACH); + sim.setSolver(THOMAS_ALGORITHM_SOLVER); - if(argc == 2){ - int numThreads = atoi(argv[1]); - sim.setNumberThreads(numThreads); - } - else{ - sim.setNumberThreads(1); - } + if (argc == 2) { + int numThreads = atoi(argv[1]); + sim.setNumberThreads(numThreads); + } else { + sim.setNumberThreads(threads[l]); + } - sim.setTimestep(0.001); - sim.setIterations(iterations[j]); - sim.setOutputCSV(CSV_OUTPUT_OFF); + sim.setTimestep(0.01); + sim.setIterations(iterations[j]); + sim.setOutputCSV(CSV_OUTPUT_OFF); - auto begin = std::chrono::high_resolution_clock::now(); - sim.run(); - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = std::chrono::duration_cast(end - begin); - myfile << milliseconds.count() << endl; - } + auto begin = std::chrono::high_resolution_clock::now(); + sim.run(); + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - + begin); + myfile << milliseconds.count() << endl; } - cout << endl; - myfile << endl; - + } + cout << endl; + myfile << endl; } myfile.close(); -} \ No newline at end of file + } +} diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp new file mode 100644 index 0000000..5dc13ef --- /dev/null +++ b/examples/profiling_speedup.cpp @@ -0,0 +1,70 @@ +#include +#include +#include +#include +#include +#include + +using namespace Eigen; +using namespace std; + +int main(int argc, char *argv[]) { + + int n[] = {2000}; + int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + int iterations[1] = {1}; + int repetition = 10; + + for (int l = 0; l < size(threads); l++) { + // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv"; + ofstream myfile; + myfile.open("speedup_1000.csv", std::ios::app); + myfile << "Number threads: " << threads[l] << endl; + + for (int i = 0; i < size(n); i++) { + cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; + // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl; + for (int j = 0; j < size(iterations); j++) { + cout << "Iterations: " << iterations[j] << endl; + // myfile << "Iterations: " << iterations[j] << endl; + for (int k = 0; k < repetition; k++) { + cout << "Wiederholung: " << k << endl; + Grid grid = Grid(n[i], n[i]); + grid.setDomain(1, 1); + + MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); + concentrations(n[i] / 2, n[i] / 2) = 1; + grid.setConcentrations(concentrations); + MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); + + Boundary bc = Boundary(grid); + + Simulation sim = Simulation(grid, bc, BTCS_APPROACH); + sim.setSolver(THOMAS_ALGORITHM_SOLVER); + + if (argc == 2) { + int numThreads = atoi(argv[1]); + sim.setNumberThreads(numThreads); + } else { + sim.setNumberThreads(threads[l]); + } + + sim.setTimestep(0.01); + sim.setIterations(iterations[j]); + sim.setOutputCSV(CSV_OUTPUT_OFF); + + auto begin = std::chrono::high_resolution_clock::now(); + sim.run(); + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - + begin); + myfile << milliseconds.count() << endl; + } + } + cout << endl; + myfile << endl; + } + myfile.close(); + } +} diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index 7b51e67..340a6fd 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -1,55 +1,52 @@ -#include #include "Eigen/Core" #include +#include using namespace std; +using namespace Eigen; int main(int argc, char *argv[]) { - int row = 50; - int col = 50; - int domain_row = 10; - int domain_col = 10; + int row = 50; + int col = 50; + int domain_row = 10; + int domain_col = 10; + // Grid + Grid grid = Grid(row, col); + grid.setDomain(domain_row, domain_col); - // Grid - Grid grid = Grid(row, col); - grid.setDomain(domain_row, domain_col); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(5, 5) = 1; + grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5,5) = 1; - grid.setConcentrations(concentrations); - - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } + MatrixXd alpha = MatrixXd::Constant(row, col, 1); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 6; j++) { + alpha(i, j) = 0.01; } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } + } + for (int i = 0; i < 5; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.001; } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } + } + for (int i = 5; i < 11; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.1; } - grid.setAlpha(alpha, alpha); + } + grid.setAlpha(alpha, alpha); + // Boundary + Boundary bc = Boundary(grid); - // Boundary - Boundary bc = Boundary(grid); + // Simulation + Simulation sim = Simulation(grid, bc, FTCS_APPROACH); + sim.setTimestep(0.001); + sim.setIterations(10000); + sim.setOutputCSV(CSV_OUTPUT_OFF); + sim.setOutputConsole(CONSOLE_OUTPUT_OFF); - - // Simulation - Simulation sim = Simulation(grid, bc, FTCS_APPROACH); - sim.setTimestep(0.001); - sim.setIterations(10000); - sim.setOutputCSV(CSV_OUTPUT_OFF); - sim.setOutputConsole(CONSOLE_OUTPUT_OFF); - - - // RUN - sim.run(); -} \ No newline at end of file + // RUN + sim.run(); +} diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index 1109a1c..b803e75 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -1,37 +1,26 @@ /** * @file Boundary.hpp - * @brief API of Boundary class, that holds all information for each boundary condition - * at the edges of the diffusion grid. - * + * @brief API of Boundary class, that holds all information for each boundary + * condition at the edges of the diffusion grid. + * */ #ifndef BOUNDARY_H_ #define BOUNDARY_H_ -#include #include "Grid.hpp" - -using namespace std; -using namespace Eigen; +#include /** - * @brief Enum defining the two implemented boundary conditions. - * + * @brief Enum defining the two implemented boundary conditions. + * */ -enum BC_TYPE { - BC_TYPE_CLOSED, - BC_TYPE_CONSTANT -}; +enum BC_TYPE { BC_TYPE_CLOSED, BC_TYPE_CONSTANT }; /** - * @brief Enum defining all 4 possible sides to a 1D and 2D grid. - * + * @brief Enum defining all 4 possible sides to a 1D and 2D grid. + * */ -enum BC_SIDE { - BC_SIDE_LEFT, - BC_SIDE_RIGHT, - BC_SIDE_TOP, - BC_SIDE_BOTTOM -}; +enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; /** * This class defines the boundary conditions of individual boundary elements. @@ -39,175 +28,184 @@ enum BC_SIDE { * The class serves as an auxiliary class for structuring the Boundary class. */ class BoundaryElement { - public: - - /** - * @brief Construct a new Boundary Element object for the closed case. - * The boundary type is here automatically set to the type - * BC_TYPE_CLOSED, where the value takes NaN. - */ - BoundaryElement(); +public: + /** + * @brief Construct a new Boundary Element object for the closed case. + * The boundary type is here automatically set to the type + * BC_TYPE_CLOSED, where the value takes -1 and does not hold any + * physical meaning. + */ + BoundaryElement(); - /** - * @brief Construct a new Boundary Element object for the constant case. - * The boundary type is automatically set to the type - * BC_TYPE_CONSTANT. - * - * @param value Value of the constant concentration to be assumed at the - * corresponding boundary element. - */ - BoundaryElement(double value); + /** + * @brief Construct a new Boundary Element object for the constant case. + * The boundary type is automatically set to the type + * BC_TYPE_CONSTANT. + * + * @param value Value of the constant concentration to be assumed at the + * corresponding boundary element. + */ + BoundaryElement(double value); - /** - * @brief Allows changing the boundary type of a corresponding - * BoundaryElement object. - * - * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or - BC_TYPE_CLOSED. - */ - void setType(BC_TYPE type); - - /** - * @brief Sets the value of a boundary condition for the constant case. - * - * @param value Concentration to be considered constant for the - * corresponding boundary element. - */ - void setValue(double value); + /** + * @brief Allows changing the boundary type of a corresponding + * BoundaryElement object. + * + * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or + BC_TYPE_CLOSED. + */ + void setType(BC_TYPE type); - /** - * @brief Return the type of the boundary condition, i.e. whether the - * boundary is considered closed or constant. - * - * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or - BC_TYPE_CONSTANT. - */ - BC_TYPE getType(); + /** + * @brief Sets the value of a boundary condition for the constant case. + * + * @param value Concentration to be considered constant for the + * corresponding boundary element. + */ + void setValue(double value); - /** - * @brief Return the concentration value for the constant boundary condition. - * - * @return double Value of the concentration. - */ - double getValue(); + /** + * @brief Return the type of the boundary condition, i.e. whether the + * boundary is considered closed or constant. + * + * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or + BC_TYPE_CONSTANT. + */ + BC_TYPE getType(); - private: - BC_TYPE type; - double value; + /** + * @brief Return the concentration value for the constant boundary condition. + * + * @return double Value of the concentration. + */ + double getValue(); + +private: + BC_TYPE type{BC_TYPE_CLOSED}; + double value{-1}; }; - /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. */ class Boundary { - public: - /** - * @brief Creates a boundary object based on the passed grid object and - * initializes the boundaries as closed. - * - * @param grid Grid object on the basis of which the simulation is to take place. - */ - Boundary(Grid grid); +public: + /** + * @brief Creates a boundary object based on the passed grid object and + * initializes the boundaries as closed. + * + * @param grid Grid object on the basis of which the simulation takes place + * and from which the dimensions (in 2D case) are taken. + */ + Boundary(Grid grid); - /** - * @brief Sets all elements of the specified boundary side to the boundary - * condition closed. - * - * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. - */ - void setBoundarySideClosed(BC_SIDE side); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition closed. + * + * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. + */ + void setBoundarySideClosed(BC_SIDE side); - /** - * @brief Sets all elements of the specified boundary side to the boundary - * condition constant. Thereby the concentration values of the - * boundaries are set to the passed value. - * - * @param side Side to be set to constant, e.g. BC_SIDE_LEFT. - * @param value Concentration to be set for all elements of the specified page. - */ - void setBoundarySideConstant(BC_SIDE side, double value); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition constant. Thereby the concentration values of the + * boundaries are set to the passed value. + * + * @param side Side to be set to constant, e.g. BC_SIDE_LEFT. + * @param value Concentration to be set for all elements of the specified + * page. + */ + void setBoundarySideConstant(BC_SIDE side, double value); - /** - * @brief Specifically sets the boundary element of the specified side - * defined by the index to the boundary condition closed. - * - * @param side Side in which an element is to be defined as closed. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - */ - void setBoundaryElementClosed(BC_SIDE side, int index); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition closed. + * + * @param side Side in which an element is to be defined as closed. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + */ + void setBoundaryElementClosed(BC_SIDE side, int index); - /** - * @brief Specifically sets the boundary element of the specified side - * defined by the index to the boundary condition constant with the - given concentration value. - * - * @param side Side in which an element is to be defined as constant. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @param value Concentration value to which the boundary element should be set. - */ - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition constant with the + given concentration value. + * + * @param side Side in which an element is to be defined as constant. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + side. + * @param value Concentration value to which the boundary element should be + set. + */ + void setBoundaryElementConstant(BC_SIDE side, int index, double value); - /** - * @brief Returns the boundary condition of a specified side as a vector - * of BoundarsElement objects. - * - * @param side Boundary side from which the boundary conditions are to be returned. - * @return vector Contains the boundary conditions as BoundaryElement objects. - */ - vector getBoundarySide(BC_SIDE side); + /** + * @brief Returns the boundary condition of a specified side as a vector + * of BoundarsElement objects. + * + * @param side Boundary side from which the boundary conditions are to be + * returned. + * @return vector Contains the boundary conditions as + * BoundaryElement objects. + */ + const std::vector getBoundarySide(BC_SIDE side); - /** - * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific - boundary is closed. - * - * @param side Boundary side for which the values are to be returned. - * @return VectorXd Vector with values as doubles. - */ - VectorXd getBoundarySideValues(BC_SIDE side); + /** + * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some + specific boundary is closed. + * + * @param side Boundary side for which the values are to be returned. + * @return VectorXd Vector with values as doubles. + */ + Eigen::VectorXd getBoundarySideValues(BC_SIDE side); - /** - * @brief Returns the boundary condition of a specified element on a given side. - * - * @param side Boundary side in which the boundary condition is located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @return BoundaryElement Boundary condition as a BoundaryElement object. - */ - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + /** + * @brief Returns the boundary condition of a specified element on a given + * side. + * + * @param side Boundary side in which the boundary condition is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + * @return BoundaryElement Boundary condition as a BoundaryElement object. + */ + BoundaryElement getBoundaryElement(BC_SIDE side, int index); - /** - * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED or - BC_TYPE_CONSTANT. - * - * @param side Boundary side in which the boundary condition type is located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @return BC_TYPE Boundary Type of the corresponding boundary condition. - */ - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + /** + * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED + or BC_TYPE_CONSTANT. + * + * @param side Boundary side in which the boundary condition type is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + side. + * @return BC_TYPE Boundary Type of the corresponding boundary condition. + */ + BC_TYPE getBoundaryElementType(BC_SIDE side, int index); - /** - * @brief Returns the concentration value of a corresponding - * BoundaryElement object if it is a constant boundary condition. - * - * @param side Boundary side in which the boundary condition value is - * located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding - * side. - * @return double Concentration of the corresponding BoundaryElement object. - */ - double getBoundaryElementValue(BC_SIDE side, int index); + /** + * @brief Returns the concentration value of a corresponding + * BoundaryElement object if it is a constant boundary condition. + * + * @param side Boundary side in which the boundary condition value is + * located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + * @return double Concentration of the corresponding BoundaryElement object. + */ + double getBoundaryElementValue(BC_SIDE side, int index); - private: - Grid grid; // Boundary is directly dependent on the dimensions of a predefined - - vector> boundaries; // Vector with Boundary Element information +private: + Grid grid; // Boundary is directly dependent on the dimensions of a predefined + + std::vector> + boundaries; // Vector with Boundary Element information }; -#endif - +#endif // BOUNDARY_H_ diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index db88429..6a6f973 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -1,171 +1,179 @@ +#ifndef GRID_H_ +#define GRID_H_ + /** * @file Grid.hpp - * @brief API of Grid class, that holds a matrix with concenctrations and a - * respective matrix/matrices of alpha coefficients. - * + * @brief API of Grid class, that holds a matrix with concenctrations and a + * respective matrix/matrices of alpha coefficients. + * */ #include #include -using namespace Eigen; - -// TODO document default values and perhaps adjust them class Grid { - public: +public: + /** + * @brief Constructs a new 1D-Grid object of a given length, which holds a + * matrix with concentrations and a respective matrix of alpha coefficients. + * The domain length is per default the same as the length. The + * concentrations are all 20 by default and the alpha coefficients are 1. + * + * @param length Length of the 1D-Grid. Must be greater than 3. + */ + Grid(int length); - /** - * @brief Constructs a new 1D-Grid object of a given length, which holds a matrix - * with concentrations and a respective matrix of alpha coefficients. - * - * @param length Length of the 1D-Grid. Must be greater than 3. - */ - Grid(int length); + /** + * @brief Constructs a new 2D-Grid object of given dimensions, which holds a + * matrix with concentrations and the respective matrices of alpha coefficient + * for each direction. The domain in x- and y-direction is per default equal + * to the col length and row length, respectively. The concentrations are all + * 20 by default across the entire grid and the alpha coefficients 1 in both + * directions. + * + * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. + * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. + */ + Grid(int row, int col); - /** - * @brief Constructs a new 2D-Grid object of given dimensions, which holds a matrix - * with concentrations and the respective matrices of alpha coefficient for - * each direction. - * - * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. - * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. - */ - Grid(int row, int col); + /** + * @brief Sets the concentrations matrix for a 1D or 2D-Grid. + * + * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix + * must have correct dimensions as defined in row and col. (Or length, in 1D + * case). + */ + void setConcentrations(const Eigen::MatrixXd &concentrations); - /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. - * - * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix must - * have correct dimensions as defined in row and col, or length, - * respectively. - */ - void setConcentrations(MatrixXd concentrations); + /** + * @brief Gets the concentrations matrix for a Grid. + * + * @return MatrixXd An Eigen3 matrix holding the concentrations and having the + * same dimensions as the grid. + */ + const Eigen::MatrixXd &getConcentrations() { return this->concentrations; } - /** - * @brief Gets the concentrations matrix for a Grid. - * - * @return MatrixXd An Eigen3 matrix holding the concentrations and having the - * same dimensions as the grid. - */ - MatrixXd getConcentrations(); + /** + * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one + * dimensional. + * + * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. + * Matrix columns must have same size as length of grid. + */ + void setAlpha(const Eigen::MatrixXd &alpha); - /** - * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one dimensional. - * - * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. Matrix - * columns must have same size as length of grid. - */ - void setAlpha(MatrixXd alpha); + /** + * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two + * dimensional. + * + * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in + * x-direction. Matrix must be of same size as the grid. + * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in + * y-direction. Matrix must be of same size as the grid. + */ + void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY); - /** - * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two dimensional. - * - * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in x-direction. - * Matrix must be of same size as the grid. - * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in y-direction. - * Matrix must be of same size as the grid. - */ - void setAlpha(MatrixXd alphaX, MatrixXd alphaY); + /** + * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one + * dimensional. + * + * @return MatrixXd A matrix with 1 row holding the alpha coefficients. + */ + const Eigen::MatrixXd &getAlpha(); - /** - * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one dimensional. - * - * @return MatrixXd A matrix with 1 row holding the alpha coefficients. - */ - MatrixXd getAlpha(); + /** + * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. + * Grid must be two dimensional. + * + * @return MatrixXd A matrix holding the alpha coefficients in x-direction. + */ + const Eigen::MatrixXd &getAlphaX(); - /** - * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. Grid must be - * two dimensional. - * - * @return MatrixXd A matrix holding the alpha coefficients in x-direction. - */ - MatrixXd getAlphaX(); + /** + * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. + * Grid must be two dimensional. + * + * @return MatrixXd A matrix holding the alpha coefficients in y-direction. + */ + const Eigen::MatrixXd &getAlphaY(); - /** - * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. Grid must be - * two dimensional. - * - * @return MatrixXd A matrix holding the alpha coefficients in y-direction. - */ - MatrixXd getAlphaY(); + /** + * @brief Gets the dimensions of the grid. + * + * @return int Dimensions, either 1 or 2. + */ + int getDim(); - /** - * @brief Gets the dimensions of the grid. - * - * @return int Dimensions, either 1 or 2. - */ - int getDim(); + /** + * @brief Gets length of 1D grid. Must be one dimensional grid. + * + * @return int Length of 1D grid. + */ + int getLength(); - /** - * @brief Gets length of 1D grid. Must be one dimensional grid. - * - * @return int Length of 1D grid. - */ - int getLength(); + /** + * @brief Gets the number of rows of the grid. + * + * @return int Number of rows. + */ + int getRow(); - /** - * @brief Gets the number of rows of the grid. - * - * @return int Number of rows. - */ - int getRow(); + /** + * @brief Gets the number of columns of the grid. + * + * @return int Number of columns. + */ + int getCol(); - /** - * @brief Gets the number of columns of the grid. - * - * @return int Number of columns. - */ - int getCol(); + /** + * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. + * + * @param domainLength A double value of the domain length. Must be positive. + */ + void setDomain(double domainLength); - /** - * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. - * - * @param domainLength A double value of the domain length. Must be positive. - */ - void setDomain(double domainLength); + /** + * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. + * + * @param domainRow A double value of the domain size in y-direction. Must be + * positive. + * @param domainCol A double value of the domain size in x-direction. Must be + * positive. + */ + void setDomain(double domainRow, double domainCol); - /** - * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. - * - * @param domainRow A double value of the domain size in y-direction. Must be positive. - * @param domainCol A double value of the domain size in x-direction. Must be positive. - */ - void setDomain(double domainRow,double domainCol); + /** + * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. + * + * @return double Delta value. + */ + double getDelta(); - /** - * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. - * - * @return double Delta value. - */ - double getDelta(); + /** + * @brief Gets the delta value in x-direction. + * + * @return double Delta value in x-direction. + */ + double getDeltaCol(); - /** - * @brief Gets the delta value in x-direction. - * - * @return double Delta value in x-direction. - */ - double getDeltaCol(); - - /** - * @brief Gets the delta value in y-direction. Must be two dimensional grid. - * - * @return double Delta value in y-direction. - */ - double getDeltaRow(); - - - private: - - int col; // number of grid columns - int row; // number of grid rows - int dim; // 1D or 2D - double domainCol; // number of domain columns - double domainRow; // number of domain rows - double deltaCol; // delta in x-direction (between columns) - double deltaRow; // delta in y-direction (between rows) - MatrixXd concentrations; // Matrix holding grid concentrations - MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction - MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction + /** + * @brief Gets the delta value in y-direction. Must be two dimensional grid. + * + * @return double Delta value in y-direction. + */ + double getDeltaRow(); +private: + int col; // number of grid columns + int row{1}; // number of grid rows + int dim; // 1D or 2D + double domainCol; // number of domain columns + double domainRow{0}; // number of domain rows + double deltaCol; // delta in x-direction (between columns) + double deltaRow{0}; // delta in y-direction (between rows) + Eigen::MatrixXd concentrations; // Matrix holding grid concentrations + Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction + Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction }; + +#endif // GRID_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 3b43212..0c09a7e 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -1,62 +1,74 @@ /** * @file Simulation.hpp - * @brief API of Simulation class, that holds all information regarding a specific simulation - * run like its timestep, number of iterations and output options. Simulation object - * also holds a predefined Grid and Boundary object. + * @brief API of Simulation class, that holds all information regarding a + * specific simulation run like its timestep, number of iterations and output + * options. Simulation object also holds a predefined Grid and Boundary object. * */ -#include "Boundary.hpp" -#include -using namespace std; +#ifndef SIMULATION_H_ +#define SIMULATION_H_ + +#include "Boundary.hpp" +#include "Grid.hpp" +#include +#include + +#ifdef _OPENMP +#include +#else +#define omp_get_num_procs() 1 +#endif /** - * @brief Enum defining the two implemented solution approaches. - * + * @brief Enum defining the two implemented solution approaches. + * */ enum APPROACH { - FTCS_APPROACH, // Forward Time-Centered Space - BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver - CRANK_NICOLSON_APPROACH + FTCS_APPROACH, // Forward Time-Centered Space + BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver + CRANK_NICOLSON_APPROACH }; /** * @brief Enum defining the Linear Equation solvers - * + * */ enum SOLVER { - EIGEN_LU_SOLVER, // EigenLU solver - THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for tridiagonal matrices + EIGEN_LU_SOLVER, // EigenLU solver + THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for + // tridiagonal matrices }; /** * @brief Enum holding different options for .csv output. - * + * */ enum CSV_OUTPUT { - CSV_OUTPUT_OFF, // do not produce csv output - CSV_OUTPUT_ON, // produce csv output with last concentration matrix - CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices - CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary conditions at beginning + CSV_OUTPUT_OFF, // do not produce csv output + CSV_OUTPUT_ON, // produce csv output with last concentration matrix + CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices + CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary + // conditions at beginning }; /** * @brief Enum holding different options for console output. - * + * */ enum CONSOLE_OUTPUT { - CONSOLE_OUTPUT_OFF, // do not print any output to console - CONSOLE_OUTPUT_ON, // print before and after concentrations to console - CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console + CONSOLE_OUTPUT_OFF, // do not print any output to console + CONSOLE_OUTPUT_ON, // print before and after concentrations to console + CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console }; /** - * @brief Enum holding different options for time measurement. - * + * @brief Enum holding different options for time measurement. + * */ enum TIME_MEASURE { - TIME_MEASURE_OFF, // do not print any time measures - TIME_MEASURE_ON // print time measure after last iteration + TIME_MEASURE_OFF, // do not print any time measures + TIME_MEASURE_ON // print time measure after last iteration }; /** @@ -66,147 +78,159 @@ enum TIME_MEASURE { * */ class Simulation { - public: - /** - * @brief Set up a simulation environment. The timestep and number of iterations - * must be set. For the BTCS approach, the Thomas algorithm is used as - * the default linear equation solver as this is faster for tridiagonal - * matrices. CSV output, console output and time measure are off by default. - * - * @param grid Valid grid object - * @param bc Valid boundary condition object - * @param approach Approach to solving the problem. Either FTCS or BTCS. - */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); +public: + /** + * @brief Set up a simulation environment. The timestep and number of + * iterations must be set. For the BTCS approach, the Thomas algorithm is used + * as the default linear equation solver as this is faster for tridiagonal + * matrices. CSV output, console output and time measure are off by + * default. Also, the number of cores is set to the maximum number of cores -1 + * by default. + * + * @param grid Valid grid object + * @param bc Valid boundary condition object + * @param approach Approach to solving the problem. Either FTCS or BTCS. + */ + Simulation(Grid &grid, Boundary &bc, APPROACH approach); - /** - * @brief Set the option to output the results to a CSV file. Off by default. - * - * - * @param csv_output Valid output option. The following options can be set - * here: - * - CSV_OUTPUT_OFF: do not produce csv output - * - CSV_OUTPUT_ON: produce csv output with last - * concentration matrix - * - CSV_OUTPUT_VERBOSE: produce csv output with all - * concentration matrices - * - CSV_OUTPUT_XTREME: produce csv output with all - * concentration matrices and simulation environment - */ - void setOutputCSV(CSV_OUTPUT csv_output); + /** + * @brief Set the option to output the results to a CSV file. Off by default. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - CSV_OUTPUT_OFF: do not produce csv output + * - CSV_OUTPUT_ON: produce csv output with last + * concentration matrix + * - CSV_OUTPUT_VERBOSE: produce csv output with all + * concentration matrices + * - CSV_OUTPUT_XTREME: produce csv output with all + * concentration matrices and simulation environment + */ + void setOutputCSV(CSV_OUTPUT csv_output); - /** - * @brief Set the options for outputting information to the console. Off by default. - * - * @param console_output Valid output option. The following options can be set - * here: - * - CONSOLE_OUTPUT_OFF: do not print any output to console - * - CONSOLE_OUTPUT_ON: print before and after concentrations to console - * - CONSOLE_OUTPUT_VERBOSE: print all concentration matrices to console - */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + /** + * @brief Set the options for outputting information to the console. Off by + * default. + * + * @param console_output Valid output option. The following options can be set + * here: + * - CONSOLE_OUTPUT_OFF: do not print any output to + * console + * - CONSOLE_OUTPUT_ON: print before and after + * concentrations to console + * - CONSOLE_OUTPUT_VERBOSE: print all concentration + * matrices to console + */ + void setOutputConsole(CONSOLE_OUTPUT console_output); - // TODO document method - /** - * @brief Set the Time Measure object. Off by default. - * - * @param time_measure - */ - void setTimeMeasure(TIME_MEASURE time_measure); + /** + * @brief Set the Time Measure option. Off by default. + * + * @param time_measure The following options are allowed: + * - TIME_MEASURE_OFF: Time of simulation is not printed + * to console + * - TIME_MEASURE_ON: Time of simulation run is printed to + * console + */ + void setTimeMeasure(TIME_MEASURE time_measure); - /** - * @brief Setting the time step for each iteration step. Time step must be - * greater than zero. Setting the timestep is required. - * - * @param timestep Valid timestep greater than zero. - */ - void setTimestep(double timestep); + /** + * @brief Setting the time step for each iteration step. Time step must be + * greater than zero. Setting the timestep is required. + * + * @param timestep Valid timestep greater than zero. + */ + void setTimestep(double timestep); - /** - * @brief Currently set time step is returned. - * - * @return double timestep - */ - double getTimestep(); + /** + * @brief Currently set time step is returned. + * + * @return double timestep + */ + double getTimestep(); - /** - * @brief Set the desired iterations to be calculated. A value greater - * than zero must be specified here. Setting iterations is required. - * - * @param iterations Number of iterations to be simulated. - */ - void setIterations(int iterations); + /** + * @brief Set the desired iterations to be calculated. A value greater + * than zero must be specified here. Setting iterations is required. + * + * @param iterations Number of iterations to be simulated. + */ + void setIterations(int iterations); - /** - * @brief Set the desired linear equation solver to be used for BTCS approach. Without effect - * in case of FTCS approach. - * - * @param solver Solver to be used. Default is Thomas Algorithm as it is more efficient for - * tridiagonal Matrices. - */ - void setSolver(SOLVER solver); + /** + * @brief Set the desired linear equation solver to be used for BTCS approach. + * Without effect in case of FTCS approach. + * + * @param solver Solver to be used. Default is Thomas Algorithm as it is more + * efficient for tridiagonal Matrices. + */ + void setSolver(SOLVER solver); - /** - * @brief Set the number of desired openMP Threads. - * - * @param num_threads Number of desired threads. Must have a value between - * 1 and the maximum available number of processors. The maximum number of - * processors is set as the default case. - */ - void setNumberThreads(int num_threads); + /** + * @brief Set the number of desired openMP Threads. + * + * @param num_threads Number of desired threads. Must have a value between + * 1 and the maximum available number of processors. The + * maximum number of processors is set as the default case during Simulation + * construction. + */ + void setNumberThreads(int num_threads); - /** - * @brief Return the currently set iterations to be calculated. - * - * @return int Number of iterations. - */ - int getIterations(); + /** + * @brief Return the currently set iterations to be calculated. + * + * @return int Number of iterations. + */ + int getIterations(); - /** - * @brief Outputs the current concentrations of the grid on the console. - * - */ - void printConcentrationsConsole(); + /** + * @brief Outputs the current concentrations of the grid on the console. + * + */ + void printConcentrationsConsole(); - // TODO move create CSVfile to TugUtils - /** - * @brief Creates a CSV file with a name containing the current simulation - * parameters. If the data name already exists, an additional counter is - * appended to the name. The name of the file is built up as follows: - * + + + +.csv - * - * @return string Filename with given simulation parameter. - */ - string createCSVfile(); + /** + * @brief Creates a CSV file with a name containing the current simulation + * parameters. If the data name already exists, an additional counter + * is appended to the name. The name of the file is built up as follows: + * + + + +.csv + * + * @return string Filename with configured simulation parameters. + */ + std::string createCSVfile(); - /** - * @brief Writes the currently calculated concentration values of the grid - * into the CSV file with the passed filename. - * - * @param filename Name of the file to which the concentration values are - * to be written. - */ - void printConcentrationsCSV(string filename); + /** + * @brief Writes the currently calculated concentration values of the grid + * into the CSV file with the passed filename. + * + * @param filename Name of the file to which the concentration values are + * to be written. + */ + void printConcentrationsCSV(const std::string &filename); - /** - * @brief Method starts the simulation process with the previously set - * parameters. - */ - void run(); + /** + * @brief Method starts the simulation process with the previously set + * parameters. + */ + void run(); - private: +private: + double timestep{-1}; + int iterations{-1}; + int innerIterations{1}; + int numThreads{omp_get_num_procs()}; + CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; + CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; + TIME_MEASURE time_measure{TIME_MEASURE_OFF}; - double timestep; - int iterations; - int innerIterations; - int numThreads; - CSV_OUTPUT csv_output; - CONSOLE_OUTPUT console_output; - TIME_MEASURE time_measure; - - Grid &grid; - Boundary &bc; - APPROACH approach; - SOLVER solver; + Grid &grid; + Boundary &bc; + APPROACH approach; + SOLVER solver{THOMAS_ALGORITHM_SOLVER}; + const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; }; + +#endif // SIMULATION_H_ diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index 6210cd8..e2f4ef5 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -286,7 +286,11 @@ ADIHetDir <- function(field, dt, iter, alpha) { return(out) } -harm <- function(x,y) 1/(1/x+1/y) +harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & z have different lengths") + 2/(1/x+1/y) +} harm(1,4) @@ -428,10 +432,10 @@ FTCS_2D <- function(field, dt, iter, alpha) { ## res[i,j] *(mean(alpha[i,j+1],alpha[i,j])+mean(alpha[i,j-1],alpha[i,j])) + ## res[i,j-1]*mean(alpha[i,j-1],alpha[i,j])) tmp[i,j] <- res[i,j] + - + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * mean(alpha[i+1,j],alpha[i,j]) - - (res[i,j]-res[i-1,j]) * mean(alpha[i-1,j],alpha[i,j])) + - + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * mean(alpha[i,j+1],alpha[i,j]) - - (res[i,j]-res[i,j-1]) * mean(alpha[i,j-1],alpha[i,j])) + + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * harm(alpha[i+1,j],alpha[i,j]) - + (res[i,j]-res[i-1,j]) * harm(alpha[i-1,j],alpha[i,j])) + + + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * harm(alpha[i,j+1],alpha[i,j]) - + (res[i,j]-res[i,j-1]) * harm(alpha[i,j-1],alpha[i,j])) } } ## swap back tmp to res for the next inner iteration diff --git a/proto/BTCS.ipynb b/scripts/hannes-philip/BTCS.ipynb similarity index 100% rename from proto/BTCS.ipynb rename to scripts/hannes-philip/BTCS.ipynb diff --git a/proto/FTCS.ipynb b/scripts/hannes-philip/FTCS.ipynb similarity index 100% rename from proto/FTCS.ipynb rename to scripts/hannes-philip/FTCS.ipynb diff --git a/src/BTCS.cpp b/src/BTCS.cpp new file mode 100644 index 0000000..3b314e7 --- /dev/null +++ b/src/BTCS.cpp @@ -0,0 +1,418 @@ +/** + * @file BTCSv2.cpp + * @brief Implementation of heterogenous BTCS (backward time-centered space) + * solution of diffusion equation in 1D and 2D space. Internally the + * alternating-direction implicit (ADI) method is used. Version 2, because + * Version 1 was an implementation for the homogeneous BTCS solution. + * + */ + +#include "Schemes.hpp" +#include "TugUtils.hpp" + +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#else +#define omp_get_thread_num() 0 +#endif + +// calculates coefficient for boundary in constant case +template +constexpr std::pair calcBoundaryCoeffConstant(T alpha_center, + T alpha_side, T sx) { + const T centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) + + 2 * alpha_center); + const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side); + + return {centerCoeff, sideCoeff}; +} + +// calculates coefficient for boundary in closed case +template +constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, + T sx) { + const T centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side); + + const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side); + + return {centerCoeff, sideCoeff}; +} + +// creates coefficient matrix for next time step from alphas in x-direction +static Eigen::SparseMatrix +createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector &bcLeft, + std::vector &bcRight, int numCols, + int rowIndex, double sx) { + + // square matrix of column^2 dimension for the coefficients + Eigen::SparseMatrix cm(numCols, numCols); + cm.reserve(Eigen::VectorXi::Constant(numCols, 3)); + + // left column + + switch (bcLeft[rowIndex].getType()) { + case BC_TYPE_CONSTANT: { + auto [centerCoeffTop, rightCoeffTop] = + calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx); + cm.insert(0, 0) = centerCoeffTop; + cm.insert(0, 1) = rightCoeffTop; + break; + } + case BC_TYPE_CLOSED: { + auto [centerCoeffTop, rightCoeffTop] = + calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx); + cm.insert(0, 0) = centerCoeffTop; + cm.insert(0, 1) = rightCoeffTop; + break; + } + default: { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Left or Top!"); + } + } + + // inner columns + int n = numCols - 1; + for (int i = 1; i < n; i++) { + cm.insert(i, i - 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)); + cm.insert(i, i) = + 1 + + sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) + + calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i))); + cm.insert(i, i + 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)); + } + + // right column + + switch (bcRight[rowIndex].getType()) { + case BC_TYPE_CONSTANT: { + auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant( + alpha(rowIndex, n), alpha(rowIndex, n - 1), sx); + cm.insert(n, n - 1) = leftCoeffBottom; + cm.insert(n, n) = centerCoeffBottom; + break; + } + case BC_TYPE_CLOSED: { + auto [centerCoeffBottom, leftCoeffBottom] = + calcBoundaryCoeffClosed(alpha(rowIndex, n), alpha(rowIndex, n - 1), sx); + cm.insert(n, n - 1) = leftCoeffBottom; + cm.insert(n, n) = centerCoeffBottom; + break; + } + default: { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } + } + + cm.makeCompressed(); // important for Eigen solver + + return cm; +} + +// calculates explicit concentration at boundary in closed case +template +constexpr T calcExplicitConcentrationsBoundaryClosed(T conc_center, + T alpha_center, + T alpha_neigbor, T sy) { + return sy * calcAlphaIntercell(alpha_center, alpha_neigbor) * conc_center + + (1 - sy * (calcAlphaIntercell(alpha_center, alpha_neigbor))) * + conc_center; +} + +// calculates explicity concentration at boundary in constant case +template +constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, + T alpha_center, + T alpha_neighbor, T sy) { + return sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center + + (1 - sy * (calcAlphaIntercell(alpha_center, alpha_center) + + 2 * alpha_center)) * + conc_center + + sy * alpha_center * conc_bc; +} + +// creates a solution vector for next time step from the current state of +// concentrations +static Eigen::VectorXd createSolutionVector( + Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX, + Eigen::MatrixXd &alphaY, std::vector &bcLeft, + std::vector &bcRight, std::vector &bcTop, + std::vector &bcBottom, int length, int rowIndex, double sx, + double sy) { + + Eigen::VectorXd sv(length); + const std::size_t numRows = concentrations.rows(); + + // inner rows + if (rowIndex > 0 && rowIndex < numRows - 1) { + for (int i = 0; i < length; i++) { + sv(i) = + sy * + calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) * + concentrations(rowIndex + 1, i) + + (1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i), + alphaY(rowIndex + 1, i)) + + calcAlphaIntercell(alphaY(rowIndex - 1, i), + alphaY(rowIndex, i)))) * + concentrations(rowIndex, i) + + sy * + calcAlphaIntercell(alphaY(rowIndex - 1, i), alphaY(rowIndex, i)) * + concentrations(rowIndex - 1, i); + } + } + + // first row + else if (rowIndex == 0) { + for (int i = 0; i < length; i++) { + switch (bcTop[i].getType()) { + case BC_TYPE_CONSTANT: { + sv(i) = calcExplicitConcentrationsBoundaryConstant( + concentrations(rowIndex, i), bcTop[i].getValue(), + alphaY(rowIndex, i), alphaY(rowIndex + 1, i), sy); + break; + } + case BC_TYPE_CLOSED: { + sv(i) = calcExplicitConcentrationsBoundaryClosed( + concentrations(rowIndex, i), alphaY(rowIndex, i), + alphaY(rowIndex + 1, i), sy); + break; + } + default: + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Left or Top!"); + } + } + } + + // last row + else if (rowIndex == numRows - 1) { + for (int i = 0; i < length; i++) { + switch (bcBottom[i].getType()) { + case BC_TYPE_CONSTANT: { + sv(i) = calcExplicitConcentrationsBoundaryConstant( + concentrations(rowIndex, i), bcBottom[i].getValue(), + alphaY(rowIndex, i), alphaY(rowIndex - 1, i), sy); + break; + } + case BC_TYPE_CLOSED: { + sv(i) = calcExplicitConcentrationsBoundaryClosed( + concentrations(rowIndex, i), alphaY(rowIndex, i), + alphaY(rowIndex - 1, i), sy); + break; + } + default: + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } + } + } + + // first column -> additional fixed concentration change from perpendicular + // dimension in constant bc case + if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue(); + } + + // last column -> additional fixed concentration change from perpendicular + // dimension in constant bc case + if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(length - 1) += + 2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue(); + } + + return sv; +} + +// solver for linear equation system; A corresponds to coefficient matrix, +// b to the solution vector +// use of EigenLU solver +static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorXd &b) { + + Eigen::SparseLU> solver; + solver.analyzePattern(A); + solver.factorize(A); + + return solver.solve(b); +} + +// solver for linear equation system; A corresponds to coefficient matrix, +// b to the solution vector +// implementation of Thomas Algorithm +static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorXd &b) { + Eigen::Index n = b.size(); + + Eigen::VectorXd a_diag(n); + Eigen::VectorXd b_diag(n); + Eigen::VectorXd c_diag(n); + Eigen::VectorXd x_vec = b; + + // Fill diagonals vectors + b_diag[0] = A.coeff(0, 0); + c_diag[0] = A.coeff(0, 1); + + for (Eigen::Index i = 1; i < n - 1; i++) { + a_diag[i] = A.coeff(i, i - 1); + b_diag[i] = A.coeff(i, i); + c_diag[i] = A.coeff(i, i + 1); + } + a_diag[n - 1] = A.coeff(n - 1, n - 2); + b_diag[n - 1] = A.coeff(n - 1, n - 1); + + // start solving - c_diag and x_vec are overwritten + n--; + c_diag[0] /= b_diag[0]; + x_vec[0] /= b_diag[0]; + + for (Eigen::Index i = 1; i < n; i++) { + c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; + x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / + (b_diag[i] - a_diag[i] * c_diag[i - 1]); + } + + x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / + (b_diag[n] - a_diag[n] * c_diag[n - 1]); + + for (Eigen::Index i = n; i-- > 0;) { + x_vec[i] -= c_diag[i] * x_vec[i + 1]; + } + + return x_vec; +} + +// BTCS solution for 1D grid +static void +BTCS_1D(Grid &grid, Boundary &bc, double timestep, + Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorXd &b)) { + int length = grid.getLength(); + double sx = timestep / (grid.getDelta() * grid.getDelta()); + + Eigen::VectorXd concentrations_t1(length); + + Eigen::SparseMatrix A; + Eigen::VectorXd b(length); + + Eigen::MatrixXd alpha = grid.getAlpha(); + std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + + Eigen::MatrixXd concentrations = grid.getConcentrations(); + int rowIndex = 0; + A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, + sx); // this is exactly same as in 2D + for (int i = 0; i < length; i++) { + b(i) = concentrations(0, i); + } + if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { + b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue(); + } + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { + b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue(); + } + + concentrations_t1 = solverFunc(A, b); + + for (int j = 0; j < length; j++) { + concentrations(0, j) = concentrations_t1(j); + } + + grid.setConcentrations(concentrations); +} + +// BTCS solution for 2D grid +static void +BTCS_2D(Grid &grid, Boundary &bc, double timestep, + Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorXd &b), + int numThreads) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); + + Eigen::MatrixXd concentrations_t1 = + Eigen::MatrixXd::Constant(rowMax, colMax, 0); + Eigen::VectorXd row_t1(colMax); + + Eigen::SparseMatrix A; + Eigen::VectorXd b; + + Eigen::MatrixXd alphaX = grid.getAlphaX(); + Eigen::MatrixXd alphaY = grid.getAlphaY(); + std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + std::vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); + std::vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + + Eigen::MatrixXd concentrations = grid.getConcentrations(); + +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < rowMax; i++) { + + A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); + b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, + bcTop, bcBottom, colMax, i, sx, sy); + + Eigen::SparseLU> solver; + + row_t1 = solverFunc(A, b); + + concentrations_t1.row(i) = row_t1; + } + + concentrations_t1.transposeInPlace(); + concentrations.transposeInPlace(); + alphaX.transposeInPlace(); + alphaY.transposeInPlace(); + +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < colMax; i++) { + // swap alphas, boundary conditions and sx/sy for column-wise calculation + A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); + b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, + bcLeft, bcRight, rowMax, i, sy, sx); + + row_t1 = solverFunc(A, b); + + concentrations.row(i) = row_t1; + } + + concentrations.transposeInPlace(); + + grid.setConcentrations(concentrations); +} + +// entry point for EigenLU solver; differentiate between 1D and 2D grid +void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); + } else { + throw_invalid_argument( + "Error: Only 1- and 2-dimensional grids are defined!"); + } +} + +// entry point for Thomas algorithm solver; differentiate 1D and 2D grid +void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, ThomasAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); + } else { + throw_invalid_argument( + "Error: Only 1- and 2-dimensional grids are defined!"); + } +} diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp deleted file mode 100644 index 12b9844..0000000 --- a/src/BTCSv2.cpp +++ /dev/null @@ -1,487 +0,0 @@ -/** - * @file BTCSv2.cpp - * @brief Implementation of heterogenous BTCS (backward time-centered space) - * solution of diffusion equation in 1D and 2D space. Internally the - * alternating-direction implicit (ADI) method is used. Version 2, because - * Version 1 was an implementation for the homogeneous BTCS solution. - * - */ - -#include "FTCS.cpp" -#include -#include - -#ifdef WRITE_THOMAS_CSV -#include -#include -#endif - -#define NUM_THREADS_BTCS 10 - -using namespace Eigen; - -// calculates coefficient for left boundary in constant case -static tuple -calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; - - centerCoeff = - 1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) + - 2 * alpha(rowIndex, 0)); - rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - - return {centerCoeff, rightCoeff}; -} - -// calculates coefficient for left boundary in closed case -static tuple -calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; - - centerCoeff = - 1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - - return {centerCoeff, rightCoeff}; -} - -// calculates coefficient for right boundary in constant case -static tuple calcRightBoundaryCoeffConstant(MatrixXd &alpha, - int rowIndex, int n, - double sx) { - double leftCoeff; - double centerCoeff; - - leftCoeff = - -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); - centerCoeff = - 1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) + - 2 * alpha(rowIndex, n)); - - return {leftCoeff, centerCoeff}; -} - -// calculates coefficient for right boundary in closed case -static tuple -calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { - double leftCoeff; - double centerCoeff; - - leftCoeff = - -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); - centerCoeff = - 1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); - - return {leftCoeff, centerCoeff}; -} - -// creates coefficient matrix for next time step from alphas in x-direction -static SparseMatrix createCoeffMatrix(MatrixXd &alpha, - vector &bcLeft, - vector &bcRight, - int numCols, int rowIndex, - double sx) { - - // square matrix of column^2 dimension for the coefficients - SparseMatrix cm(numCols, numCols); - cm.reserve(VectorXi::Constant(numCols, 3)); - - // left column - BC_TYPE type = bcLeft[rowIndex].getType(); - if (type == BC_TYPE_CONSTANT) { - auto [centerCoeffTop, rightCoeffTop] = - calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); - cm.insert(0, 0) = centerCoeffTop; - cm.insert(0, 1) = rightCoeffTop; - } else if (type == BC_TYPE_CLOSED) { - auto [centerCoeffTop, rightCoeffTop] = - calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); - cm.insert(0, 0) = centerCoeffTop; - cm.insert(0, 1) = rightCoeffTop; - } else { - throw_invalid_argument( - "Undefined Boundary Condition Type somewhere on Left or Top!"); - } - - // inner columns - int n = numCols - 1; - for (int i = 1; i < n; i++) { - cm.insert(i, i - 1) = - -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)); - cm.insert(i, i) = - 1 + - sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) + - calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i))); - cm.insert(i, i + 1) = - -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)); - } - - // right column - type = bcRight[rowIndex].getType(); - if (type == BC_TYPE_CONSTANT) { - auto [leftCoeffBottom, centerCoeffBottom] = - calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); - cm.insert(n, n - 1) = leftCoeffBottom; - cm.insert(n, n) = centerCoeffBottom; - } else if (type == BC_TYPE_CLOSED) { - auto [leftCoeffBottom, centerCoeffBottom] = - calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); - cm.insert(n, n - 1) = leftCoeffBottom; - cm.insert(n, n) = centerCoeffBottom; - } else { - throw_invalid_argument( - "Undefined Boundary Condition Type somewhere on Right or Bottom!"); - } - - cm.makeCompressed(); // important for Eigen solver - - return cm; -} - -// calculates explicity concentration at top boundary in constant case -static double calcExplicitConcentrationsTopBoundaryConstant( - MatrixXd &concentrations, MatrixXd &alpha, vector &bcTop, - int rowIndex, int i, double sy) { - double c; - - c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * - concentrations(rowIndex, i) + - (1 - - sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) + - 2 * alpha(rowIndex, i))) * - concentrations(rowIndex, i) + - sy * alpha(rowIndex, i) * bcTop[i].getValue(); - - return c; -} - -// calculates explicit concentration at top boundary in closed case -static double calcExplicitConcentrationsTopBoundaryClosed( - MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { - double c; - - c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * - concentrations(rowIndex, i) + - (1 - - sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)))) * - concentrations(rowIndex, i); - - return c; -} - -// calculates explicit concentration at bottom boundary in constant case -static double calcExplicitConcentrationsBottomBoundaryConstant( - MatrixXd &concentrations, MatrixXd &alpha, - vector &bcBottom, int rowIndex, int i, double sy) { - double c; - - c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() + - (1 - - sy * (2 * alpha(rowIndex, i) + - calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) * - concentrations(rowIndex, i) + - sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) * - concentrations(rowIndex - 1, i); - - return c; -} - -// calculates explicit concentration at bottom boundary in closed case -static double calcExplicitConcentrationsBottomBoundaryClosed( - MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { - double c; - - c = (1 - - sy * (+calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) * - concentrations(rowIndex, i) + - sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) * - concentrations(rowIndex - 1, i); - - return c; -} - -// creates a solution vector for next time step from the current state of -// concentrations -static VectorXd createSolutionVector( - MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, - vector &bcLeft, vector &bcRight, - vector &bcTop, vector &bcBottom, - int length, int rowIndex, double sx, double sy) { - - VectorXd sv(length); - int numRows = concentrations.rows(); - BC_TYPE type; - - // inner rows - if (rowIndex > 0 && rowIndex < numRows - 1) { - for (int i = 0; i < length; i++) { - sv(i) = - sy * - calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) * - concentrations(rowIndex + 1, i) + - (1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i), - alphaY(rowIndex + 1, i)) + - calcAlphaIntercell(alphaY(rowIndex - 1, i), - alphaY(rowIndex, i)))) * - concentrations(rowIndex, i) + - sy * - calcAlphaIntercell(alphaY(rowIndex - 1, i), alphaY(rowIndex, i)) * - concentrations(rowIndex - 1, i); - } - } - - // first row - else if (rowIndex == 0) { - for (int i = 0; i < length; i++) { - type = bcTop[i].getType(); - if (type == BC_TYPE_CONSTANT) { - sv(i) = calcExplicitConcentrationsTopBoundaryConstant( - concentrations, alphaY, bcTop, rowIndex, i, sy); - } else if (type == BC_TYPE_CLOSED) { - sv(i) = calcExplicitConcentrationsTopBoundaryClosed( - concentrations, alphaY, rowIndex, i, sy); - } else { - throw_invalid_argument( - "Undefined Boundary Condition Type somewhere on Left or Top!"); - } - } - } - - // last row - else if (rowIndex == numRows - 1) { - for (int i = 0; i < length; i++) { - type = bcBottom[i].getType(); - if (type == BC_TYPE_CONSTANT) { - sv(i) = calcExplicitConcentrationsBottomBoundaryConstant( - concentrations, alphaY, bcBottom, rowIndex, i, sy); - } else if (type == BC_TYPE_CLOSED) { - sv(i) = calcExplicitConcentrationsBottomBoundaryClosed( - concentrations, alphaY, rowIndex, i, sy); - } else { - throw_invalid_argument( - "Undefined Boundary Condition Type somewhere on Right or Bottom!"); - } - } - } - - // first column -> additional fixed concentration change from perpendicular - // dimension in constant bc case - if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue(); - } - - // last column -> additional fixed concentration change from perpendicular - // dimension in constant bc case - if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(length - 1) += - 2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue(); - } - - return sv; -} - -// solver for linear equation system; A corresponds to coefficient matrix, -// b to the solution vector -// use of EigenLU solver -static VectorXd EigenLUAlgorithm(SparseMatrix &A, VectorXd &b) { - - SparseLU> solver; - solver.analyzePattern(A); - solver.factorize(A); - - return solver.solve(b); -} - -// solver for linear equation system; A corresponds to coefficient matrix, -// b to the solution vector -// implementation of Thomas Algorithm -static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { - uint32_t n = b.size(); - - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b; - - // Fill diagonals vectors - b_diag[0] = A.coeff(0, 0); - c_diag[0] = A.coeff(0, 1); - - for (int i = 1; i < n - 1; i++) { - a_diag[i] = A.coeff(i, i - 1); - b_diag[i] = A.coeff(i, i); - c_diag[i] = A.coeff(i, i + 1); - } - a_diag[n - 1] = A.coeff(n - 1, n - 2); - b_diag[n - 1] = A.coeff(n - 1, n - 1); - - // HACK: write CSV to file -#ifdef WRITE_THOMAS_CSV -#include -#include - static std::uint32_t file_index = 0; - std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; - - std::ofstream out_file; - - out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); - - // print header - out_file << "Aa, Ab, Ac, b\n"; - - // iterate through all elements - for (std::size_t i = 0; i < n; i++) { - out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " - << b[i] << "\n"; - } - - out_file.close(); -#endif - - // start solving - c_diag and x_vec are overwritten - n--; - c_diag[0] /= b_diag[0]; - x_vec[0] /= b_diag[0]; - - for (int i = 1; i < n; i++) { - c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; - x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / - (b_diag[i] - a_diag[i] * c_diag[i - 1]); - } - - x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / - (b_diag[n] - a_diag[n] * c_diag[n - 1]); - - for (int i = n; i-- > 0;) { - x_vec[i] -= c_diag[i] * x_vec[i + 1]; - } - - return x_vec; -} - -// BTCS solution for 1D grid -static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, - VectorXd (*solverFunc)(SparseMatrix &A, - VectorXd &b)) { - int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); - - VectorXd concentrations_t1(length); - - SparseMatrix A; - VectorXd b(length); - - MatrixXd alpha = grid.getAlpha(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - - MatrixXd concentrations = grid.getConcentrations(); - int rowIndex = 0; - A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, - sx); // this is exactly same as in 2D - for (int i = 0; i < length; i++) { - b(i) = concentrations(0, i); - } - if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { - b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue(); - } - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { - b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue(); - } - - concentrations_t1 = solverFunc(A, b); - - for (int j = 0; j < length; j++) { - concentrations(0, j) = concentrations_t1(j); - } - - grid.setConcentrations(concentrations); -} - -// BTCS solution for 2D grid -static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, - VectorXd (*solverFunc)(SparseMatrix &A, - VectorXd &b), - int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); - VectorXd row_t1(colMax); - - SparseMatrix A; - VectorXd b; - - MatrixXd alphaX = grid.getAlphaX(); - MatrixXd alphaY = grid.getAlphaY(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); - - MatrixXd concentrations = grid.getConcentrations(); - -#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < rowMax; i++) { - - A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); - b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, - bcTop, bcBottom, colMax, i, sx, sy); - - SparseLU> solver; - - row_t1 = solverFunc(A, b); - - concentrations_t1.row(i) = row_t1; - } - - concentrations_t1.transposeInPlace(); - concentrations.transposeInPlace(); - alphaX.transposeInPlace(); - alphaY.transposeInPlace(); - -#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < colMax; i++) { - // swap alphas, boundary conditions and sx/sy for column-wise calculation - A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); - b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, - bcLeft, bcRight, rowMax, i, sy, sx); - - row_t1 = solverFunc(A, b); - - concentrations.row(i) = row_t1; - } - - concentrations.transposeInPlace(); - - grid.setConcentrations(concentrations); -} - -// entry point for EigenLU solver; differentiate between 1D and 2D grid -static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { - if (grid.getDim() == 1) { - BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); - } else if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); - } else { - throw_invalid_argument( - "Error: Only 1- and 2-dimensional grids are defined!"); - } -} - -// entry point for Thomas algorithm solver; differentiate 1D and 2D grid -static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, - int numThreads) { - if (grid.getDim() == 1) { - BTCS_1D(grid, bc, timestep, ThomasAlgorithm); - } else if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); - } else { - throw_invalid_argument( - "Error: Only 1- and 2-dimensional grids are defined!"); - } -} diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 4f9890b..def9bbd 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,164 +1,149 @@ #include "TugUtils.hpp" -#include -#include -#include + +#include #include +#include -using namespace std; +BoundaryElement::BoundaryElement() {} -BoundaryElement::BoundaryElement() { - - this->type = BC_TYPE_CLOSED; - this->value = NAN; -} +BoundaryElement::BoundaryElement(double _value) + : value(_value), type(BC_TYPE_CONSTANT) {} -BoundaryElement::BoundaryElement(double value) { - this->type = BC_TYPE_CONSTANT; - this->value = value; -} - -void BoundaryElement::setType(BC_TYPE type) { - this->type = type; -} +void BoundaryElement::setType(BC_TYPE type) { this->type = type; } void BoundaryElement::setValue(double value) { - if(value < 0){ - throw_invalid_argument("No negative concentration allowed."); - } - if(type == BC_TYPE_CLOSED){ - throw_invalid_argument( - "No constant boundary concentrations can be set for closed " - "boundaries. Please change type first."); - } - this->value = value; + if (value < 0) { + throw_invalid_argument("No negative concentration allowed."); + } + if (type == BC_TYPE_CLOSED) { + throw_invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } + this->value = value; } -BC_TYPE BoundaryElement::getType() { - return this->type; -} +BC_TYPE BoundaryElement::getType() { return this->type; } -double BoundaryElement::getValue() { - return this->value; -} +double BoundaryElement::getValue() { return this->value; } Boundary::Boundary(Grid grid) : grid(grid) { - //probably to DEBUG assignment grid + if (grid.getDim() == 1) { + this->boundaries = std::vector>( + 2); // in 1D only left and right boundary - if (grid.getDim() == 1) { - this->boundaries = vector>(2); + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (grid.getDim() == 2) { + this->boundaries = std::vector>(4); - this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); - } else if (grid.getDim() == 2) { - this->boundaries = vector>(4); - - this->boundaries[BC_SIDE_LEFT] = vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT] = vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_TOP] = vector(grid.getCol(), BoundaryElement()); - this->boundaries[BC_SIDE_BOTTOM] = vector(grid.getCol(), BoundaryElement()); - } + this->boundaries[BC_SIDE_LEFT] = + std::vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = + std::vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = + std::vector(grid.getCol(), BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = + std::vector(grid.getCol(), BoundaryElement()); + } } void Boundary::setBoundarySideClosed(BC_SIDE side) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } + } - int n; - if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { - n = grid.getRow(); - } else { - n = grid.getCol(); - } - this->boundaries[side] = vector(n, BoundaryElement()); + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? grid.getRow() : grid.getCol(); + + this->boundaries[side] = std::vector(n, BoundaryElement()); } void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } + } - int n; - if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { - n = grid.getRow(); - } else { - n = grid.getCol(); - } - this->boundaries[side] = vector(n, BoundaryElement(value)); + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? grid.getRow() : grid.getCol(); + + this->boundaries[side] = + std::vector(n, BoundaryElement(value)); } void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { - // tests whether the index really points to an element of the boundary side. - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CLOSED); + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CLOSED); } -void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) { - // tests whether the index really points to an element of the boundary side. - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CONSTANT); - this->boundaries[side][index].setValue(value); +void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, + double value) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); } -vector Boundary::getBoundarySide(BC_SIDE side) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } +const std::vector Boundary::getBoundarySide(BC_SIDE side) { + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } - return this->boundaries[side]; + } + return this->boundaries[side]; } -VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { - int length = boundaries[side].size(); - VectorXd values(length); +Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { + const std::size_t length = boundaries[side].size(); + Eigen::VectorXd values(length); - for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { - values(i) = -1; - continue; - } - values(i) = getBoundaryElementValue(side, i); + for (int i = 0; i < length; i++) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { + values(i) = -1; + continue; } + values(i) = getBoundaryElementValue(side, i); + } - return values; + return values; } BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index]; + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index]; } BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index].getType(); + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index].getType(); } double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - if(boundaries[side][index].getType() != BC_TYPE_CONSTANT){ - throw_invalid_argument("A value can only be output if it is a constant boundary condition."); - } - return this->boundaries[side][index].getValue(); + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { + throw_invalid_argument( + "A value can only be output if it is a constant boundary condition."); + } + return this->boundaries[side][index].getValue(); } - diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 04db461..5527d09 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp) +add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp) IF(TUG_NAAICE_EXAMPLE) target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) @@ -10,4 +10,4 @@ if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) target_link_libraries(tug OpenMP::OpenMP_CXX) endif() -target_include_directories(tug PUBLIC ../include) +target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 449a1c1..375fe50 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -1,435 +1,381 @@ /** * @file FTCS.cpp - * @brief Implementation of heterogenous FTCS (forward time-centered space) solution - * of diffusion equation in 1D and 2D space. - * + * @brief Implementation of heterogenous FTCS (forward time-centered space) + * solution of diffusion equation in 1D and 2D space. + * */ +#include "Schemes.hpp" #include "TugUtils.hpp" + #include -#include #include +#include + +#ifdef _OPENMP #include - -#define NUM_THREADS 6 - -using namespace std; - +#else +#define omp_get_thread_num() 0 +#endif // calculates horizontal change on one cell independent of boundary type -static double calcHorizontalChange(Grid &grid, int &row, int &col) { +static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { - double result = - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - ( - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1); - - return result; + return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col + 1) - + (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col))) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col - 1); } - // calculates vertical change on one cell independent of boundary type -static double calcVerticalChange(Grid &grid, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row+1,col) - - ( - calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row-1,col); +static inline double calcVerticalChange(Grid &grid, int &row, int &col) { - return result; + return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row + 1, col) - + (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) + + calcAlphaIntercell(grid.getAlphaY()(row - 1, col), + grid.getAlphaY()(row, col))) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaY()(row - 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row - 1, col); } - // calculates horizontal change on one cell with a constant left boundary -static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, + int &col) { - double result = - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - ( - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + 2 * grid.getAlphaX()(row,col) - ) - * grid.getConcentrations()(row,col) - + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row); - - return result; + return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col + 1) - + (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) + + 2 * grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col) + + 2 * grid.getAlphaX()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_LEFT, row); } - // calculates horizontal change on one cell with a closed left boundary -static double calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) - * (grid.getConcentrations()(row, col+1) - grid.getConcentrations()(row, col)); - - return result; -} +static inline double +calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { + return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) * + (grid.getConcentrations()(row, col + 1) - + grid.getConcentrations()(row, col)); +} // checks boundary condition type for a cell on the left edge of grid -static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, int &row, int &col) { - if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { - return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) { - return calcHorizontalChangeLeftBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } +static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { + if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { + return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); + } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) { + return calcHorizontalChangeLeftBoundaryClosed(grid, row, col); + } else { + throw_invalid_argument("Undefined Boundary Condition Type!"); + } } - // calculates horizontal change on one cell with a constant right boundary -static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, + int &col) { - double result = - 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - - ( - calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - + 2 * grid.getAlphaX()(row,col) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1); - - return result; + return 2 * grid.getAlphaX()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - + (calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) + + 2 * grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col - 1); } - // calculates horizontal change on one cell with a closed right boundary -static double calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { - - double result = - - (calcAlphaIntercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col)) - * (grid.getConcentrations()(row, col) - grid.getConcentrations()(row, col-1))); - - return result; -} +static inline double +calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { + return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) * + (grid.getConcentrations()(row, col) - + grid.getConcentrations()(row, col - 1))); +} // checks boundary condition type for a cell on the right edge of grid -static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, int &row, int &col) { - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { - return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) { - return calcHorizontalChangeRightBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } +static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { + return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); + } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) { + return calcHorizontalChangeRightBoundaryClosed(grid, row, col); + } else { + throw_invalid_argument("Undefined Boundary Condition Type!"); + } } - // calculates vertical change on one cell with a constant top boundary -static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) - * grid.getConcentrations()(row+1,col) - - ( - calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) - + 2 * grid.getAlphaY()(row, col) - ) - * grid.getConcentrations()(row, col) - + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col); +static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { - return result; + return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row + 1, col) - + (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) + + 2 * grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row, col) + + 2 * grid.getAlphaY()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_TOP, col); } - // calculates vertical change on one cell with a closed top boundary -static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getConcentrations()(row, col)) - * (grid.getConcentrations()(row+1, col) - grid.getConcentrations()(row, col)); +static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, + int &col) { - return result; + return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getConcentrations()(row, col)) * + (grid.getConcentrations()(row + 1, col) - + grid.getConcentrations()(row, col)); } - // checks boundary condition type for a cell on the top edge of grid -static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row, int &col) { - if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { - return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { - return calcVerticalChangeTopBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } +static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { + if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { + return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); + } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { + return calcVerticalChangeTopBoundaryClosed(grid, row, col); + } else { + throw_invalid_argument("Undefined Boundary Condition Type!"); + } } - // calculates vertical change on one cell with a constant bottom boundary -static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, + int &col) { - double result = - 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - - ( - calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) - + 2 * grid.getAlphaY()(row, col) - ) - * grid.getConcentrations()(row, col) - + calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) - * grid.getConcentrations()(row-1,col); - - return result; + return 2 * grid.getAlphaY()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - + (calcAlphaIntercell(grid.getAlphaY()(row, col), + grid.getAlphaY()(row - 1, col)) + + 2 * grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaY()(row, col), + grid.getAlphaY()(row - 1, col)) * + grid.getConcentrations()(row - 1, col); } - // calculates vertical change on one cell with a closed bottom boundary -static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { +static inline double +calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { - double result = - - (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) - * (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col))); - - return result; + return -(calcAlphaIntercell(grid.getAlphaY()(row, col), + grid.getAlphaY()(row - 1, col)) * + (grid.getConcentrations()(row, col) - + grid.getConcentrations()(row - 1, col))); } - // checks boundary condition type for a cell on the bottom edge of grid -static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, int &row, int &col) { - if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { - return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); - } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) { - return calcVerticalChangeBottomBoundaryClosed(grid, row, col); - } else { - throw_invalid_argument("Undefined Boundary Condition Type!"); - } +static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { + if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { + return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); + } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) { + return calcVerticalChangeBottomBoundaryClosed(grid, row, col); + } else { + throw_invalid_argument("Undefined Boundary Condition Type!"); + } } - // FTCS solution for 1D grid static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { - int colMax = grid.getCol(); - double deltaCol = grid.getDeltaCol(); + int colMax = grid.getCol(); + double deltaCol = grid.getDeltaCol(); - // matrix for concentrations at time t+1 - MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); + // matrix for concentrations at time t+1 + Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0); - // only one row in 1D case -> row constant at index 0 - int row = 0; + // only one row in 1D case -> row constant at index 0 + int row = 0; - // inner cells - // independent of boundary condition type - for (int col = 1; col < colMax-1; col++) { - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChange(grid, row, col) - ) - ; - } + // inner cells + // independent of boundary condition type + for (int col = 1; col < colMax - 1; col++) { + concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChange(grid, row, col)); + } - // left boundary; hold column constant at index 0 - int col = 0; - concentrations_t1(row, col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeLeftBoundary(grid, bc, row, col) - ) - ; + // left boundary; hold column constant at index 0 + int col = 0; + concentrations_t1(row, col) = + grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChangeLeftBoundary(grid, bc, row, col)); + // right boundary; hold column constant at max index + col = colMax - 1; + concentrations_t1(row, col) = + grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChangeRightBoundary(grid, bc, row, col)); - // right boundary; hold column constant at max index - col = colMax-1; - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeRightBoundary(grid, bc, row, col) - ) - ; - - // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); + // overwrite obsolete concentrations + grid.setConcentrations(concentrations_t1); } - // FTCS solution for 2D grid -static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - double deltaRow = grid.getDeltaRow(); - double deltaCol = grid.getDeltaCol(); +static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, + int numThreads) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double deltaRow = grid.getDeltaRow(); + double deltaCol = grid.getDeltaCol(); - // matrix for concentrations at time t+1 - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); - - // inner cells - // these are independent of the boundary condition type - // omp_set_num_threads(10); -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - for (int col = 1; col < colMax-1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChange(grid, row, col) - ) - ; - } - } - - // boundary conditions - // left without corners / looping over rows - // hold column constant at index 0 - int col = 0; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row, col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeLeftBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - ; - } - - // right without corners / looping over rows - // hold column constant at max index - col = colMax-1; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeRightBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - ; - } - - - // top without corners / looping over columns - // hold row constant at index 0 - int row = 0; -#pragma omp parallel for num_threads(numThreads) - for (int col=1; col + #include +#include -Grid::Grid(int length) { - if (length <= 3) { - throw_invalid_argument("Given grid length too small. Must be greater than 3."); - } +constexpr double MAT_INIT_VAL = 0; - this->row = 1; - this->col = length; - this->domainCol = length; // default: same size as length - this->deltaCol = double(this->domainCol)/double(this->col); // -> 1 - this->dim = 1; +Grid::Grid(int length) : col(length), domainCol(length) { + if (length <= 3) { + throw_invalid_argument( + "Given grid length too small. Must be greater than 3."); + } - this->concentrations = MatrixXd::Constant(1, col, 20); - this->alphaX = MatrixXd::Constant(1, col, 1); + this->dim = 1; + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + + this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); } -Grid::Grid(int row, int col) { - if (row <= 3 || col <= 3) { - throw_invalid_argument("Given grid dimensions too small. Must each be greater than 3."); - } +Grid::Grid(int _row, int _col) + : row(_row), col(_col), domainRow(_row), domainCol(_col) { + if (row <= 3 || col <= 3) { + throw_invalid_argument( + "Given grid dimensions too small. Must each be greater than 3."); + } - this->row = row; - this->col = col; - this->domainRow = row; // default: same size as row - this->domainCol = col; // default: same size as col - this->deltaRow = double(this->domainRow)/double(this->row); // -> 1 - this->deltaCol = double(this->domainCol)/double(this->col); // -> 1 - this->dim = 2; + this->dim = 2; + this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - this->concentrations = MatrixXd::Constant(row, col, 20); - this->alphaX = MatrixXd::Constant(row, col, 1); - this->alphaY = MatrixXd::Constant(row, col, 1); + this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); } -void Grid::setConcentrations(MatrixXd concentrations) { - if (concentrations.rows() != this->row || concentrations.cols() != this->col) { - throw_invalid_argument("Given matrix of concentrations mismatch with Grid dimensions!"); - } +void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) { + if (concentrations.rows() != this->row || + concentrations.cols() != this->col) { + throw_invalid_argument( + "Given matrix of concentrations mismatch with Grid dimensions!"); + } - this->concentrations = concentrations; + this->concentrations = concentrations; } -MatrixXd Grid::getConcentrations() { - return this->concentrations; +void Grid::setAlpha(const Eigen::MatrixXd &alpha) { + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + if (alpha.rows() != 1 || alpha.cols() != this->col) { + throw_invalid_argument( + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); + } + + this->alphaX = alpha; } -void Grid::setAlpha(MatrixXd alpha) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } +void Grid::setAlpha(const Eigen::MatrixXd &alphaX, + const Eigen::MatrixXd &alphaY) { + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, you should probably " + "use 1D setter function!"); + } + if (alphaX.rows() != this->row || alphaX.cols() != this->col) { + throw_invalid_argument("Given matrix of alpha coefficients in x-direction " + "mismatch with GRid dimensions!"); + } + if (alphaY.rows() != this->row || alphaY.cols() != this->col) { + throw_invalid_argument("Given matrix of alpha coefficients in y-direction " + "mismatch with GRid dimensions!"); + } - this->alphaX = alpha; + this->alphaX = alphaX; + this->alphaY = alphaY; } -void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in x-direction mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in y-direction mismatch with GRid dimensions!"); - } +const Eigen::MatrixXd &Grid::getAlpha() { + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use either getAlphaX() or getAlphaY()!"); + } - this->alphaX = alphaX; - this->alphaY = alphaY; + return this->alphaX; } -MatrixXd Grid::getAlpha() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use either getAlphaX() or getAlphaY()!"); - } +const Eigen::MatrixXd &Grid::getAlphaX() { + if (dim != 2) { + throw_invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } - return this->alphaX; + return this->alphaX; } -MatrixXd Grid::getAlphaX() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!"); - } +const Eigen::MatrixXd &Grid::getAlphaY() { + if (dim != 2) { + throw_invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } - return this->alphaX; + return this->alphaY; } -MatrixXd Grid::getAlphaY() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaY; -} - -int Grid::getDim() { - return dim; -} +int Grid::getDim() { return dim; } int Grid::getLength() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use getRow() or getCol()!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use getRow() or getCol()!"); + } - return col; + return col; } -int Grid::getRow() { - return row; -} +int Grid::getRow() { return row; } -int Grid::getCol() { - return col; -} +int Grid::getCol() { return col; } void Grid::setDomain(double domainLength) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probaly use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw_invalid_argument("Given domain length is not positive!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probaly " + "use the 2D domain setter!"); + } + if (domainLength <= 0) { + throw_invalid_argument("Given domain length is not positive!"); + } - this->domainCol = domainLength; - this->deltaCol = double(this->domainCol)/double(this->col); + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->col); } void Grid::setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw_invalid_argument("Given domain size is not positive!"); - } + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, you should probably " + "use the 1D domain setter!"); + } + if (domainRow <= 0 || domainCol <= 0) { + throw_invalid_argument("Given domain size is not positive!"); + } - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = double(this->domainRow)/double(this->row); - this->deltaCol = double(this->domainCol)/double(this->col); + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = double(this->domainRow) / double(this->row); + this->deltaCol = double(this->domainCol) / double(this->col); } double Grid::getDelta() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use the 2D delta getters"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use the 2D delta getters"); + } - return this->deltaCol; + return this->deltaCol; } -double Grid::getDeltaCol() { - return this->deltaCol; -} +double Grid::getDeltaCol() { return this->deltaCol; } double Grid::getDeltaRow() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, meaning there is no delta in y-direction!"); - } + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, meaning there is no " + "delta in y-direction!"); + } - return this->deltaRow; + return this->deltaRow; } diff --git a/src/Schemes.hpp b/src/Schemes.hpp new file mode 100644 index 0000000..d10e0c5 --- /dev/null +++ b/src/Schemes.hpp @@ -0,0 +1,28 @@ +/** + * @file BTCSv2.cpp + * @brief Implementation of heterogenous BTCS (backward time-centered space) + * solution of diffusion equation in 1D and 2D space. Internally the + * alternating-direction implicit (ADI) method is used. Version 2, because + * Version 1 was an implementation for the homogeneous BTCS solution. + * + */ + +#ifndef SCHEMES_H_ +#define SCHEMES_H_ + +#include "TugUtils.hpp" + +#include +#include + +// entry point; differentiate between 1D and 2D grid +extern void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads); + +// entry point for EigenLU solver; differentiate between 1D and 2D grid +extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads); + +// entry point for Thomas algorithm solver; differentiate 1D and 2D grid +extern void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, + int numThreads); + +#endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 6e7ae76..0a5a961 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -1,305 +1,326 @@ #include #include #include +#include +#include +#include #include #include + #include -#include -#include +#include "Schemes.hpp" +#include "TugUtils.hpp" -#ifndef SIMULATION_H_ -#define SIMULATION_H_ - -#include "BTCSv2.cpp" - - -using namespace std; - -Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { - - this->approach = approach; - this->solver = THOMAS_ALGORITHM_SOLVER; - this->timestep = -1; // error per default - this->iterations = -1; - this->innerIterations = 1; - this->numThreads = omp_get_num_procs(); - - this->csv_output = CSV_OUTPUT_OFF; - this->console_output = CONSOLE_OUTPUT_OFF; - this->time_measure = TIME_MEASURE_OFF; -} +Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) + : grid(_grid), bc(_bc), approach(_approach) {} void Simulation::setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid CSV output option given!"); - } + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw_invalid_argument("Invalid CSV output option given!"); + } - this->csv_output = csv_output; + this->csv_output = csv_output; } void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { - if (console_output < CONSOLE_OUTPUT_OFF && console_output > CONSOLE_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid console output option given!"); - } + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw_invalid_argument("Invalid console output option given!"); + } - this->console_output = console_output; + this->console_output = console_output; } void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { - if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { - throw_invalid_argument("Invalid time measure option given!"); - } + if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { + throw_invalid_argument("Invalid time measure option given!"); + } - this->time_measure = time_measure; + this->time_measure = time_measure; } void Simulation::setTimestep(double timestep) { - if(timestep <= 0){ - throw_invalid_argument("Timestep has to be greater than zero."); - } + if (timestep <= 0) { + throw_invalid_argument("Timestep has to be greater than zero."); + } - if (approach == FTCS_APPROACH) { + if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - double deltaRowSquare; - double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - double minDeltaSquare; - double maxAlphaX, maxAlphaY, maxAlpha; - string dim; - if (grid.getDim() == 2) { - dim = "2D"; + const double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + // will be 0 if 1D, else ... + const double deltaRowSquare = grid.getDim() != 1 + ? grid.getDeltaRow() * grid.getDeltaRow() + : deltaColSquare; + const double minDeltaSquare = + (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + double maxAlpha = std::numeric_limits::quiet_NaN(); - minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - maxAlphaX = grid.getAlphaX().maxCoeff(); - maxAlphaY = grid.getAlphaY().maxCoeff(); - maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + // determine maximum alpha + if (grid.getDim() == 2) { - } else if (grid.getDim() == 1) { - dim = "1D"; - minDeltaSquare = deltaColSquare; - maxAlpha = grid.getAlpha().maxCoeff(); - - } else { - throw_invalid_argument("Critical error: Undefined number of dimensions!"); - } + const double maxAlphaX = grid.getAlphaX().maxCoeff(); + const double maxAlphaY = grid.getAlphaY().maxCoeff(); + maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - // Courant-Friedrichs-Lewy condition - double cfl = minDeltaSquare / (4*maxAlpha); - - // stability equation from Wikipedia; might be useful if applied cfl does not work in some cases - // double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); - - cout << "FTCS_" << dim << " :: CFL condition MDL: " << cfl << endl; - cout << "FTCS_" << dim << " :: required dt=" << timestep << endl; - - if (timestep > cfl) { - - this->innerIterations = (int)ceil(timestep / cfl); - this->timestep = timestep / (double)innerIterations; - - cerr << "Warning :: Timestep was adjusted, because of stability " - "conditions. Time duration was approximately preserved by " - "adjusting internal number of iterations." - << endl; - cout << "FTCS_" << dim << " :: Required " << this->innerIterations - << " inner iterations with dt=" << this->timestep << endl; - - } else { - - this->timestep = timestep; - cout << "FTCS_" << dim << " :: No inner iterations required, dt=" << timestep << endl; - - } + } else if (grid.getDim() == 1) { + maxAlpha = grid.getAlpha().maxCoeff(); } else { - this->timestep = timestep; + throw_invalid_argument("Critical error: Undefined number of dimensions!"); } + const std::string dim = std::to_string(grid.getDim()) + "D"; + + // Courant-Friedrichs-Lewy condition + double cfl = minDeltaSquare / (4 * maxAlpha); + + // stability equation from Wikipedia; might be useful if applied cfl does + // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * + // ((1/deltaRowSquare) + (1/deltaColSquare))); + + const std::string &approachPrefix = this->approach_names[approach]; + std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl + << std::endl; + std::cout << approachPrefix << "_" << dim << " :: required dt=" << timestep + << std::endl; + + if (timestep > cfl) { + + this->innerIterations = (int)ceil(timestep / cfl); + this->timestep = timestep / (double)innerIterations; + + std::cerr << "Warning :: Timestep was adjusted, because of stability " + "conditions. Time duration was approximately preserved by " + "adjusting internal number of iterations." + << std::endl; + std::cout << approachPrefix << "_" << dim << " :: Required " + << this->innerIterations + << " inner iterations with dt=" << this->timestep << std::endl; + + } else { + + this->timestep = timestep; + std::cout << approachPrefix << "_" << dim + << " :: No inner iterations required, dt=" << timestep + << std::endl; + } + + } else { + this->timestep = timestep; + } } -double Simulation::getTimestep() { - return this->timestep; -} +double Simulation::getTimestep() { return this->timestep; } void Simulation::setIterations(int iterations) { - if(iterations <= 0){ - throw_invalid_argument("Number of iterations must be greater than zero."); - } - this->iterations = iterations; + if (iterations <= 0) { + throw_invalid_argument("Number of iterations must be greater than zero."); + } + this->iterations = iterations; } void Simulation::setSolver(SOLVER solver) { - if (this->approach == FTCS_APPROACH) { - cerr << "Warning: Solver was set, but FTCS approach initialized. Setting the solver " - "is thus without effect." - << endl; - } + if (this->approach == FTCS_APPROACH) { + std::cerr + << "Warning: Solver was set, but FTCS approach initialized. Setting " + "the solver " + "is thus without effect." + << std::endl; + } - this->solver = solver; + this->solver = solver; } -void Simulation::setNumberThreads(int numThreads){ - if(numThreads>0 && numThreads<=omp_get_num_procs()){ - this->numThreads=numThreads; - } - else{ - int maxThreadNumber = omp_get_num_procs(); - - throw_invalid_argument("Number of threads exceeds the number of processor cores or is less than 1."); - } +void Simulation::setNumberThreads(int numThreads) { + if (numThreads > 0 && numThreads <= omp_get_num_procs()) { + this->numThreads = numThreads; + } else { + int maxThreadNumber = omp_get_num_procs(); + std::string outputMessage = + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ") or is less than 1."; + + throw_invalid_argument(outputMessage); + } } -int Simulation::getIterations() { - return this->iterations; -} +int Simulation::getIterations() { return this->iterations; } void Simulation::printConcentrationsConsole() { - cout << grid.getConcentrations() << endl; - cout << endl; + std::cout << grid.getConcentrations() << std::endl; + std::cout << std::endl; } -string Simulation::createCSVfile() { - ofstream file; - int appendIdent = 0; - string appendIdentString; +std::string Simulation::createCSVfile() { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; - string approachString = (approach == 0) ? "FTCS" : "BTCS"; - string row = to_string(grid.getRow()); - string col = to_string(grid.getCol()); - string numIterations = to_string(iterations); + // string approachString = (approach == 0) ? "FTCS" : "BTCS"; + const std::string &approachString = this->approach_names[approach]; + std::string row = std::to_string(grid.getRow()); + std::string col = std::to_string(grid.getCol()); + std::string numIterations = std::to_string(iterations); - string filename = approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; + std::string filename = + approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; - while (filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = to_string(appendIdent); - filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv"; - } + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = approachString + "_" + row + "_" + col + "_" + numIterations + + "-" + appendIdentString + ".csv"; + } - file.open(filename); - if (!file) { - exit(1); - } + file.open(filename); + if (!file) { + exit(1); + } - // adds lines at the beginning of verbose output csv that represent the boundary conditions and their values - // -1 in case of closed boundary - if (csv_output == CSV_OUTPUT_XTREME) { - IOFormat one_row(StreamPrecision, DontAlignCols, "", " "); - file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) << endl; // boundary left - file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) << endl; // boundary right - file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) << endl; // boundary top - file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) << endl; // boundary bottom - file << endl << endl; - } + // adds lines at the beginning of verbose output csv that represent the + // boundary conditions and their values -1 in case of closed boundary + if (csv_output == CSV_OUTPUT_XTREME) { + Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", + " "); + file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) + << std::endl; // boundary left + file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) + << std::endl; // boundary right + file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) + << std::endl; // boundary top + file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) + << std::endl; // boundary bottom + file << std::endl << std::endl; + } - file.close(); + file.close(); - return filename; + return filename; } -void Simulation::printConcentrationsCSV(string filename) { - ofstream file; +void Simulation::printConcentrationsCSV(const std::string &filename) { + std::ofstream file; - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } - IOFormat do_not_align(StreamPrecision, DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << endl; - file << endl << endl; - file.close(); + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); } void Simulation::run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } + if (this->timestep == -1) { + throw_invalid_argument("Timestep is not set!"); + } + if (this->iterations == -1) { + throw_invalid_argument("Number of iterations are not set!"); + } - string filename; - if (this->console_output > CONSOLE_OUTPUT_OFF) { + std::string filename; + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + filename = createCSVfile(); + } + + auto begin = std::chrono::high_resolution_clock::now(); + + if (approach == FTCS_APPROACH) { // FTCS case + + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - filename = createCSVfile(); - } - - auto begin = std::chrono::high_resolution_clock::now(); - - if (approach == FTCS_APPROACH) { // FTCS case - - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - - // if (i % (iterations * innerIterations / 100) == 0) { - // double percentage = (double)i / ((double)iterations * (double)innerIterations) * 100; - // if ((int)percentage % 10 == 0) { - // cout << "Progress: " << percentage << "%" << endl; - // } - // } - } - - } else if (approach == BTCS_APPROACH) { // BTCS case - - if (solver == EIGEN_LU_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); - - } - } else if (solver == THOMAS_ALGORITHM_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - - } - } - - } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - - // TODO - - } - - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = std::chrono::duration_cast(end - begin); - - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { printConcentrationsCSV(filename); - } - if (this->time_measure > TIME_MEASURE_OFF) { - string approachString = (approach == 0) ? "FTCS" : "BTCS"; - string dimString = (grid.getDim() == 1) ? "-1D" : "-2D"; - cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << endl; - } - -} + } -#endif + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + + // if (i % (iterations * innerIterations / 100) == 0) { + // double percentage = (double)i / ((double)iterations * + // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { + // cout << "Progress: " << percentage << "%" << endl; + // } + // } + } + + } else if (approach == BTCS_APPROACH) { // BTCS case + + if (solver == EIGEN_LU_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); + } + } else if (solver == THOMAS_ALGORITHM_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + } + } + + } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case + + double beta = 0.5; + + // TODO this implementation is very inefficient! + // a separate implementation that sets up a specific tridiagonal matrix for + // Crank-Nicolson would be better + Eigen::MatrixXd concentrations; + Eigen::MatrixXd concentrationsFTCS; + Eigen::MatrixXd concentrationsResult; + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + concentrations = grid.getConcentrations(); + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsFTCS = grid.getConcentrations(); + grid.setConcentrations(concentrations); + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsResult = + beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations(); + grid.setConcentrations(concentrationsResult); + } + } + + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - begin); + + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + printConcentrationsCSV(filename); + } + if (this->time_measure > TIME_MEASURE_OFF) { + const std::string &approachString = this->approach_names[approach]; + const std::string dimString = std::to_string(grid.getDim()) + "D"; + std::cout << approachString << dimString << ":: run() finished in " + << milliseconds.count() << "ms" << std::endl; + } +} diff --git a/src/TugUtils.hpp b/src/TugUtils.hpp index 6b56f06..e6a2741 100644 --- a/src/TugUtils.hpp +++ b/src/TugUtils.hpp @@ -1,5 +1,5 @@ -#ifndef BTCSUTILS_H_ -#define BTCSUTILS_H_ +#ifndef TUGUTILS_H_ +#define TUGUTILS_H_ #include #include @@ -23,13 +23,15 @@ start); \ duration.count(); \ }) -#endif // BTCSUTILS_H_ // calculates arithmetic or harmonic mean of alpha between two cells -static double calcAlphaIntercell(double &alpha1, double &alpha2, bool useHarmonic = true) { - if (useHarmonic) { - return double(2) / ((double(1)/alpha1) + (double(1)/alpha2)); - } else { - return 0.5 * (alpha1 + alpha2); - } -} \ No newline at end of file +template +constexpr T calcAlphaIntercell(T alpha1, T alpha2, + bool useHarmonic = true) { + if (useHarmonic) { + return double(2) / ((double(1) / alpha1) + (double(1) / alpha2)); + } else { + return 0.5 * (alpha1 + alpha2); + } +} +#endif // TUGUTILS_H_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5210b95..19c4daf 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,15 +11,15 @@ if(NOT DOCTEST_LIB) FetchContent_MakeAvailable(DocTest) endif() -add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp) +add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) target_link_libraries(testTug doctest tug) # get relative path of the CSV file -get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH CACHE) +get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH) # set relative path in header file configure_file(testSimulation.hpp.in testSimulation.hpp) # include test directory with generated header file from above -target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}") +target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src") add_custom_target( check diff --git a/test/TestUtils.cpp b/test/TestUtils.cpp deleted file mode 100644 index ef0d3db..0000000 --- a/test/TestUtils.cpp +++ /dev/null @@ -1,47 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -using namespace std; -using namespace Eigen; - -MatrixXd CSV2Eigen(string file2Convert) { - - vector matrixEntries; - - ifstream matrixDataFile(file2Convert); - if (matrixDataFile.fail()) { - throw invalid_argument("File probably non-existent!"); - } - - string matrixRowString; - string matrixEntry; - int matrixRowNumber = 0; - - while(getline(matrixDataFile, matrixRowString)){ - stringstream matrixRowStringStream(matrixRowString); - while(getline(matrixRowStringStream, matrixEntry, ' ')){ - matrixEntries.push_back(stod(matrixEntry)); - } - if (matrixRowString.length() > 1) { - matrixRowNumber++; - } - } - - return Map>(matrixEntries.data(), matrixRowNumber, matrixEntries.size() / matrixRowNumber); -} - -bool checkSimilarity(MatrixXd a, MatrixXd b, double precision=1e-5) { - return a.isApprox(b, precision); -} - -bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) { - - MatrixXd diff = a - b; - double maxCoeff = diff.maxCoeff(); - return abs(maxCoeff) < maxDiff; -} \ No newline at end of file diff --git a/test/TestUtils.hpp b/test/TestUtils.hpp new file mode 100644 index 0000000..e6375de --- /dev/null +++ b/test/TestUtils.hpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include +#include +#include + +inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) { + + std::vector matrixEntries; + + std::ifstream matrixDataFile(file2Convert); + if (matrixDataFile.fail()) { + throw std::invalid_argument("File probably non-existent!"); + } + + std::string matrixRowString; + std::string matrixEntry; + int matrixRowNumber = 0; + + while (getline(matrixDataFile, matrixRowString)) { + std::stringstream matrixRowStringStream(matrixRowString); + while (getline(matrixRowStringStream, matrixEntry, ' ')) { + matrixEntries.push_back(stod(matrixEntry)); + } + if (matrixRowString.length() > 1) { + matrixRowNumber++; + } + } + + return Eigen::Map< + Eigen::Matrix>( + matrixEntries.data(), matrixRowNumber, + matrixEntries.size() / matrixRowNumber); +} + +inline bool checkSimilarity(Eigen::MatrixXd a, Eigen::MatrixXd b, + double precision = 1e-5) { + return a.isApprox(b, precision); +} + +inline bool checkSimilarityV2(Eigen::MatrixXd a, Eigen::MatrixXd b, + double maxDiff) { + + Eigen::MatrixXd diff = a - b; + double maxCoeff = diff.maxCoeff(); + return abs(maxCoeff) < maxDiff; +} diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 47ca1fa..1d2bff2 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -1,68 +1,76 @@ -#include #include -#include -#include -#include #include +#include +#include +#include +#include -TEST_CASE("BoundaryElement"){ +using namespace std; - SUBCASE("Closed case"){ - BoundaryElement boundaryElementClosed = BoundaryElement(); - CHECK_NOTHROW(BoundaryElement()); - CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); - CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN)); - CHECK_THROWS(boundaryElementClosed.setValue(0.2)); - } +TEST_CASE("BoundaryElement") { - SUBCASE("Constant case"){ - BoundaryElement boundaryElementConstant = BoundaryElement(0.1); - CHECK_NOTHROW(BoundaryElement(0.1)); - CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); - CHECK_EQ(boundaryElementConstant.getValue(), 0.1); - CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); - CHECK_EQ(boundaryElementConstant.getValue(), 0.2); - } + SUBCASE("Closed case") { + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); + CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); + CHECK_EQ(boundaryElementClosed.getValue(), -1); + CHECK_THROWS(boundaryElementClosed.setValue(0.2)); + } + + SUBCASE("Constant case") { + BoundaryElement boundaryElementConstant = BoundaryElement(0.1); + CHECK_NOTHROW(BoundaryElement(0.1)); + CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); + CHECK_EQ(boundaryElementConstant.getValue(), 0.1); + CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); + CHECK_EQ(boundaryElementConstant.getValue(), 0.2); + } } -TEST_CASE("Boundary Class"){ - Grid grid1D = Grid(10); - Grid grid2D = Grid(10, 12); - Boundary boundary1D = Boundary(grid1D); - Boundary boundary2D = Boundary(grid2D); - vector boundary1DVector(1, BoundaryElement(1.0)); +TEST_CASE("Boundary Class") { + Grid grid1D = Grid(10); + Grid grid2D = Grid(10, 12); + Boundary boundary1D = Boundary(grid1D); + Boundary boundary2D = Boundary(grid2D); + vector boundary1DVector(1, BoundaryElement(1.0)); - SUBCASE("Boundaries 1D case"){ - CHECK_NOTHROW(Boundary boundary(grid1D)); - CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); - CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); - CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); - CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); - CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); - CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); - CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); - CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - } + SUBCASE("Boundaries 1D case") { + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CLOSED); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CONSTANT); + CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), + boundary1DVector[0].getType()); + } - SUBCASE("Boundaries 2D case"){ - CHECK_NOTHROW(Boundary boundary(grid1D)); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); - CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); - CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); - CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); - CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); - CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); - CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - } -} \ No newline at end of file + SUBCASE("Boundaries 2D case") { + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CLOSED); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CONSTANT); + CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), + boundary1DVector[0].getType()); + } +} diff --git a/test/testFTCS.cpp b/test/testFTCS.cpp index d77999a..209ac8a 100644 --- a/test/testFTCS.cpp +++ b/test/testFTCS.cpp @@ -1,19 +1,20 @@ +#include + #include -#include <../src/FTCS.cpp> #include TEST_CASE("Maths") { - SUBCASE("mean between two alphas") { - double alpha1 = 10; - double alpha2 = 20; - double average = 15; - double harmonicMean = double(2) / ((double(1)/alpha1)+(double(1)/alpha2)); + SUBCASE("mean between two alphas") { + double alpha1 = 10; + double alpha2 = 20; + double average = 15; + double harmonicMean = + double(2) / ((double(1) / alpha1) + (double(1) / alpha2)); - // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) - harmonicMean); - // CHECK(difference < std::numeric_limits::epsilon()); - CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean); - CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average); - } - - -} \ No newline at end of file + // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) - + // harmonicMean); CHECK(difference < + // std::numeric_limits::epsilon()); + CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean); + CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average); + } +} diff --git a/test/testGrid.cpp b/test/testGrid.cpp index 4963976..a1f2886 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,251 +1,253 @@ +#include #include #include -#include + +using namespace Eigen; +using namespace std; + TEST_CASE("1D Grid, too small length") { - int l = 2; - CHECK_THROWS(Grid(l)); + int l = 2; + CHECK_THROWS(Grid(l)); } TEST_CASE("2D Grid, too small side") { - int r = 2; - int c = 4; - CHECK_THROWS(Grid(r, c)); + int r = 2; + int c = 4; + CHECK_THROWS(Grid(r, c)); - r = 4; - c = 2; - CHECK_THROWS(Grid(r, c)); + r = 4; + c = 2; + CHECK_THROWS(Grid(r, c)); } TEST_CASE("1D Grid") { - int l = 12; - Grid grid(l); + int l = 12; + Grid grid(l); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 1); - CHECK_EQ(grid.getLength(), l); - CHECK_EQ(grid.getCol(), l); - CHECK_EQ(grid.getRow(), 1); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 1); + CHECK_EQ(grid.getLength(), l); + CHECK_EQ(grid.getCol(), l); + CHECK_EQ(grid.getRow(), 1); - CHECK_EQ(grid.getConcentrations().rows(), 1); - CHECK_EQ(grid.getConcentrations().cols(), l); - CHECK_EQ(grid.getAlpha().rows(), 1); - CHECK_EQ(grid.getAlpha().cols(), l); - CHECK_EQ(grid.getDeltaCol(), 1); + CHECK_EQ(grid.getConcentrations().rows(), 1); + CHECK_EQ(grid.getConcentrations().cols(), l); + CHECK_EQ(grid.getAlpha().rows(), 1); + CHECK_EQ(grid.getAlpha().cols(), l); + CHECK_EQ(grid.getDeltaCol(), 1); - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); - CHECK_THROWS(grid.getDeltaRow()); - } + CHECK_THROWS(grid.getAlphaX()); + CHECK_THROWS(grid.getAlphaY()); + CHECK_THROWS(grid.getDeltaRow()); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(1, l, 3); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - SUBCASE("setting alpha") { - // correct alpha matrix - MatrixXd alpha = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setAlpha(alpha)); + SUBCASE("setting alpha") { + // correct alpha matrix + MatrixXd alpha = MatrixXd::Constant(1, l, 3); + CHECK_NOTHROW(grid.setAlpha(alpha)); - CHECK_THROWS(grid.setAlpha(alpha, alpha)); + CHECK_THROWS(grid.setAlpha(alpha, alpha)); - grid.setAlpha(alpha); - CHECK_EQ(grid.getAlpha(), alpha); - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); + grid.setAlpha(alpha); + CHECK_EQ(grid.getAlpha(), alpha); + CHECK_THROWS(grid.getAlphaX()); + CHECK_THROWS(grid.getAlphaY()); - // false alpha matrix - MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); - CHECK_THROWS(grid.setAlpha(wAlpha)); - } + // false alpha matrix + MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); + CHECK_THROWS(grid.setAlpha(wAlpha)); + } - SUBCASE("setting domain") { - int d = 8; - // set 1D domain - CHECK_NOTHROW(grid.setDomain(d)); + SUBCASE("setting domain") { + int d = 8; + // set 1D domain + CHECK_NOTHROW(grid.setDomain(d)); - // set 2D domain - CHECK_THROWS(grid.setDomain(d, d)); + // set 2D domain + CHECK_THROWS(grid.setDomain(d, d)); - grid.setDomain(d); - CHECK_EQ(grid.getDeltaCol(), double(d)/double(l)); - CHECK_THROWS(grid.getDeltaRow()); + grid.setDomain(d); + CHECK_EQ(grid.getDeltaCol(), double(d) / double(l)); + CHECK_THROWS(grid.getDeltaRow()); - // set too small domain - d = 0; - CHECK_THROWS(grid.setDomain(d)); - d = -2; - CHECK_THROWS(grid.setDomain(d)); - } + // set too small domain + d = 0; + CHECK_THROWS(grid.setDomain(d)); + d = -2; + CHECK_THROWS(grid.setDomain(d)); + } } TEST_CASE("2D Grid quadratic") { - int rc = 12; - Grid grid(rc, rc); + int rc = 12; + Grid grid(rc, rc); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); - CHECK_EQ(grid.getCol(), rc); - CHECK_EQ(grid.getRow(), rc); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 2); + CHECK_THROWS(grid.getLength()); + CHECK_EQ(grid.getCol(), rc); + CHECK_EQ(grid.getRow(), rc); - CHECK_EQ(grid.getConcentrations().rows(), rc); - CHECK_EQ(grid.getConcentrations().cols(), rc); - CHECK_THROWS(grid.getAlpha()); + CHECK_EQ(grid.getConcentrations().rows(), rc); + CHECK_EQ(grid.getConcentrations().cols(), rc); + CHECK_THROWS(grid.getAlpha()); - CHECK_EQ(grid.getAlphaX().rows(), rc); - CHECK_EQ(grid.getAlphaX().cols(), rc); - CHECK_EQ(grid.getAlphaY().rows(), rc); - CHECK_EQ(grid.getAlphaY().cols(), rc); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } + CHECK_EQ(grid.getAlphaX().rows(), rc); + CHECK_EQ(grid.getAlphaX().cols(), rc); + CHECK_EQ(grid.getAlphaY().rows(), rc); + CHECK_EQ(grid.getAlphaY().cols(), rc); + CHECK_EQ(grid.getDeltaRow(), 1); + CHECK_EQ(grid.getDeltaCol(), 1); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(rc, rc+3, 1); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc+3, rc, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc+2, rc+4, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + SUBCASE("setting alphas") { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); + MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); + CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); - MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); + CHECK_THROWS(grid.setAlpha(alphax)); - CHECK_THROWS(grid.setAlpha(alphax)); + grid.setAlpha(alphax, alphay); + CHECK_EQ(grid.getAlphaX(), alphax); + CHECK_EQ(grid.getAlphaY(), alphay); + CHECK_THROWS(grid.getAlpha()); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); + // false alpha matrices + alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + } - // false alpha matrices - alphax = MatrixXd::Constant(rc+3, rc+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(rc+2, rc+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - } + SUBCASE("setting domain") { + int dr = 8; + int dc = 9; - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; + // set 1D domain + CHECK_THROWS(grid.setDomain(dr)); - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); + // set 2D domain + CHECK_NOTHROW(grid.setDomain(dr, dc)); - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); + grid.setDomain(dr, dc); + CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); + CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc)/double(rc)); - CHECK_EQ(grid.getDeltaRow(), double(dr)/double(rc)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); - } + // set too small domain + dr = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = 8; + dc = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = -2; + CHECK_THROWS(grid.setDomain(dr, dc)); + } } TEST_CASE("2D Grid non-quadratic") { - int r = 12; - int c = 15; - Grid grid(r, c); + int r = 12; + int c = 15; + Grid grid(r, c); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); - CHECK_EQ(grid.getCol(), c); - CHECK_EQ(grid.getRow(), r); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 2); + CHECK_THROWS(grid.getLength()); + CHECK_EQ(grid.getCol(), c); + CHECK_EQ(grid.getRow(), r); - CHECK_EQ(grid.getConcentrations().rows(), r); - CHECK_EQ(grid.getConcentrations().cols(), c); - CHECK_THROWS(grid.getAlpha()); + CHECK_EQ(grid.getConcentrations().rows(), r); + CHECK_EQ(grid.getConcentrations().cols(), c); + CHECK_THROWS(grid.getAlpha()); - CHECK_EQ(grid.getAlphaX().rows(), r); - CHECK_EQ(grid.getAlphaX().cols(), c); - CHECK_EQ(grid.getAlphaY().rows(), r); - CHECK_EQ(grid.getAlphaY().cols(), c); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } + CHECK_EQ(grid.getAlphaX().rows(), r); + CHECK_EQ(grid.getAlphaX().cols(), c); + CHECK_EQ(grid.getAlphaY().rows(), r); + CHECK_EQ(grid.getAlphaY().cols(), c); + CHECK_EQ(grid.getDeltaRow(), 1); + CHECK_EQ(grid.getDeltaCol(), 1); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(r, c, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(r, c, 2); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(r + 3, c, 3); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(r, c+3, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r+3, c, 3); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r+2, c+4, 2); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + SUBCASE("setting alphas") { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(r, c, 2); + MatrixXd alphay = MatrixXd::Constant(r, c, 4); + CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(r, c, 2); - MatrixXd alphay = MatrixXd::Constant(r, c, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); + CHECK_THROWS(grid.setAlpha(alphax)); - CHECK_THROWS(grid.setAlpha(alphax)); + grid.setAlpha(alphax, alphay); + CHECK_EQ(grid.getAlphaX(), alphax); + CHECK_EQ(grid.getAlphaY(), alphay); + CHECK_THROWS(grid.getAlpha()); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); + // false alpha matrices + alphax = MatrixXd::Constant(r + 3, c + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + alphay = MatrixXd::Constant(r + 2, c + 1, 5); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + } - // false alpha matrices - alphax = MatrixXd::Constant(r+3, c+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(r+2, c+1, 5); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - } + SUBCASE("setting domain") { + int dr = 8; + int dc = 9; - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; + // set 1D domain + CHECK_THROWS(grid.setDomain(dr)); - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); + // set 2D domain + CHECK_NOTHROW(grid.setDomain(dr, dc)); - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); + grid.setDomain(dr, dc); + CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); + CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc)/double(c)); - CHECK_EQ(grid.getDeltaRow(), double(dr)/double(r)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = -1; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); - } -} \ No newline at end of file + // set too small domain + dr = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = 8; + dc = -1; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = -2; + CHECK_THROWS(grid.setDomain(dr, dc)); + } +} diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index c70a475..004799a 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -1,106 +1,108 @@ -#include +#include "TestUtils.hpp" + #include -#include -#include "TestUtils.cpp" +#include #include +#include // include the configured header file #include -static Grid setupSimulation(APPROACH approach, double timestep, int iterations) { - int row = 11; - int col = 11; - int domain_row = 10; - int domain_col = 10; +using namespace Eigen; +using namespace std; +static Grid setupSimulation(APPROACH approach, double timestep, + int iterations) { + int row = 11; + int col = 11; + int domain_row = 10; + int domain_col = 10; - // Grid - Grid grid = Grid(row, col); - grid.setDomain(domain_row, domain_col); + // Grid + Grid grid = Grid(row, col); + grid.setDomain(domain_row, domain_col); - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5,5) = 1; - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(5, 5) = 1; + grid.setConcentrations(concentrations); - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } + MatrixXd alpha = MatrixXd::Constant(row, col, 1); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 6; j++) { + alpha(i, j) = 0.01; } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } + } + for (int i = 0; i < 5; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.001; } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } + } + for (int i = 5; i < 11; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.1; } - grid.setAlpha(alpha, alpha); + } + grid.setAlpha(alpha, alpha); + // Boundary + Boundary bc = Boundary(grid); - // Boundary - Boundary bc = Boundary(grid); + // Simulation + Simulation sim = Simulation(grid, bc, approach); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); - - // Simulation - Simulation sim = Simulation(grid, bc, approach); - sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - // RUN - return grid; - + // RUN + return grid; } TEST_CASE("equality to reference matrix with FTCS") { - // set string from the header file - string test_path = testSimulationCSVDir; - MatrixXd reference = CSV2Eigen(test_path); - Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); - CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); + // set string from the header file + string test_path = testSimulationCSVDir; + MatrixXd reference = CSV2Eigen(test_path); + cout << "FTCS Test: " << endl; + Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); + cout << endl; + CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); } TEST_CASE("equality to reference matrix with BTCS") { - // set string from the header file - string test_path = testSimulationCSVDir; - MatrixXd reference = CSV2Eigen(test_path); - Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); - CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); + // set string from the header file + string test_path = testSimulationCSVDir; + MatrixXd reference = CSV2Eigen(test_path); + cout << "BTCS Test: " << endl; + Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); + cout << endl; + CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); } -TEST_CASE("Initialize environment"){ - int rc = 12; - Grid grid(rc, rc); - Boundary boundary(grid); +TEST_CASE("Initialize environment") { + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); - CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); + CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); } -TEST_CASE("Simulation environment"){ - int rc = 12; - Grid grid(rc, rc); - Boundary boundary(grid); - Simulation sim(grid, boundary, FTCS_APPROACH); +TEST_CASE("Simulation environment") { + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); + Simulation sim(grid, boundary, FTCS_APPROACH); - SUBCASE("default paremeters"){ - CHECK_EQ(sim.getIterations(), -1); - } + SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); } - SUBCASE("set iterations"){ - CHECK_NOTHROW(sim.setIterations(2000)); - CHECK_EQ(sim.getIterations(), 2000); - CHECK_THROWS(sim.setIterations(-300)); - } + SUBCASE("set iterations") { + CHECK_NOTHROW(sim.setIterations(2000)); + CHECK_EQ(sim.getIterations(), 2000); + CHECK_THROWS(sim.setIterations(-300)); + } - SUBCASE("set timestep"){ - CHECK_NOTHROW(sim.setTimestep(0.1)); - CHECK_EQ(sim.getTimestep(), 0.1); - CHECK_THROWS(sim.setTimestep(-0.3)); - } + SUBCASE("set timestep") { + CHECK_NOTHROW(sim.setTimestep(0.1)); + CHECK_EQ(sim.getTimestep(), 0.1); + CHECK_THROWS(sim.setTimestep(-0.3)); + } } -