max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 09:28:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'main' into naaice

Browse Source
This commit is contained in:
Max Lübke 2023-09-15 08:15:19 +02:00
commit 46f9cef3a9
42 changed files with 3023 additions and 2827 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitlab-ci.yml
View File

@ -43,9 +43,9 @@ pages:
lint:
stage: static_analyze
before_script:
- apk add clang-extra-tools
- apk add clang-extra-tools openmp-dev
script:
- mkdir lint && cd lint
- cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*, modernize-*" -DTUG_ENABLE_TESTING=OFF ..
- cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*" -DTUG_ENABLE_TESTING=OFF ..
- make tug
when: manual

16
CMakeLists.txt
View File

@ -33,7 +33,15 @@ endif()
option(TUG_ENABLE_TESTING
"Run tests after succesfull compilation"
ON)
OFF)
option(TUG_HANNESPHILIPP_EXAMPLES
"Compile example applications"
OFF)
option(TUG_NAAICE_EXAMPLE
"Enables NAAICE examples with export of CSV files"
OFF)
add_subdirectory(src)
@ -41,9 +49,9 @@ if(TUG_ENABLE_TESTING)
add_subdirectory(test)
endif()
add_subdirectory(examples)
option(TUG_NAAICE_EXAMPLE "Enables NAAICE examples with export of CSV files" OFF)
if(TUG_HANNESPHILIPP_EXAMPLES)
add_subdirectory(examples)
endif()
if (TUG_NAAICE_EXAMPLE)
add_subdirectory(naaice)

339
LICENSE Normal file
View File

@ -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.

5
README.org
View File

@ -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?

1
cmake-build-debug/_deps/doctest-src

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

4
docs_sphinx/_build/html/_images/activity_diagram_run.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 15 KiB

4
docs_sphinx/_build/html/_images/class_diagram.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 23 KiB

224
docs_sphinx/_build/vis.html
View File

@ -1,224 +0,0 @@
<!DOCTYPE html>
<div id="container"></div>
<form>
<label for="c_file">Select an output file: </label>
<input type="file" id="c_file" name="c_file" accept="text/csv" />
</form>
<script type="module">
import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
// Declare the chart dimensions and margins.
const width = 600;
const height = 600;
const svgMargin = { top: 50, right: 50, bottom: 50, left: 50 };
var gridWidth = width - svgMargin.left - svgMargin.right;
var gridHeight = height - svgMargin.top - svgMargin.bottom;
const gap = 1;
var state = 0;
var numIterations = 0;
var cellXcount = 0;
var cellYcount = 0;
var concentrations = [];
// Create the outer SVG container.
const svg = d3.create("svg")
.attr("width", width)
.attr("height", height);
svg.append("rect")
.attr("x", 0)
.attr("y", 0)
.attr("width", width)
.attr("height", height)
.attr("fill-opacity", 0)
.attr("stroke-width", 4)
.attr("stroke", "black");
var grid = svg.append("g")
.attr("class", "grid")
.attr("transform", `translate(${svgMargin.left}, ${svgMargin.top})`);
// color scale
var colorScale = d3.scaleLinear()
.range(["#d0d0ec", "#7f0000"])
.domain([1, 2000]);
function getColor(c, max) {
if (c == -1) {
return "#000000";
}
colorScale.domain([0, max]);
return colorScale(c);
}
function calcMaxConcentrationInIteration(state) {
var maxRow = concentrations[state].map(function(row){ return Math.max.apply(Math, row); });
var maxC = Math.max.apply(null, maxRow);
return maxC;
}
// Load data from file
async function createVisualFromFile(csv_input) {
// console.log(csv_input);
// var data = await d3.text(csv_input).then(function(csv){
// return d3.dsvFormat(" ").parseRows(csv);
// })
var data = d3.dsvFormat(" ").parseRows(csv_input);
// console.log(data);
var leftBoundary = data[0];
var rightBoundary = data[1];
var topBoundary = data[2];
var bottomBoundary = data[3];
cellXcount = topBoundary.length;
cellYcount = leftBoundary.length;
concentrations = []; // reset concentrations
grid.selectAll("rect").remove();
svg.selectAll("rect").remove();
state = 0;
console.log(topBoundary);
var cellWidth = gridWidth / cellXcount;
var cellHeight = gridHeight / cellYcount;
var iteration = [];
numIterations = (data.length - 6) / (cellYcount + 2);
console.log(numIterations);
for (let i = 0; i < numIterations; i++) {
iteration = [];
for (let j = 0; j < cellYcount; j++) {
iteration.push(data[i * (cellYcount + 2) + 6 + j])
}
concentrations.push(iteration);
}
console.log(concentrations);
var maxC = calcMaxConcentrationInIteration(state);
// Create grid
for (let i = 0; i < cellYcount; i++) {
for (let j = 0; j < cellXcount; j++) {
grid.append("rect")
.attr("x", i * cellWidth + gap/2)
.attr("y", j * cellHeight + gap/2)
.attr("width", cellWidth - gap)
.attr("height", cellHeight - gap)
.attr("fill", getColor(concentrations[state][i][j], maxC));
}
}
// Create Boundaries
// left and right
for (let j = 0; j < cellYcount; j++) {
svg.append("rect")
.attr("x", svgMargin.left - 10)
.attr("y", svgMargin.top + j * cellHeight + gap/2)
.attr("width", 7)
.attr("height", cellHeight - gap)
.attr("fill", getColor(leftBoundary[j], maxC));
svg.append("rect")
.attr("x", width - svgMargin.right + 3)
.attr("y", svgMargin.top + j * cellHeight + gap/2)
.attr("width", 7)
.attr("height", cellHeight - gap)
.attr("fill", getColor(rightBoundary[j], maxC));
}
// top and bottom
for (let i = 0; i < cellXcount; i++) {
svg.append("rect")
.attr("x", svgMargin.left + i * cellWidth + gap/2)
.attr("y", svgMargin.top - 10)
.attr("width", cellWidth - gap)
.attr("height", 7)
.attr("fill", getColor(topBoundary[i], maxC));
svg.append("rect")
.attr("x", svgMargin.left + i * cellWidth + gap/2)
.attr("y", height - svgMargin.bottom + 3)
.attr("width", cellWidth - gap)
.attr("height", 7)
.attr("fill", getColor(bottomBoundary[i], maxC));
}
}
function updateGrid(new_state) {
var maxC = calcMaxConcentrationInIteration(new_state);
console.log(maxC);
grid.selectAll("rect")
.attr("fill", function (d,i) {
var row = Math.floor(i/20);
var col = i%20;
return getColor(concentrations[new_state][row][col], maxC);
})
}
// key events for changing grid iteration state
addEventListener("keydown", (event) => {
if (event.isComposing || event.keyCode === 81) {
if (state > 0) {
state -= 1;
updateGrid(state);
}
}
if (event.isComposing || event.keyCode === 69) {
if (state < numIterations-1) {
state += 1;
updateGrid(state);
}
}
if (event.isComposing || event.keyCode === 65) {
if (state > 9) {
state -= 10;
} else {
state = 0;
}
updateGrid(state);
}
if (event.isComposing || event.keyCode === 68) {
if (state < numIterations-10) {
state += 10;
} else {
state = numIterations-1;
}
updateGrid(state);
}
});
async function logFile() {
var file = document.getElementById("c_file").files[0];
const reader = new FileReader();
reader.readAsText(file);
console.log(reader.result);
}
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", logFile);
document.getElementById("c_file").addEventListener("change", (evt) => {
// No files, do nothing.
if (!evt.target.files) {
return;
}
reader(evt.target.files[0], (err, res) => {
console.log(res); // Base64 `data:image/...` String result.
createVisualFromFile(res);
});
});
// Append the SVG element.
container.append(svg.node());
</script>

2
docs_sphinx/conf.py
View File

@ -22,7 +22,7 @@ subprocess.call('doxygen Doxyfile.in', shell=True)
# -- Project information -----------------------------------------------------
project = 'TUG'
copyright = 'MIT'
copyright = 'GPL2'
author = 'Philipp Ungrund, Hannes Signer'

15
examples/BTCS_1D_proto_example.cpp
View File

@ -1,5 +1,8 @@
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
int main(int argc, char *argv[]) {
// **************
// **** GRID ****
@ -9,12 +12,11 @@ int main(int argc, char *argv[]) {
int cells = 20;
Grid grid = Grid(cells);
MatrixXd concentrations = MatrixXd::Constant(1,20,0);
concentrations(0,0) = 2000;
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 ****
// ******************
@ -24,13 +26,13 @@ int main(int argc, char *argv[]) {
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
Simulation simulation =
Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
@ -38,7 +40,8 @@ int main(int argc, char *argv[]) {
// set the number of iterations
simulation.setIterations(100);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****

23
examples/BTCS_2D_proto_example.cpp
View File

@ -1,5 +1,8 @@
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
int main(int argc, char *argv[]) {
// EASY_PROFILER_ENABLE;
// profiler::startListen();
@ -10,25 +13,23 @@ int main(int argc, char *argv[]) {
// create a grid with a 20 x 20 field
int row = 40;
int col = 50;
Grid grid = Grid(row,col);
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;
// 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 ****
// ******************
@ -44,7 +45,6 @@ int main(int argc, char *argv[]) {
// 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
@ -55,13 +55,13 @@ int main(int argc, char *argv[]) {
// 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
Simulation simulation =
Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
@ -69,7 +69,8 @@ int main(int argc, char *argv[]) {
// set the number of iterations
simulation.setIterations(300);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_XTREME);
// **** RUN SIMULATION ****

2
examples/CMakeLists.txt
View File

@ -2,6 +2,7 @@ add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp)
add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp)
add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
add_executable(profiling_openmp profiling_openmp.cpp)
@ -9,6 +10,7 @@ target_link_libraries(FTCS_1D_proto_example tug)
target_link_libraries(FTCS_2D_proto_example tug)
target_link_libraries(BTCS_1D_proto_example tug)
target_link_libraries(BTCS_2D_proto_example tug)
target_link_libraries(CRNI_2D_proto_example tug)
target_link_libraries(reference-FTCS_2D_closed tug)
target_link_libraries(profiling_openmp tug)

27
examples/CRNI_2D_proto_example.cpp Normal file
View File

@ -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();
}

13
examples/FTCS_1D_proto_example.cpp
View File

@ -1,6 +1,9 @@
#include "tug/Boundary.hpp"
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
int main(int argc, char *argv[]) {
// **************
// **** GRID ****
@ -10,10 +13,9 @@ int main(int argc, char *argv[]) {
int cells = 20;
Grid grid = Grid(cells);
MatrixXd concentrations = MatrixXd::Constant(1,20,20);
MatrixXd concentrations = MatrixXd::Constant(1, 20, 20);
grid.setConcentrations(concentrations);
// ******************
// **** BOUNDARY ****
// ******************
@ -23,13 +25,13 @@ int main(int argc, char *argv[]) {
bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
Simulation simulation =
Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// (optional) set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
@ -37,7 +39,8 @@ int main(int argc, char *argv[]) {
// (optional) set the number of iterations
simulation.setIterations(100);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_OFF);
// **** RUN SIMULATION ****

33
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -1,14 +1,18 @@
/**
* @file FTCS_2D_proto_closed_mdl.cpp
* @author Hannes Signer, Philipp Ungrund, MDL
* @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary condition; optional command line argument: number of cols and rows
* @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary
* condition; optional command line argument: number of cols and rows
*
*/
#include <Eigen/Eigen>
#include <cstdlib>
#include <iostream>
#include <tug/Simulation.hpp>
using namespace Eigen;
int main(int argc, char *argv[]) {
int row = 64;
@ -17,7 +21,7 @@ int main(int argc, char *argv[]) {
// no cmd line argument, take col=row=64
row = atoi(argv[1]);
}
int col=row;
int col = row;
std::cout << "Nrow =" << row << std::endl;
// **************
@ -25,19 +29,20 @@ int main(int argc, char *argv[]) {
// **************
// create a grid with a 20 x 20 field
int n2 = row/2-1;
Grid grid = Grid(row,col);
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;
// 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.:
@ -45,7 +50,6 @@ int main(int argc, char *argv[]) {
MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value
grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
@ -59,13 +63,13 @@ int main(int argc, char *argv[]) {
bc.setBoundarySideClosed(BC_SIDE_TOP);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
Simulation simulation =
Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(10000); // timestep
@ -73,7 +77,8 @@ int main(int argc, char *argv[]) {
// set the number of iterations
simulation.setIterations(100);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****

29
examples/FTCS_2D_proto_example.cpp
View File

@ -1,16 +1,18 @@
/**
* @file FTCS_2D_proto_example.cpp
* @author Hannes Signer, Philipp Ungrund
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach
* and constant boundary condition
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS
* approach and constant boundary condition
*
*/
#include <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();
@ -21,25 +23,23 @@ int main(int argc, char *argv[]) {
// create a grid with a 20 x 20 field
int row = 20;
int col = 20;
Grid grid = Grid(row,col);
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;
// 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 ****
// ******************
@ -51,7 +51,6 @@ int main(int argc, char *argv[]) {
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
@ -62,13 +61,13 @@ int main(int argc, char *argv[]) {
// 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
Simulation simulation =
Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
@ -76,10 +75,10 @@ int main(int argc, char *argv[]) {
// set the number of iterations
simulation.setIterations(10000);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****
// run the simulation

32
examples/FTCS_2D_proto_example_mdl.cpp
View File

@ -1,13 +1,16 @@
/**
* @file FTCS_2D_proto_example.cpp
* @author Hannes Signer, Philipp Ungrund
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach
* and constant boundary condition
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS
* approach and constant boundary condition
*
*/
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
int main(int argc, char *argv[]) {
// **************
@ -17,19 +20,20 @@ int main(int argc, char *argv[]) {
// 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);
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;
// 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.:
@ -37,7 +41,6 @@ int main(int argc, char *argv[]) {
MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value
grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
@ -51,13 +54,13 @@ int main(int argc, char *argv[]) {
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
Simulation simulation =
Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// (optional) set the timestep of the simulation
simulation.setTimestep(1000); // timestep
@ -65,7 +68,8 @@ int main(int argc, char *argv[]) {
// (optional) set the number of iterations
simulation.setIterations(5);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_OFF);
// **** RUN SIMULATION ****

51
examples/profiling_openmp.cpp
View File

@ -1,31 +1,39 @@
#include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
#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[4] = {100, 500, 1000, 2000};
int threads[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
int iterations[1] = {5};
int repetition = 1;
int n[] = {2000};
int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
int iterations[1] = {1};
int repetition = 10;
for (int l = 0; l < size(threads); l++) {
// string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
ofstream myfile;
myfile.open("testLarge.csv");
myfile.open("speedup_1000.csv", std::ios::app);
myfile << "Number threads: " << threads[l] << endl;
for (int i = 0; i < size(n); i++){
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++){
// 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++){
// 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;
concentrations(n[i] / 2, n[i] / 2) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
@ -34,28 +42,29 @@ int main(int argc, char *argv[]) {
Simulation sim = Simulation(grid, bc, BTCS_APPROACH);
sim.setSolver(THOMAS_ALGORITHM_SOLVER);
if(argc == 2){
if (argc == 2) {
int numThreads = atoi(argv[1]);
sim.setNumberThreads(numThreads);
}
else{
sim.setNumberThreads(1);
} else {
sim.setNumberThreads(threads[l]);
}
sim.setTimestep(0.001);
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);
auto milliseconds =
std::chrono::duration_cast<std::chrono::milliseconds>(end -
begin);
myfile << milliseconds.count() << endl;
}
}
cout << endl;
myfile << endl;
}
myfile.close();
}
}

70
examples/profiling_speedup.cpp Normal file
View File

@ -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();
}
}

9
examples/reference-FTCS_2D_closed.cpp
View File

@ -1,8 +1,9 @@
#include <tug/Simulation.hpp>
#include "Eigen/Core"
#include <iostream>
#include <tug/Simulation.hpp>
using namespace std;
using namespace Eigen;
int main(int argc, char *argv[]) {
int row = 50;
@ -10,13 +11,12 @@ int main(int argc, char *argv[]) {
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;
concentrations(5, 5) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
@ -37,11 +37,9 @@ int main(int argc, char *argv[]) {
}
grid.setAlpha(alpha, alpha);
// Boundary
Boundary bc = Boundary(grid);
// Simulation
Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
sim.setTimestep(0.001);
@ -49,7 +47,6 @@ int main(int argc, char *argv[]) {
sim.setOutputCSV(CSV_OUTPUT_OFF);
sim.setOutputConsole(CONSOLE_OUTPUT_OFF);
// RUN
sim.run();
}

86
include/tug/Boundary.hpp
View File

@ -1,37 +1,26 @@
/**
* @file Boundary.hpp
* @brief API of Boundary class, that holds all information for each boundary condition
* at the edges of the diffusion grid.
* @brief API of Boundary class, that holds all information for each boundary
* condition at the edges of the diffusion grid.
*
*/
#ifndef BOUNDARY_H_
#define BOUNDARY_H_
#include <cstddef>
#include "Grid.hpp"
using namespace std;
using namespace Eigen;
#include <cstddef>
/**
* @brief Enum defining the two implemented boundary conditions.
*
*/
enum BC_TYPE {
BC_TYPE_CLOSED,
BC_TYPE_CONSTANT
};
enum BC_TYPE { BC_TYPE_CLOSED, BC_TYPE_CONSTANT };
/**
* @brief Enum defining all 4 possible sides to a 1D and 2D grid.
*
*/
enum BC_SIDE {
BC_SIDE_LEFT,
BC_SIDE_RIGHT,
BC_SIDE_TOP,
BC_SIDE_BOTTOM
};
enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM };
/**
* This class defines the boundary conditions of individual boundary elements.
@ -39,12 +28,12 @@ enum BC_SIDE {
* The class serves as an auxiliary class for structuring the Boundary class.
*/
class BoundaryElement {
public:
public:
/**
* @brief Construct a new Boundary Element object for the closed case.
* The boundary type is here automatically set to the type
* BC_TYPE_CLOSED, where the value takes NaN.
* BC_TYPE_CLOSED, where the value takes -1 and does not hold any
* physical meaning.
*/
BoundaryElement();
@ -91,23 +80,23 @@ class BoundaryElement {
*/
double getValue();
private:
BC_TYPE type;
double value;
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:
public:
/**
* @brief Creates a boundary object based on the passed grid object and
* initializes the boundaries as closed.
*
* @param grid Grid object on the basis of which the simulation is to take place.
* @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);
@ -125,7 +114,8 @@ class Boundary {
* 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.
* @param value Concentration to be set for all elements of the specified
* page.
*/
void setBoundarySideConstant(BC_SIDE side, double value);
@ -135,7 +125,8 @@ class Boundary {
*
* @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.
* boundary side. Must index an element of the corresponding
* side.
*/
void setBoundaryElementClosed(BC_SIDE side, int index);
@ -146,8 +137,10 @@ class Boundary {
*
* @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.
* 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);
@ -155,37 +148,42 @@ class Boundary {
* @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.
* @param side Boundary side from which the boundary conditions are to be
* returned.
* @return vector<BoundaryElement> Contains the boundary conditions as
* BoundaryElement objects.
*/
vector<BoundaryElement> getBoundarySide(BC_SIDE side);
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.
* @brief Get thes Boundary Side Values as a vector. Value is -1 in case some
specific boundary is closed.
*
* @param side Boundary side for which the values are to be returned.
* @return VectorXd Vector with values as doubles.
*/
VectorXd getBoundarySideValues(BC_SIDE side);
Eigen::VectorXd getBoundarySideValues(BC_SIDE side);
/**
* @brief Returns the boundary condition of a specified element on a given 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.
* 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.
* @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.
* 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);
@ -203,11 +201,11 @@ class Boundary {
*/
double getBoundaryElementValue(BC_SIDE side, int index);
private:
private:
Grid grid; // Boundary is directly dependent on the dimensions of a predefined
vector<vector<BoundaryElement>> boundaries; // Vector with Boundary Element information
std::vector<std::vector<BoundaryElement>>
boundaries; // Vector with Boundary Element information
};
#endif
#endif // BOUNDARY_H_

100
include/tug/Grid.hpp
View File

@ -1,3 +1,6 @@
#ifndef GRID_H_
#define GRID_H_
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
@ -8,24 +11,25 @@
#include <Eigen/Core>
#include <Eigen/Sparse>
using namespace Eigen;
// TODO document default values and perhaps adjust them
class Grid {
public:
public:
/**
* @brief Constructs a new 1D-Grid object of a given length, which holds a matrix
* with concentrations and a respective matrix of alpha coefficients.
* @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.
* @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.
@ -35,11 +39,11 @@ class Grid {
/**
* @brief Sets the concentrations matrix for a 1D or 2D-Grid.
*
* @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix must
* have correct dimensions as defined in row and col, or length,
* respectively.
* @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(MatrixXd concentrations);
void setConcentrations(const Eigen::MatrixXd &concentrations);
/**
* @brief Gets the concentrations matrix for a Grid.
@ -47,48 +51,51 @@ class Grid {
* @return MatrixXd An Eigen3 matrix holding the concentrations and having the
* same dimensions as the grid.
*/
MatrixXd getConcentrations();
const Eigen::MatrixXd &getConcentrations() { return this->concentrations; }
/**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one dimensional.
* @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.
* @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
* Matrix columns must have same size as length of grid.
*/
void setAlpha(MatrixXd alpha);
void setAlpha(const Eigen::MatrixXd &alpha);
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two dimensional.
* @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.
* @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in
* x-direction. Matrix must be of same size as the grid.
* @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid.
*/
void setAlpha(MatrixXd alphaX, MatrixXd alphaY);
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.
* @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @return MatrixXd A matrix with 1 row holding the alpha coefficients.
*/
MatrixXd getAlpha();
const Eigen::MatrixXd &getAlpha();
/**
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. Grid must be
* two dimensional.
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return MatrixXd A matrix holding the alpha coefficients in x-direction.
*/
MatrixXd getAlphaX();
const Eigen::MatrixXd &getAlphaX();
/**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. Grid must be
* two dimensional.
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return MatrixXd A matrix holding the alpha coefficients in y-direction.
*/
MatrixXd getAlphaY();
const Eigen::MatrixXd &getAlphaY();
/**
* @brief Gets the dimensions of the grid.
@ -128,10 +135,12 @@ class Grid {
/**
* @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.
* @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);
void setDomain(double domainRow, double domainCol);
/**
* @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
@ -154,18 +163,17 @@ class Grid {
*/
double getDeltaRow();
private:
private:
int col; // number of grid columns
int row; // number of grid rows
int row{1}; // number of grid rows
int dim; // 1D or 2D
double domainCol; // number of domain columns
double domainRow; // number of domain rows
double domainRow{0}; // number of domain rows
double deltaCol; // delta in x-direction (between columns)
double deltaRow; // delta in y-direction (between rows)
MatrixXd concentrations; // Matrix holding grid concentrations
MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
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_

110
include/tug/Simulation.hpp
View File

@ -1,14 +1,24 @@
/**
* @file Simulation.hpp
* @brief API of Simulation class, that holds all information regarding a specific simulation
* run like its timestep, number of iterations and output options. Simulation object
* also holds a predefined Grid and Boundary object.
* @brief API of Simulation class, that holds all information regarding a
* specific simulation run like its timestep, number of iterations and output
* options. Simulation object also holds a predefined Grid and Boundary object.
*
*/
#include "Boundary.hpp"
#include <ios>
using namespace std;
#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.
@ -26,7 +36,8 @@ enum APPROACH {
*/
enum SOLVER {
EIGEN_LU_SOLVER, // EigenLU solver
THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for tridiagonal matrices
THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for
// tridiagonal matrices
};
/**
@ -37,7 +48,8 @@ enum CSV_OUTPUT {
CSV_OUTPUT_OFF, // do not produce csv output
CSV_OUTPUT_ON, // produce csv output with last concentration matrix
CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices
CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary conditions at beginning
CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary
// conditions at beginning
};
/**
@ -66,12 +78,14 @@ enum TIME_MEASURE {
*
*/
class Simulation {
public:
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.
* @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
@ -96,21 +110,28 @@ class Simulation {
void setOutputCSV(CSV_OUTPUT csv_output);
/**
* @brief Set the options for outputting information to the console. Off by default.
* @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
* - CONSOLE_OUTPUT_OFF: do not print any output to
* console
* - CONSOLE_OUTPUT_ON: print before and after
* concentrations to console
* - CONSOLE_OUTPUT_VERBOSE: print all concentration
* matrices to console
*/
void setOutputConsole(CONSOLE_OUTPUT console_output);
// TODO document method
/**
* @brief Set the Time Measure object. Off by default.
* @brief Set the Time Measure option. Off by default.
*
* @param time_measure
* @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);
@ -138,11 +159,11 @@ class Simulation {
void setIterations(int iterations);
/**
* @brief Set the desired linear equation solver to be used for BTCS approach. Without effect
* in case of FTCS approach.
* @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.
* @param solver Solver to be used. Default is Thomas Algorithm as it is more
* efficient for tridiagonal Matrices.
*/
void setSolver(SOLVER solver);
@ -150,8 +171,9 @@ class Simulation {
* @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.
* 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);
@ -168,16 +190,16 @@ class Simulation {
*/
void printConcentrationsConsole();
// TODO move create CSVfile to TugUtils
/**
* @brief Creates a CSV file with a name containing the current simulation
* parameters. If the data name already exists, an additional counter is
* appended to the name. The name of the file is built up as follows:
* <approach> + <number rows> + <number columns> + <number of iterations>+<counter>.csv
* 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 given simulation parameter.
* @return string Filename with configured simulation parameters.
*/
string createCSVfile();
std::string createCSVfile();
/**
* @brief Writes the currently calculated concentration values of the grid
@ -186,7 +208,7 @@ class Simulation {
* @param filename Name of the file to which the concentration values are
* to be written.
*/
void printConcentrationsCSV(string filename);
void printConcentrationsCSV(const std::string &filename);
/**
* @brief Method starts the simulation process with the previously set
@ -194,19 +216,21 @@ class Simulation {
*/
void run();
private:
double timestep;
int iterations;
int innerIterations;
int numThreads;
CSV_OUTPUT csv_output;
CONSOLE_OUTPUT console_output;
TIME_MEASURE time_measure;
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;
SOLVER solver{THOMAS_ALGORITHM_SOLVER};
const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
};
#endif // SIMULATION_H_

14
scripts/Adi2D_Reference.R
View File

@ -286,7 +286,11 @@ ADIHetDir <- function(field, dt, iter, alpha) {
return(out)
}
harm <- function(x,y) 1/(1/x+1/y)
harm <- function(x,y) {
if (length(x) != 1 || length(y) != 1)
stop("x & z have different lengths")
2/(1/x+1/y)
}
harm(1,4)
@ -428,10 +432,10 @@ FTCS_2D <- function(field, dt, iter, alpha) {
## res[i,j] *(mean(alpha[i,j+1],alpha[i,j])+mean(alpha[i,j-1],alpha[i,j])) +
## res[i,j-1]*mean(alpha[i,j-1],alpha[i,j]))
tmp[i,j] <- res[i,j] +
+ tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * mean(alpha[i+1,j],alpha[i,j]) -
(res[i,j]-res[i-1,j]) * mean(alpha[i-1,j],alpha[i,j])) +
+ tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * mean(alpha[i,j+1],alpha[i,j]) -
(res[i,j]-res[i,j-1]) * mean(alpha[i,j-1],alpha[i,j]))
+ tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * harm(alpha[i+1,j],alpha[i,j]) -
(res[i,j]-res[i-1,j]) * harm(alpha[i-1,j],alpha[i,j])) +
+ tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * harm(alpha[i,j+1],alpha[i,j]) -
(res[i,j]-res[i,j-1]) * harm(alpha[i,j-1],alpha[i,j]))
}
}
## swap back tmp to res for the next inner iteration

0
proto/BTCS.ipynb → scripts/hannes-philip/BTCS.ipynb
View File

0
proto/FTCS.ipynb → scripts/hannes-philip/FTCS.ipynb
View File

418
src/BTCS.cpp Normal file
View File

@ -0,0 +1,418 @@
/**
* @file BTCSv2.cpp
* @brief Implementation of heterogenous BTCS (backward time-centered space)
* solution of diffusion equation in 1D and 2D space. Internally the
* alternating-direction implicit (ADI) method is used. Version 2, because
* Version 1 was an implementation for the homogeneous BTCS solution.
*
*/
#include "Schemes.hpp"
#include "TugUtils.hpp"
#include <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
// calculates coefficient for boundary in constant case
template <class T>
constexpr std::pair<T, T> calcBoundaryCoeffConstant(T alpha_center,
T alpha_side, T sx) {
const T centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) +
2 * alpha_center);
const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
return {centerCoeff, sideCoeff};
}
// calculates coefficient for boundary in closed case
template <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);
const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
return {centerCoeff, sideCoeff};
}
// creates coefficient matrix for next time step from alphas in x-direction
static Eigen::SparseMatrix<double>
createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
std::vector<BoundaryElement> &bcRight, int numCols,
int rowIndex, double sx) {
// square matrix of column^2 dimension for the coefficients
Eigen::SparseMatrix<double> cm(numCols, numCols);
cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
// left column
switch (bcLeft[rowIndex].getType()) {
case BC_TYPE_CONSTANT: {
auto [centerCoeffTop, rightCoeffTop] =
calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
break;
}
case BC_TYPE_CLOSED: {
auto [centerCoeffTop, rightCoeffTop] =
calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
break;
}
default: {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
}
}
// inner columns
int n = numCols - 1;
for (int i = 1; i < n; i++) {
cm.insert(i, i - 1) =
-sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
cm.insert(i, i) =
1 +
sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) +
calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)));
cm.insert(i, i + 1) =
-sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1));
}
// right column
switch (bcRight[rowIndex].getType()) {
case BC_TYPE_CONSTANT: {
auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
break;
}
case BC_TYPE_CLOSED: {
auto [centerCoeffBottom, leftCoeffBottom] =
calcBoundaryCoeffClosed(alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
break;
}
default: {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
}
}
cm.makeCompressed(); // important for Eigen solver
return cm;
}
// calculates explicit concentration at boundary in closed case
template <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);
}
}
// first row
else if (rowIndex == 0) {
for (int i = 0; i < length; i++) {
switch (bcTop[i].getType()) {
case BC_TYPE_CONSTANT: {
sv(i) = calcExplicitConcentrationsBoundaryConstant(
concentrations(rowIndex, i), bcTop[i].getValue(),
alphaY(rowIndex, i), alphaY(rowIndex + 1, i), sy);
break;
}
case BC_TYPE_CLOSED: {
sv(i) = calcExplicitConcentrationsBoundaryClosed(
concentrations(rowIndex, i), alphaY(rowIndex, i),
alphaY(rowIndex + 1, i), sy);
break;
}
default:
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
}
}
}
// last row
else if (rowIndex == numRows - 1) {
for (int i = 0; i < length; i++) {
switch (bcBottom[i].getType()) {
case BC_TYPE_CONSTANT: {
sv(i) = calcExplicitConcentrationsBoundaryConstant(
concentrations(rowIndex, i), bcBottom[i].getValue(),
alphaY(rowIndex, i), alphaY(rowIndex - 1, i), sy);
break;
}
case BC_TYPE_CLOSED: {
sv(i) = calcExplicitConcentrationsBoundaryClosed(
concentrations(rowIndex, i), alphaY(rowIndex, i),
alphaY(rowIndex - 1, i), sy);
break;
}
default:
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
}
}
}
// first column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
}
// last column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(length - 1) +=
2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
}
return sv;
}
// solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector
// use of EigenLU solver
static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<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!");
}
}

487
src/BTCSv2.cpp
View File

@ -1,487 +0,0 @@
/**
* @file BTCSv2.cpp
* @brief Implementation of heterogenous BTCS (backward time-centered space)
* solution of diffusion equation in 1D and 2D space. Internally the
* alternating-direction implicit (ADI) method is used. Version 2, because
* Version 1 was an implementation for the homogeneous BTCS solution.
*
*/
#include "FTCS.cpp"
#include <omp.h>
#include <tug/Boundary.hpp>
#ifdef WRITE_THOMAS_CSV
#include <fstream>
#include <string>
#endif
#define NUM_THREADS_BTCS 10
using namespace Eigen;
// calculates coefficient for left boundary in constant case
static tuple<double, double>
calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff;
double rightCoeff;
centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) +
2 * alpha(rowIndex, 0));
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff};
}
// calculates coefficient for left boundary in closed case
static tuple<double, double>
calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff;
double rightCoeff;
centerCoeff =
1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff};
}
// calculates coefficient for right boundary in constant case
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha,
int rowIndex, int n,
double sx) {
double leftCoeff;
double centerCoeff;
leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) +
2 * alpha(rowIndex, n));
return {leftCoeff, centerCoeff};
}
// calculates coefficient for right boundary in closed case
static tuple<double, double>
calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
double leftCoeff;
double centerCoeff;
leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
centerCoeff =
1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
return {leftCoeff, centerCoeff};
}
// creates coefficient matrix for next time step from alphas in x-direction
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha,
vector<BoundaryElement> &bcLeft,
vector<BoundaryElement> &bcRight,
int numCols, int rowIndex,
double sx) {
// square matrix of column^2 dimension for the coefficients
SparseMatrix<double> cm(numCols, numCols);
cm.reserve(VectorXi::Constant(numCols, 3));
// left column
BC_TYPE type = bcLeft[rowIndex].getType();
if (type == BC_TYPE_CONSTANT) {
auto [centerCoeffTop, rightCoeffTop] =
calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx);
cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
} else if (type == BC_TYPE_CLOSED) {
auto [centerCoeffTop, rightCoeffTop] =
calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx);
cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop;
} else {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
}
// inner columns
int n = numCols - 1;
for (int i = 1; i < n; i++) {
cm.insert(i, i - 1) =
-sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
cm.insert(i, i) =
1 +
sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) +
calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)));
cm.insert(i, i + 1) =
-sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1));
}
// right column
type = bcRight[rowIndex].getType();
if (type == BC_TYPE_CONSTANT) {
auto [leftCoeffBottom, centerCoeffBottom] =
calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx);
cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
} else if (type == BC_TYPE_CLOSED) {
auto [leftCoeffBottom, centerCoeffBottom] =
calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx);
cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom;
} else {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
}
cm.makeCompressed(); // important for Eigen solver
return cm;
}
// calculates explicity concentration at top boundary in constant case
static double calcExplicitConcentrationsTopBoundaryConstant(
MatrixXd &concentrations, MatrixXd &alpha, vector<BoundaryElement> &bcTop,
int rowIndex, int i, double sy) {
double c;
c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
concentrations(rowIndex, i) +
(1 -
sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) +
2 * alpha(rowIndex, i))) *
concentrations(rowIndex, i) +
sy * alpha(rowIndex, i) * bcTop[i].getValue();
return c;
}
// calculates explicit concentration at top boundary in closed case
static double calcExplicitConcentrationsTopBoundaryClosed(
MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
double c;
c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
concentrations(rowIndex, i) +
(1 -
sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)))) *
concentrations(rowIndex, i);
return c;
}
// calculates explicit concentration at bottom boundary in constant case
static double calcExplicitConcentrationsBottomBoundaryConstant(
MatrixXd &concentrations, MatrixXd &alpha,
vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
double c;
c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() +
(1 -
sy * (2 * alpha(rowIndex, i) +
calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) *
concentrations(rowIndex, i) +
sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) *
concentrations(rowIndex - 1, i);
return c;
}
// calculates explicit concentration at bottom boundary in closed case
static double calcExplicitConcentrationsBottomBoundaryClosed(
MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) {
double c;
c = (1 -
sy * (+calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) *
concentrations(rowIndex, i) +
sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) *
concentrations(rowIndex - 1, i);
return c;
}
// creates a solution vector for next time step from the current state of
// concentrations
static VectorXd createSolutionVector(
MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom,
int length, int rowIndex, double sx, double sy) {
VectorXd sv(length);
int numRows = concentrations.rows();
BC_TYPE type;
// inner rows
if (rowIndex > 0 && rowIndex < numRows - 1) {
for (int i = 0; i < length; i++) {
sv(i) =
sy *
calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) *
concentrations(rowIndex + 1, i) +
(1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i),
alphaY(rowIndex + 1, i)) +
calcAlphaIntercell(alphaY(rowIndex - 1, i),
alphaY(rowIndex, i)))) *
concentrations(rowIndex, i) +
sy *
calcAlphaIntercell(alphaY(rowIndex - 1, i), alphaY(rowIndex, i)) *
concentrations(rowIndex - 1, i);
}
}
// first row
else if (rowIndex == 0) {
for (int i = 0; i < length; i++) {
type = bcTop[i].getType();
if (type == BC_TYPE_CONSTANT) {
sv(i) = calcExplicitConcentrationsTopBoundaryConstant(
concentrations, alphaY, bcTop, rowIndex, i, sy);
} else if (type == BC_TYPE_CLOSED) {
sv(i) = calcExplicitConcentrationsTopBoundaryClosed(
concentrations, alphaY, rowIndex, i, sy);
} else {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Left or Top!");
}
}
}
// last row
else if (rowIndex == numRows - 1) {
for (int i = 0; i < length; i++) {
type = bcBottom[i].getType();
if (type == BC_TYPE_CONSTANT) {
sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(
concentrations, alphaY, bcBottom, rowIndex, i, sy);
} else if (type == BC_TYPE_CLOSED) {
sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(
concentrations, alphaY, rowIndex, i, sy);
} else {
throw_invalid_argument(
"Undefined Boundary Condition Type somewhere on Right or Bottom!");
}
}
}
// first column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
}
// last column -> additional fixed concentration change from perpendicular
// dimension in constant bc case
if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
sv(length - 1) +=
2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
}
return sv;
}
// solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector
// use of EigenLU solver
static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
SparseLU<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 VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
uint32_t n = b.size();
Eigen::VectorXd a_diag(n);
Eigen::VectorXd b_diag(n);
Eigen::VectorXd c_diag(n);
Eigen::VectorXd x_vec = b;
// Fill diagonals vectors
b_diag[0] = A.coeff(0, 0);
c_diag[0] = A.coeff(0, 1);
for (int i = 1; i < n - 1; i++) {
a_diag[i] = A.coeff(i, i - 1);
b_diag[i] = A.coeff(i, i);
c_diag[i] = A.coeff(i, i + 1);
}
a_diag[n - 1] = A.coeff(n - 1, n - 2);
b_diag[n - 1] = A.coeff(n - 1, n - 1);
// HACK: write CSV to file
#ifdef WRITE_THOMAS_CSV
#include <fstream>
#include <string>
static std::uint32_t file_index = 0;
std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv";
std::ofstream out_file;
out_file.open(file_name, std::ofstream::trunc | std::ofstream::out);
// print header
out_file << "Aa, Ab, Ac, b\n";
// iterate through all elements
for (std::size_t i = 0; i < n; i++) {
out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", "
<< b[i] << "\n";
}
out_file.close();
#endif
// start solving - c_diag and x_vec are overwritten
n--;
c_diag[0] /= b_diag[0];
x_vec[0] /= b_diag[0];
for (int i = 1; i < n; i++) {
c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1];
x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) /
(b_diag[i] - a_diag[i] * c_diag[i - 1]);
}
x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) /
(b_diag[n] - a_diag[n] * c_diag[n - 1]);
for (int i = n; i-- > 0;) {
x_vec[i] -= c_diag[i] * x_vec[i + 1];
}
return x_vec;
}
// BTCS solution for 1D grid
static void BTCS_1D(Grid &grid, Boundary &bc, double timestep,
VectorXd (*solverFunc)(SparseMatrix<double> &A,
VectorXd &b)) {
int length = grid.getLength();
double sx = timestep / (grid.getDelta() * grid.getDelta());
VectorXd concentrations_t1(length);
SparseMatrix<double> A;
VectorXd b(length);
MatrixXd alpha = grid.getAlpha();
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
MatrixXd concentrations = grid.getConcentrations();
int rowIndex = 0;
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
sx); // this is exactly same as in 2D
for (int i = 0; i < length; i++) {
b(i) = concentrations(0, i);
}
if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue();
}
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
}
concentrations_t1 = solverFunc(A, b);
for (int j = 0; j < length; j++) {
concentrations(0, j) = concentrations_t1(j);
}
grid.setConcentrations(concentrations);
}
// BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double timestep,
VectorXd (*solverFunc)(SparseMatrix<double> &A,
VectorXd &b),
int numThreads) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
VectorXd row_t1(colMax);
SparseMatrix<double> A;
VectorXd b;
MatrixXd alphaX = grid.getAlphaX();
MatrixXd alphaY = grid.getAlphaY();
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
MatrixXd concentrations = grid.getConcentrations();
#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
for (int i = 0; i < rowMax; i++) {
A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy);
SparseLU<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
static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
} else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
} else {
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
}
}
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid
static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep,
int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
} else {
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
}
}

117
src/Boundary.cpp
View File

@ -1,31 +1,21 @@
#include "TugUtils.hpp"
#include <iostream>
#include <omp.h>
#include <tug/Boundary.hpp>
#include <cstddef>
#include <stdexcept>
#include <tug/Boundary.hpp>
using namespace std;
BoundaryElement::BoundaryElement() {}
BoundaryElement::BoundaryElement() {
BoundaryElement::BoundaryElement(double _value)
: value(_value), type(BC_TYPE_CONSTANT) {}
this->type = BC_TYPE_CLOSED;
this->value = NAN;
}
BoundaryElement::BoundaryElement(double value) {
this->type = BC_TYPE_CONSTANT;
this->value = value;
}
void BoundaryElement::setType(BC_TYPE type) {
this->type = type;
}
void BoundaryElement::setType(BC_TYPE type) { this->type = type; }
void BoundaryElement::setValue(double value) {
if(value < 0){
if (value < 0) {
throw_invalid_argument("No negative concentration allowed.");
}
if(type == BC_TYPE_CLOSED){
if (type == BC_TYPE_CLOSED) {
throw_invalid_argument(
"No constant boundary concentrations can be set for closed "
"boundaries. Please change type first.");
@ -33,88 +23,83 @@ void BoundaryElement::setValue(double value) {
this->value = value;
}
BC_TYPE BoundaryElement::getType() {
return this->type;
}
BC_TYPE BoundaryElement::getType() { return this->type; }
double BoundaryElement::getValue() {
return this->value;
}
double BoundaryElement::getValue() { return this->value; }
Boundary::Boundary(Grid grid) : grid(grid) {
//probably to DEBUG assignment grid
if (grid.getDim() == 1) {
this->boundaries = vector<vector<BoundaryElement>>(2);
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 = vector<vector<BoundaryElement>>(4);
this->boundaries = std::vector<std::vector<BoundaryElement>>(4);
this->boundaries[BC_SIDE_LEFT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_RIGHT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_TOP] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());
this->boundaries[BC_SIDE_BOTTOM] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());
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)){
if (grid.getDim() == 1) {
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument(
"For the one-dimensional case, only the BC_SIDE_LEFT and "
"BC_SIDE_RIGHT borders exist.");
}
}
int n;
if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) {
n = grid.getRow();
} else {
n = grid.getCol();
}
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement());
const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
const int n = is_vertical ? grid.getRow() : grid.getCol();
this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement());
}
void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
if(grid.getDim() == 1){
if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){
if (grid.getDim() == 1) {
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument(
"For the one-dimensional case, only the BC_SIDE_LEFT and "
"BC_SIDE_RIGHT borders exist.");
}
}
int n;
if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) {
n = grid.getRow();
} else {
n = grid.getCol();
}
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement(value));
const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
const int n = is_vertical ? grid.getRow() : grid.getCol();
this->boundaries[side] =
std::vector<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){
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) {
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){
if ((boundaries[side].size() < index) || index < 0) {
throw_invalid_argument("Index is selected either too large or too small.");
}
this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
this->boundaries[side][index].setValue(value);
}
vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
if(grid.getDim() == 1){
if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){
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.");
@ -123,9 +108,9 @@ vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
return this->boundaries[side];
}
VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
int length = boundaries[side].size();
VectorXd values(length);
Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
const std::size_t length = boundaries[side].size();
Eigen::VectorXd values(length);
for (int i = 0; i < length; i++) {
if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {
@ -139,26 +124,26 @@ VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
}
BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) {
if((boundaries[side].size() < index) || index < 0){
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){
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){
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.");
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();
}

4
src/CMakeLists.txt
View File

@ -1,4 +1,4 @@
add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp)
add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp)
IF(TUG_NAAICE_EXAMPLE)
target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV)
@ -10,4 +10,4 @@ if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND)
target_link_libraries(tug OpenMP::OpenMP_CXX)
endif()
target_include_directories(tug PUBLIC ../include)
target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include)

456
src/FTCS.cpp
View File

@ -1,87 +1,85 @@
/**
* @file FTCS.cpp
* @brief Implementation of heterogenous FTCS (forward time-centered space) solution
* of diffusion equation in 1D and 2D space.
* @brief Implementation of heterogenous FTCS (forward time-centered space)
* solution of diffusion equation in 1D and 2D space.
*
*/
#include "Schemes.hpp"
#include "TugUtils.hpp"
#include <cstddef>
#include <tug/Boundary.hpp>
#include <iostream>
#include <tug/Boundary.hpp>
#ifdef _OPENMP
#include <omp.h>
#define NUM_THREADS 6
using namespace std;
#else
#define omp_get_thread_num() 0
#endif
// calculates horizontal change on one cell independent of boundary type
static double calcHorizontalChange(Grid &grid, int &row, int &col) {
static inline double calcHorizontalChange(Grid &grid, int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col+1)
- (
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
)
* grid.getConcentrations()(row,col)
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col-1);
return result;
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col))) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1);
}
// calculates vertical change on one cell independent of boundary type
static double calcVerticalChange(Grid &grid, int &row, int &col) {
static inline double calcVerticalChange(Grid &grid, int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row+1,col)
- (
calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+ calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
)
* grid.getConcentrations()(row,col)
+ calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row-1,col);
return result;
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col))) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row - 1, col);
}
// calculates horizontal change on one cell with a constant left boundary
static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) {
static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid,
Boundary &bc,
int &row,
int &col) {
double result =
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col+1)
- (
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+ 2 * grid.getAlphaX()(row,col)
)
* grid.getConcentrations()(row,col)
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
return result;
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) +
2 * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col) +
2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
}
// calculates horizontal change on one cell with a closed left boundary
static double calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
static inline double
calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col))
* (grid.getConcentrations()(row, col+1) - grid.getConcentrations()(row, col));
return result;
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col + 1) -
grid.getConcentrations()(row, col));
}
// checks boundary condition type for a cell on the left edge of grid
static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, int &row, int &col) {
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) {
@ -91,37 +89,36 @@ static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, int &ro
}
}
// calculates horizontal change on one cell with a constant right boundary
static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) {
static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
Boundary &bc,
int &row,
int &col) {
double result =
2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row)
- (
calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+ 2 * grid.getAlphaX()(row,col)
)
* grid.getConcentrations()(row,col)
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
* grid.getConcentrations()(row,col-1);
return result;
return 2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) +
2 * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1);
}
// calculates horizontal change on one cell with a closed right boundary
static double calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
static inline double
calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
double result =
- (calcAlphaIntercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row, col-1)));
return result;
return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col) -
grid.getConcentrations()(row, col - 1)));
}
// checks boundary condition type for a cell on the right edge of grid
static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, int &row, int &col) {
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) {
@ -131,37 +128,35 @@ static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, int &r
}
}
// calculates vertical change on one cell with a constant top boundary
static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) {
static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid,
Boundary &bc,
int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col))
* grid.getConcentrations()(row+1,col)
- (
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col))
+ 2 * grid.getAlphaY()(row, col)
)
* grid.getConcentrations()(row, col)
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
return result;
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
2 * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row, col) +
2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_TOP, col);
}
// calculates vertical change on one cell with a closed top boundary
static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, int &col) {
static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getConcentrations()(row, col))
* (grid.getConcentrations()(row+1, col) - grid.getConcentrations()(row, col));
return result;
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getConcentrations()(row, col)) *
(grid.getConcentrations()(row + 1, col) -
grid.getConcentrations()(row, col));
}
// checks boundary condition type for a cell on the top edge of grid
static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row, int &col) {
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) {
@ -171,37 +166,36 @@ static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row,
}
}
// calculates vertical change on one cell with a constant bottom boundary
static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) {
static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid,
Boundary &bc,
int &row,
int &col) {
double result =
2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col)
- (
calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
+ 2 * grid.getAlphaY()(row, col)
)
* grid.getConcentrations()(row, col)
+ calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
* grid.getConcentrations()(row-1,col);
return result;
return 2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) +
2 * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) *
grid.getConcentrations()(row - 1, col);
}
// calculates vertical change on one cell with a closed bottom boundary
static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
static inline double
calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
double result =
- (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col)));
return result;
return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) *
(grid.getConcentrations()(row, col) -
grid.getConcentrations()(row - 1, col)));
}
// checks boundary condition type for a cell on the bottom edge of grid
static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, int &row, int &col) {
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) {
@ -211,79 +205,66 @@ static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, int &ro
}
}
// 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
MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0);
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)
)
;
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)
)
;
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)
)
;
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) {
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
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
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)
)
;
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));
}
}
@ -292,144 +273,109 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep, int numThreads)
// 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)
)
;
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;
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)
)
;
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)
)
;
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;
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)
)
;
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)
)
;
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)
)
;
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;
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)
)
;
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)
)
;
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
static void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
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!");
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
}
}

123
src/Grid.cpp
View File

@ -1,137 +1,140 @@
#include "TugUtils.hpp"
#include <tug/Grid.hpp>
#include <iostream>
#include <tug/Grid.hpp>
Grid::Grid(int length) {
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.");
throw_invalid_argument(
"Given grid length too small. Must be greater than 3.");
}
this->row = 1;
this->col = length;
this->domainCol = length; // default: same size as length
this->deltaCol = double(this->domainCol)/double(this->col); // -> 1
this->dim = 1;
this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
this->concentrations = MatrixXd::Constant(1, col, 20);
this->alphaX = MatrixXd::Constant(1, 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) {
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.");
throw_invalid_argument(
"Given grid dimensions too small. Must each be greater than 3.");
}
this->row = row;
this->col = col;
this->domainRow = row; // default: same size as row
this->domainCol = col; // default: same size as col
this->deltaRow = double(this->domainRow)/double(this->row); // -> 1
this->deltaCol = double(this->domainCol)/double(this->col); // -> 1
this->dim = 2;
this->deltaRow = double(this->domainRow) / double(this->row); // -> 1
this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
this->concentrations = MatrixXd::Constant(row, col, 20);
this->alphaX = MatrixXd::Constant(row, col, 1);
this->alphaY = MatrixXd::Constant(row, col, 1);
this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
}
void Grid::setConcentrations(MatrixXd concentrations) {
if (concentrations.rows() != this->row || concentrations.cols() != this->col) {
throw_invalid_argument("Given matrix of concentrations mismatch with Grid dimensions!");
void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) {
if (concentrations.rows() != this->row ||
concentrations.cols() != this->col) {
throw_invalid_argument(
"Given matrix of concentrations mismatch with Grid dimensions!");
}
this->concentrations = concentrations;
}
MatrixXd Grid::getConcentrations() {
return this->concentrations;
}
void Grid::setAlpha(MatrixXd alpha) {
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!");
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!");
throw_invalid_argument(
"Given matrix of alpha coefficients mismatch with Grid dimensions!");
}
this->alphaX = alpha;
}
void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) {
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!");
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!");
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!");
throw_invalid_argument("Given matrix of alpha coefficients in y-direction "
"mismatch with GRid dimensions!");
}
this->alphaX = alphaX;
this->alphaY = alphaY;
}
MatrixXd Grid::getAlpha() {
const Eigen::MatrixXd &Grid::getAlpha() {
if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probably use either getAlphaX() or getAlphaY()!");
throw_invalid_argument("Grid is not one dimensional, you should probably "
"use either getAlphaX() or getAlphaY()!");
}
return this->alphaX;
}
MatrixXd Grid::getAlphaX() {
const Eigen::MatrixXd &Grid::getAlphaX() {
if (dim != 2) {
throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!");
throw_invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!");
}
return this->alphaX;
}
MatrixXd Grid::getAlphaY() {
const Eigen::MatrixXd &Grid::getAlphaY() {
if (dim != 2) {
throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!");
throw_invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!");
}
return this->alphaY;
}
int Grid::getDim() {
return dim;
}
int Grid::getDim() { return dim; }
int Grid::getLength() {
if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probably use getRow() or getCol()!");
throw_invalid_argument("Grid is not one dimensional, you should probably "
"use getRow() or getCol()!");
}
return col;
}
int Grid::getRow() {
return row;
}
int Grid::getRow() { return row; }
int Grid::getCol() {
return col;
}
int Grid::getCol() { return col; }
void Grid::setDomain(double domainLength) {
if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probaly use the 2D domain setter!");
throw_invalid_argument("Grid is not one dimensional, you should probaly "
"use the 2D domain setter!");
}
if (domainLength <= 0) {
throw_invalid_argument("Given domain length is not positive!");
}
this->domainCol = domainLength;
this->deltaCol = double(this->domainCol)/double(this->col);
this->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!");
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!");
@ -139,25 +142,25 @@ void Grid::setDomain(double domainRow, double domainCol) {
this->domainRow = domainRow;
this->domainCol = domainCol;
this->deltaRow = double(this->domainRow)/double(this->row);
this->deltaCol = double(this->domainCol)/double(this->col);
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");
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::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!");
throw_invalid_argument("Grid is not two dimensional, meaning there is no "
"delta in y-direction!");
}
return this->deltaRow;

28
src/Schemes.hpp Normal file
View File

@ -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_

231
src/Simulation.cpp
View File

@ -1,34 +1,19 @@
#include <cmath>
#include <cstddef>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <tug/Simulation.hpp>
#include <omp.h>
#include <fstream>
#include "Schemes.hpp"
#include "TugUtils.hpp"
#ifndef SIMULATION_H_
#define SIMULATION_H_
#include "BTCSv2.cpp"
using namespace std;
Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) {
this->approach = approach;
this->solver = THOMAS_ALGORITHM_SOLVER;
this->timestep = -1; // error per default
this->iterations = -1;
this->innerIterations = 1;
this->numThreads = omp_get_num_procs();
this->csv_output = CSV_OUTPUT_OFF;
this->console_output = CONSOLE_OUTPUT_OFF;
this->time_measure = TIME_MEASURE_OFF;
}
Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach)
: grid(_grid), bc(_bc), approach(_approach) {}
void Simulation::setOutputCSV(CSV_OUTPUT csv_output) {
if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
@ -39,7 +24,8 @@ void Simulation::setOutputCSV(CSV_OUTPUT csv_output) {
}
void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) {
if (console_output < CONSOLE_OUTPUT_OFF && console_output > CONSOLE_OUTPUT_VERBOSE) {
if (console_output < CONSOLE_OUTPUT_OFF &&
console_output > CONSOLE_OUTPUT_VERBOSE) {
throw_invalid_argument("Invalid console output option given!");
}
@ -55,76 +41,81 @@ void Simulation::setTimeMeasure(TIME_MEASURE time_measure) {
}
void Simulation::setTimestep(double timestep) {
if(timestep <= 0){
if (timestep <= 0) {
throw_invalid_argument("Timestep has to be greater than zero.");
}
if (approach == FTCS_APPROACH) {
if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) {
double deltaRowSquare;
double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
double minDeltaSquare;
double maxAlphaX, maxAlphaY, maxAlpha;
string dim;
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) {
dim = "2D";
deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare;
maxAlphaX = grid.getAlphaX().maxCoeff();
maxAlphaY = grid.getAlphaY().maxCoeff();
const double maxAlphaX = grid.getAlphaX().maxCoeff();
const double maxAlphaY = grid.getAlphaY().maxCoeff();
maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
} else if (grid.getDim() == 1) {
dim = "1D";
minDeltaSquare = deltaColSquare;
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);
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)));
// stability equation from Wikipedia; might be useful if applied cfl does
// not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha *
// ((1/deltaRowSquare) + (1/deltaColSquare)));
cout << "FTCS_" << dim << " :: CFL condition MDL: " << cfl << endl;
cout << "FTCS_" << dim << " :: required dt=" << timestep << endl;
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;
cerr << "Warning :: Timestep was adjusted, because of stability "
std::cerr << "Warning :: Timestep was adjusted, because of stability "
"conditions. Time duration was approximately preserved by "
"adjusting internal number of iterations."
<< endl;
cout << "FTCS_" << dim << " :: Required " << this->innerIterations
<< " inner iterations with dt=" << this->timestep << endl;
<< std::endl;
std::cout << approachPrefix << "_" << dim << " :: Required "
<< this->innerIterations
<< " inner iterations with dt=" << this->timestep << std::endl;
} else {
this->timestep = timestep;
cout << "FTCS_" << dim << " :: No inner iterations required, dt=" << timestep << endl;
std::cout << approachPrefix << "_" << dim
<< " :: No inner iterations required, dt=" << timestep
<< std::endl;
}
} else {
this->timestep = timestep;
}
}
double Simulation::getTimestep() {
return this->timestep;
}
double Simulation::getTimestep() { return this->timestep; }
void Simulation::setIterations(int iterations) {
if(iterations <= 0){
if (iterations <= 0) {
throw_invalid_argument("Number of iterations must be greater than zero.");
}
this->iterations = iterations;
@ -132,50 +123,55 @@ void Simulation::setIterations(int iterations) {
void Simulation::setSolver(SOLVER solver) {
if (this->approach == FTCS_APPROACH) {
cerr << "Warning: Solver was set, but FTCS approach initialized. Setting the solver "
std::cerr
<< "Warning: Solver was set, but FTCS approach initialized. Setting "
"the solver "
"is thus without effect."
<< endl;
<< std::endl;
}
this->solver = solver;
}
void Simulation::setNumberThreads(int numThreads){
if(numThreads>0 && numThreads<=omp_get_num_procs()){
this->numThreads=numThreads;
}
else{
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("Number of threads exceeds the number of processor cores or is less than 1.");
throw_invalid_argument(outputMessage);
}
}
int Simulation::getIterations() {
return this->iterations;
}
int Simulation::getIterations() { return this->iterations; }
void Simulation::printConcentrationsConsole() {
cout << grid.getConcentrations() << endl;
cout << endl;
std::cout << grid.getConcentrations() << std::endl;
std::cout << std::endl;
}
string Simulation::createCSVfile() {
ofstream file;
std::string Simulation::createCSVfile() {
std::ofstream file;
int appendIdent = 0;
string appendIdentString;
std::string appendIdentString;
string approachString = (approach == 0) ? "FTCS" : "BTCS";
string row = to_string(grid.getRow());
string col = to_string(grid.getCol());
string numIterations = to_string(iterations);
// string approachString = (approach == 0) ? "FTCS" : "BTCS";
const std::string &approachString = this->approach_names[approach];
std::string row = std::to_string(grid.getRow());
std::string col = std::to_string(grid.getCol());
std::string numIterations = std::to_string(iterations);
string filename = approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
std::string filename =
approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
while (filesystem::exists(filename)) {
while (std::filesystem::exists(filename)) {
appendIdent += 1;
appendIdentString = to_string(appendIdent);
filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv";
appendIdentString = std::to_string(appendIdent);
filename = approachString + "_" + row + "_" + col + "_" + numIterations +
"-" + appendIdentString + ".csv";
}
file.open(filename);
@ -183,15 +179,20 @@ string Simulation::createCSVfile() {
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
// adds lines at the beginning of verbose output csv that represent the
// boundary conditions and their values -1 in case of closed boundary
if (csv_output == CSV_OUTPUT_XTREME) {
IOFormat one_row(StreamPrecision, DontAlignCols, "", " ");
file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) << endl; // boundary left
file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) << endl; // boundary right
file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) << endl; // boundary top
file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) << endl; // boundary bottom
file << endl << endl;
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();
@ -199,17 +200,17 @@ string Simulation::createCSVfile() {
return filename;
}
void Simulation::printConcentrationsCSV(string filename) {
ofstream file;
void Simulation::printConcentrationsCSV(const std::string &filename) {
std::ofstream file;
file.open(filename, std::ios_base::app);
if (!file) {
exit(1);
}
IOFormat do_not_align(StreamPrecision, DontAlignCols);
file << grid.getConcentrations().format(do_not_align) << endl;
file << endl << endl;
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();
}
@ -221,7 +222,7 @@ void Simulation::run() {
throw_invalid_argument("Number of iterations are not set!");
}
string filename;
std::string filename;
if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole();
}
@ -244,8 +245,8 @@ void Simulation::run() {
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) {
// double percentage = (double)i / ((double)iterations *
// (double)innerIterations) * 100; if ((int)percentage % 10 == 0) {
// cout << "Progress: " << percentage << "%" << endl;
// }
// }
@ -263,7 +264,6 @@ void Simulation::run() {
}
BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads);
}
} else if (solver == THOMAS_ALGORITHM_SOLVER) {
for (int i = 0; i < iterations; i++) {
@ -275,18 +275,41 @@ void Simulation::run() {
}
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
}
}
} else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
// TODO
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);
auto milliseconds =
std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole();
@ -295,11 +318,9 @@ void Simulation::run() {
printConcentrationsCSV(filename);
}
if (this->time_measure > TIME_MEASURE_OFF) {
string approachString = (approach == 0) ? "FTCS" : "BTCS";
string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << endl;
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;
}
}
#endif

12
src/TugUtils.hpp
View File

@ -1,5 +1,5 @@
#ifndef BTCSUTILS_H_
#define BTCSUTILS_H_
#ifndef TUGUTILS_H_
#define TUGUTILS_H_
#include <chrono>
#include <stdexcept>
@ -23,13 +23,15 @@
start); \
duration.count(); \
})
#endif // BTCSUTILS_H_
// calculates arithmetic or harmonic mean of alpha between two cells
static double calcAlphaIntercell(double &alpha1, double &alpha2, bool useHarmonic = true) {
template <typename T>
constexpr T calcAlphaIntercell(T alpha1, T alpha2,
bool useHarmonic = true) {
if (useHarmonic) {
return double(2) / ((double(1)/alpha1) + (double(1)/alpha2));
return double(2) / ((double(1) / alpha1) + (double(1) / alpha2));
} else {
return 0.5 * (alpha1 + alpha2);
}
}
#endif // TUGUTILS_H_

6
test/CMakeLists.txt
View File

@ -11,15 +11,15 @@ if(NOT DOCTEST_LIB)
FetchContent_MakeAvailable(DocTest)
endif()
add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp)
add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp)
target_link_libraries(testTug doctest tug)
# get relative path of the CSV file
get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH CACHE)
get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH)
# set relative path in header file
configure_file(testSimulation.hpp.in testSimulation.hpp)
# include test directory with generated header file from above
target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}")
target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src")
add_custom_target(
check

47
test/TestUtils.cpp
View File

@ -1,47 +0,0 @@
#include <ios>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <fstream>
#include <sstream>
#include <stdexcept>
using namespace std;
using namespace Eigen;
MatrixXd CSV2Eigen(string file2Convert) {
vector<double> matrixEntries;
ifstream matrixDataFile(file2Convert);
if (matrixDataFile.fail()) {
throw invalid_argument("File probably non-existent!");
}
string matrixRowString;
string matrixEntry;
int matrixRowNumber = 0;
while(getline(matrixDataFile, matrixRowString)){
stringstream matrixRowStringStream(matrixRowString);
while(getline(matrixRowStringStream, matrixEntry, ' ')){
matrixEntries.push_back(stod(matrixEntry));
}
if (matrixRowString.length() > 1) {
matrixRowNumber++;
}
}
return Map<Matrix<double, Dynamic, Dynamic, RowMajor>>(matrixEntries.data(), matrixRowNumber, matrixEntries.size() / matrixRowNumber);
}
bool checkSimilarity(MatrixXd a, MatrixXd b, double precision=1e-5) {
return a.isApprox(b, precision);
}
bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) {
MatrixXd diff = a - b;
double maxCoeff = diff.maxCoeff();
return abs(maxCoeff) < maxDiff;
}

49
test/TestUtils.hpp Normal file
View File

@ -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;
}

42
test/testBoundary.cpp
View File

@ -1,21 +1,23 @@
#include <stdio.h>
#include <doctest/doctest.h>
#include <tug/Boundary.hpp>
#include <string>
#include <typeinfo>
#include <iostream>
#include <stdio.h>
#include <string>
#include <tug/Boundary.hpp>
#include <typeinfo>
TEST_CASE("BoundaryElement"){
using namespace std;
SUBCASE("Closed case"){
TEST_CASE("BoundaryElement") {
SUBCASE("Closed case") {
BoundaryElement boundaryElementClosed = BoundaryElement();
CHECK_NOTHROW(BoundaryElement());
CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN));
CHECK_EQ(boundaryElementClosed.getValue(), -1);
CHECK_THROWS(boundaryElementClosed.setValue(0.2));
}
SUBCASE("Constant case"){
SUBCASE("Constant case") {
BoundaryElement boundaryElementConstant = BoundaryElement(0.1);
CHECK_NOTHROW(BoundaryElement(0.1));
CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT);
@ -25,18 +27,19 @@ TEST_CASE("BoundaryElement"){
}
}
TEST_CASE("Boundary Class"){
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"){
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_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));
@ -44,17 +47,20 @@ TEST_CASE("Boundary Class"){
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());
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"){
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_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));
@ -62,7 +68,9 @@ TEST_CASE("Boundary Class"){
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());
CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
BC_TYPE_CONSTANT);
CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
boundary1DVector[0].getType());
}
}

13
test/testFTCS.cpp
View File

@ -1,5 +1,6 @@
#include <TugUtils.hpp>
#include <doctest/doctest.h>
#include <../src/FTCS.cpp>
#include <limits>
TEST_CASE("Maths") {
@ -7,13 +8,13 @@ TEST_CASE("Maths") {
double alpha1 = 10;
double alpha2 = 20;
double average = 15;
double harmonicMean = double(2) / ((double(1)/alpha1)+(double(1)/alpha2));
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());
// 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);
}
}

38
test/testGrid.cpp
View File

@ -1,6 +1,10 @@
#include <Eigen/Core>
#include <doctest/doctest.h>
#include <tug/Grid.hpp>
#include <Eigen/Core>
using namespace Eigen;
using namespace std;
TEST_CASE("1D Grid, too small length") {
int l = 2;
@ -74,7 +78,7 @@ TEST_CASE("1D Grid") {
CHECK_THROWS(grid.setDomain(d, d));
grid.setDomain(d);
CHECK_EQ(grid.getDeltaCol(), double(d)/double(l));
CHECK_EQ(grid.getDeltaCol(), double(d) / double(l));
CHECK_THROWS(grid.getDeltaRow());
// set too small domain
@ -112,13 +116,12 @@ TEST_CASE("2D Grid quadratic") {
MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2);
CHECK_NOTHROW(grid.setConcentrations(concentrations));
// false concentrations matrix
MatrixXd wConcentrations = MatrixXd::Constant(rc, rc+3, 1);
MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(rc+3, rc, 4);
wConcentrations = MatrixXd::Constant(rc + 3, rc, 4);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(rc+2, rc+4, 6);
wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
@ -136,9 +139,9 @@ TEST_CASE("2D Grid quadratic") {
CHECK_THROWS(grid.getAlpha());
// false alpha matrices
alphax = MatrixXd::Constant(rc+3, rc+1, 3);
alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
alphay = MatrixXd::Constant(rc+2, rc+1, 3);
alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
@ -153,8 +156,8 @@ TEST_CASE("2D Grid quadratic") {
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));
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
// set too small domain
dr = 0;
@ -195,13 +198,12 @@ TEST_CASE("2D Grid non-quadratic") {
MatrixXd concentrations = MatrixXd::Constant(r, c, 2);
CHECK_NOTHROW(grid.setConcentrations(concentrations));
// false concentrations matrix
MatrixXd wConcentrations = MatrixXd::Constant(r, c+3, 6);
MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(r+3, c, 3);
wConcentrations = MatrixXd::Constant(r + 3, c, 3);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(r+2, c+4, 2);
wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
@ -219,9 +221,9 @@ TEST_CASE("2D Grid non-quadratic") {
CHECK_THROWS(grid.getAlpha());
// false alpha matrices
alphax = MatrixXd::Constant(r+3, c+1, 3);
alphax = MatrixXd::Constant(r + 3, c + 1, 3);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
alphay = MatrixXd::Constant(r+2, c+1, 5);
alphay = MatrixXd::Constant(r + 2, c + 1, 5);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
@ -236,8 +238,8 @@ TEST_CASE("2D Grid non-quadratic") {
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));
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
// set too small domain
dr = 0;

38
test/testSimulation.cpp
View File

@ -1,25 +1,29 @@
#include <stdio.h>
#include "TestUtils.hpp"
#include <doctest/doctest.h>
#include <tug/Simulation.hpp>
#include "TestUtils.cpp"
#include <stdio.h>
#include <string>
#include <tug/Simulation.hpp>
// include the configured header file
#include <testSimulation.hpp>
static Grid setupSimulation(APPROACH approach, double timestep, int iterations) {
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;
concentrations(5, 5) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
@ -40,28 +44,27 @@ static Grid setupSimulation(APPROACH approach, double timestep, int iterations)
}
grid.setAlpha(alpha, alpha);
// Boundary
Boundary bc = Boundary(grid);
// Simulation
Simulation sim = Simulation(grid, bc, approach);
sim.setOutputConsole(CONSOLE_OUTPUT_ON);
// 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);
}
@ -69,11 +72,13 @@ 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"){
TEST_CASE("Initialize environment") {
int rc = 12;
Grid grid(rc, rc);
Boundary boundary(grid);
@ -81,26 +86,23 @@ TEST_CASE("Initialize environment"){
CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH));
}
TEST_CASE("Simulation environment"){
TEST_CASE("Simulation environment") {
int rc = 12;
Grid grid(rc, rc);
Boundary boundary(grid);
Simulation sim(grid, boundary, FTCS_APPROACH);
SUBCASE("default paremeters"){
CHECK_EQ(sim.getIterations(), -1);
}
SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); }
SUBCASE("set iterations"){
SUBCASE("set iterations") {
CHECK_NOTHROW(sim.setIterations(2000));
CHECK_EQ(sim.getIterations(), 2000);
CHECK_THROWS(sim.setIterations(-300));
}
SUBCASE("set timestep"){
SUBCASE("set timestep") {
CHECK_NOTHROW(sim.setTimestep(0.1));
CHECK_EQ(sim.getTimestep(), 0.1);
CHECK_THROWS(sim.setTimestep(-0.3));
}
}