BREAKING CHANGE: reworked API

BREAKING CHANGE: added heterogeneous diffusion

BREAKING CHANGE: added FTCS scheme

See merge request naaice/tug!16

.gitignore

@@ -9,3 +9,4 @@ compile_commands.json
 .ccls-cache/
 /iwyu/
 .Rhistory
+.vscode/

.gitlab-ci.yml

@@ -4,6 +4,7 @@ stages:
     - build
     - test
     - static_analyze
+    - doc
 
 build_release:
   stage: build
@@ -15,18 +16,36 @@ build_release:
       - mkdir build && cd build
       - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON ..
       - make -j$(nproc)
+      - cp ../test/FTCS_11_11_7000.csv test/
 
 test:
     stage: test
     script:
         - ./build/test/testTug
 
+pages:
+  stage: doc
+  image: python:slim
+  before_script:
+    - apt-get update && apt-get install --no-install-recommends -y graphviz imagemagick doxygen make
+    - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme
+    - mkdir public
+  script:
+    - pushd docs_sphinx
+    - make html
+    - popd && mv docs_sphinx/_build/html/* public/
+  artifacts:
+    paths:
+      - public
+  rules:
+    - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH
+
 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

CMakeLists.txt

@@ -7,6 +7,8 @@ set(CMAKE_CXX_STANDARD 17)
 
 find_package(Eigen3 REQUIRED NO_MODULE)
 find_package(OpenMP)
+# find_package(easy_profiler)
+# option(EASY_OPTION_LOG "Verbose easy_profiler" 1)
 
 ## SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -O2 -mfma")
 option(TUG_USE_OPENMP "Compile with OpenMP support" ON)
@@ -31,10 +33,18 @@ endif()
 
 option(TUG_ENABLE_TESTING
   "Run tests after succesfull compilation"
-  ON)
+  OFF)
+
+option(TUG_ENABLE_EXAMPLES
+  "Compile example applications"
+  OFF)
 
 add_subdirectory(src)
 
 if(TUG_ENABLE_TESTING)
   add_subdirectory(test)
 endif()
+
+if(TUG_ENABLE_EXAMPLES)
+  add_subdirectory(examples)
+endif()

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.
+
+  <signature of Ty Coon>, 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.

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(<your_target> 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?

doc/ADI_scheme.org

@@ -1,6 +1,7 @@
 #+TITLE: Finite Difference Schemes for the numerical solution of heterogeneous diffusion equation in 2D
 #+LaTeX_CLASS_OPTIONS: [a4paper,10pt]
 #+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{charter}
 #+LATEX_HEADER: \usepackage{amsmath, systeme, cancel, xcolor}
 #+OPTIONS: toc:nil
 
@@ -487,10 +488,39 @@ C_{i,j}^{t+1} = & C^t_{i,j} +\\
 \end{aligned}
 \end{equation}
 
-The Courant-Friedrichs-Lewy stability criterion for this scheme reads:
+The Courant-Friedrichs-Lewy stability criterion (cfr Lee, 2017) for
+this scheme reads:
+
+#+NAME: eqn:CFL2DFTCS_Lee
+\begin{equation}
+\Delta t \leq \frac{1}{2 \max(\alpha_{i,j})} \cdot \frac{1}{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}
+\end{equation}
+
+
+Note that other derivations for the CFL condition are found in
+literature. For example, the sources cited by [[https://en.wikipedia.org/wiki/FTCS_scheme][Wikipedia solution]] give:
+
+#+NAME: eqn:CFL2DFTCS_wiki
+\begin{equation}
+\displaystyle \Delta t\leq {\frac {1}{4 \max(\alpha) \left({\frac {1}{\Delta x^{2}}}+{\frac {1}{\Delta y^{2}}}\right)}}
+\end{equation}
+
+We can produce a more restrictive condition than equation
+[[eqn:CFL2DFTCS_Lee]] by considering the min of the $\Delta x$ and $\Delta
+y$:
+
 #+NAME: eqn:CFL2DFTCS
 \begin{equation}
-\Delta t \leq \frac{(\Delta x^2, \Delta y^2)}{2 \max(\alpha_{i,j})}
+\Delta t \leq \frac{\min(\Delta x, \Delta y)^2}{4 \max(\alpha_{i,j})}
+\end{equation}
+
+In practice for the implementation it is advantageous to specify an
+optional parameter $C$, $C \in [0, 1]$ so that the user can restrict
+the "inner time stepping":
+
+#+NAME: eqn:CFL2DFTCS_impl
+\begin{equation}
+\Delta t \leq C \cdot \frac{\min(\Delta x, \Delta y)^2}{4 \max(\alpha_{i,j})}
 \end{equation}
 
 ** Boundary conditions
@@ -536,4 +566,3 @@ left derivative becomes zero and we are left with:
 &\displaystyle =\frac{\alpha_{i-1/2}}{\Delta x^2} (C_{i-1} - C_i)
 \end{aligned}
 \end{equation}
-

doc/ValidationHetDiff.org

@@ -0,0 +1,115 @@
+#+TITLE: Validation Examples for 2D Heterogeneous Diffusion
+#+AUTHOR: MDL <delucia@gfz-potsdam.de>
+#+DATE: 2023-08-26
+#+STARTUP: inlineimages
+#+LATEX_CLASS_OPTIONS: [a4paper,9pt]
+#+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{amsmath, systeme}
+#+LATEX_HEADER: \usepackage{graphicx}
+#+LATEX_HEADER: \usepackage{charter}
+#+OPTIONS: toc:nil
+
+
+* Simple setup using deSolve/ReacTran
+
+- Square of side 10
+- Discretization: 11 \times 11 cells
+- All boundaries closed
+- Initial state: 0 everywhere, 1 in the center (6,6)
+- Time step: 1 s, 10 iterations
+
+The matrix of spatially variable diffusion coefficients \alpha is
+constant in 4 quadrants:
+
+\alpha_{x,y} = 
+
+| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 |
+| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 |
+| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 |
+| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 |
+| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+|    1 |    1 |    1 |    1 |    1 |    1 |   0.1 |   0.1 |   0.1 |   0.1 |   0.1 |
+
+The relevant part of the R script used to produce these results is
+presented in listing 1; the whole script is at [[file:scripts/HetDiff.R]].
+A visualization of the output of the reference simulation is given in
+figure [[fig:1][1]].
+
+Note: all results from this script are stored in the =outc= matrix by
+the =deSolve= function. I stored a different version into
+[[file:../scripts/gold/HetDiff1.csv]]: this file contains 11 columns (one
+for each time step including initial conditions) and 121 rows, one for
+each domain element, with no headers.
+
+#+caption: Result of ReacTran/deSolve solution of the above problem at 4
+#+name: fig:1
+[[./images/deSolve_AlphaHet1.png]]
+
+
+#+name: lst:1
+#+begin_src R :language R :frame single :caption Listing 1, generate reference simulation using R packages deSolve/ReacTran :captionpos b :label lst:1
+  library(ReacTran)
+  library(deSolve)
+
+  ## harmonic mean
+  harm <- function(x,y) {
+    if (length(x) != 1 || length(y) != 1)
+    stop("x & y have different lengths")
+    2/(1/x+1/y)
+  }
+
+  N     <- 11          # number of grid cells
+  ini   <- 1           # initial value at x=0
+  N2    <- ceiling(N/2)
+  L     <- 10          # domain side
+
+  ## Define diff.coeff per cell, in  4 quadrants
+  alphas <- matrix(0, N, N) 
+  alphas[1:N2, 1:N2] <- 1
+  alphas[1:N2, seq(N2+1,N)] <- 0.1
+  alphas[seq(N2+1,N), 1:N2] <- 0.01
+  alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001
+
+  cmpharm <- function(x) {
+    y <- c(0, x, 0)
+    ret <- numeric(length(x)+1)
+    for (i in seq(2, length(y))) {
+      ret[i-1] <- harm(y[i], y[i-1])
+    }
+    ret
+  }
+
+  ## Construction of the 2D grid
+  x.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  y.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  grid2D <- setup.grid.2D(x.grid, y.grid)
+  dx <- dy <- L/N
+
+  D.grid <- list()
+  ## Diffusion coefs on x-interfaces
+  D.grid$x.int <- apply(alphas, 1, cmpharm)
+  ## Diffusion coefs on y-interfaces
+  D.grid$y.int <- t(apply(alphas, 2, cmpharm))
+
+  # The model
+  Diff2Dc <- function(t, y, parms) {
+    CONC  <- matrix(nrow = N, ncol = N, data = y)
+    dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC
+    return(list(dCONC))
+  }
+
+  ## initial condition: 0 everywhere, except in central point
+  y <- matrix(nrow = N, ncol = N, data = 0)
+  y[N2, N2] <- ini  # initial concentration in the central point...
+
+  ## solve for 10 time units
+  times <- 0:10
+  outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL,
+		 dim = c(N, N), lrw = 1860000)
+#+end_src
+

docs_sphinx/Boundary.rst

@@ -0,0 +1,4 @@
+Boundary
+========
+
+.. doxygenclass:: Boundary

docs_sphinx/Grid.rst

@@ -0,0 +1,4 @@
+Grid
+====
+
+.. doxygenclass:: Grid

docs_sphinx/Makefile

@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

docs_sphinx/Simulation.rst

@@ -0,0 +1,4 @@
+Simulation
+==========
+
+.. doxygenclass:: Simulation

docs_sphinx/conf.py

@@ -0,0 +1,97 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# This file only contains a selection of the most common options. For a full
+# list see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+# import os
+# import sys
+# sys.path.insert(0, os.path.abspath('.'))
+from sphinx.builders.html import StandaloneHTMLBuilder
+import subprocess, os
+
+# Doxygen
+subprocess.call('doxygen Doxyfile.in', shell=True)
+
+# -- Project information -----------------------------------------------------
+
+project = 'TUG'
+copyright = 'GPL2'
+author = 'Philipp Ungrund, Hannes Signer'
+
+
+# -- General configuration ---------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.intersphinx',
+    'sphinx.ext.autosectionlabel',
+    'sphinx.ext.todo',
+    'sphinx.ext.coverage',
+    'sphinx.ext.mathjax',
+    'sphinx.ext.ifconfig',
+    'sphinx.ext.viewcode',
+    'sphinx.ext.inheritance_diagram',
+    'breathe'
+]
+
+html_baseurl = "/index.html"
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+highlight_language = 'c++'
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'sphinx_rtd_theme'
+html_theme_options = {
+    'canonical_url': '',
+    'analytics_id': '',  #  Provided by Google in your dashboard
+    'display_version': True,
+    'prev_next_buttons_location': 'bottom',
+    'style_external_links': False,
+    
+    'logo_only': False,
+
+    # Toc options
+    'collapse_navigation': True,
+    'sticky_navigation': True,
+    'navigation_depth': 4,
+    'includehidden': True,
+    'titles_only': False
+}
+# html_logo = ''
+# github_url = ''
+# html_baseurl = ''
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+html_logo = "images/tug_logo.svg"
+
+# -- Breathe configuration -------------------------------------------------
+
+breathe_projects = {
+	"Tug": "_build/xml/"
+}
+breathe_default_project = "Tug"
+breathe_default_members = ('members', 'undoc-members')

docs_sphinx/developer.rst

@@ -0,0 +1,38 @@
+Developer Guide
+===============
+
+=========================
+Class Diagram of user API
+=========================
+
+The following graphic shows the class diagram of the user API. The FTCS and 
+BTCS functionalities are externally outsourced and not visible to the user.
+
+.. image:: images/class_diagram.svg
+    :width: 2000
+    :alt: Class diagram for the user API
+
+====================================================
+Activity Diagram for run routine in simulation class
+====================================================
+
+The following activity diagram represents the actions when the run method is called within the simulation class. 
+For better distinction, the activities of the calculation methods FTCS and BTCS are shown in two separate activity diagrams.
+
+.. image:: images/activity_diagram_run.svg
+    :width: 2000
+    :alt: Activity diagram for the run method in the simulation class
+
+
+**Activity Diagram for FTCS method**
+
+.. image:: images/activity_diagram_FTCS.svg
+    :width: 400
+    :alt: Activity diagram for the FTCS method
+
+
+**Activity Diagram for BTCS method**
+
+.. image:: images/activity_diagram_BTCS.svg
+    :width: 400
+    :alt: Activity diagram for the BTCS method

docs_sphinx/developper.rst

@@ -0,0 +1,2 @@
+Developper API
+==============

docs_sphinx/examples.rst

@@ -0,0 +1,62 @@
+Examples
+========
+
+At this point, some typical commented examples are presented to illustrate how Tug works.
+In general, each simulation is divided into three blocks:
+    - the initialization of the grid, which is to be simulated
+    - the setting of the boundary conditions
+    - the setting of the simulation parameters and the start of the simulation
+
+Two dimensional grid with constant boundaries and FTCS method
+-------------------------------------------------------------
+**Initialization of the grid**
+
+For example, the initalization of a grid with 20 by 20 cells and a domain size (physical extent of the grid) of 
+also 20 by 20 length units can be done as follows. The setting of the domain is optional here and is set to the 
+same size as the number of cells in the standard case. As seen in the code, the cells of the grid are set to an 
+initial value of 0 and only in the upper left corner (0,0) the starting concentration is set to the value 20.
+
+.. code-block:: cpp
+    
+    int row = 20
+    int col = 20;
+    Grid grid = Grid(row,col);
+    grid.setDomain(row, col);
+    MatrixXd concentrations = MatrixXd::Constant(row,col,0);
+    concentrations(0,0) = 20;
+    grid.setConcentrations(concentrations);
+
+**Setting of the boundary conditions**
+
+First, a boundary class is created and then the corresponding boundary conditions are set. In this case, all four sides
+of the grid are set as constant edges with a concentration of 0.
+
+.. code-block:: cpp
+
+    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);
+
+**Setting of the simulation parameters and simulation start**
+In the last block, a simulation class is created and the objects of the grid and the boundary conditions are passed. The solution
+method is also specified (either FCTS or BTCS). Furthermore, the desired time step and the number of iterations are set. The penultimate
+parameter specifies the output of the simulated results in a CSV file. In the present case, the result of each iteration step is written
+one below the other into the corresponding CSV file.
+
+.. code-block:: cpp
+    
+    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH);
+    simulation.setTimestep(0.1);
+    simulation.setIterations(1000);
+    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    simulation.run();
+
+
+
+
+
+
+Setting special boundary conditions on individual cells
+-------------------------------------------------------

docs_sphinx/index.rst

@@ -0,0 +1,51 @@
+.. Tug documentation master file, created by
+   sphinx-quickstart on Mon Aug 14 11:30:23 2023.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+Welcome to Tug's documentation!
+===============================
+
+Welcome to the documentation of the TUG project, a simulation program
+for solving transport equations in one- and two-dimensional uniform
+grids using cell centered finite differences.
+
+---------
+Diffusion
+---------
+
+TUG can solve diffusion problems with heterogeneous and anisotropic
+diffusion coefficients. The partial differential equation expressing
+diffusion reads:
+
+.. math::
+   \frac{\partial C}{\partial t} =  \nabla \cdot \left[ \mathbf{\alpha} \nabla C \right]
+
+In 2D, the equation reads:
+   
+.. math::
+   \frac{\partial C}{\partial t} =  \frac{\partial}{\partial x}\left[ \alpha_x \frac{\partial C}{\partial x}\right] + \frac{\partial}{\partial y}\left[ \alpha_y \frac{\partial C}{\partial y}\right]
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`search`
+
+Table of Contents
+^^^^^^^^^^^^^^^^^
+.. toctree::    
+
+    :maxdepth: 2
+    self
+    installation
+    theory
+    user
+    developer
+    examples
+    visualization

docs_sphinx/installation.rst

@@ -0,0 +1,3 @@
+Installation
+============
+

docs_sphinx/requirements.txt

@@ -0,0 +1,2 @@
+breathe
+sphinx-rtd-theme

docs_sphinx/theory.rst

@@ -0,0 +1,15 @@
+Theoretical Foundations
+=======================
+
+=====================
+The Diffusion Problem
+=====================
+
+================
+Numerical Solver
+================
+
+**Backward Time-Centered Space (BTCS) Method**
+
+
+**Forward Time-Centered Space (BTCS) Method**

docs_sphinx/user.rst

@@ -0,0 +1,9 @@
+User API
+========
+
+.. toctree::    
+  
+    :maxdepth: 2
+    Grid
+    Boundary
+    Simulation

docs_sphinx/visualization.rst

@@ -0,0 +1,3 @@
+Visualization
+=============
+

examples/BTCS_1D_proto_example.cpp

@@ -0,0 +1,51 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+
+int main(int argc, char *argv[]) {
+  // **************
+  // **** GRID ****
+  // **************
+
+  // 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);
+
+  // ******************
+  // **** BOUNDARY ****
+  // ******************
+
+  // 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
+
+  // set the timestep of the simulation
+  simulation.setTimestep(0.1); // timestep
+
+  // set the number of iterations
+  simulation.setIterations(100);
+
+  // 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();
+}

examples/BTCS_2D_proto_example.cpp

@@ -0,0 +1,85 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+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);
+
+  // (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 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 ****
+  // ******************
+
+  // 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 ****
+  // ************************
+
+  // 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
+
+  // set the number of iterations
+  simulation.setIterations(300);
+
+  // 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();
+}

examples/CMakeLists.txt

@@ -0,0 +1,20 @@
+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)
+
+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)
+
+add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
+add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+target_link_libraries(FTCS_2D_proto_closed_mdl tug)
+target_link_libraries(FTCS_2D_proto_example_mdl tug)

examples/CRNI_2D_proto_example.cpp

@@ -0,0 +1,27 @@
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+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();
+}

examples/FTCS_1D_proto_example.cpp

@@ -0,0 +1,50 @@
+#include "tug/Boundary.hpp"
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+
+int main(int argc, char *argv[]) {
+  // **************
+  // **** GRID ****
+  // **************
+
+  // create a linear grid with 20 cells
+  int cells = 20;
+  Grid grid = Grid(cells);
+
+  MatrixXd concentrations = MatrixXd::Constant(1, 20, 20);
+  grid.setConcentrations(concentrations);
+
+  // ******************
+  // **** BOUNDARY ****
+  // ******************
+
+  // 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
+
+  // (optional) set the timestep of the simulation
+  simulation.setTimestep(0.1); // timestep
+
+  // (optional) set the number of iterations
+  simulation.setIterations(100);
+
+  // (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();
+}

examples/FTCS_2D_proto_closed_mdl.cpp

@@ -0,0 +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
+ *
+ */
+
+#include <Eigen/Eigen>
+#include <cstdlib>
+#include <iostream>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+
+int main(int argc, char *argv[]) {
+
+  int row = 64;
+
+  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 ****
+  // **************
+
+  // 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 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 ****
+  // ******************
+
+  // 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);
+
+  // ************************
+  // **** SIMULATION ENV ****
+  // ************************
+
+  // set up a simulation environment
+  Simulation simulation =
+      Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
+
+  // set the timestep of the simulation
+  simulation.setTimestep(10000); // timestep
+
+  // set the number of iterations
+  simulation.setIterations(100);
+
+  // (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();
+
+  return 0;
+}

examples/FTCS_2D_proto_example.cpp

@@ -0,0 +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
+ *
+ */
+
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+using namespace Eigen;
+// #include <easy/profiler.h>
+// #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);
+
+  // (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 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 ****
+  // ******************
+
+  // 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 ****
+  // ************************
+
+  // 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
+
+  // set the number of iterations
+  simulation.setIterations(10000);
+
+  // 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();
+}

examples/FTCS_2D_proto_example_mdl.cpp

@@ -0,0 +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
+ *
+ */
+
+#include <Eigen/Eigen>
+#include <tug/Simulation.hpp>
+
+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);
+
+  // (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 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 ****
+  // ******************
+
+  // 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);
+
+  // ************************
+  // **** SIMULATION ENV ****
+  // ************************
+
+  // 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 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);
+
+  // **** RUN SIMULATION ****
+
+  // run the simulation
+  simulation.run();
+
+  return 0;
+}

examples/Rcpp-BTCS-1d.cpp

@@ -1,59 +0,0 @@
-#include "../include/diffusion/BTCSDiffusion.hpp"
-#include "../include/diffusion/BoundaryCondition.hpp"
-#include <algorithm> // for copy, max
-#include <cmath>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-#include <Rcpp.h>
-
-using namespace std;
-using namespace Diffusion;
-using namespace Rcpp;
-//using namespace Eigen;
-
-// [[Rcpp::depends(RcppEigen)]]
-// [[Rcpp::plugins("cpp11")]]
-// [[Rcpp::export]]
-std::vector<double> & diff1D(int n,
-                             double length,
-                             std::vector<double> & field,
-                             std::vector<double> & alpha,
-                             double timestep,
-                             double bc_left,
-                             double bc_right,
-                             int iterations) {
-    // dimension of grid
-    int dim = 1;
-    
-    // create input + diffusion coefficients for each grid cell
-    // std::vector<double> alpha(n, 1 * pow(10, -1));
-    // std::vector<double> field(n, 1 * std::pow(10, -6));
-    std::vector<boundary_condition> bc(n, {0,0});
-    
-    // create instance of diffusion module
-    BTCSDiffusion diffu(dim);
-    
-    diffu.setXDimensions(length, n);
-    
-    // set the boundary condition for the left ghost cell to dirichlet
-    bc[0] = {Diffusion::BC_CONSTANT, bc_left};
-    bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
-    
-    // set timestep for simulation to 1 second
-    diffu.setTimestep(timestep);
-    
-    //cout << setprecision(12);
-    
-    // loop 100 times
-    // output is currently generated by the method itself
-    for (int i = 0; i < iterations; i++) {
-        diffu.simulate(field.data(), alpha.data(), bc.data());
-    }
-
-    // for (auto & cell : field) {
-    //     Rcout << cell << "\n";
-    // }
-    
-    return(field);
-}

examples/Rcpp-BTCS-2d.cpp

@@ -1,60 +0,0 @@
-#include "../include/diffusion/BTCSDiffusion.hpp"
-#include "../include/diffusion/BoundaryCondition.hpp"
-#include <algorithm> // for copy, max
-#include <cmath>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-#include <Rcpp.h>
-
-using namespace std;
-using namespace Diffusion;
-using namespace Rcpp;
-//using namespace Eigen;
-
-// [[Rcpp::depends(RcppEigen)]]
-// [[Rcpp::plugins("cpp11")]]
-// [[Rcpp::export]]
-std::vector<double> & diff2D(int nx,
-                             int ny,
-                             double lenx,
-                             double leny,
-                             std::vector<double> & field,
-                             std::vector<double> & alpha,
-                             double timestep,
-                             int iterations)
-{
-    // problem dimensionality
-    int dim = 2;
-    // total number of grid cells
-    int n = nx*ny;
-    
-    std::vector<boundary_condition> bc(n, {0,0});
-    
-    // create instance of diffusion module
-    BTCSDiffusion diffu(dim);
-    
-    diffu.setXDimensions(lenx, nx);
-    diffu.setXDimensions(leny, ny);
-    
-    // set the boundary condition for the left ghost cell to dirichlet
-    // bc[0] = {Diffusion::BC_CONSTANT, bc_left};
-    // bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
-    
-    // set timestep for simulation to 1 second
-    diffu.setTimestep(timestep);
-    
-    //cout << setprecision(12);
-    
-    // loop 100 times
-    // output is currently generated by the method itself
-    for (int i = 0; i < iterations; i++) {
-        diffu.simulate(field.data(), alpha.data(), bc.data());
-    }
-
-    // for (auto & cell : field) {
-    //     Rcout << cell << "\n";
-    // }
-    
-    return(field);
-}

examples/Rcpp-interface.R

@@ -1,159 +0,0 @@
-## Time-stamp: "Last modified 2022-03-16 14:01:11 delucia"
-library(Rcpp)
-library(RcppEigen)
-library(ReacTran)
-library(deSolve)
-
-options(width=110)
-
-setwd("app")
-
-## This creates the "diff1D" function with our BTCSdiffusion code
-sourceCpp("Rcpp-BTCS-1d.cpp")
-
-### FTCS explicit (same name)
-sourceCpp("RcppFTCS.cpp")
-
-
-
-## Grid 101
-## Set initial conditions
-N <- 1001
-D.coeff <- 1E-3
-C0 <- 1             ## Initial concentration (mg/L)
-X0 <- 0             ## Location of initial concentration (m)
-## Yini <- c(C0, rep(0,N-1))
-
-## Ode1d solution
-xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N)
-x <- xgrid$x.mid
-Diffusion <- function (t, Y, parms){
-  tran <- tran.1D(C = Y, C.up = 0, C.down = 0, D = parms$D, dx = xgrid)
-  return(list(tran$dC))
-}
-
-
-## gaussian pulse as initial condition
-sigma <- 0.02
-Yini <- 0.5*exp(-0.5*((x-1/2.0)**2)/sigma**2)
-
-## plot(x, Yini, type="l")
-
-parms1 <- list(D=D.coeff)
-# 1 timestep, 10 s
-times <- seq(from = 0, to = 1, by = 0.1)
-
-system.time({
-    out1 <- ode.1D(y = Yini, times = times, func = Diffusion,
-                  parms = parms1, dimens = N)[11,-1]
-})
-
-## Now with BTCS 
-alpha <- rep(D.coeff, N)
-
-system.time({
-    out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 0.1, 0, 0, iterations = 10)
-})
-
-## plot(out1, out2)
-## abline(0,1)
-
-## matplot(cbind(out1,out2),type="l", col=c("black","red"),lty="solid", lwd=2,
-##         xlab="grid element", ylab="Concentration", las=1)
-## legend("topright", c("ReacTran ode1D", "BTCS 1d"), text.col=c("black","red"), bty = "n")
-
-
-
-
-
-system.time({
-    out3 <- RcppFTCS(n=N, length=1, field=Yini, alpha=1E-3, bc_left = 0, bc_right = 0, timestep = 1)
-})
-
-## Poor man's 
-mm <- colMeans(rbind(out2,out3))
-
-
-
-matplot(cbind(Yini,out1, out2, out3, mm),type="l", col=c("grey","black","red","blue","green4"), lty="solid", lwd=2,
-        xlab="grid element", ylab="Concentration", las=1)
-legend("topright", c("init","ReacTran ode1D", "BTCS 1d", "FTCS", "poor man's CN"), text.col=c("grey","black","red","blue","green4"), bty = "n")
-
-
-sum(Yini)
-sum(out1)
-sum(out2)
-sum(out3)
-sum(mm)
-
-## Yini <- 0.2*sin(pi/0.1*x)+0.2
-## plot(Yini)
-
-## plot(out3)
-
-Fun <- function(dx) {
-    tmp <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = floor(1/dx))
-    sqrt(sum({out1-tmp}^2))
-}
-
-reso <- optimise(f=Fun, interval=c(1E-5, 1E-1), maximum = FALSE)
-
-
-dx <- 0.0006038284
-floor(1/dx)
-
-1/dx
-
-system.time({
-    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1656)
-})
-
-matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
-        xlab="grid element", ylab="Concentration", las=1)
-legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
-
-
-dx <- 0.05
-
-system.time({
-    out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1/dx)
-})
-
-matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2,
-        xlab="grid element", ylab="Concentration", las=1)
-legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n")
-
-Matplot
-
-
-
-
-## This creates the "diff1D" function with our BTCSdiffusion code
-sourceCpp("Rcpp-BTCS-2d.cpp")
-
-n <- 256
-a2d <- rep(1E-3, n^2)
-
-init2d <- readRDS("gs1.rds")
-
-ll <- {init2d - min(init2d)}/diff(range(init2d))
-
-system.time({
-    res1 <- diff2D(nx=N, ny=N, lenx=1, leny=1, field=ll, alpha=a2d, timestep = 0.1, iterations = 10)
-})
-
-hist(ll,32)
-
-
-
-
-
-
-
-
-
-
-
-
-
-

examples/Rcpp-interface.cpp

@@ -1,59 +0,0 @@
-#include "../include/diffusion/BTCSDiffusion.hpp"
-#include "../include/diffusion/BoundaryCondition.hpp"
-#include <algorithm> // for copy, max
-#include <cmath>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-#include <Rcpp.h>
-
-using namespace std;
-using namespace Diffusion;
-using namespace Rcpp;
-//using namespace Eigen;
-
-// [[Rcpp::depends(RcppEigen)]]
-// [[Rcpp::plugins("cpp11")]]
-// [[Rcpp::export]]
-std::vector<double> & diff1D(int n,
-                             double length,
-                             std::vector<double> & field,
-                             std::vector<double> & alpha,
-                             double timestep,
-                             double bc_left,
-                             double bc_right,
-                             int iterations) {
-    // dimension of grid
-    int dim = 1;
-    
-    // create input + diffusion coefficients for each grid cell
-    // std::vector<double> alpha(n, 1 * pow(10, -1));
-    // std::vector<double> field(n, 1 * std::pow(10, -6));
-    std::vector<boundary_condition> bc(n, {0,0});
-    
-    // create instance of diffusion module
-    BTCSDiffusion diffu(dim);
-    
-    diffu.setXDimensions(length, n);
-    
-    // set the boundary condition for the left ghost cell to dirichlet
-    bc[0] = {Diffusion::BC_CONSTANT, bc_left};
-    bc[n-1] = {Diffusion::BC_CONSTANT, bc_right};
-    
-    // set timestep for simulation to 1 second
-    diffu.setTimestep(timestep);
-    
-    //cout << setprecision(12);
-    
-    // loop 100 times
-    // output is currently generated by the method itself
-    for (int i = 0; i < iterations; i++) {
-        diffu.simulate(field.data(), alpha.data(), bc.data());
-    }
-
-    // for (auto & cell : field) {
-    //     Rcout << cell << "\n";
-    // }
-    
-    return(field);
-}

examples/RcppFTCS.cpp

@@ -1,40 +0,0 @@
-//  Time-stamp: "Last modified 2022-03-15 18:15:39 delucia"
-#include <Rcpp.h>
-#include <iostream> // for std
-#include <vector>   // for vector
-
-using namespace std;
-using namespace Rcpp;
-
-// [[Rcpp::plugins("cpp11")]]
-// [[Rcpp::export]]
-NumericVector RcppFTCS(int n,
-                       double length,    
-                       NumericVector & field,
-                       double alpha,
-                       double bc_left,
-                       double bc_right,
-                       double timestep)
-{
-    // dimension of grid
-    NumericVector ext (clone(field));
-    double dx = length / ((double) n - 1.);
-    double dt = 0.25*dx*dx/alpha;
-
-    
-    double afac = alpha*dt/dx/dx;
-    int iter = (int) (timestep/dt);
-
-    Rcout << "dt: " << dt << "; inner iterations: " << iter << endl;
-
-
-    for (int it = 0; it < iter; it++){
-        for (int i = 1; i < ext.size()-1; i++) {
-            ext[i] = (1. - 2*afac)*ext[i] + afac*(ext[i+1]+ext[i-1]);
-        }
-        ext[0] = bc_left;
-        ext[n-1] = bc_right;
-    }
-    return(ext);
-}
-

examples/_deps/doctest-src

@@ -0,0 +1 @@
+Subproject commit b7c21ec5ceeadb4951b00396fc1e4642dd347e5f

examples/main_1D.cpp

@@ -1,55 +0,0 @@
-#include "diffusion/BTCSBoundaryCondition.hpp"
-#include <diffusion/BTCSDiffusion.hpp>
-
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-using namespace std;
-
-using namespace Diffusion;
-
-int main(int argc, char *argv[]) {
-
-  // dimension of grid
-  int dim = 1;
-
-  int n = 20;
-
-  // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n, 1e-1);
-  std::vector<double> field(n, 1e-6);
-
-  BTCSBoundaryCondition bc(n);
-
-  // create instance of diffusion module
-  BTCSDiffusion diffu(dim);
-
-  diffu.setXDimensions(1, n);
-
-  // set the boundary condition for the left ghost cell to dirichlet
-  bc(BC_SIDE_LEFT) = {BC_TYPE_CONSTANT, 5e-6};
-  // bc[0] = {Diffusion::BC_CONSTANT, 5e-6};
-  // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT,
-  //                            5. * std::pow(10, -6));
-
-  // set timestep for simulation to 1 second
-  diffu.setTimestep(1.);
-
-  cout << setprecision(12);
-
-  // loop 100 times
-  // output is currently generated by the method itself
-  for (int i = 0; i < 100; i++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
-
-    cout << "Iteration: " << i << "\n\n";
-
-    for (auto &cell : field) {
-      cout << cell << "\n";
-    }
-
-    cout << "\n" << endl;
-  }
-
-  return 0;
-}

examples/main_2D.cpp

@@ -1,58 +0,0 @@
-#include "diffusion/BTCSBoundaryCondition.hpp"
-#include <diffusion/BTCSDiffusion.hpp>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-using namespace std;
-
-using namespace Diffusion;
-
-int main(int argc, char *argv[]) {
-
-  // dimension of grid
-  int dim = 2;
-
-  int n = 5;
-  int m = 5;
-
-  // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n * m, 1e-6);
-  std::vector<double> field(n * m, 0);
-  BTCSBoundaryCondition bc(n, m);
-
-  // create instance of diffusion module
-  BTCSDiffusion diffu(dim);
-
-  diffu.setXDimensions(n, n);
-  diffu.setYDimensions(m, m);
-
-  // set inlet to higher concentration without setting bc
-  field[12] = 1;
-
-  // set timestep for simulation to 1 second
-  diffu.setTimestep(1);
-
-  cout << setprecision(12);
-
-  for (int t = 0; t < 1000; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
-
-    cout << "Iteration: " << t << "\n\n";
-
-    double sum = 0;
-
-    // iterate through rows
-    for (int i = 0; i < m; i++) {
-      // iterate through columns
-      for (int j = 0; j < n; j++) {
-        cout << left << std::setw(20) << field[i * n + j];
-        sum += field[i * n + j];
-      }
-      cout << "\n";
-    }
-
-    cout << "sum: " << sum << "\n" << endl;
-  }
-
-  return 0;
-}

examples/main_2D_const.cpp

@@ -1,59 +0,0 @@
-#include "diffusion/BTCSBoundaryCondition.hpp"
-#include <diffusion/BTCSDiffusion.hpp>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-using namespace std;
-
-using namespace Diffusion;
-
-int main(int argc, char *argv[]) {
-
-  // dimension of grid
-  int dim = 2;
-
-  int n = 5;
-  int m = 5;
-
-  // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n * m, 1e-1);
-  std::vector<double> field(n * m, 1e-6);
-  BTCSBoundaryCondition bc(n, m);
-
-  // create instance of diffusion module
-  BTCSDiffusion diffu(dim);
-
-  diffu.setXDimensions(1, n);
-  diffu.setYDimensions(1, m);
-
-  boundary_condition input = {Diffusion::BC_TYPE_CONSTANT, 5e-6};
-
-  bc.setSide(BC_SIDE_LEFT, input);
-
-  // for (int i = 1; i <= n; i++) {
-  //   bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5e-6};
-  // }
-  // set timestep for simulation to 1 second
-  diffu.setTimestep(1.);
-
-  cout << setprecision(12);
-
-  for (int t = 0; t < 10; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
-
-    cout << "Iteration: " << t << "\n\n";
-
-    // iterate through rows
-    for (int i = 0; i < m; i++) {
-      // iterate through columns
-      for (int j = 0; j < n; j++) {
-        cout << left << std::setw(20) << field[i * n + j];
-      }
-      cout << "\n";
-    }
-
-    cout << "\n" << endl;
-  }
-
-  return 0;
-}

examples/main_2D_mdl.cpp

@@ -1,67 +0,0 @@
-#include "diffusion/BTCSBoundaryCondition.hpp"
-#include <diffusion/BTCSDiffusion.hpp>
-#include <iomanip>
-#include <iostream> // for std
-#include <vector>   // for vector
-using namespace std;
-
-using namespace Diffusion;
-
-int main(int argc, char *argv[]) {
-
-  // dimension of grid
-  int dim = 2;
-
-  int n = 501;
-  int m = 501;
-
-  // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n * m, 1e-3);
-  std::vector<double> field(n * m, 0.);
-  BTCSBoundaryCondition bc(n, m);
-
-  field[125500] = 1;
-
-  // create instance of diffusion module
-  BTCSDiffusion diffu(dim);
-
-  diffu.setXDimensions(1., n);
-  diffu.setYDimensions(1., m);
-
-  // set the boundary condition for the left ghost cell to dirichlet
-  // diffu.setBoundaryCondition(250, 250, BTCSDiffusion::BC_CONSTANT, 1);
-  // for (int d=0; d<5;d++){
-  //     diffu.setBoundaryCondition(d, 0, BC_CONSTANT, .1);
-  // }
-  // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1);
-  // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1);
-
-  // set timestep for simulation to 1 second
-  diffu.setTimestep(1.);
-
-  cout << setprecision(7);
-
-  // First we output the initial state
-  cout << 0;
-
-  for (int i = 0; i < m * n; i++) {
-    cout << "," << field[i];
-  }
-  cout << endl;
-
-  // Now we simulate and output 8 steps à 1 sec
-  for (int t = 1; t < 6; t++) {
-    double time = diffu.simulate(field.data(), alpha.data(), bc);
-
-    cerr << "time elapsed: " << time << " seconds" << endl;
-
-    cout << t;
-
-    for (int i = 0; i < m * n; i++) {
-      cout << "," << field[i];
-    }
-    cout << endl;
-  }
-
-  return 0;
-}

examples/profiling_openmp.cpp

@@ -0,0 +1,70 @@
+#include <Eigen/Eigen>
+#include <chrono>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <tug/Simulation.hpp>
+
+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<std::chrono::milliseconds>(end -
+                                                                    begin);
+          myfile << milliseconds.count() << endl;
+        }
+      }
+      cout << endl;
+      myfile << endl;
+    }
+    myfile.close();
+  }
+}

examples/profiling_speedup.cpp

@@ -0,0 +1,70 @@
+#include <Eigen/Eigen>
+#include <chrono>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <tug/Simulation.hpp>
+
+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<std::chrono::milliseconds>(end -
+                                                                    begin);
+          myfile << milliseconds.count() << endl;
+        }
+      }
+      cout << endl;
+      myfile << endl;
+    }
+    myfile.close();
+  }
+}

examples/reference-FTCS_2D_closed.cpp

@@ -0,0 +1,52 @@
+#include "Eigen/Core"
+#include <iostream>
+#include <tug/Simulation.hpp>
+
+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;
+
+  // 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 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 = 5; i < 11; i++) {
+    for (int j = 6; j < 11; j++) {
+      alpha(i, j) = 0.1;
+    }
+  }
+  grid.setAlpha(alpha, alpha);
+
+  // 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);
+
+  // RUN
+  sim.run();
+}

include/tug/Boundary.hpp

@@ -0,0 +1,211 @@
+/**
+ * @file Boundary.hpp
+ * @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 "Grid.hpp"
+#include <cstddef>
+
+/**
+ * @brief Enum defining the two implemented boundary conditions.
+ *
+ */
+enum BC_TYPE { BC_TYPE_CLOSED, BC_TYPE_CONSTANT };
+
+/**
+ * @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 };
+
+/**
+ * This class defines the boundary conditions of individual boundary elements.
+ * These can be flexibly used and combined later in other classes.
+ * 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 -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 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 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 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 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 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 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<BoundaryElement> Contains the boundary conditions as
+   * BoundaryElement objects.
+   */
+  const std::vector<BoundaryElement> 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.
+   */
+  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 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);
+
+private:
+  Grid grid; // Boundary is directly dependent on the dimensions of a predefined
+
+  std::vector<std::vector<BoundaryElement>>
+      boundaries; // Vector with Boundary Element information
+};
+
+#endif // BOUNDARY_H_

include/tug/BoundaryCondition.hpp

@@ -1,278 +0,0 @@
-#ifndef BOUNDARYCONDITION_H_
-#define BOUNDARYCONDITION_H_
-
-#include <array>
-#include <map>
-#include <stdexcept>
-#include <stdint.h>
-#include <vector>
-
-typedef uint8_t bctype;
-
-namespace tug {
-namespace bc {
-
-enum {
-  BC_TYPE_CLOSED,   /**< Defines a closed/Neumann boundary condition. */
-  BC_TYPE_FLUX,     /**< Defines a flux/Cauchy boundary condition. */
-  BC_TYPE_CONSTANT, /**< Defines a constant/Dirichlet boundary condition. */
-  BC_UNSET          /**< Indicates undefined boundary condition*/
-};
-
-enum {
-  BC_SIDE_LEFT,  /**< Defines boundary conditions for the left side of the grid.
-                  */
-  BC_SIDE_RIGHT, /**< Defines boundary conditions for the right side of the
-                    grid. */
-  BC_SIDE_TOP,   /**< Defines boundary conditions for the top of the grid. */
-  BC_SIDE_BOTTOM, /**< Defines boundary conditions for the bottom of the grid.
-                   */
-  BC_INNER
-};
-
-/**
- * Defines the boundary condition type and according value.
- */
-typedef struct boundary_condition_s {
-  bctype type;  /**< Type of the boundary condition */
-  double value; /**< Value of the boundary condition. Either a concrete value of
-                   concentration for BC_TYPE_CONSTANT or gradient when type is
-                   BC_TYPE_FLUX. Unused if BC_TYPE_CLOSED.*/
-} boundary_condition;
-
-/**
- * Represents both boundary conditions of a row/column.
- */
-typedef std::array<boundary_condition, 2> bc_tuple;
-typedef std::vector<boundary_condition> bc_vec;
-
-/**
- * Class to define the boundary condition of a grid.
- */
-class BoundaryCondition {
-public:
-  /**
-   *  Creates a new instance with two elements. Used when defining boundary
-   *  conditions of 1D grids.
-   *
-   *  \param x Number of grid cells in x-direction
-   */
-  BoundaryCondition(int x);
-
-  /**
-   * Creates a new instance with 4 * max(x,y) elements. Used to describe the
-   * boundary conditions for 2D grids. NOTE: On non-squared grids there are more
-   * elements than needed and no exception is thrown if some index exceeding
-   * grid limits.
-   *
-   * \param x Number of grid cells in x-direction
-   * \param y Number of grid cells in y-direction
-   *
-   */
-  BoundaryCondition(int x, int y);
-
-  /**
-   * Sets the boundary condition for a specific side of the grid.
-   *
-   * \param side Side for which the given boundary condition should be set.
-   * \param input_bc Instance of struct boundary_condition with desired boundary
-   * condition.
-   *
-   * \throws std::invalid_argument Indicates wrong dimensions of the grid.
-   * \throws std::out_of_range Indicates a out of range value for side.
-   */
-  void setSide(uint8_t side, boundary_condition &input_bc);
-
-  /**
-   * Sets the boundary condition for a specific side of the grid.
-   *
-   * \param side Side for which the given boundary condition should be set.
-   * \param input_bc Vector of boundary conditions for specific side.
-   *
-   * \throws std::invalid_argument Indicates wrong dimensions of the grid.
-   * \throws std::out_of_range Indicates a out of range value for side or
-   * invalid size of input vector.
-   */
-  void setSide(uint8_t side, std::vector<boundary_condition> &input_bc);
-
-  /**
-   * Returns a vector of boundary conditions of given side. Can be used to set
-   * custom boundary conditions and set back via setSide() with vector input.
-   *
-   * \param side Side which boundary conditions should be returned
-   *
-   * \returns Vector of boundary conditions
-   *
-   * \throws std::invalid_argument If given dimension is less or equal to 1.
-   * \throws std::out_of_range Indicates a out of range value for side.
-   */
-  auto getSide(uint8_t side) -> std::vector<boundary_condition>;
-
-  /**
-   * Get both boundary conditions of a given row (left and right).
-   *
-   * \param i Index of row
-   *
-   * \returns Left and right boundary values as an array (defined as data
-   * type bc_tuple).
-   */
-  auto row_boundary(uint32_t i) const -> bc_tuple;
-
-  /**
-   * Get both boundary conditions of a given column (top and bottom).
-   *
-   * \param i Index of column
-   *
-   * \returns Top and bottom boundary values as an array (defined as data
-   * type bc_tuple).
-   */
-  auto col_boundary(uint32_t i) const -> bc_tuple;
-
-  /**
-   * Sets one cell of the inner grid to a given boundary condition.
-   *
-   * \param bc Boundary condition to be set.
-   * \param x Index of grid cell in x direction.
-   * \param y Index of grid cell in y direction.
-   */
-  void setInnerBC(boundary_condition bc, int x, int y);
-
-  /**
-   * Unsets a previously set inner boundary condition.
-   *
-   * \param x Index of grid cell in x direction.
-   * \param y Index of grid cell in y direction.
-   */
-  void unsetInnerBC(int x, int y);
-
-  /**
-   * Returns the current boundary condition set for specific inner grid cell.
-   *
-   * \param x Index of grid cell in x direction.
-   * \param y Index of grid cell in y direction.
-   */
-  auto getInnerBC(int x, int y) -> boundary_condition;
-
-  /**
-   * Get a row of field and its inner boundary conditions.
-   *
-   * \param i Index of the row starting at 0.
-   *
-   * \returns Row of the inner boundary conditions of the field.
-   */
-  auto getInnerRow(uint32_t i) const -> bc_vec;
-
-  /**
-   * Get a column of field and its inner boundary conditions.
-   *
-   * \param i Index of the column starting at 0.
-   *
-   * \returns Column of the inner boundary conditions of the field.
-   */
-  auto getInnerCol(uint32_t i) const -> bc_vec;
-
-  /**
-   * Create an instance of boundary_condition data type. Can be seen as a helper
-   * function.
-   *
-   * \param type Type of the boundary condition.
-   * \param value According value of condition.
-   *
-   * \returns Instance of boundary_condition
-   */
-  static boundary_condition returnBoundaryCondition(bctype type, double value) {
-    return {type, value};
-  }
-
-private:
-  std::vector<boundary_condition> bc_internal;
-
-  std::map<uint32_t, boundary_condition> inner_cells;
-
-  uint8_t dim;
-
-  std::array<uint32_t, 2> sizes;
-  uint32_t maxsize;
-  uint32_t maxindex;
-
-  enum { X_DIM, Y_DIM };
-
-public:
-  /**
-   * Returns the left/right boundary condition for 1D grid.
-   *
-   * \param side Side of the boundary condition to get.
-   *
-   * \returns Boundary condition
-   */
-  boundary_condition operator()(uint8_t side) const {
-    if (dim != 1) {
-      throw std::invalid_argument(
-          "Only 1D grid support 1 parameter in operator");
-    }
-    if (side > 1) {
-      throw std::out_of_range("1D index out of range");
-    }
-    return bc_internal[side];
-  }
-
-  /**
-   * Returns the boundary condition of a given side for 2D grids.
-   *
-   * \param side Side of the boundary condition to get.
-   * \param i Index of the boundary condition.
-   *
-   * \returns Boundary condition
-   */
-  boundary_condition operator()(uint8_t side, uint32_t i) const {
-    if (dim != 2) {
-      throw std::invalid_argument(
-          "Only 2D grids support 2 parameters in operator");
-    }
-    if (side > 3) {
-      throw std::out_of_range("2D index out of range");
-    }
-    return bc_internal[side * maxsize + i];
-  }
-
-  /**
-   * Returns the left/right boundary condition for 1D grid.
-   *
-   * \param side Side of the boundary condition to get.
-   *
-   * \returns Boundary condition
-   */
-  boundary_condition &operator()(uint8_t side) {
-    if (dim != 1) {
-      throw std::invalid_argument(
-          "Only 1D grid support 1 parameter in operator");
-    }
-    if (side > 1) {
-      throw std::out_of_range("1D index out of range");
-    }
-    return bc_internal[side];
-  }
-
-  /**
-   * Returns the boundary condition of a given side for 2D grids.
-   *
-   * \param side Side of the boundary condition to get.
-   * \param i Index of the boundary condition.
-   *
-   * \returns Boundary condition
-   */
-  boundary_condition &operator()(uint8_t side, uint32_t i) {
-    if (dim != 2) {
-      throw std::invalid_argument("Explicit setting of bc value with 2 "
-                                  "parameters is only supported for 2D grids");
-    }
-    if (side > 3) {
-      throw std::out_of_range("2D index out of range");
-    }
-    return bc_internal[side * maxsize + i];
-  }
-};
-
-} // namespace bc
-} // namespace tug
-#endif // BOUNDARYCONDITION_H_

include/tug/Diffusion.hpp

@@ -1,138 +0,0 @@
-#ifndef DIFFUSION_H_
-#define DIFFUSION_H_
-
-#include "BoundaryCondition.hpp"
-#include "Solver.hpp"
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-#include <array>
-#include <vector>
-
-namespace tug {
-namespace diffusion {
-
-constexpr uint8_t MAX_ARR_SIZE = 3;
-
-/**
- * Defines grid dimensions and boundary conditions.
- */
-typedef struct tug_grid_s {
-  std::array<uint32_t, MAX_ARR_SIZE>
-      grid_cells; /**< Count of grid cells in each of the 3 directions.*/
-  std::array<double, MAX_ARR_SIZE>
-      domain_size; /**< Domain sizes in each of the 3 directions.*/
-  bc::BoundaryCondition *bc = NULL; /**< Boundary conditions for the grid.*/
-} TugGrid;
-
-/**
- * Besides containing the grid structure it holds also information about the
- * desired time step to simulate and which solver to use.
- */
-typedef struct tug_input_s {
-  double time_step; /**< Time step which should be simulated by diffusion.*/
-  Eigen::VectorXd (*solver)(const Eigen::SparseMatrix<double> &,
-                            const Eigen::VectorXd &) =
-      tug::solver::ThomasAlgorithm; /**< Solver function to use.*/
-  TugGrid grid;                     /**< Grid specification.*/
-
-  /**
-   * Set the desired time step for diffusion simulation.
-   *
-   * \param dt Time step in seconds.
-   */
-  void setTimestep(double dt) { time_step = dt; }
-
-  /**
-   * Set the count of grid cells in each dimension.
-   *
-   * \param x Count of grid cells in x direction.
-   * \param y Count of grid cells in y direction.
-   * \param z Count of grid cells in z direction.
-   */
-  void setGridCellN(uint32_t x, uint32_t y = 0, uint32_t z = 0) {
-    grid.grid_cells[0] = x;
-    grid.grid_cells[1] = y;
-    grid.grid_cells[2] = z;
-  }
-
-  /**
-   * Set the domain size of the grid in each direction.
-
-   * \param Domain size in x direction.
-   * \param Domain size in y direction.
-   * \param Domain size in z direction.
-   */
-  void setDomainSize(double x, double y = 0, double z = 0) {
-    grid.domain_size[0] = x;
-    grid.domain_size[1] = y;
-    grid.domain_size[2] = z;
-  }
-
-  /**
-   * Set boundary conditions for grid instance.
-   *
-   * \param bc Boundary conditions to be set.
-   */
-  void setBoundaryCondition(bc::BoundaryCondition &bc) { grid.bc = &bc; }
-
-  /**
-   * Retrieve the set boundary condition from grid instance.
-   *
-   * \return Boundary condition object if boundary conditions were set,
-   * otherwise NULL.
-   */
-  auto getBoundaryCondition() -> bc::BoundaryCondition { return *(grid.bc); }
-
-  /**
-   * Set the solver function.
-   *
-   * \param f_in Pointer to function which takes a sparse matrix and a vector as
-   * input and returns another vector.
-   */
-  void
-  setSolverFunction(Eigen::VectorXd (*f_in)(const Eigen::SparseMatrix<double> &,
-                                            const Eigen::VectorXd &)) {
-    solver = f_in;
-  }
-} TugInput;
-
-/**
- * Solving 1D diffusion problems with Backward Time Centred Space scheme.
- *
- * \param input_param Object with information for the simulation e.g. grid
- * dimensions, time step ...
- *
- * \param field Pointer to continious memory holding the concentrations for each
- * grid cell of the grid. It doesn't matter if stored in row (most likely) or
- * column major.
- *
- * \param alpha Pointer to continious memory holding the alpha for each grid
- * cell of the grid. (NOTE: only constant alphas are supported currently)
- *
- * \return Runtime in seconds
- */
-auto BTCS_1D(const TugInput &input_param, double *field, const double *alpha)
-    -> double;
-
-/**
- * Solving 2D diffusion problems with Alternating-direction implicit method.
- *
- * \param input_param Object with information for the simulation e.g. grid
- * dimensions, time step ...
- *
- * \param field Pointer to continious memory holding the concentrations for each
- * grid cell of the grid. It doesn't matter if stored in row (most likely) or
- * column major.
- *
- * \param alpha Pointer to continious memory holding the alpha for each grid
- * cell of the grid. (NOTE: only constant alphas are supported currently)
- *
- * \return Runtime in seconds
- */
-auto ADI_2D(const TugInput &input_param, double *field, const double *alpha)
-    -> double;
-
-} // namespace diffusion
-} // namespace tug
-
-#endif // DIFFUSION_H_

include/tug/Grid.hpp

@@ -0,0 +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.
+ *
+ */
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+class Grid {
+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 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 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 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 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 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 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 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 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 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 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 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 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 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_

include/tug/Simulation.hpp

@@ -0,0 +1,236 @@
+/**
+ * @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.
+ *
+ */
+
+#ifndef SIMULATION_H_
+#define SIMULATION_H_
+
+#include "Boundary.hpp"
+#include "Grid.hpp"
+#include <string>
+#include <vector>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_num_procs() 1
+#endif
+
+/**
+ * @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
+};
+
+/**
+ * @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
+};
+
+/**
+ * @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
+};
+
+/**
+ * @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
+};
+
+/**
+ * @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
+};
+
+/**
+ * @brief The class forms the interface for performing the diffusion simulations
+ * and contains all the methods for controlling the desired parameters, such as
+ * time step, number of simulations, etc.
+ *
+ */
+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. 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 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 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 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 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 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 Outputs the current concentrations of the grid on the console.
+   *
+   */
+  void printConcentrationsConsole();
+
+  /**
+   * @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:
+   *        <approach> + <number rows> + <number columns> + <number of
+   * iterations>+<counter>.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(const std::string &filename);
+
+  /**
+   * @brief Method starts the simulation process with the previously set
+   *        parameters.
+   */
+  void run();
+
+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};
+
+  Grid &grid;
+  Boundary &bc;
+  APPROACH approach;
+  SOLVER solver{THOMAS_ALGORITHM_SOLVER};
+
+  const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
+};
+
+#endif // SIMULATION_H_

include/tug/Solver.hpp

@@ -1,41 +0,0 @@
-#ifndef SOLVER_H_
-#define SOLVER_H_
-
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-
-namespace tug {
-namespace solver {
-
-/**
- * Solving linear equation system with LU-decomposition implemented in Eigen
- * library.
- *
- * \param A_matrix The A matrix represented as a sparse matrix using Eigen
- * library.
- *
- * \param b_vector Right hand side vector of the linear equation system.
- *
- * \return Solution represented as vector.
- */
-auto EigenLU(const Eigen::SparseMatrix<double> &A_matrix,
-             const Eigen::VectorXd &b_vector) -> Eigen::VectorXd;
-
-/**
- * Solving linear equation system with brutal implementation of the Thomas
- * algorithm.
- *
- * \param A_matrix The A matrix represented as a sparse matrix using Eigen
- * library.
- *
- * \param b_vector Right hand side vector of the linear equation system.
- *
- * \return Solution represented as vector.
- */
-auto ThomasAlgorithm(const Eigen::SparseMatrix<double> &A_matrix,
-                     const Eigen::VectorXd &b_vector) -> Eigen::VectorXd;
-
-} // namespace solver
-} // namespace tug
-
-#endif // SOLVER_H_

scripts/Adi2D_Reference.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-07-20 15:37:25 delucia"
+## Time-stamp: "Last modified 2023-07-31 14:26:49 delucia"
 
 ## Brutal implementation of 2D ADI scheme
 ## Square NxN grid with dx=dy=1

scripts/HetDiff.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-07-21 11:23:08 delucia"
+## Time-stamp: "Last modified 2023-07-31 16:28:48 delucia"
 
 library(ReacTran)
 library(deSolve)
@@ -99,15 +99,25 @@ times <- 0:10
 outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL,
                 dim = c(N, N), lrw = 1860000)
 
-outtimes <- c(1, 4, 7, 10)
+outtimes <- c(0, 4, 7, 10)
 
-cairo_pdf("deSolve_AlphaHet1.pdf", family="serif", width=12, height=12)
+## NB: assuming current working dir is "tug"
+cairo_pdf("doc/images/deSolve_AlphaHet1.pdf", family="serif", width=12, height=12)
 image(outc, ask = FALSE, mfrow = c(2, 2), main = paste("time", outtimes),
-      legend = TRUE, add.contour = FALSE, subset = time %in% outtimes, xlab="",ylab="", axes=FALSE)
+      legend = TRUE, add.contour = FALSE, subset = time %in% outtimes,
+      xlab="",ylab="", axes=FALSE, asp=1)
 dev.off()
 
-outc
+## outc is a matrix with 11 rows and 122 columns (first column is
+## simulation time);
+str(outc)
 
+## extract only the results and transpose the matrix for storage
+ret <- data.matrix(t(outc[ , -1]))
+rownames(ret) <- NULL
+
+## NB: assuming current working dir is "tug"
+data.table::fwrite(ret, file="scripts/gold/HetDiff1.csv", col.names=FALSE)
 
 
 

scripts/gold/HetDiff1.csv

@@ -0,0 +1,121 @@
+0,1.15723009391899e-05,0.000573422585977779,0.00271483227696389,0.00592017488450261,0.00919024372055321,0.0119992110381186,0.01421611640169,0.0158871549414673,0.0171115136326332,0.0179901833258047
+0,4.32496365970161e-05,0.00110003761455522,0.0038588402926536,0.00723956931728921,0.0103665385368182,0.0129275133140941,0.0149018727052112,0.0163722917415511,0.0174425542133075,0.0182067473313329
+0,0.000156410551955808,0.0022501860978233,0.00596846775043427,0.00951979828485079,0.0123439077799255,0.0144650737976258,0.0160242190667615,0.0171550743320483,0.0179654704632021,0.0185369538227271
+0,0.000449900131743341,0.00397802650798475,0.0085351147388942,0.0120552559726064,0.0144482796517419,0.0160547847343675,0.0171511707942764,0.0179099936829702,0.0184375120234885,0.0188001645647094
+0,0.00096247784605253,0.00575161551352107,0.0106416083625011,0.0139078265679269,0.0158706072121663,0.0170493964826529,0.0177834857154376,0.0182584391833927,0.018573228563985,0.0187811039627322
+0,0.00139148060908994,0.00653728609085073,0.0112226963812925,0.0141895931925161,0.0159026822227563,0.0168925530607027,0.0174883152516615,0.017865325174865,0.0181135647391488,0.0182789062492789
+0,7.88755181655377e-05,0.000890855516524495,0.0025489601065516,0.00451888907818989,0.00640234215109739,0.00803400181298965,0.00938380805959316,0.0104782803692741,0.0113602330285538,0.0120715601783272
+0,1.96835046894113e-06,5.20032828635902e-05,0.000248131836356188,0.000630746867714403,0.00117274437938983,0.00182064447664627,0.00252211668991391,0.00323643607872941,0.00393564702867779,0.00460249317726503
+0,4.05263620284678e-08,2.41553600349545e-06,1.87665958722685e-05,6.74081501157623e-05,0.000163264089823103,0.000313023273063137,0.000515631508985307,0.000764816320228696,0.00105167332970099,0.00136658044925295
+0,7.12888460576857e-10,9.36408593803959e-08,1.16671799220012e-06,5.86484773654181e-06,1.84017858631692e-05,4.34716444970184e-05,8.51739231661631e-05,0.000146375097555955,0.000228524220562364,0.000331765311847696
+0,1.11266213313442e-11,3.21569209939452e-09,6.45997582539074e-08,4.58674913815976e-07,1.88637915450085e-06,5.56437226690733e-06,1.31531030829188e-05,2.6577288329037e-05,4.78110315422556e-05,7.86848274922097e-05
+0,4.35985871981009e-05,0.0011173013784839,0.00393465640259938,0.00738648425859891,0.0105643563174128,0.0131485926879987,0.0151246465122267,0.016583690239707,0.0176362908990393,0.0183809352240671
+0,0.000162940275870536,0.00214415472789055,0.0055971341202626,0.00904266600518215,0.0119317484324307,0.01418453738493,0.0158744170951516,0.0171103507405317,0.0179967639568099,0.0186199630282087
+0,0.000589272728343116,0.00438867743664915,0.00866893626959719,0.0119151946783314,0.0142428604635392,0.0159134847420158,0.0171142649822806,0.0179721981758363,0.0185775454158148,0.0189955063245207
+0,0.00169526769334991,0.00776694628989879,0.0124247154337811,0.0151382671118746,0.0167367107498905,0.0177366662921586,0.0183940071716407,0.0188368236363926,0.0191342578330509,0.0193277227804878
+0,0.00362839452566188,0.0112526737744286,0.0155521246084951,0.0175590393037948,0.0184985843185778,0.0189571544442947,0.0191923705608016,0.0193168273546316,0.0193791769704849,0.0194019806796609
+0,0.00525177970919126,0.0128429513741672,0.0165189505188843,0.0180778870279157,0.0187197075826292,0.018968965691842,0.0190518322037241,0.0190652772723224,0.0190482483782045,0.019016041284815
+0,0.000385153961868949,0.00229577940499986,0.00480460165448024,0.00711204718132809,0.0089909576548431,0.0104472162269959,0.0115548772926463,0.0123946614477843,0.0130350658282232,0.013528343285464
+0,1.16858258297074e-05,0.000164584714174317,0.000570116218846185,0.00119271144355875,0.00194723110061389,0.00275621955267251,0.00356441433332025,0.00433765719695555,0.00505760032168857,0.00571654514201105
+0,2.81208380577636e-07,8.97423573740115e-06,5.03422264992399e-05,0.000147427646174569,0.000310274427544454,0.000536832367026033,0.000817602280027168,0.00113989443932459,0.00149069663261676,0.00185828031564749
+0,5.63399964199174e-09,3.96964588417007e-07,3.55564115271452e-06,1.44716592579598e-05,3.91537025779409e-05,8.28348829677685e-05,0.000148977185800308,0.000239043494747069,0.00035275053097983,0.000488513445484224
+0,9.84398346768558e-11,1.5293568296117e-08,2.20334810416761e-07,1.26108145384204e-06,4.44865566281344e-06,1.16871019830672e-05,2.52195633384727e-05,4.73253788486059e-05,8.00591501276681e-05,0.000125072519296284
+0,0.00015939225040812,0.00233357034425003,0.00625562530069953,0.0100135877171516,0.0129676882444301,0.0151374162693059,0.0166874294020243,0.0177764984876221,0.0185308441084237,0.0190433426894296
+0,0.000595680647430986,0.0044801608514742,0.00890743169921897,0.0122772025886856,0.0146734676860174,0.0163633416618241,0.0175504470842378,0.0183769265178609,0.018943746718196,0.019322547574371
+0,0.00215448195385171,0.0091767592756406,0.0138199838635887,0.0162232592499362,0.0175799863177681,0.0184333811027871,0.0190005234819769,0.0193808295698239,0.0196287328140203,0.01977995018463
+0,0.00619992143295609,0.0162622253524985,0.0198667054266907,0.0207089785376823,0.0207820327860086,0.0206824193920932,0.0205605370395911,0.0204471756023412,0.0203412141522796,0.0202389042951594
+0,0.0132779136133337,0.0236231555894036,0.0250044710521133,0.0242144479160946,0.0231938485225262,0.0223380815152074,0.0216793603038312,0.0211797341820616,0.020794392809717,0.0204895366674834
+0,0.0192489696423083,0.0271184370488408,0.026840217635384,0.0252838506050641,0.023849962885927,0.0227245048675836,0.0218698389880889,0.0212219274277141,0.0207240679736698,0.020334308600127
+0,0.00188379770767293,0.00657190925731914,0.0104892667525962,0.0130700680836906,0.0146349533533907,0.015535980568487,0.0160269602381333,0.0162728756742906,0.0163767130755721,0.016400710298742
+0,7.01829175866307e-05,0.000581093706531354,0.00153370940121776,0.00269094081762158,0.00386803514667116,0.00496316154582596,0.00593261862345271,0.00676642929553236,0.0074719703393928,0.00806413175855109
+0,1.97847089603411e-06,3.69925955384929e-05,0.000157406497985584,0.000384860592255597,0.000709846476804252,0.0011079382186538,0.00155170108970297,0.00201681216519113,0.00248426641823041,0.00294058581228933
+0,4.51078087575571e-08,1.85104501980116e-06,1.25057013282593e-05,4.22892094350572e-05,9.98395815944757e-05,0.000189787308802236,0.000312690644970973,0.000466034643916367,0.000645463690913993,0.000845828528229199
+0,8.80252798763183e-10,7.92215439962928e-08,8.56765592551646e-07,4.05803525996784e-06,1.2449180396882e-05,2.9295973763362e-05,5.7746988001425e-05,0.0001003733295813,0.000158928915289496,0.00023429575744269
+0,0.000466368704740309,0.00424791180746615,0.00926944539909274,0.0131713216643424,0.0157625616947689,0.0174133722552654,0.0184577152890545,0.0191156877417616,0.0195250042740468,0.0197700479556243
+0,0.00174300015594418,0.00815930380362134,0.0132116161354557,0.016172754442502,0.0178696343095053,0.0188636958720877,0.0194550603754561,0.0198042367684176,0.0200009359043926,0.0200979478510297
+0,0.00630533087106852,0.0167264504542955,0.0205352907522532,0.021433296617797,0.0214908287288404,0.021342901494693,0.0211590163834043,0.0209808421912509,0.020813382900871,0.0206559438303523
+0,0.0181516683021334,0.0296877506877057,0.0296188711707772,0.0274991387847508,0.0255697170579367,0.0241218925670229,0.0230708035427124,0.0223012794968431,0.0217225548818433,0.0212745536462067
+0,0.0389067518217964,0.0432743493327086,0.037525495598965,0.0324527348417642,0.0288543692365208,0.0263663031569143,0.0246230583185604,0.0233728213888174,0.0224524193870472,0.0217568092848378
+0,0.0565396239160925,0.0500897261990562,0.0408387347867424,0.0344846493850116,0.0302543240032215,0.0273643802632421,0.0253287379325492,0.0238539136707014,0.0227584403765214,0.02192583980509
+0,0.00794755170946868,0.0172628249510397,0.0219912124473273,0.0237906919799591,0.0240591222673087,0.0235982544067714,0.0228369786001951,0.021996032718532,0.0211823304839058,0.0204423856134907
+0,0.000373684524467729,0.00190283496658901,0.00395283626547245,0.00595433814031928,0.00766198835968331,0.0090150942497423,0.0100370905098692,0.0107818312810106,0.0113083452714233,0.0116697843722264
+0,1.24794246503553e-05,0.000141270114233274,0.000466648605475989,0.000969170133300241,0.00158727378005724,0.00225743530947124,0.00292973311195242,0.00357075743939607,0.00416127927601125,0.00469265749237697
+0,3.25056776407433e-07,7.95461921640221e-06,4.12089493515841e-05,0.000117250753288648,0.00024402042254145,0.000420314286347587,0.000638908995931823,0.000889809192557268,0.00116261453581859,0.00144792147140295
+0,7.0854379238238e-09,3.75524392933622e-07,3.08321489667331e-06,1.21925554449326e-05,3.27749649029249e-05,6.95563221553643e-05,0.000126018794389096,0.000204052811537285,0.000304018569398211,0.000425019521621274
+0,0.00102836224521133,0.0064274670666241,0.0121491816508084,0.015968397765568,0.0181492003638804,0.019314880428558,0.0199102032969297,0.0201928751968812,0.0203044557297043,0.020320071184116
+0,0.00384385388740204,0.0123515070854559,0.0173313777131353,0.0196324133235912,0.0206081184457373,0.0209612687545928,0.0210258040457677,0.0209599611581841,0.0208374482363417,0.0206927731773186
+0,0.0139084634226096,0.0253420008838306,0.0269858208570455,0.0260863431776468,0.0248662406090429,0.023805495300228,0.0229587141516884,0.0222941985053244,0.0217679373130175,0.021344966099159
+0,0.0400577489259572,0.0450598366852064,0.0390556903095381,0.0336304680408348,0.0297597604093353,0.0270807237205804,0.0252035706273644,0.0238573134643687,0.0228662562144913,0.0221177566182138
+0,0.0859596622787531,0.0659645822190159,0.0498440895724295,0.0400638374763904,0.0339433683852068,0.0299353070092985,0.0272039578032914,0.0252772077773997,0.0238777953326208,0.0228341669415018
+0,0.125398084679699,0.0772427043755313,0.0551767906604444,0.0434193651923798,0.0363278914606691,0.0317025734216055,0.0285266548704347,0.0262622056865752,0.0246012510386757,0.0233529211513904
+0,0.0279719026476077,0.0399912836159061,0.041824217737992,0.0399960375533372,0.0370306079722165,0.0339334669014817,0.0310904350363235,0.0286238301128732,0.0265419722346029,0.0248094228214593
+0,0.0016878543861307,0.00546802879678612,0.00908467643641864,0.0118685267101036,0.0137867836918316,0.014992473171871,0.0156668769769553,0.0159679687584234,0.0160191464240754,0.0159109388085067
+0,6.66210210366172e-05,0.000466311941721689,0.0012074539008488,0.00214399617569523,0.00313678410870533,0.0040909928004654,0.0049535861902635,0.00570178710489221,0.00633199242057907,0.00685150908457396
+0,1.96581789972705e-06,2.90655193631781e-05,0.000116181708737399,0.000279461298722284,0.000515297156287373,0.00080891264649635,0.00114182790936584,0.00149639884138984,0.0018579271762422,0.00221530254986257
+0,4.73511582533572e-08,1.48929389156945e-06,9.32003158799015e-06,3.08850626482769e-05,7.30754947692817e-05,0.000140629698309681,0.000235669226477032,0.00035796975388455,0.000505546298234415,0.000675276938000184
+0,0.0015836697501133,0.00789553113761796,0.0138353209275783,0.0174900423844659,0.0193829171959464,0.0202511267409907,0.0205839676229735,0.0206522039767607,0.0205960670144616,0.0204840229205406
+0,0.00592024795026366,0.015178056357401,0.0197491887189647,0.0215219374447417,0.0220321907378727,0.0220033162511111,0.0217646356223314,0.0214641136669779,0.0211632832330552,0.0208849932589192
+0,0.0214261698026196,0.0311627748380295,0.030789710824671,0.0286485315817415,0.0266434190236869,0.0250513323715184,0.023828703885917,0.0228921858496036,0.0221672538845208,0.0215984913078681
+0,0.0617381544584604,0.0554953966664163,0.0446769294750762,0.0370604778930011,0.0320153173963507,0.0286236282907804,0.0262786287048439,0.0246098514264461,0.0233900725628177,0.0224762767353041
+0,0.132662095447386,0.0815715457115294,0.0573620962224714,0.0444656436513204,0.0367970084501995,0.0318900160341415,0.0285856373696197,0.0262711290074283,0.0245996244789251,0.0233597987792184
+1,0.194600594053136,0.0967188378926377,0.0645033312825157,0.0489815255136541,0.0400071807557426,0.0342746856436233,0.0303863403322649,0.0276362424095596,0.0256323186910777,0.0241347191227446
+0,0.072774737081717,0.0731608818750184,0.0649513569830106,0.0562393519385766,0.048657668424965,0.0424321620057059,0.0374285709142598,0.0334382577022264,0.0302593069735787,0.0277213581035032
+0,0.00546229967617301,0.0117769944796152,0.0160960402630076,0.0186595966553889,0.0199729770538671,0.020458608892436,0.0204223397121343,0.0200740130964532,0.0195534451842866,0.0189514347157146
+0,0.000246612705267533,0.00110781620007196,0.00231152127659001,0.00359442519681962,0.00480227348450595,0.00586145140404475,0.00674755083035105,0.00746338966694999,0.00802538295172145,0.00845539933096027
+0,8.01024172201437e-06,7.40436516662321e-05,0.000234915899100426,0.000490064251651924,0.000819680663206465,0.00119845837778001,0.00160264998868245,0.00201295787243647,0.00241514785548912,0.00279964947569617
+0,2.07955309779492e-07,4.0104486961948e-06,1.96923787162901e-05,5.61707039451707e-05,0.00011992118154755,0.000214121267381263,0.000338952769617425,0.000492335832038914,0.000670725912883758,0.000869802668989746
+0,9.15481631874651e-06,0.000117693874573294,0.000377019455101549,0.000741577814721888,0.00115788487944848,0.00159137882061179,0.00202342112731787,0.00244488663055897,0.00285171452271575,0.00324239719313011
+0,4.46105775530757e-05,0.000303984477552903,0.000717806961735601,0.00118873254936734,0.00166857393355206,0.00213667568842936,0.00258503807646796,0.00301133151201256,0.00341571030049681,0.00379926644315538
+0,0.000217726325331875,0.000873770408536327,0.00159096186837063,0.00224649712918517,0.00282908355295164,0.00334893383349626,0.00381769528011674,0.00424489367607894,0.00463786438443504,0.00500213301923337
+0,0.000918933831995098,0.00232435924877404,0.00343741337091249,0.00429315823548186,0.00497145196445075,0.00552787579177755,0.00599805857236328,0.00640506369246164,0.00676421367084299,0.00708600072229698
+0,0.00331122058413453,0.00568447123669356,0.00712540944284069,0.00809397304876091,0.00879084024979697,0.0093162979589995,0.00972657995173228,0.0100559006360597,0.010326180081839,0.0105521545670149
+0,0.0100043670600405,0.0127097661456827,0.0139937515689219,0.0147091809677894,0.0151257626062355,0.0153636070507909,0.0154869837094123,0.0155343002312855,0.0155299222027964,0.0154901124905561
+0,0.000136637116798085,0.000338905367123057,0.000531352383440882,0.000703792701127588,0.000856724436487964,0.000992875271716,0.00111515974854185,0.00122613935565325,0.00132792886232481,0.00142222359205319
+0,5.41197942669333e-06,2.66179151421015e-05,6.06837139232639e-05,0.00010289135389313,0.00014958545120669,0.00019826296968678,0.000247306756529671,0.000295722537183644,0.000342937360699378,0.000388655525576309
+0,1.72186371215188e-07,1.70696816511905e-06,5.79640506738781e-06,1.29195485396102e-05,2.30515604056451e-05,3.58992272238309e-05,5.10570165270148e-05,6.80982556623907e-05,8.66227559703461e-05,0.000106279145950884
+0,4.32364439009495e-09,8.64589599304522e-08,4.3955011662034e-07,1.29702112369114e-06,2.8633412989769e-06,5.28573284839916e-06,8.65071350621842e-06,1.29921291564899e-05,1.83031949023414e-05,2.45484723579094e-05
+0,9.12457420765116e-11,3.74315585927427e-09,2.9038454685416e-08,1.15682654773554e-07,3.22205377431422e-07,7.18623792678699e-07,1.37865133447057e-06,2.37328126021504e-06,3.76631554161963e-06,5.61170605601181e-06
+0,2.08011169409904e-08,6.36518490608328e-07,3.48383663648875e-06,1.01031069547314e-05,2.13412103139757e-05,3.74893709216357e-05,5.85171433417359e-05,8.42328798186352e-05,0.000114372709309901,0.000148647128283624
+0,1.23164013985779e-07,2.01619957543725e-06,8.06746622712326e-06,1.9408951787608e-05,3.62807655737012e-05,5.85218420212134e-05,8.58109894680849e-05,0.000117776723517294,0.000154044674846032,0.00019425785338382
+0,7.38055239000408e-07,7.13661815273176e-06,2.1957201331878e-05,4.48051934222451e-05,7.46425998373469e-05,0.000110486321507166,0.000151514246571192,0.000197051618892773,0.000246540311467432,0.00029951243350813
+0,3.94931915614762e-06,2.37664477362685e-05,5.85078377488893e-05,0.000104356634872865,0.000158457684921854,0.000218839976363591,0.000284118381194624,0.000353286447661602,0.00042558888273829,0.000500441217237315
+0,1.87536920881212e-05,7.42297836614388e-05,0.000150953066786234,0.000240330072584709,0.000337698410142427,0.000440229993907094,0.00054607757879758,0.000653971107402446,0.000763006862111349,0.000872522857886958
+0,7.67060564087528e-05,0.000211583511011325,0.000364290214247424,0.000523711941940329,0.000684842914248833,0.000844968149576811,0.00100249413331233,0.00115645811127351,0.00130628569191617,0.00145165008061455
+0,1.29198147797745e-07,7.26337646539414e-07,1.8789774469963e-06,3.58894576036159e-06,5.83743184885366e-06,8.59992729697426e-06,1.1850973027563e-05,1.55658536023551e-05,1.97211620171109e-05,2.42949568717564e-05
+0,1.89355863862761e-09,2.0149347443802e-08,7.32356468455249e-08,1.74272463156743e-07,3.31196104192811e-07,5.481872642976e-07,8.2692218882057e-07,1.16749109225162e-06,1.56902754100022e-06,2.03013177037545e-06
+0,4.63537374431423e-11,9.81049734780514e-10,5.2662674195959e-09,1.63744082880102e-08,3.80213424268806e-08,7.37013316252193e-08,1.26462911129735e-07,1.98832947513683e-07,2.92825213625396e-07,4.09988785684554e-07
+0,9.52762309519183e-13,4.03064693534399e-11,3.21954008431541e-10,1.31921085909251e-09,3.77652149951806e-09,8.65283546407388e-09,1.70473483728495e-08,3.01291907541136e-08,4.90801983841489e-08,7.50526597099419e-08
+0,1.69833113269482e-14,1.45989004470562e-12,1.76611013922202e-11,9.70786436099117e-11,3.4867494868619e-10,9.60502789076322e-10,2.20870660779088e-09,4.45777595611339e-09,8.15402734448277e-09,1.38148719622932e-08
+0,3.8915196963406e-11,2.71147134577767e-09,2.44901018232562e-08,1.01911652264688e-07,2.85170970090303e-07,6.29937137214719e-07,1.19197488155735e-06,2.02429502148315e-06,3.17588061786517e-06,4.69121110232959e-06
+0,2.69185075542407e-10,1.00640042548647e-08,6.5998653123948e-08,2.25483644613014e-07,5.52286794898488e-07,1.10885707763952e-06,1.95314930874714e-06,3.13753072469243e-06,4.70863852449393e-06,6.70761014036735e-06
+0,1.89155326006871e-09,4.15998208034506e-08,2.08604035348676e-07,6.01200476492204e-07,1.30548633543253e-06,2.39315754080736e-06,3.9233691392589e-06,5.9449170992954e-06,8.49806448732054e-06,1.16159780561656e-05
+0,1.20839005460873e-08,1.63063271190061e-07,6.45781047622635e-07,1.61034741128995e-06,3.16276730648069e-06,5.37778778885908e-06,8.30838167561027e-06,1.19919722977472e-05,1.64544997795969e-05,2.17132051523398e-05
+0,6.98490407129936e-08,6.0448889346693e-07,1.94055108474347e-06,4.2597649069398e-06,7.66233240335793e-06,1.22016422946016e-05,1.79016972891742e-05,2.476703497099e-05,3.27887852389906e-05,4.19486886713683e-05
+0,3.50437914406735e-07,2.02949943667367e-06,5.38075813382185e-06,1.04873870799375e-05,1.73426483579829e-05,2.58972599259919e-05,3.60796128969124e-05,4.78063686693464e-05,6.09885319751288e-05,7.55352438512354e-05
+0,2.52830836713648e-10,3.03681967057117e-09,1.23485455469166e-08,3.26088143094103e-08,6.82703234078585e-08,1.23666594902017e-07,2.02940810635938e-07,3.10010651321603e-07,4.48552679684649e-07,6.21997431145188e-07
+0,5.62803310698392e-13,1.27640446069762e-11,7.32343235940447e-11,2.42833323537223e-10,6.0008786533641e-10,1.23564347077632e-09,2.24827721329559e-09,3.74216806129839e-09,5.82509155188699e-09,8.60722506603327e-09
+0,1.01387223716053e-14,4.50063177569099e-13,3.76367908055613e-12,1.6114550692278e-11,4.81228068266938e-11,1.1484654272007e-10,2.3535141406549e-10,4.32106599374873e-10,7.30344360289206e-10,1.15746018360108e-09
+0,1.76594647095188e-16,1.5591711552107e-14,1.93446159551681e-13,1.08935554130944e-12,4.00524857658289e-12,1.12878863819327e-11,2.65435340079881e-11,5.47630083633411e-11,1.02366893517743e-10,1.77190911910303e-10
+0,2.72938929463367e-18,4.86658533444932e-16,9.10403238648642e-15,6.85236595138709e-14,3.15126909211129e-13,1.06497145870696e-12,2.91648358432319e-12,6.85836462119616e-12,1.4373187057492e-11,2.75299294530364e-11
+0,6.21517803140922e-14,9.59223157644613e-12,1.40071205391153e-10,8.23460517316941e-10,3.01588393968386e-09,8.29925719783238e-09,1.88974373890384e-08,3.76435597362723e-08,6.79236252804073e-08,1.13611867015149e-07
+0,4.89626154105836e-13,4.05815676062707e-11,4.27904064325323e-10,2.04932427478431e-09,6.51652855169281e-09,1.61738279532313e-08,3.40432642093828e-08,6.37469204528233e-08,1.09424110270049e-07,1.75651555512571e-07
+0,3.92394931909311e-12,1.90205223022698e-10,1.52475376523412e-09,6.12909774497683e-09,1.7202341627476e-08,3.88279099539024e-08,7.57896022082591e-08,1.33417356513312e-07,2.17461256827561e-07,3.33987690067168e-07
+0,2.8933367514824e-11,8.49634390328452e-10,5.32093344063423e-09,1.83496200775791e-08,4.62765938242325e-08,9.63901823877104e-08,1.76587279811171e-07,2.95134919493546e-07,4.60499444520641e-07,6.81219273509691e-07
+0,1.95394555097804e-10,3.60599662518081e-09,1.80378231091857e-08,5.41614780599373e-08,1.24057163023099e-07,2.4042528230335e-07,4.1611636054975e-07,6.638317374418e-07,9.95928543382885e-07,1.42428849426674e-06
+0,1.14722473747572e-09,1.37465796551839e-08,5.57323434828026e-08,1.46670400626854e-07,3.05907787720263e-07,5.51856057412464e-07,9.01663873101398e-07,1.37107045554111e-06,1.97435903337969e-06,2.72436530869665e-06
+0,6.01331068747146e-13,1.48666863208301e-11,9.22909151658988e-11,3.28914376451997e-10,8.6846657139908e-10,1.90065440926915e-09,3.65824807559149e-09,6.41362850871594e-09,1.04751195906183e-08,1.61832994397799e-08
+0,2.47209629421286e-16,1.20861421473939e-14,1.10842575020953e-13,5.18453630267307e-13,1.68538118789818e-12,4.36405036350932e-12,9.67329876973588e-12,1.91546745620659e-11,3.48213414272979e-11,5.92006053141907e-11
+0,1.89201707214526e-18,1.74461832628427e-16,2.25617452264583e-15,1.32210769025309e-14,5.05114773691188e-14,1.47735746305447e-13,3.60121598845708e-13,7.69378961925591e-13,1.48783448922978e-12,2.66185680931199e-12
+0,2.82938853517197e-20,5.166160471115e-18,9.88324205334377e-17,7.5994903425817e-16,3.56758223375618e-15,1.23000563514818e-14,3.43481298564478e-14,8.23311838199757e-14,1.75814766445246e-13,3.43039404178284e-13
+0,3.8648935971775e-22,1.41879457944149e-19,4.08084894427787e-18,4.18430606959238e-17,2.4523734952935e-16,1.01226704556591e-15,3.28784556445063e-15,8.97275887121034e-15,2.14634799098414e-14,4.63096283150802e-14
+0,8.69726640954885e-17,2.92087981800005e-14,6.81319014247032e-13,5.60977017733479e-12,2.67241290013728e-11,9.12044259353017e-11,2.49116450278525e-10,5.80809267726238e-10,1.20365994453422e-09,2.27799582449115e-09
+0,7.66381761027786e-16,1.38259280480535e-13,2.31866103567873e-12,1.54608384966719e-11,6.35500210965569e-11,1.94448402919345e-10,4.88264062451473e-10,1.06479889696641e-09,2.08990883233503e-09,3.78080528306937e-09
+0,6.87383174408679e-15,7.21096117919272e-13,9.14517774138443e-12,5.09570157318565e-11,1.84185808276166e-10,5.10862690038548e-10,1.1861963221612e-09,2.42559699182021e-09,4.50988933263573e-09,7.78908854017501e-09
+0,5.72300830893883e-14,3.59756333500574e-12,3.53239285399593e-11,1.6766501693894e-10,5.4153652527201e-10,1.38002614943133e-09,2.99692285950551e-09,5.80167590670959e-09,1.03017611203301e-08,1.71032265873292e-08
+0,4.4015162593732e-13,1.71041885299831e-11,1.32555658504776e-10,5.42962213883109e-10,1.58183574169237e-09,3.73024467366084e-09,7.61934524058875e-09,1.40282403312057e-08,2.38787854849993e-08,3.82284985005463e-08
+0,2.94015789946406e-12,7.24491356536544e-11,4.48309024211327e-10,1.59267035200818e-09,4.19217068372523e-09,9.14631653977701e-09,1.75503380261564e-08,3.06757992440402e-08,4.99505561914961e-08,7.69389888117878e-08
+0,1.25520662515885e-15,6.34202350482921e-14,5.98701785935928e-13,2.87279812512801e-12,9.55263541481702e-12,2.52380922965864e-11,5.69558750757137e-11,1.14607040894637e-10,2.11363987866348e-10,3.6401786096932e-10
+0,2.69293069884163e-19,2.76242819230984e-17,3.94910273521441e-16,2.5438021672209e-15,1.0629681164291e-14,3.38496697614216e-14,8.94675031954185e-14,2.06465258826618e-13,4.29763614969243e-13,8.24924625699432e-13
+0,3.40393634256781e-22,6.53632617070714e-20,1.31310268914671e-18,1.05893998255463e-17,5.20829887103865e-17,1.87956724398352e-16,5.48951609054929e-16,1.3751099157628e-15,3.06669713938502e-15,6.2447989892164e-15
+0,4.00278935798164e-24,1.50166583042231e-21,4.40866620105182e-20,4.60970828253607e-19,2.75310436445269e-18,1.15735040672427e-17,3.82669184362989e-17,1.06268988633135e-16,2.58592708871164e-16,5.67426085674954e-16
+0,4.89860423986417e-26,3.68055944030794e-23,1.6212071379999e-21,2.25639400694407e-20,1.6800782168471e-19,8.44489572108405e-19,3.24448865148392e-18,1.02486142861079e-17,2.79143092738667e-17,6.76923776023426e-17

src/BTCS.cpp

@@ -1,327 +1,418 @@
-#include <array>
+/**
-#include <iostream>
+ * @file BTCSv2.cpp
-#include <tug/BoundaryCondition.hpp>
+ * @brief Implementation of heterogenous BTCS (backward time-centered space)
-#include <tug/Diffusion.hpp>
+ * solution of diffusion equation in 1D and 2D space. Internally the
-#include <tug/Solver.hpp>
+ * alternating-direction implicit (ADI) method is used. Version 2, because
-
+ * Version 1 was an implementation for the homogeneous BTCS solution.
-#include <Eigen/Dense>
+ *
-#include <Eigen/Sparse>
+ */
-#include <Eigen/src/Core/Matrix.h>
-#include <chrono>
-#include <vector>
 
+#include "Schemes.hpp"
 #include "TugUtils.hpp"
 
+#include <Eigen/src/Core/util/Meta.h>
+#include <cstddef>
+#include <tug/Boundary.hpp>
+#include <tug/Grid.hpp>
+#include <utility>
+
 #ifdef _OPENMP
 #include <omp.h>
 #else
 #define omp_get_thread_num() 0
 #endif
 
-inline auto
+// calculates coefficient for boundary in constant case
-init_delta(const std::array<double, tug::diffusion::MAX_ARR_SIZE> &domain_size,
+template <class T>
-           const std::array<uint32_t, tug::diffusion::MAX_ARR_SIZE> &grid_cells,
+constexpr std::pair<T, T> calcBoundaryCoeffConstant(T alpha_center,
-           const uint8_t dim) -> std::vector<double> {
+                                                    T alpha_side, T sx) {
-  std::vector<double> out(dim);
+  const T centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) +
-  for (uint8_t i = 0; i < dim; i++) {
+                                  2 * alpha_center);
-    out[i] = (double)(domain_size.at(i) / grid_cells.at(i));
+  const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
-  }
+
-  return out;
+  return {centerCoeff, sideCoeff};
 }
 
-namespace {
+// calculates coefficient for boundary in closed case
-enum { GRID_1D = 1, GRID_2D, GRID_3D };
+template <class T>
+constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
+                                                  T sx) {
+  const T centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side);
 
-constexpr int BTCS_MAX_DEP_PER_CELL = 3;
+  const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
-constexpr int BTCS_2D_DT_SIZE = 2;
 
-using DMatrixRowMajor =
+  return {centerCoeff, sideCoeff};
-    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-using DVectorRowMajor =
-    Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>;
-
-inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c,
-                          double neighbor_alpha) -> double {
-  double val = 0;
-
-  if (bc.type == tug::bc::BC_TYPE_CLOSED) {
-    val = neighbor_c;
-  } else if (bc.type == tug::bc::BC_TYPE_FLUX) {
-    // TODO
-    //  val = bc[index].value;
-  } else {
-    // TODO: implement error handling here. Type was set to wrong value.
-  }
-
-  return val;
 }
 
-auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
+// creates coefficient matrix for next time step from alphas in x-direction
-                  const tug::bc::BoundaryCondition &bc, bool transposed,
+static Eigen::SparseMatrix<double>
-                  double time_step, double dx) -> DMatrixRowMajor {
+createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
+                  std::vector<BoundaryElement> &bcRight, int numCols,
+                  int rowIndex, double sx) {
 
-  uint8_t upper = (transposed ? tug::bc::BC_SIDE_LEFT : tug::bc::BC_SIDE_TOP);
+  // square matrix of column^2 dimension for the coefficients
-  uint8_t lower =
+  Eigen::SparseMatrix<double> cm(numCols, numCols);
-      (transposed ? tug::bc::BC_SIDE_RIGHT : tug::bc::BC_SIDE_BOTTOM);
+  cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
 
-  int n_rows = c.rows();
+  // left column
-  int n_cols = c.cols();
 
-  DMatrixRowMajor d_ortho(n_rows, n_cols);
+  switch (bcLeft[rowIndex].getType()) {
-
+  case BC_TYPE_CONSTANT: {
-  std::array<double, 3> y_values{};
+    auto [centerCoeffTop, rightCoeffTop] =
-
+        calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
-  // first, iterate over first row
+    cm.insert(0, 0) = centerCoeffTop;
-  for (int j = 0; j < n_cols; j++) {
+    cm.insert(0, 1) = rightCoeffTop;
-    tug::bc::boundary_condition tmp_bc = bc(upper, j);
+    break;
-    double sy = (time_step * alpha(0, j)) / (dx * dx);
+  }
-
+  case BC_TYPE_CLOSED: {
-    y_values[0] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT
+    auto [centerCoeffTop, rightCoeffTop] =
-                       ? tmp_bc.value
+        calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
-                       : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)));
+    cm.insert(0, 0) = centerCoeffTop;
-    y_values[1] = c(0, j);
+    cm.insert(0, 1) = rightCoeffTop;
-    y_values[2] = c(1, j);
+    break;
-
+  }
-    d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]);
+  default: {
+    throw_invalid_argument(
+        "Undefined Boundary Condition Type somewhere on Left or Top!");
+  }
   }
 
-  // then iterate over inlet
+  // inner columns
-  for (int i = 1; i < n_rows - 1; i++) {
+  int n = numCols - 1;
-    for (int j = 0; j < n_cols; j++) {
+  for (int i = 1; i < n; i++) {
-      double sy = (time_step * alpha(i, j)) / (dx * dx);
+    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));
+  }
 
-      y_values[0] = c(i - 1, j);
+  // right column
-      y_values[1] = c(i, j);
-      y_values[2] = c(i + 1, j);
 
-      d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]);
+  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 <typename T>
+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 <typename T>
+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<BoundaryElement> &bcLeft,
+    std::vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcTop,
+    std::vector<BoundaryElement> &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);
     }
   }
 
-  int end = n_rows - 1;
+  // first row
-
+  else if (rowIndex == 0) {
-  // and finally over last row
+    for (int i = 0; i < length; i++) {
-  for (int j = 0; j < n_cols; j++) {
+      switch (bcTop[i].getType()) {
-    tug::bc::boundary_condition tmp_bc = bc(lower, j);
+      case BC_TYPE_CONSTANT: {
-    double sy = (time_step * alpha(end, j)) / (dx * dx);
+        sv(i) = calcExplicitConcentrationsBoundaryConstant(
-
+            concentrations(rowIndex, i), bcTop[i].getValue(),
-    y_values[0] = c(end - 1, j);
+            alphaY(rowIndex, i), alphaY(rowIndex + 1, i), sy);
-    y_values[1] = c(end, j);
+        break;
-    y_values[2] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT
-                       ? tmp_bc.value
-                       : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)));
-
-    d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]);
-  }
-
-  return d_ortho;
-}
-
-auto fillMatrixFromRow(const DVectorRowMajor &alpha,
-                       const tug::bc::bc_vec &bc_inner, int size, double dx,
-                       double time_step) -> Eigen::SparseMatrix<double> {
-
-  Eigen::SparseMatrix<double> A_matrix(size + 2, size + 2);
-
-  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL));
-
-  double sx = 0;
-
-  int A_size = A_matrix.cols();
-
-  A_matrix.insert(0, 0) = 1;
-
-  if (bc_inner[0].type != tug::bc::BC_UNSET) {
-    if (bc_inner[0].type != tug::bc::BC_TYPE_CONSTANT) {
-      throw_invalid_argument("Inner boundary conditions with other type than "
-                             "BC_TYPE_CONSTANT are currently not supported.");
-    }
-    A_matrix.insert(1, 1) = 1;
-  } else {
-    sx = (alpha[0] * time_step) / (dx * dx);
-    A_matrix.insert(1, 1) = -1. - 3. * sx;
-    A_matrix.insert(1, 0) = 2. * sx;
-    A_matrix.insert(1, 2) = sx;
-  }
-
-  for (int j = 2, k = j - 1; k < size - 1; j++, k++) {
-    if (bc_inner[k].type != tug::bc::BC_UNSET) {
-      if (bc_inner[k].type != tug::bc::BC_TYPE_CONSTANT) {
-        throw_invalid_argument("Inner boundary conditions with other type than "
-                               "BC_TYPE_CONSTANT are currently not supported.");
       }
-      A_matrix.insert(j, j) = 1;
+      case BC_TYPE_CLOSED: {
-      continue;
+        sv(i) = calcExplicitConcentrationsBoundaryClosed(
-    }
+            concentrations(rowIndex, i), alphaY(rowIndex, i),
-    sx = (alpha[k] * time_step) / (dx * dx);
+            alphaY(rowIndex + 1, i), sy);
-
+        break;
-    A_matrix.insert(j, j) = -1. - 2. * sx;
+      }
-    A_matrix.insert(j, (j - 1)) = sx;
+      default:
-    A_matrix.insert(j, (j + 1)) = sx;
+        throw_invalid_argument(
-  }
+            "Undefined Boundary Condition Type somewhere on Left or Top!");
-
-  if (bc_inner[size - 1].type != tug::bc::BC_UNSET) {
-    if (bc_inner[size - 1].type != tug::bc::BC_TYPE_CONSTANT) {
-      throw_invalid_argument("Inner boundary conditions with other type than "
-                             "BC_TYPE_CONSTANT are currently not supported.");
-    }
-    A_matrix.insert(A_size - 2, A_size - 2) = 1;
-  } else {
-    sx = (alpha[size - 1] * time_step) / (dx * dx);
-    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
-    A_matrix.insert(A_size - 2, A_size - 3) = sx;
-    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
-  }
-
-  A_matrix.insert(A_size - 1, A_size - 1) = 1;
-
-  return A_matrix;
-}
-
-auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha,
-                       const tug::bc::bc_tuple &bc,
-                       const tug::bc::bc_vec &bc_inner,
-                       const DVectorRowMajor &d_ortho, int size, double dx,
-                       double time_step) -> Eigen::VectorXd {
-
-  Eigen::VectorXd b_vector(size + 2);
-
-  tug::bc::boundary_condition left = bc[0];
-  tug::bc::boundary_condition right = bc[1];
-
-  bool left_constant = (left.type == tug::bc::BC_TYPE_CONSTANT);
-  bool right_constant = (right.type == tug::bc::BC_TYPE_CONSTANT);
-
-  int b_size = b_vector.size();
-
-  for (int j = 0; j < size; j++) {
-    if (bc_inner[j].type != tug::bc::BC_UNSET) {
-      if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT) {
-        throw_invalid_argument("Inner boundary conditions with other type than "
-                               "BC_TYPE_CONSTANT are currently not supported.");
       }
-      b_vector[j + 1] = bc_inner[j].value;
-      continue;
     }
-    b_vector[j + 1] = -c[j] + d_ortho[j];
   }
 
-  // this is not correct currently.We will fix this when we are able to define
+  // last row
-  // FLUX boundary conditions
+  else if (rowIndex == numRows - 1) {
-  b_vector[0] =
+    for (int i = 0; i < length; i++) {
-      (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
+      switch (bcBottom[i].getType()) {
-
+      case BC_TYPE_CONSTANT: {
-  b_vector[b_size - 1] =
+        sv(i) = calcExplicitConcentrationsBoundaryConstant(
-      (right_constant ? right.value
+            concentrations(rowIndex, i), bcBottom[i].getValue(),
-                      : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
+            alphaY(rowIndex, i), alphaY(rowIndex - 1, i), sy);
-
+        break;
-  return b_vector;
+      }
-}
+      case BC_TYPE_CLOSED: {
-
+        sv(i) = calcExplicitConcentrationsBoundaryClosed(
-auto setupBTCSAndSolve(
+            concentrations(rowIndex, i), alphaY(rowIndex, i),
-    DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost,
+            alphaY(rowIndex - 1, i), sy);
-    const tug::bc::bc_vec &bc_inner, const DVectorRowMajor &alpha, double dx,
+        break;
-    double time_step, int size, const DVectorRowMajor &d_ortho,
+      }
-    Eigen::VectorXd (*solver)(const Eigen::SparseMatrix<double> &,
+      default:
-                              const Eigen::VectorXd &)) -> DVectorRowMajor {
+        throw_invalid_argument(
-
+            "Undefined Boundary Condition Type somewhere on Right or Bottom!");
-  const Eigen::SparseMatrix<double> A_matrix =
+      }
-      fillMatrixFromRow(alpha, bc_inner, size, dx, time_step);
+    }
-
-  const Eigen::VectorXd b_vector = fillVectorFromRow(
-      c, alpha, bc_ghost, bc_inner, d_ortho, size, dx, time_step);
-
-  // solving of the LEQ
-  Eigen::VectorXd x_vector = solver(A_matrix, b_vector);
-
-  DVectorRowMajor out_vector = x_vector.segment(1, size);
-
-  return out_vector;
-}
-
-} // namespace
-  //
-auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param,
-                             double *field, const double *alpha) -> double {
-
-  auto start = time_marker();
-
-  uint32_t size = input_param.grid.grid_cells[0];
-
-  auto deltas = init_delta(input_param.grid.domain_size,
-                           input_param.grid.grid_cells, GRID_1D);
-  double dx = deltas[0];
-
-  double time_step = input_param.time_step;
-
-  const tug::bc::BoundaryCondition bc =
-      (input_param.grid.bc != nullptr ? *input_param.grid.bc
-                                      : tug::bc::BoundaryCondition(size));
-
-  Eigen::Map<DVectorRowMajor> c_in(field, size);
-  Eigen::Map<const DVectorRowMajor> alpha_in(alpha, size);
-
-  DVectorRowMajor input_field = c_in.row(0);
-
-  DVectorRowMajor output = setupBTCSAndSolve(
-      input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha_in, dx,
-      time_step, size, Eigen::VectorXd::Constant(size, 0), input_param.solver);
-
-  c_in.row(0) << output;
-
-  auto end = time_marker();
-
-  return diff_time(start, end);
-}
-
-auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param,
-                            double *field, const double *alpha) -> double {
-
-  auto start = time_marker();
-
-  uint32_t n_cols = input_param.grid.grid_cells[0];
-  uint32_t n_rows = input_param.grid.grid_cells[1];
-
-  auto deltas = init_delta(input_param.grid.domain_size,
-                           input_param.grid.grid_cells, GRID_2D);
-  double dx = deltas[0];
-  double dy = deltas[1];
-
-  double local_dt = input_param.time_step / BTCS_2D_DT_SIZE;
-
-  tug::bc::BoundaryCondition bc =
-      (input_param.grid.bc != nullptr
-           ? *input_param.grid.bc
-           : tug::bc::BoundaryCondition(n_cols, n_rows));
-
-  Eigen::Map<DMatrixRowMajor> c_in(field, n_rows, n_cols);
-  Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, n_rows, n_cols);
-
-  DMatrixRowMajor d_ortho =
-      calc_d_ortho(c_in, alpha_in, bc, false, local_dt, dx);
-
-#pragma omp parallel for schedule(dynamic)
-  for (int i = 0; i < n_rows; i++) {
-    DVectorRowMajor input_field = c_in.row(i);
-    DVectorRowMajor output = setupBTCSAndSolve(
-        input_field, bc.row_boundary(i), bc.getInnerRow(i), alpha_in.row(i), dx,
-        local_dt, n_cols, d_ortho.row(i), input_param.solver);
-    c_in.row(i) << output;
   }
 
-  d_ortho = calc_d_ortho(c_in.transpose(), alpha_in.transpose(), bc, true,
+  // first column -> additional fixed concentration change from perpendicular
-                         local_dt, dy);
+  // dimension in constant bc case
-
+  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
-#pragma omp parallel for schedule(dynamic)
+    sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
-  for (int i = 0; i < n_cols; i++) {
-    DVectorRowMajor input_field = c_in.col(i);
-    DVectorRowMajor output = setupBTCSAndSolve(
-        input_field, bc.col_boundary(i), bc.getInnerCol(i), alpha_in.col(i), dy,
-        local_dt, n_rows, d_ortho.row(i), input_param.solver);
-    c_in.col(i) << output.transpose();
   }
 
-  auto end = time_marker();
+  // 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 diff_time(start, end);
+  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<double> &A,
+                                        Eigen::VectorXd &b) {
+
+  Eigen::SparseLU<Eigen::SparseMatrix<double>> 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<double> &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<double> &A,
+                                      Eigen::VectorXd &b)) {
+  int length = grid.getLength();
+  double sx = timestep / (grid.getDelta() * grid.getDelta());
+
+  Eigen::VectorXd concentrations_t1(length);
+
+  Eigen::SparseMatrix<double> A;
+  Eigen::VectorXd b(length);
+
+  Eigen::MatrixXd alpha = grid.getAlpha();
+  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  std::vector<BoundaryElement> 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<double> &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<double> A;
+  Eigen::VectorXd b;
+
+  Eigen::MatrixXd alphaX = grid.getAlphaX();
+  Eigen::MatrixXd alphaY = grid.getAlphaY();
+  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+  std::vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
+  std::vector<BoundaryElement> 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<Eigen::SparseMatrix<double>> 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!");
+  }
 }

src/Boundary.cpp

@@ -0,0 +1,149 @@
+#include "TugUtils.hpp"
+
+#include <cstddef>
+#include <stdexcept>
+#include <tug/Boundary.hpp>
+
+BoundaryElement::BoundaryElement() {}
+
+BoundaryElement::BoundaryElement(double _value)
+    : value(_value), type(BC_TYPE_CONSTANT) {}
+
+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;
+}
+
+BC_TYPE BoundaryElement::getType() { return this->type; }
+
+double BoundaryElement::getValue() { return this->value; }
+
+Boundary::Boundary(Grid grid) : grid(grid) {
+  if (grid.getDim() == 1) {
+    this->boundaries = std::vector<std::vector<BoundaryElement>>(
+        2); // in 1D only left and right boundary
+
+    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<std::vector<BoundaryElement>>(4);
+
+    this->boundaries[BC_SIDE_LEFT] =
+        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
+    this->boundaries[BC_SIDE_RIGHT] =
+        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
+    this->boundaries[BC_SIDE_TOP] =
+        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
+    this->boundaries[BC_SIDE_BOTTOM] =
+        std::vector<BoundaryElement>(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.");
+    }
+  }
+
+  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<BoundaryElement>(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.");
+    }
+  }
+
+  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<BoundaryElement>(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);
+}
+
+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);
+}
+
+const std::vector<BoundaryElement> 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];
+}
+
+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);
+  }
+
+  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];
+}
+
+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();
+}
+
+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();
+}

src/BoundaryCondition.cpp

@@ -1,208 +0,0 @@
-#include <algorithm>
-#include <tug/BoundaryCondition.hpp>
-#include <vector>
-
-#include "TugUtils.hpp"
-
-constexpr uint8_t DIM_1D = 2;
-constexpr uint8_t DIM_2D = 4;
-
-tug::bc::BoundaryCondition::BoundaryCondition(int x) {
-  this->bc_internal.resize(DIM_1D, {0, 0});
-  this->dim = 1;
-  // this value is actually unused
-  this->maxsize = 1;
-
-  this->sizes[X_DIM] = x;
-  this->sizes[Y_DIM] = 1;
-
-  this->maxindex = x - 1;
-}
-
-tug::bc::BoundaryCondition::BoundaryCondition(int x, int y) {
-  this->maxsize = (x >= y ? x : y);
-  this->bc_internal.resize(DIM_2D * maxsize, {0, 0});
-  this->dim = 2;
-
-  this->sizes[X_DIM] = x;
-  this->sizes[Y_DIM] = y;
-
-  this->maxindex = (x * y) - 1;
-}
-
-void tug::bc::BoundaryCondition::setSide(
-    uint8_t side, tug::bc::boundary_condition &input_bc) {
-  if (this->dim == 1) {
-    throw_invalid_argument("setSide requires at least a 2D grid");
-  }
-  if (side > 3) {
-    throw_out_of_range("Invalid range for 2D grid");
-  }
-
-  uint32_t size =
-      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-           ? this->sizes[Y_DIM]
-           : this->sizes[X_DIM]);
-
-  for (uint32_t i = 0; i < size; i++) {
-    this->bc_internal[side * maxsize + i] = input_bc;
-  }
-}
-
-void tug::bc::BoundaryCondition::setSide(
-    uint8_t side, std::vector<tug::bc::boundary_condition> &input_bc) {
-  if (this->dim == 1) {
-    throw_invalid_argument("setSide requires at least a 2D grid");
-  }
-  if (side > 3) {
-    throw_out_of_range("Invalid range for 2D grid");
-  }
-
-  uint32_t size =
-      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-           ? this->sizes[Y_DIM]
-           : this->sizes[X_DIM]);
-
-  if (input_bc.size() > size) {
-    throw_out_of_range("Input vector is greater than maximum excpected value");
-  }
-
-  for (int i = 0; i < size; i++) {
-    bc_internal[this->maxsize * side + i] = input_bc[i];
-  }
-}
-
-auto tug::bc::BoundaryCondition::getSide(uint8_t side)
-    -> std::vector<tug::bc::boundary_condition> {
-  if (this->dim == 1) {
-    throw_invalid_argument("getSide requires at least a 2D grid");
-  }
-  if (side > 3) {
-    throw_out_of_range("Invalid range for 2D grid");
-  }
-
-  uint32_t size =
-      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-           ? this->sizes[Y_DIM]
-           : this->sizes[X_DIM]);
-
-  std::vector<tug::bc::boundary_condition> out(size);
-
-  for (int i = 0; i < size; i++) {
-    out[i] = this->bc_internal[this->maxsize * side + i];
-  }
-
-  return out;
-}
-
-auto tug::bc::BoundaryCondition::row_boundary(uint32_t i) const
-    -> tug::bc::bc_tuple {
-  if (i >= this->sizes[Y_DIM]) {
-    throw_out_of_range("Index out of range");
-  }
-
-  return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i],
-          this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]};
-}
-
-auto tug::bc::BoundaryCondition::col_boundary(uint32_t i) const
-    -> tug::bc::bc_tuple {
-  if (this->dim == 1) {
-    throw_invalid_argument("Access of column requires at least 2D grid");
-  }
-  if (i >= this->sizes[X_DIM]) {
-    throw_out_of_range("Index out of range");
-  }
-
-  return {this->bc_internal[BC_SIDE_TOP * this->maxsize + i],
-          this->bc_internal[BC_SIDE_BOTTOM * this->maxsize + i]};
-}
-
-auto tug::bc::BoundaryCondition::getInnerRow(uint32_t i) const -> bc_vec {
-  if (i >= this->sizes[Y_DIM]) {
-    throw_out_of_range("Index is out of range");
-  }
-
-  bc_vec row(this->sizes[X_DIM], {tug::bc::BC_UNSET, 0});
-
-  if (this->inner_cells.empty()) {
-    return row;
-  }
-
-  uint32_t index_min = i * this->sizes[X_DIM];
-  uint32_t index_max = ((i + 1) * this->sizes[X_DIM]) - 1;
-
-  for (auto const &cell : this->inner_cells) {
-    if (cell.first < index_min) {
-      continue;
-    }
-    if (cell.first > index_max) {
-      break;
-    }
-
-    row[cell.first - index_min] = cell.second;
-  }
-
-  return row;
-}
-
-auto tug::bc::BoundaryCondition::getInnerCol(uint32_t i) const -> bc_vec {
-  if (this->dim != 2) {
-    throw_invalid_argument("getInnerCol is only applicable for 2D grids");
-  }
-  if (i >= this->sizes[X_DIM]) {
-    throw_out_of_range("Index is out of range");
-  }
-
-  bc_vec col(this->sizes[Y_DIM], {tug::bc::BC_UNSET, 0});
-
-  if (this->inner_cells.empty()) {
-    return col;
-  }
-
-  for (auto const &cell : this->inner_cells) {
-    if (cell.first % this->sizes[X_DIM] == i) {
-      col[cell.first / this->sizes[X_DIM]] = cell.second;
-    }
-  }
-
-  return col;
-}
-
-void tug::bc::BoundaryCondition::setInnerBC(boundary_condition bc, int x,
-                                            int y = 0) {
-  if (x >= this->sizes[X_DIM] || y >= this->sizes[Y_DIM]) {
-    throw_out_of_range("One input parameter is out of range");
-  }
-  uint32_t index = y * this->sizes[X_DIM] + x;
-  auto it = this->inner_cells.find(index);
-
-  if (it != this->inner_cells.end()) {
-    it->second = bc;
-    return;
-  }
-
-  this->inner_cells.insert({index, bc});
-}
-
-void tug::bc::BoundaryCondition::unsetInnerBC(int x, int y) {
-  uint32_t index = y * this->sizes[X_DIM] + x;
-  this->inner_cells.erase(index);
-}
-
-auto tug::bc::BoundaryCondition::getInnerBC(int x, int y = 0)
-    -> boundary_condition {
-  if (x >= this->sizes[X_DIM] || y >= this->sizes[Y_DIM]) {
-    throw_out_of_range("One input parameter is out of range");
-  }
-
-  uint32_t index = y * this->sizes[X_DIM] + x;
-
-  auto it = this->inner_cells.find(index);
-
-  if (it != this->inner_cells.end()) {
-    return it->second;
-  }
-
-  return {tug::bc::BC_UNSET, 0};
-}

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp)
+add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp)
 
 target_link_libraries(tug Eigen3::Eigen)
 
@@ -6,4 +6,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)

src/FTCS.cpp

@@ -0,0 +1,381 @@
+/**
+ * @file FTCS.cpp
+ * @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 <cstddef>
+#include <iostream>
+#include <tug/Boundary.hpp>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_thread_num() 0
+#endif
+
+// calculates horizontal change on one cell independent of boundary type
+static inline double calcHorizontalChange(Grid &grid, int &row, int &col) {
+
+  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 inline double calcVerticalChange(Grid &grid, int &row, int &col) {
+
+  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 inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid,
+                                                              Boundary &bc,
+                                                              int &row,
+                                                              int &col) {
+
+  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 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 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 inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
+                                                               Boundary &bc,
+                                                               int &row,
+                                                               int &col) {
+
+  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 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 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 inline double calcVerticalChangeTopBoundaryConstant(Grid &grid,
+                                                           Boundary &bc,
+                                                           int &row, int &col) {
+
+  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 inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
+                                                         int &col) {
+
+  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 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 inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid,
+                                                              Boundary &bc,
+                                                              int &row,
+                                                              int &col) {
+
+  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 inline double
+calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
+
+  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 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 &timestep) {
+  int colMax = grid.getCol();
+  double deltaCol = grid.getDeltaCol();
+
+  // 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;
+
+  // 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));
+
+  // 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);
+}
+
+// FTCS solution for 2D grid
+static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
+                    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
+  Eigen::MatrixXd concentrations_t1 =
+      Eigen::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 < colMax - 1; col++) {
+    concentrations_t1(row, col) =
+        grid.getConcentrations()(row, col) +
+        timestep / (deltaRow * deltaRow) *
+            (calcVerticalChangeTopBoundary(grid, bc, row, col)) +
+        timestep / (deltaCol * deltaCol) *
+            (calcHorizontalChange(grid, row, col));
+  }
+
+  // bottom without corners / looping over columns
+  // hold row constant at max index
+  row = rowMax - 1;
+#pragma omp parallel for num_threads(numThreads)
+  for (int col = 1; col < colMax - 1; col++) {
+    concentrations_t1(row, col) =
+        grid.getConcentrations()(row, col) +
+        timestep / (deltaRow * deltaRow) *
+            (calcVerticalChangeBottomBoundary(grid, bc, row, col)) +
+        timestep / (deltaCol * deltaCol) *
+            (calcHorizontalChange(grid, row, col));
+  }
+
+  // corner top left
+  // hold row and column constant at 0
+  row = 0;
+  col = 0;
+  concentrations_t1(row, col) =
+      grid.getConcentrations()(row, col) +
+      timestep / (deltaCol * deltaCol) *
+          (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
+      timestep / (deltaRow * deltaRow) *
+          (calcVerticalChangeTopBoundary(grid, bc, row, col));
+
+  // corner top right
+  // hold row constant at 0 and column constant at max index
+  row = 0;
+  col = colMax - 1;
+  concentrations_t1(row, col) =
+      grid.getConcentrations()(row, col) +
+      timestep / (deltaCol * deltaCol) *
+          (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
+      timestep / (deltaRow * deltaRow) *
+          (calcVerticalChangeTopBoundary(grid, bc, row, col));
+
+  // corner bottom left
+  // hold row constant at max index and column constant at 0
+  row = rowMax - 1;
+  col = 0;
+  concentrations_t1(row, col) =
+      grid.getConcentrations()(row, col) +
+      timestep / (deltaCol * deltaCol) *
+          (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
+      timestep / (deltaRow * deltaRow) *
+          (calcVerticalChangeBottomBoundary(grid, bc, row, col));
+
+  // corner bottom right
+  // hold row and column constant at max index
+  row = rowMax - 1;
+  col = colMax - 1;
+  concentrations_t1(row, col) =
+      grid.getConcentrations()(row, col) +
+      timestep / (deltaCol * deltaCol) *
+          (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
+      timestep / (deltaRow * deltaRow) *
+          (calcVerticalChangeBottomBoundary(grid, bc, row, col));
+
+  // overwrite obsolete concentrations
+  grid.setConcentrations(concentrations_t1);
+  // }
+}
+
+// entry point; differentiate between 1D and 2D grid
+void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
+  if (grid.getDim() == 1) {
+    FTCS_1D(grid, bc, timestep);
+  } else if (grid.getDim() == 2) {
+    FTCS_2D(grid, bc, timestep, numThreads);
+  } else {
+    throw_invalid_argument(
+        "Error: Only 1- and 2-dimensional grids are defined!");
+  }
+}

src/Grid.cpp

@@ -0,0 +1,167 @@
+#include "TugUtils.hpp"
+
+#include <iostream>
+#include <tug/Grid.hpp>
+
+constexpr double MAT_INIT_VAL = 0;
+
+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->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)
+    : 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->dim = 2;
+  this->deltaRow = double(this->domainRow) / double(this->row); // -> 1
+  this->deltaCol = double(this->domainCol) / double(this->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(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;
+}
+
+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(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 = alphaX;
+  this->alphaY = alphaY;
+}
+
+const Eigen::MatrixXd &Grid::getAlpha() {
+  if (dim != 1) {
+    throw_invalid_argument("Grid is not one dimensional, you should probably "
+                           "use either getAlphaX() or getAlphaY()!");
+  }
+
+  return this->alphaX;
+}
+
+const Eigen::MatrixXd &Grid::getAlphaX() {
+  if (dim != 2) {
+    throw_invalid_argument(
+        "Grid is not two dimensional, you should probably use getAlpha()!");
+  }
+
+  return this->alphaX;
+}
+
+const Eigen::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::getLength() {
+  if (dim != 1) {
+    throw_invalid_argument("Grid is not one dimensional, you should probably "
+                           "use getRow() or getCol()!");
+  }
+
+  return col;
+}
+
+int Grid::getRow() { return row; }
+
+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!");
+  }
+
+  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!");
+  }
+
+  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");
+  }
+
+  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!");
+  }
+
+  return this->deltaRow;
+}

src/README.md

@@ -0,0 +1,36 @@
+<h1>src-Directory</h1>
+This is the src-directory that holds the source code to the implementation of the TUG framework.
+
+<pre>
+src/  
+  ├── CMakeFiles/
+  ├── Boundary.cpp  
+  ├── BTCSv2.cpp
+  ├── CMakeLists.txt  
+  ├── FTCS.cpp
+  ├── Grid.cpp 
+  ├── README.md
+  ├── Simulation.cpp
+  └── TugUtils.hpp 
+</pre>
+
+**src/** Directory with the source code.
+
+**CMakeFiles/** Automatically generated directory by CMake.
+
+**Boundary.cpp** Implementation of Boundary class, that holds the boundary conditions.
+
+**BTCSv2.cpp** Implementation of BTCS solution to heterogeneous diffusion in 1D and 2D. 
+
+**CMakeLists.txt** CMakeLists for this directory.
+
+**FTCS.cpp** Implementation of FTCS solution to heterogeneous diffusion in 1D and 2D. 
+
+**Grid.cpp** Implementation of Grid class, that holds all of the concentrations alpha coefficient in x- and y-direction.
+
+**README.md** <i>This</i> file.
+
+**Simulation.cpp** Implementation of Simulation class, that holds all of the information for a specific simulation run, as well as the Boundary and Grid objects. 
+
+**TugUtils.hpp** Helper functions for other source files. 
+

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 <tug/Boundary.hpp>
+#include <tug/Grid.hpp>
+
+// entry point; differentiate between 1D and 2D grid
+extern void FTCS(Grid &grid, Boundary &bc, double &timestep, 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_

src/Simulation.cpp

@@ -0,0 +1,326 @@
+#include <cmath>
+#include <cstddef>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <string>
+
+#include <tug/Simulation.hpp>
+
+#include "Schemes.hpp"
+#include "TugUtils.hpp"
+
+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!");
+  }
+
+  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!");
+  }
+
+  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!");
+  }
+
+  this->time_measure = time_measure;
+}
+
+void Simulation::setTimestep(double timestep) {
+  if (timestep <= 0) {
+    throw_invalid_argument("Timestep has to be greater than zero.");
+  }
+
+  if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) {
+
+    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;
+
+    double maxAlpha = std::numeric_limits<double>::quiet_NaN();
+
+    // determine maximum alpha
+    if (grid.getDim() == 2) {
+
+      const double maxAlphaX = grid.getAlphaX().maxCoeff();
+      const double maxAlphaY = grid.getAlphaY().maxCoeff();
+      maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
+
+    } else if (grid.getDim() == 1) {
+      maxAlpha = grid.getAlpha().maxCoeff();
+
+    } else {
+      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; }
+
+void Simulation::setIterations(int 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) {
+    std::cerr
+        << "Warning: Solver was set, but FTCS approach initialized. Setting "
+           "the solver "
+           "is thus without effect."
+        << std::endl;
+  }
+
+  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();
+    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; }
+
+void Simulation::printConcentrationsConsole() {
+  std::cout << grid.getConcentrations() << std::endl;
+  std::cout << std::endl;
+}
+
+std::string Simulation::createCSVfile() {
+  std::ofstream file;
+  int appendIdent = 0;
+  std::string appendIdentString;
+
+  // 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);
+
+  std::string filename =
+      approachString + "_" + row + "_" + col + "_" + numIterations + ".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);
+  }
+
+  // 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();
+
+  return filename;
+}
+
+void Simulation::printConcentrationsCSV(const std::string &filename) {
+  std::ofstream file;
+
+  file.open(filename, std::ios_base::app);
+  if (!file) {
+    exit(1);
+  }
+
+  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!");
+  }
+
+  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 (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
+
+    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<std::chrono::milliseconds>(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;
+  }
+}

src/Solver.cpp

@@ -1,61 +0,0 @@
-#include <tug/Solver.hpp>
-
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-
-auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> &A_matrix,
-                          const Eigen::VectorXd &b_vector) -> Eigen::VectorXd {
-
-  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
-      solver;
-  solver.analyzePattern(A_matrix);
-
-  solver.factorize(A_matrix);
-
-  Eigen::VectorXd x_vec = solver.solve(b_vector);
-
-  return solver.solve(b_vector);
-}
-
-auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix<double> &A_matrix,
-                                  const Eigen::VectorXd &b_vector)
-    -> Eigen::VectorXd {
-  uint32_t n = b_vector.size();
-
-  Eigen::VectorXd a_diag(n);
-  Eigen::VectorXd b_diag(n);
-  Eigen::VectorXd c_diag(n);
-  Eigen::VectorXd x_vec = b_vector;
-
-  // Fill diagonals vectors
-  b_diag[0] = A_matrix.coeff(0, 0);
-  c_diag[0] = A_matrix.coeff(0, 1);
-
-  for (int i = 1; i < n - 1; i++) {
-    a_diag[i] = A_matrix.coeff(i, i - 1);
-    b_diag[i] = A_matrix.coeff(i, i);
-    c_diag[i] = A_matrix.coeff(i, i + 1);
-  }
-  a_diag[n - 1] = A_matrix.coeff(n - 1, n - 2);
-  b_diag[n - 1] = A_matrix.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 (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;
-}

src/TugUtils.hpp

@@ -1,5 +1,5 @@
-#ifndef BTCSUTILS_H_
+#ifndef TUGUTILS_H_
-#define BTCSUTILS_H_
+#define TUGUTILS_H_
 
 #include <chrono>
 #include <stdexcept>
@@ -23,4 +23,15 @@
                                                                   start);      \
     duration.count();                                                          \
   })
-#endif // BTCSUTILS_H_
+
+// calculates arithmetic or harmonic mean of alpha between two cells
+template <typename T>
+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_

test/CMakeLists.txt

@@ -11,9 +11,16 @@ if(NOT DOCTEST_LIB)
   FetchContent_MakeAvailable(DocTest)
 endif()
 
-add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.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)
+# 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}" "${PROJECT_SOURCE_DIR}/src")
+
 add_custom_target(
   check
   COMMAND $<TARGET_FILE:testTug>

test/FTCS_11_11_7000.csv

@@ -0,0 +1,13 @@
+0 0 0 0 0 0 0 0 0 0 0
+1.88664e-08 3.39962e-08 7.57021e-08 1.76412e-07 4.15752e-07 9.00973e-07 3.65403e-09 9.6579e-12 3.59442e-13 3.42591e-14 3.27595e-15
+1.19102e-06 1.95195e-06 3.92165e-06 8.30575e-06 1.78976e-05 3.60742e-05 2.02843e-07 2.24659e-09 2.35085e-10 2.64988e-11 2.90933e-12
+5.85009e-05 8.57948e-05 0.000151499 0.000284105 0.00054607 0.00100251 1.18494e-05 8.26706e-07 1.26394e-07 1.70309e-08 2.20525e-09
+0.00202345 0.00258511 0.00381783 0.00599829 0.00972689 0.0154873 0.0011152 0.000247309 5.10506e-05 8.64727e-06 1.37747e-06
+0.0205848 0.0217651 0.0238282 0.0262762 0.0285812 0.0303808 0.0374255 0.0204234 0.00674813 0.00160264 0.000338852
+0.0199112 0.0210265 0.0229587 0.0252019 0.0272007 0.0285225 0.0310896 0.0156681 0.00495399 0.00114176 0.000235573
+0.0184589 0.0194561 0.0211596 0.0230704 0.0246215 0.0253266 0.0228374 0.010038 0.00292986 0.000638799 0.000125942
+0.0166888 0.0175517 0.0190015 0.0205611 0.0216793 0.0218694 0.0160278 0.00593307 0.00155164 0.000312583 5.76964e-05
+0.0151262 0.0158758 0.0171155 0.0183949 0.019193 0.0190523 0.0115557 0.00356455 0.000817461 0.000148892 2.51893e-05
+0.0142177 0.0149034 0.0160255 0.0171522 0.0177843 0.0174891 0.0093846 0.0025221 0.000515469 8.51039e-05 1.31328e-05
+
+

test/TestUtils.hpp

@@ -0,0 +1,49 @@
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <fstream>
+#include <ios>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+
+inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
+
+  std::vector<double> 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<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
+      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;
+}

test/testBoundary.cpp

@@ -0,0 +1,76 @@
+#include <doctest/doctest.h>
+#include <iostream>
+#include <stdio.h>
+#include <string>
+#include <tug/Boundary.hpp>
+#include <typeinfo>
+
+using namespace std;
+
+TEST_CASE("BoundaryElement") {
+
+  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<BoundaryElement> 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 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());
+  }
+}

test/testBoundaryCondition.cpp

@@ -1,196 +0,0 @@
-#include <doctest/doctest.h>
-#include <tug/BoundaryCondition.hpp>
-
-using namespace tug::bc;
-
-#define BC_CONST_VALUE 1e-5
-
-TEST_CASE("1D Boundary Condition") {
-
-  BoundaryCondition bc(5);
-  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
-
-  SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT).value, 0); }
-
-  SUBCASE("invalid get") {
-    CHECK_THROWS(bc(BC_SIDE_TOP));
-    CHECK_THROWS(bc(BC_SIDE_LEFT, 1));
-  }
-
-  SUBCASE("valid set") {
-    CHECK_NOTHROW(bc(BC_SIDE_LEFT) = bc_set);
-    CHECK_EQ(bc(BC_SIDE_LEFT).value, bc_set.value);
-    CHECK_EQ(bc(BC_SIDE_LEFT).type, bc_set.type);
-  }
-
-  SUBCASE("invalid set") {
-    CHECK_THROWS(bc(BC_SIDE_TOP) = bc_set);
-    CHECK_THROWS(bc(BC_SIDE_LEFT, 0) = bc_set);
-  }
-
-  SUBCASE("valid row getter") {
-    bc(BC_SIDE_LEFT) = bc_set;
-    bc_tuple tup = bc.row_boundary(0);
-
-    CHECK_EQ(tup[0].value, BC_CONST_VALUE);
-    CHECK_EQ(tup[1].value, 0);
-  }
-
-  SUBCASE("invalid row and col getter") {
-    CHECK_THROWS(bc.row_boundary(1));
-    CHECK_THROWS(bc.col_boundary(0));
-  }
-}
-
-TEST_CASE("2D Boundary Condition") {
-
-  BoundaryCondition bc(5, 5);
-  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
-
-  SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, 0); }
-
-  SUBCASE("invalid get") {
-    CHECK_THROWS(bc(5, 0));
-    CHECK_THROWS(bc(BC_SIDE_LEFT));
-  }
-
-  SUBCASE("valid set") {
-    CHECK_NOTHROW(bc(BC_SIDE_LEFT, 0) = bc_set);
-    CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, bc_set.value);
-    CHECK_EQ(bc(BC_SIDE_LEFT, 0).type, bc_set.type);
-  }
-
-  SUBCASE("invalid set") {
-    CHECK_THROWS(bc(BC_SIDE_LEFT) = bc_set);
-    CHECK_THROWS(bc(5, 0) = bc_set);
-  }
-
-  SUBCASE("call of setSide") {
-    CHECK_NOTHROW(bc.setSide(BC_SIDE_BOTTOM, bc_set));
-    CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).value, bc_set.value);
-    CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).type, bc_set.type);
-  }
-
-  SUBCASE("get and set of side") {
-    std::vector<boundary_condition> bc_vec;
-    CHECK_NOTHROW(bc_vec = bc.getSide(BC_SIDE_BOTTOM));
-    bc_vec[3] = {BC_TYPE_CONSTANT, 1e-5};
-    CHECK_NOTHROW(bc.setSide(BC_SIDE_BOTTOM, bc_vec));
-    CHECK_EQ(bc(BC_SIDE_BOTTOM, 3).type, BC_TYPE_CONSTANT);
-    CHECK_EQ(bc(BC_SIDE_BOTTOM, 3).value, 1e-5);
-
-    CHECK_EQ(bc(BC_SIDE_BOTTOM, 2).value, 0);
-  }
-}
-
-TEST_CASE("Boundary Condition helpers") {
-  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
-
-  SUBCASE("return boundary condition skeleton") {
-    boundary_condition bc_test =
-        BoundaryCondition::returnBoundaryCondition(bc_set.type, bc_set.value);
-    CHECK_EQ(bc_test.value, bc_set.value);
-    CHECK_EQ(bc_test.type, bc_set.type);
-  }
-}
-
-TEST_CASE("1D special inner grid cells") {
-  BoundaryCondition bc(5);
-  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
-
-  SUBCASE("valid set") { CHECK_NOTHROW(bc.setInnerBC(bc_set, 0, 0)); }
-
-  SUBCASE("valid get") {
-    bc.setInnerBC(bc_set, 0, 0);
-    CHECK_EQ(bc.getInnerBC(0, 0).type, bc_set.type);
-  }
-
-  SUBCASE("invalid get") {
-    CHECK_EQ(bc.getInnerBC(1, 0).type, BC_UNSET);
-    CHECK_THROWS(bc.getInnerBC(0, 1));
-  }
-
-  SUBCASE("invalid set") { CHECK_THROWS(bc.setInnerBC(bc_set, 0, 1)); }
-
-  SUBCASE("valid row getter") {
-    bc.setInnerBC(bc_set, 1, 0);
-
-    bc_vec ret = bc.getInnerRow(0);
-
-    CHECK_EQ(ret[0].type, BC_UNSET);
-    CHECK_EQ(ret[1].type, bc_set.type);
-  }
-
-  SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(1)); }
-
-  SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(0)); }
-}
-
-TEST_CASE("2D special inner grid cells") {
-  BoundaryCondition bc(5, 5);
-  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
-
-  SUBCASE("valid set") { CHECK_NOTHROW(bc.setInnerBC(bc_set, 0, 0)); }
-
-  SUBCASE("valid get") {
-    bc.setInnerBC(bc_set, 0, 0);
-    CHECK_EQ(bc.getInnerBC(0, 0).type, bc_set.type);
-  }
-
-  SUBCASE("invalid get") {
-    CHECK_EQ(bc.getInnerBC(1, 0).type, BC_UNSET);
-    CHECK_THROWS(bc.getInnerBC(5, 5));
-  }
-
-  SUBCASE("invalid set") { CHECK_THROWS(bc.setInnerBC(bc_set, 5, 5)); }
-
-  SUBCASE("valid row getter") {
-    bc.setInnerBC(bc_set, 0, 0);
-    bc.setInnerBC(bc_set, 1, 1);
-
-    bc_vec ret = bc.getInnerRow(0);
-
-    CHECK_EQ(ret[0].type, bc_set.type);
-    CHECK_EQ(ret[1].type, BC_UNSET);
-
-    ret = bc.getInnerRow(1);
-
-    CHECK_EQ(ret[0].type, BC_UNSET);
-    CHECK_EQ(ret[1].type, bc_set.type);
-  }
-
-  SUBCASE("valid col getter") {
-    bc.setInnerBC(bc_set, 0, 1);
-    bc.setInnerBC(bc_set, 1, 0);
-
-    bc_vec ret = bc.getInnerCol(0);
-
-    CHECK_EQ(ret[0].type, BC_UNSET);
-    CHECK_EQ(ret[1].type, bc_set.type);
-
-    ret = bc.getInnerCol(1);
-
-    CHECK_EQ(ret[0].type, bc_set.type);
-    CHECK_EQ(ret[1].type, BC_UNSET);
-  }
-
-  SUBCASE("unset boundary condition") {
-    bc.setInnerBC(bc_set, 0, 0);
-    bc.setInnerBC(bc_set, 1, 1);
-    bc.unsetInnerBC(1, 1);
-
-    bc_vec ret = bc.getInnerRow(1);
-
-    CHECK_EQ(ret[0].type, BC_UNSET);
-    CHECK_EQ(ret[1].type, BC_UNSET);
-
-    ret = bc.getInnerCol(1);
-
-    CHECK_EQ(ret[0].type, BC_UNSET);
-    CHECK_EQ(ret[1].type, BC_UNSET);
-  }
-
-  SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(5)); }
-
-  SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(5)); }
-}

test/testDiffusion.cpp

@@ -1,170 +0,0 @@
-#include <doctest/doctest.h>
-#include <tug/BoundaryCondition.hpp>
-#include <tug/Diffusion.hpp>
-#include <tug/Solver.hpp>
-#include <vector>
-
-using namespace tug::bc;
-using namespace tug::diffusion;
-
-#define DIMENSION 2
-#define N 51
-#define M 51
-#define MID 1300
-
-static std::vector<double> alpha(N *M, 1e-3);
-
-static TugInput setupDiffu() {
-  TugInput diffu;
-
-  diffu.setTimestep(1);
-  diffu.setGridCellN(N, M);
-  diffu.setDomainSize(N, M);
-
-  return diffu;
-}
-
-TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") {
-  std::vector<double> field(N * M, 0);
-
-  field[MID] = 1;
-
-  BoundaryCondition bc(N, M);
-
-  TugInput diffu = setupDiffu();
-
-  uint32_t iterations = 1000;
-  double sum = 0;
-
-  for (int t = 0; t < iterations; t++) {
-    ADI_2D(diffu, field.data(), alpha.data());
-  }
-
-  // iterate through rows
-  for (int i = 0; i < M; i++) {
-    // iterate through columns
-    for (int j = 0; j < N; j++) {
-      sum += field[i * N + j];
-    }
-  }
-  CAPTURE(sum);
-  // epsilon of 1e-8
-  CHECK(sum == doctest::Approx(1).epsilon(1e-6));
-}
-
-TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") {
-  std::vector<double> field(N * M, 0);
-
-  field[MID] = 1;
-
-  BoundaryCondition bc(N, M);
-
-  boundary_condition input = {BC_TYPE_CONSTANT, 0};
-
-  bc.setSide(BC_SIDE_LEFT, input);
-  bc.setSide(BC_SIDE_RIGHT, input);
-  bc.setSide(BC_SIDE_TOP, input);
-  bc.setSide(BC_SIDE_BOTTOM, input);
-
-  TugInput diffu = setupDiffu();
-  diffu.setBoundaryCondition(bc);
-
-  uint32_t max_iterations = 20000;
-  bool reached = false;
-
-  int t = 0;
-
-  for (t = 0; t < max_iterations; t++) {
-    ADI_2D(diffu, field.data(), alpha.data());
-
-    if (field[N * M - 1] > 1e-15) {
-      reached = true;
-      break;
-    }
-  }
-
-  if (!reached) {
-    CAPTURE(field[N * M - 1]);
-    FAIL_CHECK(
-        "Concentration did not reach boundaries after count of iterations: ",
-        t);
-  }
-}
-
-TEST_CASE(
-    "constant top and bottom (1 and 0) - left and right closed - 0 inlet") {
-  std::vector<double> field(N * M, 0);
-
-  BoundaryCondition bc(N, M);
-
-  boundary_condition top =
-      BoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 1);
-  boundary_condition bottom =
-      BoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 0);
-
-  bc.setSide(BC_SIDE_TOP, top);
-  bc.setSide(BC_SIDE_BOTTOM, bottom);
-
-  TugInput diffu = setupDiffu();
-  diffu.setBoundaryCondition(bc);
-
-  uint32_t max_iterations = 100;
-
-  for (int t = 0; t < max_iterations; t++) {
-    ADI_2D(diffu, field.data(), alpha.data());
-  }
-
-  for (int i = 0; i < N; i++) {
-    double above = field[i];
-    for (int j = 1; j < M; j++) {
-      double curr = field[j * N + i];
-      if (curr > above) {
-        CAPTURE(curr);
-        CAPTURE(above);
-        FAIL("Concentration below is greater than above @ cell ", j * N + i);
-      }
-    }
-  }
-}
-
-TEST_CASE("2D closed boundaries, 1 constant cell in the middle") {
-  std::vector<double> field(N * M, 0);
-  double val = 1e-2;
-
-  BoundaryCondition bc(N, M);
-
-  field[MID] = val;
-  bc.setInnerBC({BC_TYPE_CONSTANT, val}, N / 2, M / 2);
-
-  TugInput diffu = setupDiffu();
-  diffu.setBoundaryCondition(bc);
-
-  uint32_t max_iterations = 100;
-
-  double sum = val;
-
-  for (int t = 0; t < max_iterations; t++) {
-    ADI_2D(diffu, field.data(), alpha.data());
-
-    CHECK_EQ(field[MID], val);
-
-    double new_sum = .0;
-
-    for (int i = 0; i < M; i++) {
-      // iterate through columns
-      for (int j = 0; j < N; j++) {
-        new_sum += field[i * N + j];
-      }
-    }
-
-    if (sum > new_sum) {
-      CAPTURE(sum);
-      CAPTURE(new_sum);
-      FAIL("Sum of concentrations is equal or lower than to previous iteration "
-           "after iteration ",
-           t + 1);
-    }
-
-    sum = new_sum;
-  }
-}

test/testFTCS.cpp

@@ -0,0 +1,20 @@
+#include <TugUtils.hpp>
+
+#include <doctest/doctest.h>
+#include <limits>
+
+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));
+
+    // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) -
+    // harmonicMean); CHECK(difference <
+    // std::numeric_limits<double>::epsilon());
+    CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean);
+    CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average);
+  }
+}

test/testGrid.cpp

@@ -0,0 +1,253 @@
+#include <Eigen/Core>
+#include <doctest/doctest.h>
+#include <tug/Grid.hpp>
+
+using namespace Eigen;
+using namespace std;
+
+
+TEST_CASE("1D Grid, too small length") {
+  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));
+
+  r = 4;
+  c = 2;
+  CHECK_THROWS(Grid(r, c));
+}
+
+TEST_CASE("1D Grid") {
+  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);
+
+    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());
+  }
+
+  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));
+  }
+
+  SUBCASE("setting alpha") {
+    // correct alpha matrix
+    MatrixXd alpha = MatrixXd::Constant(1, l, 3);
+    CHECK_NOTHROW(grid.setAlpha(alpha));
+
+    CHECK_THROWS(grid.setAlpha(alpha, alpha));
+
+    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));
+  }
+
+  SUBCASE("setting domain") {
+    int d = 8;
+    // set 1D domain
+    CHECK_NOTHROW(grid.setDomain(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());
+
+    // 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);
+
+  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.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));
+
+    // 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));
+
+    CHECK_THROWS(grid.setAlpha(alphax));
+
+    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));
+  }
+
+  SUBCASE("setting domain") {
+    int dr = 8;
+    int dc = 9;
+
+    // set 1D domain
+    CHECK_THROWS(grid.setDomain(dr));
+
+    // 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));
+
+    // 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);
+
+  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.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));
+
+    // 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));
+
+    CHECK_THROWS(grid.setAlpha(alphax));
+
+    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));
+  }
+
+  SUBCASE("setting domain") {
+    int dr = 8;
+    int dc = 9;
+
+    // set 1D domain
+    CHECK_THROWS(grid.setDomain(dr));
+
+    // 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));
+
+    // 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));
+  }
+}

test/testSimulation.cpp

@@ -0,0 +1,108 @@
+#include "TestUtils.hpp"
+
+#include <doctest/doctest.h>
+#include <stdio.h>
+#include <string>
+#include <tug/Simulation.hpp>
+
+// include the configured header file
+#include <testSimulation.hpp>
+
+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);
+
+  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;
+    }
+  }
+  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;
+    }
+  }
+  grid.setAlpha(alpha, alpha);
+
+  // 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();
+
+  // 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);
+  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);
+  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);
+
+  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);
+
+  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 timestep") {
+    CHECK_NOTHROW(sim.setTimestep(0.1));
+    CHECK_EQ(sim.getTimestep(), 0.1);
+    CHECK_THROWS(sim.setTimestep(-0.3));
+  }
+}

test/testSimulation.hpp.in

@@ -0,0 +1,7 @@
+#ifndef TESTSIMULATION_H_
+#define TESTSIMULATION_H_
+
+// CSV file needed for validation
+const char *testSimulationCSVDir = "@testSimulationCSV@";
+
+#endif // TESTSIMULATION_H_

utils/output-visualization.html

@@ -0,0 +1,543 @@
+<!DOCTYPE html>
+<html lang="en-US">
+<head>
+    <!-- <link rel="stylesheet" href="style.css"> -->
+    <style>
+        body {
+            background-color: white;
+        }
+
+        :root {
+            --max-width: 1440px;
+        }
+
+        .main-wrapper {
+            display: grid;
+            grid-template-areas: "sidebar main";
+            grid-template-columns: minmax(0,1fr) minmax(0,2.5fr);
+            grid-column-gap: 50px;
+
+            margin: 0 auto;
+            max-width: var(--max-width);
+        }
+
+        #state-info {
+            background-color: lightgray;
+            width: 120px;
+            text-align: center;
+            border: 2px solid black;
+            border-radius: 5px;
+            font-family: Arial, Helvetica, sans-serif;
+            font-size: 25px;
+            margin-bottom: 20px;
+        }
+    </style>
+</head>
+<body>
+<div id="root">
+<div class="main-wrapper">
+
+    <div class="sidebar-container">
+        <aside id="sidebar">
+        <div class="menu"> 
+            <form>
+                <label for="c_file">File (XTREME-format): </label>
+                <input type="file" id="c_file" name="c_file" accept="text/csv" />
+            </form>
+        </div>
+        <div class="info">
+            <p>
+                Use the keys <i>q</i> and <i>e</i> for increasing and decreasing the state by <b>1</b>.
+                The keys <i>a</i> and <i>d</i> for increasing and decresing the state by <b>10</b>.
+            </p>
+            <div id="state-info">
+                State: -
+            </div>
+            <div id="legend-info">
+
+            </div>
+        </div>
+        </aside>
+    </div>
+
+    <div class="content-container">
+        <main id="content" class="main-content">
+            
+        </main>
+    </div>
+
+</div>
+</div>
+</body>
+</html>
+
+<script type="module">
+
+    import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
+
+    // layout information; constants
+    const width = document.getElementById('content').offsetWidth;
+    const height = document.getElementById('content').offsetWidth;
+    const svgMargin = { top: 10, right: 10, bottom: 10, left: 10 };
+    const boundary_gap = 1;
+    const gridWidth = width - svgMargin.left - svgMargin.right - 2*boundary_gap;
+    const gridHeight = height - svgMargin.top - svgMargin.bottom - 2*boundary_gap;
+
+    // simulation information
+    var state = 0; // iteration state
+    var iteration_count = 0;
+    var number_rows = 0;
+    var number_cols = 0;
+    var max_concentration = 0;
+    var data = [];
+    var leftb, rightb, topb, bottomb;
+
+    // test data; to be deleted later on
+    var test_data = [{"row": 0, "col": 0, "value": 1},
+                     {"row": 0, "col": 1, "value": 1},
+                     {"row": 0, "col": 2, "value": 1},
+                     {"row": 0, "col": 3, "value": 1},
+                     {"row": 1, "col": 0, "value": 1},
+                     {"row": 1, "col": 1, "value": 1},
+                     {"row": 1, "col": 2, "value": 1},
+                     {"row": 1, "col": 3, "value": 1},
+                     {"row": 2, "col": 0, "value": 1},
+                     {"row": 2, "col": 1, "value": 1},
+                     {"row": 2, "col": 2, "value": 1},
+                     {"row": 2, "col": 3, "value": 1},
+                     {"row": 3, "col": 0, "value": 1},
+                     {"row": 3, "col": 1, "value": 1},
+                     {"row": 3, "col": 2, "value": 1},
+                     {"row": 3, "col": 3, "value": 1}];
+
+
+    //////////////////////////////////////////
+    ////////
+    //////////////////////////////////////////
+    // Copyright 2021, Observable Inc.
+    // Released under the ISC license.
+    // https://observablehq.com/@d3/color-legend
+    function Legend(color, {
+    title,
+    tickSize = 6,
+    width = 320, 
+    height = 44 + tickSize,
+    marginTop = 18,
+    marginRight = 0,
+    marginBottom = 16 + tickSize,
+    marginLeft = 0,
+    ticks = width / 64,
+    tickFormat,
+    tickValues
+    } = {}) {
+
+    function ramp(color, n = 256) {
+        const canvas = document.createElement("canvas");
+        canvas.width = n;
+        canvas.height = 1;
+        const context = canvas.getContext("2d");
+        for (let i = 0; i < n; ++i) {
+        context.fillStyle = color(i / (n - 1));
+        context.fillRect(i, 0, 1, 1);
+        }
+        return canvas;
+    }
+
+    const svg = d3.create("svg")
+        .attr("width", width)
+        .attr("height", height)
+        .attr("viewBox", [0, 0, width, height])
+        .style("overflow", "visible")
+        .style("display", "block");
+
+    let tickAdjust = g => g.selectAll(".tick line").attr("y1", marginTop + marginBottom - height);
+    let x;
+
+    // Continuous
+    if (color.interpolate) {
+        const n = Math.min(color.domain().length, color.range().length);
+
+        x = color.copy().rangeRound(d3.quantize(d3.interpolate(marginLeft, width - marginRight), n));
+
+        svg.append("image")
+            .attr("x", marginLeft)
+            .attr("y", marginTop)
+            .attr("width", width - marginLeft - marginRight)
+            .attr("height", height - marginTop - marginBottom)
+            .attr("preserveAspectRatio", "none")
+            .attr("xlink:href", ramp(color.copy().domain(d3.quantize(d3.interpolate(0, 1), n))).toDataURL());
+    }
+
+    // Sequential
+    else if (color.interpolator) {
+        x = Object.assign(color.copy()
+            .interpolator(d3.interpolateRound(marginLeft, width - marginRight)),
+            {range() { return [marginLeft, width - marginRight]; }});
+
+        svg.append("image")
+            .attr("x", marginLeft)
+            .attr("y", marginTop)
+            .attr("width", width - marginLeft - marginRight)
+            .attr("height", height - marginTop - marginBottom)
+            .attr("preserveAspectRatio", "none")
+            .attr("xlink:href", ramp(color.interpolator()).toDataURL());
+
+        // scaleSequentialQuantile doesn’t implement ticks or tickFormat.
+        if (!x.ticks) {
+        if (tickValues === undefined) {
+            const n = Math.round(ticks + 1);
+            tickValues = d3.range(n).map(i => d3.quantile(color.domain(), i / (n - 1)));
+        }
+        if (typeof tickFormat !== "function") {
+            tickFormat = d3.format(tickFormat === undefined ? ",f" : tickFormat);
+        }
+        }
+    }
+
+    // Threshold
+    else if (color.invertExtent) {
+        const thresholds
+            = color.thresholds ? color.thresholds() // scaleQuantize
+            : color.quantiles ? color.quantiles() // scaleQuantile
+            : color.domain(); // scaleThreshold
+
+        const thresholdFormat
+            = tickFormat === undefined ? d => d
+            : typeof tickFormat === "string" ? d3.format(tickFormat)
+            : tickFormat;
+
+        x = d3.scaleLinear()
+            .domain([-1, color.range().length - 1])
+            .rangeRound([marginLeft, width - marginRight]);
+
+        svg.append("g")
+        .selectAll("rect")
+        .data(color.range())
+        .join("rect")
+            .attr("x", (d, i) => x(i - 1))
+            .attr("y", marginTop)
+            .attr("width", (d, i) => x(i) - x(i - 1))
+            .attr("height", height - marginTop - marginBottom)
+            .attr("fill", d => d);
+
+        tickValues = d3.range(thresholds.length);
+        tickFormat = i => thresholdFormat(thresholds[i], i);
+    }
+
+    // Ordinal
+    else {
+        x = d3.scaleBand()
+            .domain(color.domain())
+            .rangeRound([marginLeft, width - marginRight]);
+
+        svg.append("g")
+        .selectAll("rect")
+        .data(color.domain())
+        .join("rect")
+            .attr("x", x)
+            .attr("y", marginTop)
+            .attr("width", Math.max(0, x.bandwidth() - 1))
+            .attr("height", height - marginTop - marginBottom)
+            .attr("fill", color);
+
+        tickAdjust = () => {};
+    }
+
+    svg.append("g")
+        .attr("transform", `translate(0,${height - marginBottom})`)
+        .call(d3.axisBottom(x)
+            .ticks(ticks, typeof tickFormat === "string" ? tickFormat : undefined)
+            .tickFormat(typeof tickFormat === "function" ? tickFormat : undefined)
+            .tickSize(tickSize)
+            .tickValues(tickValues))
+        .call(tickAdjust)
+        .call(g => g.select(".domain").remove())
+        .call(g => g.append("text")
+            .attr("x", marginLeft)
+            .attr("y", marginTop + marginBottom - height - 6)
+            .attr("fill", "currentColor")
+            .attr("text-anchor", "start")
+            .attr("font-weight", "bold")
+            .attr("class", "title")
+            .text(title));
+
+    return svg.node();
+    }
+    //////////////////////////////////////////
+    ////////
+    //////////////////////////////////////////
+
+    
+
+    // max concentration
+    function calc_max_concentration(state) {
+        // var maxRow = data[state].map(function(row){ return Math.max.apply(Math, row); });
+        // var maxC = Math.max.apply(null, maxRow);
+        var maxC = Math.max(...data[state].map(d => d.value))
+        return maxC;
+    }
+
+    // scales
+    var col, row, color;
+    function create_n_update_scales() {
+        col = d3.scaleBand()
+            .range([0, gridWidth])
+            .domain([...Array(number_cols).keys()])
+            .padding(0.05);
+        row = d3.scaleBand()
+            .range([0, gridHeight])
+            .domain([...Array(number_rows).keys()])
+            .padding(0.05);
+        color = d3.scaleLinear()
+            .range(["#d0d0ec", "#7f0000"])
+            .domain([0, max_concentration]);
+    }
+
+    // grid
+    function draw_n_update_grid(grid_data) {
+        document.getElementById('state-info').innerHTML = `State: ${state}`;
+        // Grid rects
+        grid.selectAll()
+            .data(grid_data, (d,i) => d.row+':'+d.col)
+            .enter()
+            .append("rect")
+                .attr("x", (d,i) => col(d.col))
+                .attr("y", (d,i) => row(d.row))
+                .attr("width", col.bandwidth())
+                .attr("height", row.bandwidth())
+                .attr("fill", (d,i) => color(d.value))
+            .on("mouseover", mouseover)
+            .on("mousemove", mousemove)
+            .on("mouseleave", mouseleave);
+    }
+
+    // boundaries
+    function draw_boundaries() {
+        left_side.selectAll("rect")
+            .data(leftb)
+            .join("rect")
+            .attr("x", (d,i) => 0)
+            .attr("y", (d,i) => row(i))
+            .attr("width", svgMargin.left)
+            .attr("height", row.bandwidth())
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        right_side.selectAll("rect")
+            .data(rightb)
+            .join("rect")
+            .attr("x", (d,i) => 0)
+            .attr("y", (d,i) => row(i))
+            .attr("width", svgMargin.right)
+            .attr("height", row.bandwidth())
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        top_side.selectAll("rect")
+            .data(topb)
+            .join("rect")
+            .attr("x", (d,i) => col(i))
+            .attr("y", (d,i) => 0)
+            .attr("width", col.bandwidth())
+            .attr("height", svgMargin.top)
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        bottom_side.selectAll("rect")
+            .data(bottomb)
+            .join("rect")
+            .attr("x", (d,i) => col(i))
+            .attr("y", (d,i) => 0)
+            .attr("width", col.bandwidth())
+            .attr("height", svgMargin.bottom)
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+    }
+
+    function draw_legend() {
+        d3.select("#legend-info").selectAll("svg").remove();
+
+        var legend = Legend(d3.scaleLinear([0, max_concentration], ["#d0d0ec", "#7f0000"]), {
+            title: "Concentration",
+            tickFormat: ".2f",
+            tickValues: [0, max_concentration/4, max_concentration/2, max_concentration/4*3, max_concentration]
+        });
+
+        d3.select("#legend-info")
+            .node()
+            .appendChild(legend);
+    }
+
+    // tooltip events
+    var mouseover = function(event, d) {
+        tooltip.style("opacity", 1);
+    }
+    var mousemove = function(event, d) {
+        const tooltip_text = d3.select('.tooltip-text');
+        const tooltip_box = d3.select('.tooltip-box');
+
+        tooltip_text.text(`Value: ${d.value}`);
+        var box_width = tooltip_text.node().getComputedTextLength();
+        tooltip_box.attr("width", box_width+10);
+
+        const [x, y] = d3.pointer(event);
+
+        tooltip
+          .attr('transform', `translate(${x+svgMargin.left+50}, ${y+svgMargin.top})`);
+    }
+    var mouseleave = function(event, d) {
+        tooltip.style("opacity", 0)
+    }
+
+
+    // key events for changing grid iteration state
+    addEventListener("keydown", (event) => {
+    if (event.isComposing || event.keyCode === 81) {
+        if (state > 0) {
+            state -= 1;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_n_update_grid(data[state]);
+        }
+    }
+    if (event.isComposing || event.keyCode === 69) {
+        if (state < iteration_count-1) {
+            state += 1;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_n_update_grid(data[state]);
+        }
+    }
+    if (event.isComposing || event.keyCode === 65) {
+        if (state > 9) {
+            state -= 10;
+        } else {
+            state = 0;
+        }
+        max_concentration = calc_max_concentration(state);
+        create_n_update_scales();
+        draw_legend();
+        draw_n_update_grid(data[state]);
+    }
+    if (event.isComposing || event.keyCode === 68) {
+        if (state < iteration_count-10) {
+            state += 10;
+        } else {
+            state = iteration_count-1;
+        }
+        max_concentration = calc_max_concentration(state);
+        create_n_update_scales();
+        draw_legend();
+        draw_n_update_grid(data[state]);
+    }
+    });
+
+
+    // read data from file when selecting new file
+    async function create_data_struct(csv_text) {
+        var row_data = d3.dsvFormat(" ").parseRows(csv_text);
+
+        leftb = row_data[0];
+        rightb = row_data[1];
+        topb = row_data[2];
+        bottomb = row_data[3];
+
+        number_rows = row_data[0].length;
+        number_cols = row_data[2].length;
+        iteration_count = (row_data.length - 6) / (number_rows + 2);
+        console.log(iteration_count);
+        var iteration = []; // loop temp variable
+        var entry = {}; // loop temp variable
+        for (let i = 0; i < iteration_count; i++) {
+            iteration = []; // reset temp variable
+            for (let j = 0; j < number_rows; j++) {
+                for (let k = 0; k < number_cols; k++) {
+                    entry = {}; // reset temp variable
+                    entry.row = j;
+                    entry.col = k;
+                    entry.value = row_data[6 + i*(number_rows+2) + j][k];
+                    iteration.push(entry);
+                }
+            }
+            data.push(iteration);
+        }
+    }
+    function reader(file, callback) {
+        const fr = new FileReader();
+        fr.onload = () => callback(null, fr.result);
+        fr.onerror = (err) => callback(err);
+        fr.readAsText(file)
+    }
+    document.getElementById("c_file").addEventListener("change", (evt) => {
+        if (!evt.target.files) { // no file, do nothing
+            return;
+        }
+        reader(evt.target.files[0], (err, res) => {
+            create_data_struct(res);
+            state = 0;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_boundaries();
+            draw_n_update_grid(data[0]);
+        });
+    });
+
+
+    const svg = d3.create("svg")
+        .attr("width", width)
+        .attr("height", height);
+
+    var left_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(0, ${svgMargin.top+boundary_gap})`);
+    var right_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${width-svgMargin.right}, ${svgMargin.top+boundary_gap})`);
+    var top_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${svgMargin.left+boundary_gap}, 0)`);
+    var bottom_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${svgMargin.left+boundary_gap}, ${height-svgMargin.bottom})`);
+
+    var grid = svg.append("g")
+        .attr("class", "grid")
+        .attr("transform", `translate(${svgMargin.left + boundary_gap}, ${svgMargin.top + boundary_gap})`);
+        
+    // tooltip 
+    var tooltip = svg.append("g")
+        .attr("class", "tooltip-area")
+        .style("opacity", 0);
+    var tooltip_box = tooltip.append("rect")
+        .attr("class", "tooltip-box")
+        .attr("height", 25)
+        .attr("transform", `translate(-5,-18)`)
+        .attr("fill", "white")
+        .attr("stroke", "black");
+    tooltip.append("text")
+        .attr("class", "tooltip-text")
+        .text("Value: ");
+    
+    content.append(svg.node());
+
+</script>