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

Merge branch 'ml/restructure-api' into 'main'

Browse Source 
Draft: Refactor API

See merge request naaice/tug!33
This commit is contained in:
Max Lübke 2024-11-29 05:31:42 +01:00
commit a9d58cc2a3
30 changed files with 658 additions and 1398 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitlab-ci.yml
View File

@ -14,7 +14,7 @@ build_release:
expire_in: 100s expire_in: 100s
script: script:
- mkdir build && cd build - mkdir build && cd build
- cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON .. - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON ..
- make -j$(nproc) - make -j$(nproc)
- cp ../test/FTCS_11_11_7000.csv test/ - cp ../test/FTCS_11_11_7000.csv test/

4
docs_sphinx/Simulation.rst → docs_sphinx/Diffusion.rst
View File

@ -1,4 +1,4 @@
Simulation Diffusion
========== ==========
.. doxygenenum:: tug::APPROACH .. doxygenenum:: tug::APPROACH
@ -7,4 +7,4 @@ Simulation
.. doxygenenum:: tug::CONSOLE_OUTPUT .. doxygenenum:: tug::CONSOLE_OUTPUT
.. doxygenenum:: tug::TIME_MEASURE .. doxygenenum:: tug::TIME_MEASURE
.. doxygenclass:: tug::Simulation .. doxygenclass:: tug::Diffusion

2
docs_sphinx/user.rst
View File

@ -5,4 +5,4 @@ User API
Grid Grid
Boundary Boundary
Simulation Diffusion

51
examples/BTCS_1D_proto_example.cpp
View File

@ -1,51 +0,0 @@
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace tug;
int main(int argc, char *argv[]) {
// **************
// **** GRID ****
// **************
// create a linear grid with 20 cells
int cells = 20;
Grid64 grid(cells);
MatrixXd concentrations = MatrixXd::Constant(1, 20, 0);
concentrations(0, 0) = 2000;
// TODO add option to set concentrations with a vector in 1D case
grid.setConcentrations(concentrations);
// ******************
// **** BOUNDARY ****
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation = Simulation(grid, bc); // grid,boundary
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
// set the number of iterations
simulation.setIterations(100);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****
// run the simulation
simulation.run();
}

8
examples/BTCS_2D_proto_example.cpp
View File

@ -1,5 +1,5 @@
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <tug/Simulation.hpp> #include <tug/Diffusion.hpp>
using namespace Eigen; using namespace Eigen;
using namespace tug; using namespace tug;
@ -14,7 +14,6 @@ int main(int argc, char *argv[]) {
// create a grid with a 20 x 20 field // create a grid with a 20 x 20 field
int row = 40; int row = 40;
int col = 50; int col = 50;
Grid64 grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
// grid.setDomain(20, 20); // grid.setDomain(20, 20);
@ -24,7 +23,7 @@ int main(int argc, char *argv[]) {
// #row,#col,value grid.setConcentrations(concentrations); // #row,#col,value grid.setConcentrations(concentrations);
MatrixXd concentrations = MatrixXd::Constant(row, col, 0); MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(10, 10) = 2000; concentrations(10, 10) = 2000;
grid.setConcentrations(concentrations); UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.: // (optional) set alphas of the grid, e.g.:
// MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value
@ -61,8 +60,7 @@ int main(int argc, char *argv[]) {
// ************************ // ************************
// set up a simulation environment // set up a simulation environment
Simulation simulation = Diffusion simulation(grid, bc); // grid,boundary
Simulation(grid, bc); // grid,boundary
// set the timestep of the simulation // set the timestep of the simulation
simulation.setTimestep(0.1); // timestep simulation.setTimestep(0.1); // timestep

17
examples/CMakeLists.txt
View File

@ -1,20 +1,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(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp)
add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
add_executable(profiling_openmp profiling_openmp.cpp)
target_link_libraries(FTCS_1D_proto_example tug)
target_link_libraries(FTCS_2D_proto_example tug)
target_link_libraries(BTCS_1D_proto_example tug)
target_link_libraries(BTCS_2D_proto_example tug)
target_link_libraries(CRNI_2D_proto_example tug)
target_link_libraries(reference-FTCS_2D_closed tug)
target_link_libraries(profiling_openmp tug)
add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
target_link_libraries(BTCS_2D_proto_example tug)
target_link_libraries(FTCS_2D_proto_closed_mdl tug) target_link_libraries(FTCS_2D_proto_closed_mdl tug)
target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(FTCS_2D_proto_example_mdl tug)

29
examples/CRNI_2D_proto_example.cpp
View File

@ -1,29 +0,0 @@
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace tug;
int main(int argc, char *argv[]) {
int row = 20;
int col = 20;
Grid64 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<double, tug::CRANK_NICOLSON_APPROACH>(grid, bc);
simulation.setTimestep(0.1);
simulation.setIterations(50);
simulation.setOutputCSV(CSV_OUTPUT_XTREME);
simulation.run();
}

51
examples/FTCS_1D_proto_example.cpp
View File

@ -1,51 +0,0 @@
#include "tug/Boundary.hpp"
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace tug;
int main(int argc, char *argv[]) {
// **************
// **** GRID ****
// **************
// create a linear grid with 20 cells
int cells = 20;
Grid64 grid(cells);
MatrixXd concentrations = MatrixXd::Constant(1, 20, 20);
grid.setConcentrations(concentrations);
// ******************
// **** BOUNDARY ****
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation =
Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
// (optional) set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
// (optional) set the number of iterations
simulation.setIterations(100);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_OFF);
// **** RUN SIMULATION ****
// run the simulation
simulation.run();
}

9
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -9,7 +9,7 @@
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <cstdlib> #include <cstdlib>
#include <iostream> #include <iostream>
#include <tug/Simulation.hpp> #include <tug/Diffusion.hpp>
using namespace Eigen; using namespace Eigen;
using namespace tug; using namespace tug;
@ -31,7 +31,6 @@ int main(int argc, char *argv[]) {
// create a grid with a 20 x 20 field // create a grid with a 20 x 20 field
int n2 = row / 2 - 1; int n2 = row / 2 - 1;
Grid64 grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
// grid.setDomain(20, 20); // grid.setDomain(20, 20);
@ -44,7 +43,7 @@ int main(int argc, char *argv[]) {
concentrations(n2, n2 + 1) = 1; concentrations(n2, n2 + 1) = 1;
concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2) = 1;
concentrations(n2 + 1, n2 + 1) = 1; concentrations(n2 + 1, n2 + 1) = 1;
grid.setConcentrations(concentrations); UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.: // (optional) set alphas of the grid, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value
@ -69,8 +68,8 @@ int main(int argc, char *argv[]) {
// ************************ // ************************
// set up a simulation environment // set up a simulation environment
Simulation simulation = Diffusion<double, FTCS_APPROACH> simulation(
Simulation<double, FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach grid, bc); // grid,boundary,simulation-approach
// set the timestep of the simulation // set the timestep of the simulation
simulation.setTimestep(10000); // timestep simulation.setTimestep(10000); // timestep

92
examples/FTCS_2D_proto_example.cpp
View File

@ -1,92 +0,0 @@
/**
* @file FTCS_2D_proto_example.cpp
* @author Hannes Signer, Philipp Ungrund
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS
* approach and constant boundary condition
*
*/
#include <Eigen/Eigen>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace tug;
// #include <easy/profiler.h>
// #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true);
int main(int argc, char *argv[]) {
// EASY_PROFILER_ENABLE;
// profiler::startListen();
// **************
// **** GRID ****
// **************
// profiler::startListen();
// create a grid with a 20 x 20 field
int row = 20;
int col = 20;
Grid64 grid(row, col);
// (optional) set the domain, e.g.:
// grid.setDomain(20, 20);
// (optional) set the concentrations, e.g.:
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); //
// #row,#col,value grid.setConcentrations(concentrations);
MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(0, 0) = 1999;
grid.setConcentrations(concentrations);
// (optional) set alphas of the grid, e.g.:
// MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value
// MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value
// grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
// (optional) set boundary condition values for one side, e.g.:
// VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
// bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values
// VectorXd bc_zero_values = VectorXd::Constant(20,0);
// bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
// bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
// VectorXd bc_front_values = VectorXd::Constant(20,2000);
// bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
// bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation =
Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
// set the number of iterations
simulation.setIterations(10000);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
// CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****
// run the simulation
// EASY_BLOCK("SIMULATION")
simulation.run();
// EASY_END_BLOCK;
// profiler::dumpBlocksToFile("test_profile.prof");
// profiler::stopListen();
}

9
examples/FTCS_2D_proto_example_mdl.cpp
View File

@ -7,7 +7,7 @@
*/ */
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <tug/Simulation.hpp> #include <tug/Diffusion.hpp>
using namespace Eigen; using namespace Eigen;
using namespace tug; using namespace tug;
@ -22,7 +22,6 @@ int main(int argc, char *argv[]) {
int row = 64; int row = 64;
int col = 64; int col = 64;
int n2 = row / 2 - 1; int n2 = row / 2 - 1;
Grid64 grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
// grid.setDomain(20, 20); // grid.setDomain(20, 20);
@ -35,7 +34,7 @@ int main(int argc, char *argv[]) {
concentrations(n2, n2 + 1) = 1; concentrations(n2, n2 + 1) = 1;
concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2) = 1;
concentrations(n2 + 1, n2 + 1) = 1; concentrations(n2 + 1, n2 + 1) = 1;
grid.setConcentrations(concentrations); UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.: // (optional) set alphas of the grid, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value
@ -60,8 +59,8 @@ int main(int argc, char *argv[]) {
// ************************ // ************************
// set up a simulation environment // set up a simulation environment
Simulation simulation = Diffusion<double, tug::FTCS_APPROACH> simulation(
Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach grid, bc); // grid,boundary,simulation-approach
// (optional) set the timestep of the simulation // (optional) set the timestep of the simulation
simulation.setTimestep(1000); // timestep simulation.setTimestep(1000); // timestep

70
examples/profiling_openmp.cpp
View File

@ -1,70 +0,0 @@
#include <Eigen/Eigen>
#include <chrono>
#include <fstream>
#include <iostream>
#include <string>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace std;
using namespace tug;
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;
Grid64 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);
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();
}
}

70
examples/profiling_speedup.cpp
View File

@ -1,70 +0,0 @@
#include <Eigen/Eigen>
#include <chrono>
#include <fstream>
#include <iostream>
#include <string>
#include <tug/Simulation.hpp>
using namespace Eigen;
using namespace std;
using namespace tug;
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;
Grid64 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);
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();
}
}

53
examples/reference-FTCS_2D_closed.cpp
View File

@ -1,53 +0,0 @@
#include "Eigen/Core"
#include <iostream>
#include <tug/Simulation.hpp>
using namespace std;
using namespace Eigen;
using namespace tug;
int main(int argc, char *argv[]) {
int row = 50;
int col = 50;
int domain_row = 10;
int domain_col = 10;
// Grid
Grid64 grid(row, col);
grid.setDomain(domain_row, domain_col);
MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(5, 5) = 1;
grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
for (int i = 0; i < 5; i++) {
for (int j = 0; j < 6; j++) {
alpha(i, j) = 0.01;
}
}
for (int i = 0; i < 5; i++) {
for (int j = 6; j < 11; j++) {
alpha(i, j) = 0.001;
}
}
for (int i = 5; i < 11; i++) {
for (int j = 6; j < 11; j++) {
alpha(i, j) = 0.1;
}
}
grid.setAlpha(alpha, alpha);
// Boundary
Boundary bc = Boundary(grid);
// Simulation
Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
sim.setTimestep(0.001);
sim.setIterations(10000);
sim.setOutputCSV(CSV_OUTPUT_OFF);
sim.setOutputConsole(CONSOLE_OUTPUT_OFF);
// RUN
sim.run();
}

134
include/tug/Boundary.hpp
View File

@ -7,7 +7,8 @@
#ifndef BOUNDARY_H_ #ifndef BOUNDARY_H_
#define BOUNDARY_H_ #define BOUNDARY_H_
#include "Grid.hpp" #include "UniformGrid.hpp"
#include "tug/Core/TugUtils.hpp"
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
@ -114,7 +115,7 @@ public:
* *
* @param length Length of the grid * @param length Length of the grid
*/ */
Boundary(std::uint32_t length) : Boundary(Grid<T>(length)){}; Boundary(std::uint32_t length) : Boundary(UnfiormGrid<T>(length)){};
/** /**
* @brief Creates a boundary object for a 2D grid * @brief Creates a boundary object for a 2D grid
@ -123,7 +124,7 @@ public:
* @param cols Number of columns of the grid * @param cols Number of columns of the grid
*/ */
Boundary(std::uint32_t rows, std::uint32_t cols) Boundary(std::uint32_t rows, std::uint32_t cols)
: Boundary(Grid<T>(rows, cols)){}; : Boundary(UnfiormGrid<T>(rows, cols)){};
/** /**
* @brief Creates a boundary object based on the passed grid object and * @brief Creates a boundary object based on the passed grid object and
@ -132,7 +133,7 @@ public:
* @param grid Grid object on the basis of which the simulation takes place * @param grid Grid object on the basis of which the simulation takes place
* and from which the dimensions (in 2D case) are taken. * and from which the dimensions (in 2D case) are taken.
*/ */
Boundary(const Grid<T> &grid) Boundary(const UnfiormGrid<T> &grid)
: dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
if (this->dim == 1) { if (this->dim == 1) {
this->boundaries = std::vector<std::vector<BoundaryElement<T>>>( this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
@ -161,13 +162,11 @@ public:
* @param side Side to be set to closed, e.g. BC_SIDE_LEFT. * @param side Side to be set to closed, e.g. BC_SIDE_LEFT.
*/ */
void setBoundarySideClosed(BC_SIDE side) { void setBoundarySideClosed(BC_SIDE side) {
if (this->dim == 1) { tug_assert((this->dim > 1) ||
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
throw std::invalid_argument( "For the "
"For the one-dimensional case, only the BC_SIDE_LEFT and " "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
"BC_SIDE_RIGHT borders exist."); "borders exist.");
}
}
const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
const int n = is_vertical ? this->rows : this->cols; const int n = is_vertical ? this->rows : this->cols;
@ -186,13 +185,11 @@ public:
* page. * page.
*/ */
void setBoundarySideConstant(BC_SIDE side, double value) { void setBoundarySideConstant(BC_SIDE side, double value) {
if (this->dim == 1) { tug_assert((this->dim > 1) ||
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
throw std::invalid_argument( "For the "
"For the one-dimensional case, only the BC_SIDE_LEFT and " "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
"BC_SIDE_RIGHT borders exist."); "borders exist.");
}
}
const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
const int n = is_vertical ? this->rows : this->cols; const int n = is_vertical ? this->rows : this->cols;
@ -212,10 +209,9 @@ public:
*/ */
void setBoundaryElemenClosed(BC_SIDE side, int index) { void setBoundaryElemenClosed(BC_SIDE side, int index) {
// tests whether the index really points to an element of the boundary side. // tests whether the index really points to an element of the boundary side.
if ((boundaries[side].size() < index) || index < 0) { tug_assert(boundaries[side].size() > index && index >= 0,
throw std::invalid_argument( "Index is selected either too large or too small.");
"Index is selected either too large or too small.");
}
this->boundaries[side][index].setType(BC_TYPE_CLOSED); this->boundaries[side][index].setType(BC_TYPE_CLOSED);
} }
@ -233,10 +229,8 @@ public:
*/ */
void setBoundaryElementConstant(BC_SIDE side, int index, double value) { void setBoundaryElementConstant(BC_SIDE side, int index, double value) {
// tests whether the index really points to an element of the boundary side. // tests whether the index really points to an element of the boundary side.
if ((boundaries[side].size() < index) || index < 0) { tug_assert(boundaries[side].size() > index && index >= 0,
throw std::invalid_argument( "Index is selected either too large or too small.");
"Index is selected either too large or too small.");
}
this->boundaries[side][index].setType(BC_TYPE_CONSTANT); this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
this->boundaries[side][index].setValue(value); this->boundaries[side][index].setValue(value);
} }
@ -251,13 +245,11 @@ public:
* BoundaryElement<T> objects. * BoundaryElement<T> objects.
*/ */
const std::vector<BoundaryElement<T>> &getBoundarySide(BC_SIDE side) const { const std::vector<BoundaryElement<T>> &getBoundarySide(BC_SIDE side) const {
if (this->dim == 1) { tug_assert((this->dim > 1) ||
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
throw std::invalid_argument( "For the "
"For the one-dimensional trap, only the BC_SIDE_LEFT and " "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
"BC_SIDE_RIGHT borders exist."); "borders exist.");
}
}
return this->boundaries[side]; return this->boundaries[side];
} }
@ -295,10 +287,8 @@ public:
* object. * object.
*/ */
BoundaryElement<T> getBoundaryElement(BC_SIDE side, int index) const { BoundaryElement<T> getBoundaryElement(BC_SIDE side, int index) const {
if ((boundaries[side].size() < index) || index < 0) { tug_assert(boundaries[side].size() > index && index >= 0,
throw std::invalid_argument( "Index is selected either too large or too small.");
"Index is selected either too large or too small.");
}
return this->boundaries[side][index]; return this->boundaries[side][index];
} }
@ -313,10 +303,8 @@ public:
* @return Boundary Type of the corresponding boundary condition. * @return Boundary Type of the corresponding boundary condition.
*/ */
BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const {
if ((boundaries[side].size() < index) || index < 0) { tug_assert(boundaries[side].size() > index && index >= 0,
throw std::invalid_argument( "Index is selected either too large or too small.");
"Index is selected either too large or too small.");
}
return this->boundaries[side][index].getType(); return this->boundaries[side][index].getType();
} }
@ -333,14 +321,12 @@ public:
* object. * object.
*/ */
T getBoundaryElementValue(BC_SIDE side, int index) const { T getBoundaryElementValue(BC_SIDE side, int index) const {
if ((boundaries[side].size() < index) || index < 0) { tug_assert(boundaries[side].size() > index && index >= 0,
throw std::invalid_argument( "Index is selected either too large or too small.");
"Index is selected either too large or too small."); tug_assert(
} boundaries[side][index].getType() == BC_TYPE_CONSTANT,
if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { "A value can only be output if it is a constant boundary condition.");
throw std::invalid_argument(
"A value can only be output if it is a constant boundary condition.");
}
return this->boundaries[side][index].getValue(); return this->boundaries[side][index].getValue();
} }
@ -382,13 +368,8 @@ public:
* @param value Value of the inner constant boundary condition * @param value Value of the inner constant boundary condition
*/ */
void setInnerBoundary(std::uint32_t index, T value) { void setInnerBoundary(std::uint32_t index, T value) {
if (this->dim != 1) { tug_assert(this->dim == 1, "This function is only available for 1D grids.");
throw std::invalid_argument( tug_assert(index < this->cols, "Index is out of bounds.");
"This function is only available for 1D grids.");
}
if (index >= this->cols) {
throw std::invalid_argument("Index is out of bounds.");
}
this->inner_boundary[std::make_pair(0, index)] = value; this->inner_boundary[std::make_pair(0, index)] = value;
} }
@ -401,13 +382,8 @@ public:
* @param value Value of the inner constant boundary condition * @param value Value of the inner constant boundary condition
*/ */
void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) { void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) {
if (this->dim != 2) { tug_assert(this->dim == 2, "This function is only available for 2D grids.");
throw std::invalid_argument( tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
"This function is only available for 2D grids.");
}
if (row >= this->rows || col >= this->cols) {
throw std::invalid_argument("Index is out of bounds.");
}
this->inner_boundary[std::make_pair(row, col)] = value; this->inner_boundary[std::make_pair(row, col)] = value;
} }
@ -420,13 +396,8 @@ public:
* set or not) and value of the inner constant boundary condition * set or not) and value of the inner constant boundary condition
*/ */
std::pair<bool, T> getInnerBoundary(std::uint32_t index) const { std::pair<bool, T> getInnerBoundary(std::uint32_t index) const {
if (this->dim != 1) { tug_assert(this->dim == 1, "This function is only available for 1D grids.");
throw std::invalid_argument( tug_assert(index < this->cols, "Index is out of bounds.");
"This function is only available for 1D grids.");
}
if (index >= this->cols) {
throw std::invalid_argument("Index is out of bounds.");
}
auto it = this->inner_boundary.find(std::make_pair(0, index)); auto it = this->inner_boundary.find(std::make_pair(0, index));
if (it == this->inner_boundary.end()) { if (it == this->inner_boundary.end()) {
@ -445,13 +416,8 @@ public:
*/ */
std::pair<bool, T> getInnerBoundary(std::uint32_t row, std::pair<bool, T> getInnerBoundary(std::uint32_t row,
std::uint32_t col) const { std::uint32_t col) const {
if (this->dim != 2) { tug_assert(this->dim == 2, "This function is only available for 2D grids.");
throw std::invalid_argument( tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
"This function is only available for 2D grids.");
}
if (row >= this->rows || col >= this->cols) {
throw std::invalid_argument("Index is out of bounds.");
}
auto it = this->inner_boundary.find(std::make_pair(row, col)); auto it = this->inner_boundary.find(std::make_pair(row, col));
if (it == this->inner_boundary.end()) { if (it == this->inner_boundary.end()) {
@ -471,9 +437,7 @@ public:
* condition * condition
*/ */
std::vector<std::pair<bool, T>> getInnerBoundaryRow(std::uint32_t row) const { std::vector<std::pair<bool, T>> getInnerBoundaryRow(std::uint32_t row) const {
if (row >= this->rows) { tug_assert(row < this->rows, "Index is out of bounds.");
throw std::invalid_argument("Index is out of bounds.");
}
if (inner_boundary.empty()) { if (inner_boundary.empty()) {
return std::vector<std::pair<bool, T>>(this->cols, return std::vector<std::pair<bool, T>>(this->cols,
@ -499,14 +463,8 @@ public:
* condition * condition
*/ */
std::vector<std::pair<bool, T>> getInnerBoundaryCol(std::uint32_t col) const { std::vector<std::pair<bool, T>> getInnerBoundaryCol(std::uint32_t col) const {
if (this->dim != 2) { tug_assert(this->dim == 2, "This function is only available for 2D grids.");
throw std::invalid_argument( tug_assert(col < this->cols, "Index is out of bounds.");
"This function is only available for 2D grids.");
}
if (col >= this->cols) {
throw std::invalid_argument("Index is out of bounds.");
}
if (inner_boundary.empty()) { if (inner_boundary.empty()) {
return std::vector<std::pair<bool, T>>(this->rows, return std::vector<std::pair<bool, T>>(this->rows,

135
include/tug/Core/BaseSimulation.hpp Normal file
View File

@ -0,0 +1,135 @@
#pragma once
#include <stdexcept>
namespace tug {
/**
* @brief Enum holding different options for .csv output.
*
*/
enum CSV_OUTPUT {
CSV_OUTPUT_OFF, /*!< do not produce csv output */
CSV_OUTPUT_ON, /*!< produce csv output with last concentration matrix */
CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */
CSV_OUTPUT_XTREME /*!< csv output like VERBOSE but additional boundary
conditions at beginning */
};
/**
* @brief Enum holding different options for console output.
*
*/
enum CONSOLE_OUTPUT {
CONSOLE_OUTPUT_OFF, /*!< do not print any output to console */
CONSOLE_OUTPUT_ON, /*!< print before and after concentrations to console */
CONSOLE_OUTPUT_VERBOSE /*!< print all concentration matrices to console */
};
/**
* @brief Enum holding different options for time measurement.
*
*/
enum TIME_MEASURE {
TIME_MEASURE_OFF, /*!< do not print any time measures */
TIME_MEASURE_ON /*!< print time measure after last iteration */
};
class BaseSimulation {
protected:
CSV_OUTPUT csv_output{CSV_OUTPUT_OFF};
CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF};
TIME_MEASURE time_measure{TIME_MEASURE_OFF};
int iterations{1};
public:
/**
* @brief Set the option to output the results to a CSV file. Off by default.
*
*
* @param csv_output Valid output option. The following options can be set
* here:
* - CSV_OUTPUT_OFF: do not produce csv output
* - CSV_OUTPUT_ON: produce csv output with last
* concentration matrix
* - CSV_OUTPUT_VERBOSE: produce csv output with all
* concentration matrices
* - CSV_OUTPUT_XTREME: produce csv output with all
* concentration matrices and simulation environment
*/
void setOutputCSV(CSV_OUTPUT csv_output) {
if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
throw std::invalid_argument("Invalid CSV output option given!");
}
this->csv_output = csv_output;
}
/**
* @brief Set the options for outputting information to the console. Off by
* default.
*
* @param console_output Valid output option. The following options can be set
* here:
* - CONSOLE_OUTPUT_OFF: do not print any output to
* console
* - CONSOLE_OUTPUT_ON: print before and after
* concentrations to console
* - CONSOLE_OUTPUT_VERBOSE: print all concentration
* matrices to console
*/
void setOutputConsole(CONSOLE_OUTPUT console_output) {
if (console_output < CONSOLE_OUTPUT_OFF &&
console_output > CONSOLE_OUTPUT_VERBOSE) {
throw std::invalid_argument("Invalid console output option given!");
}
this->console_output = console_output;
}
/**
* @brief Set the Time Measure option. Off by default.
*
* @param time_measure The following options are allowed:
* - TIME_MEASURE_OFF: Time of simulation is not printed
* to console
* - TIME_MEASURE_ON: Time of simulation run is printed to
* console
*/
void setTimeMeasure(TIME_MEASURE time_measure) {
if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
throw std::invalid_argument("Invalid time measure option given!");
}
this->time_measure = time_measure;
}
/**
* @brief Set the desired iterations to be calculated. A value greater
* than zero must be specified here. Setting iterations is required.
*
* @param iterations Number of iterations to be simulated.
*/
void setIterations(int iterations) {
if (iterations <= 0) {
throw std::invalid_argument(
"Number of iterations must be greater than zero.");
}
this->iterations = iterations;
}
/**
* @brief Return the currently set iterations to be calculated.
*
* @return int Number of iterations.
*/
int getIterations() const { return this->iterations; }
/**
* @brief Method starts the simulation process with the previously set
* parameters.
*/
virtual void run() = 0;
};
} // namespace tug

45
include/tug/Core/Matrix.hpp
View File

@ -1,6 +1,8 @@
#pragma once #pragma once
#include <Eigen/Core> #include <Eigen/Core>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
namespace tug { namespace tug {
/** /**
@ -13,9 +15,48 @@ namespace tug {
* *
* @tparam T The type of the matrix elements. * @tparam T The type of the matrix elements.
*/ */
// template <typename T>
// using RowMajMat =
// Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
template <typename T> template <typename T>
using RowMajMat = class RowMajMat
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; : public Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> {
protected:
std::uint8_t dim;
public:
RowMajMat(Eigen::Index rows, Eigen::Index cols, T initial_value) : dim(2) {
*this = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
Eigen::RowMajor>::Constant(rows, cols, initial_value);
};
RowMajMat(Eigen::Index n_cells, T initial_value) : dim(1) {
*this = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
Eigen::RowMajor>::Constant(1, n_cells, initial_value);
};
/**
* @brief Gets the number of rows of the grid.
*
* @return Number of rows.
*/
Eigen::Index getRow() const { return this->rows(); }
/**
* @brief Gets the number of columns of the grid.
*
* @return Number of columns.
*/
Eigen::Index getCol() const { return this->cols(); }
/**
* @brief Gets the dimensions of the grid.
*
* @return Dimensions, either 1 or 2.
*/
int getDim() const { return this->dim; }
};
template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>; template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>;
} // namespace tug } // namespace tug

20
include/tug/Core/BTCS.hpp → include/tug/Core/Numeric/BTCS.hpp
View File

@ -10,8 +10,8 @@
#ifndef BTCS_H_ #ifndef BTCS_H_
#define BTCS_H_ #define BTCS_H_
#include "Matrix.hpp" #include "../Matrix.hpp"
#include "TugUtils.hpp" #include "../TugUtils.hpp"
#include <cstddef> #include <cstddef>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
@ -351,18 +351,18 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
// BTCS solution for 1D grid // BTCS solution for 1D grid
template <class T> template <class T>
static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep, static void BTCS_1D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A, Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b)) { Eigen::VectorX<T> &b)) {
int length = grid.getLength(); int length = grid.getCol();
T sx = timestep / (grid.getDelta() * grid.getDelta()); T sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol());
Eigen::VectorX<T> concentrations_t1(length); Eigen::VectorX<T> concentrations_t1(length);
Eigen::SparseMatrix<T> A; Eigen::SparseMatrix<T> A;
Eigen::VectorX<T> b(length); Eigen::VectorX<T> b(length);
const auto &alpha = grid.getAlpha(); const auto &alpha = grid.getAlphaX();
const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
@ -396,7 +396,7 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// BTCS solution for 2D grid // BTCS solution for 2D grid
template <class T> template <class T>
static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep, static void BTCS_2D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A, Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b), Eigen::VectorX<T> &b),
int numThreads) { int numThreads) {
@ -449,7 +449,8 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// entry point for EigenLU solver; differentiate between 1D and 2D grid // entry point for EigenLU solver; differentiate between 1D and 2D grid
template <class T> template <class T>
void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) { void BTCS_LU(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {
@ -462,7 +463,8 @@ void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
template <class T> template <class T>
void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) { void BTCS_Thomas(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm); BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {

90
include/tug/Core/FTCS.hpp → include/tug/Core/Numeric/FTCS.hpp
View File

@ -8,10 +8,10 @@
#ifndef FTCS_H_ #ifndef FTCS_H_
#define FTCS_H_ #define FTCS_H_
#include "TugUtils.hpp" #include "../TugUtils.hpp"
#include "tug/Core/Matrix.hpp"
#include <cstddef> #include <cstring>
#include <iostream>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#ifdef _OPENMP #ifdef _OPENMP
@ -24,7 +24,7 @@ namespace tug {
// calculates horizontal change on one cell independent of boundary type // calculates horizontal change on one cell independent of boundary type
template <class T> template <class T>
static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) { static inline T calcHorizontalChange(UnfiormGrid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
@ -41,7 +41,7 @@ static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) {
// calculates vertical change on one cell independent of boundary type // calculates vertical change on one cell independent of boundary type
template <class T> template <class T>
static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) { static inline T calcVerticalChange(UnfiormGrid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) * grid.getAlphaY()(row, col)) *
@ -58,7 +58,7 @@ static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
// calculates horizontal change on one cell with a constant left boundary // calculates horizontal change on one cell with a constant left boundary
template <class T> template <class T>
static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid, static inline T calcHorizontalChangeLeftBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc, Boundary<T> &bc,
int &row, int &col) { int &row, int &col) {
@ -75,8 +75,8 @@ static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
// calculates horizontal change on one cell with a closed left boundary // calculates horizontal change on one cell with a closed left boundary
template <class T> template <class T>
static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row, static inline T calcHorizontalChangeLeftBoundaryClosed(UnfiormGrid<T> &grid,
int &col) { int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
@ -86,8 +86,9 @@ static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the left edge of grid // checks boundary condition type for a cell on the left edge of grid
template <class T> template <class T>
static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc, static inline T calcHorizontalChangeLeftBoundary(UnfiormGrid<T> &grid,
int &row, int &col) { Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) { } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) {
@ -99,7 +100,7 @@ static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
// calculates horizontal change on one cell with a constant right boundary // calculates horizontal change on one cell with a constant right boundary
template <class T> template <class T>
static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid, static inline T calcHorizontalChangeRightBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc, Boundary<T> &bc,
int &row, int &col) { int &row, int &col) {
@ -116,8 +117,8 @@ static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
// calculates horizontal change on one cell with a closed right boundary // calculates horizontal change on one cell with a closed right boundary
template <class T> template <class T>
static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row, static inline T calcHorizontalChangeRightBoundaryClosed(UnfiormGrid<T> &grid,
int &col) { int &row, int &col) {
return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) * grid.getAlphaX()(row, col)) *
@ -127,7 +128,7 @@ static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the right edge of grid // checks boundary condition type for a cell on the right edge of grid
template <class T> template <class T>
static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid, static inline T calcHorizontalChangeRightBoundary(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row, Boundary<T> &bc, int &row,
int &col) { int &col) {
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) {
@ -141,7 +142,7 @@ static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
// calculates vertical change on one cell with a constant top boundary // calculates vertical change on one cell with a constant top boundary
template <class T> template <class T>
static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid, static inline T calcVerticalChangeTopBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row, Boundary<T> &bc, int &row,
int &col) { int &col) {
@ -158,8 +159,8 @@ static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
// calculates vertical change on one cell with a closed top boundary // calculates vertical change on one cell with a closed top boundary
template <class T> template <class T>
static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row, static inline T calcVerticalChangeTopBoundaryClosed(UnfiormGrid<T> &grid,
int &col) { int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) * grid.getAlphaY()(row, col)) *
@ -169,8 +170,9 @@ static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the top edge of grid // checks boundary condition type for a cell on the top edge of grid
template <class T> template <class T>
static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc, static inline T calcVerticalChangeTopBoundary(UnfiormGrid<T> &grid,
int &row, int &col) { Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
@ -182,7 +184,7 @@ static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc,
// calculates vertical change on one cell with a constant bottom boundary // calculates vertical change on one cell with a constant bottom boundary
template <class T> template <class T>
static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid, static inline T calcVerticalChangeBottomBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc, Boundary<T> &bc,
int &row, int &col) { int &row, int &col) {
@ -199,8 +201,8 @@ static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
// calculates vertical change on one cell with a closed bottom boundary // calculates vertical change on one cell with a closed bottom boundary
template <class T> template <class T>
static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row, static inline T calcVerticalChangeBottomBoundaryClosed(UnfiormGrid<T> &grid,
int &col) { int &row, int &col) {
return -(calcAlphaIntercell(grid.getAlphaY()(row, col), return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) * grid.getAlphaY()(row - 1, col)) *
@ -210,8 +212,9 @@ static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the bottom edge of grid // checks boundary condition type for a cell on the bottom edge of grid
template <class T> template <class T>
static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc, static inline T calcVerticalChangeBottomBoundary(UnfiormGrid<T> &grid,
int &row, int &col) { Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) { } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
@ -223,10 +226,11 @@ static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc,
// FTCS solution for 1D grid // FTCS solution for 1D grid
template <class T> template <class T>
static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) { static void FTCS_1D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep) {
int colMax = grid.getCol(); int colMax = grid.getCol();
T deltaCol = grid.getDeltaCol(); T deltaCol = grid.getDeltaCol();
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0); RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
@ -236,7 +240,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
// inner cells // inner cells
// independent of boundary condition type // independent of boundary condition type
for (int col = 1; col < colMax - 1; col++) { for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col) + concentrations_t1(row, col) = concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col)); (calcHorizontalChange(grid, row, col));
} }
@ -244,30 +248,33 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
// left boundary; hold column constant at index 0 // left boundary; hold column constant at index 0
int col = 0; int col = 0;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)); (calcHorizontalChangeLeftBoundary(grid, bc, row, col));
// right boundary; hold column constant at max index // right boundary; hold column constant at max index
col = colMax - 1; col = colMax - 1;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)); (calcHorizontalChangeRightBoundary(grid, bc, row, col));
// overwrite obsolete concentrations // overwrite obsolete concentrations
grid.setConcentrations(concentrations_t1); std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
colMax * sizeof(T));
} }
// FTCS solution for 2D grid // FTCS solution for 2D grid
template <class T> template <class T>
static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep, static void FTCS_2D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) { int numThreads) {
int rowMax = grid.getRow(); int rowMax = grid.getRow();
int colMax = grid.getCol(); int colMax = grid.getCol();
T deltaRow = grid.getDeltaRow(); T deltaRow = grid.getDeltaRow();
T deltaCol = grid.getDeltaCol(); T deltaCol = grid.getDeltaCol();
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0); RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
@ -277,7 +284,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
#pragma omp parallel for num_threads(numThreads) #pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) { for (int row = 1; row < rowMax - 1; row++) {
for (int col = 1; col < colMax - 1; col++) { for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col) + concentrations_t1(row, col) = concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
(calcVerticalChange(grid, row, col)) + (calcVerticalChange(grid, row, col)) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
@ -292,7 +299,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
#pragma omp parallel for num_threads(numThreads) #pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) { for (int row = 1; row < rowMax - 1; row++) {
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
@ -304,7 +311,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
#pragma omp parallel for num_threads(numThreads) #pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) { for (int row = 1; row < rowMax - 1; row++) {
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) + (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
@ -316,7 +323,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
#pragma omp parallel for num_threads(numThreads) #pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) { for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col)) + (calcVerticalChangeTopBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
@ -329,7 +336,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
#pragma omp parallel for num_threads(numThreads) #pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) { for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col)) + (calcVerticalChangeBottomBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
@ -341,7 +348,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
row = 0; row = 0;
col = 0; col = 0;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
@ -352,7 +359,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
row = 0; row = 0;
col = colMax - 1; col = colMax - 1;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) + (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
@ -363,7 +370,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
row = rowMax - 1; row = rowMax - 1;
col = 0; col = 0;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
@ -374,20 +381,21 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
row = rowMax - 1; row = rowMax - 1;
col = colMax - 1; col = colMax - 1;
concentrations_t1(row, col) = concentrations_t1(row, col) =
grid.getConcentrations()(row, col) + concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) * timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) + (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col)); (calcVerticalChangeBottomBoundary(grid, bc, row, col));
// overwrite obsolete concentrations // overwrite obsolete concentrations
grid.setConcentrations(concentrations_t1); std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
rowMax * colMax * sizeof(T));
// } // }
} }
// entry point; differentiate between 1D and 2D grid // entry point; differentiate between 1D and 2D grid
template <class T> template <class T>
void FTCS(Grid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) { void FTCS(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
FTCS_1D(grid, bc, timestep); FTCS_1D(grid, bc, timestep);
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {

10
include/tug/Core/TugUtils.hpp
View File

@ -1,9 +1,6 @@
#ifndef TUGUTILS_H_ #pragma once
#define TUGUTILS_H_
#include <chrono> #include <cassert>
#include <stdexcept>
#include <string>
#define throw_invalid_argument(msg) \ #define throw_invalid_argument(msg) \
throw std::invalid_argument(std::string(__FILE__) + ":" + \ throw std::invalid_argument(std::string(__FILE__) + ":" + \
@ -24,6 +21,8 @@
duration.count(); \ duration.count(); \
}) })
#define tug_assert(expr, msg) assert((expr) && msg)
// calculates arithmetic or harmonic mean of alpha between two cells // calculates arithmetic or harmonic mean of alpha between two cells
template <typename T> template <typename T>
constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) { constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) {
@ -38,4 +37,3 @@ constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) {
return 0.5 * (alpha1 + alpha2); return 0.5 * (alpha1 + alpha2);
} }
} }
#endif // TUGUTILS_H_

175
include/tug/Simulation.hpp → include/tug/Diffusion.hpp
View File

@ -1,16 +1,15 @@
/** /**
* @file Simulation.hpp * @file Diffusion.hpp
* @brief API of Simulation class, that holds all information regarding a * @brief API of Diffusion class, that holds all information regarding a
* specific simulation run like its timestep, number of iterations and output * specific simulation run like its timestep, number of iterations and output
* options. Simulation object also holds a predefined Grid and Boundary object. * options. Diffusion object also holds a predefined Grid and Boundary object.
* *
*/ */
#ifndef SIMULATION_H_ #pragma once
#define SIMULATION_H_
#include "Boundary.hpp" #include "Boundary.hpp"
#include "Grid.hpp" #include "UniformGrid.hpp"
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <fstream> #include <fstream>
@ -19,9 +18,10 @@
#include <string> #include <string>
#include <vector> #include <vector>
#include "Core/BTCS.hpp" #include "Core/Numeric/BTCS.hpp"
#include "Core/FTCS.hpp" #include "Core/Numeric/FTCS.hpp"
#include "Core/TugUtils.hpp" #include "Core/TugUtils.hpp"
#include "tug/Core/BaseSimulation.hpp"
#ifdef _OPENMP #ifdef _OPENMP
#include <omp.h> #include <omp.h>
@ -51,37 +51,6 @@ enum SOLVER {
tridiagonal matrices */ tridiagonal matrices */
}; };
/**
* @brief Enum holding different options for .csv output.
*
*/
enum CSV_OUTPUT {
CSV_OUTPUT_OFF, /*!< do not produce csv output */
CSV_OUTPUT_ON, /*!< produce csv output with last concentration matrix */
CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */
CSV_OUTPUT_XTREME /*!< csv output like VERBOSE but additional boundary
conditions at beginning */
};
/**
* @brief Enum holding different options for console output.
*
*/
enum CONSOLE_OUTPUT {
CONSOLE_OUTPUT_OFF, /*!< do not print any output to console */
CONSOLE_OUTPUT_ON, /*!< print before and after concentrations to console */
CONSOLE_OUTPUT_VERBOSE /*!< print all concentration matrices to console */
};
/**
* @brief Enum holding different options for time measurement.
*
*/
enum TIME_MEASURE {
TIME_MEASURE_OFF, /*!< do not print any time measures */
TIME_MEASURE_ON /*!< print time measure after last iteration */
};
/** /**
* @brief The class forms the interface for performing the diffusion simulations * @brief The class forms the interface for performing the diffusion simulations
* and contains all the methods for controlling the desired parameters, such as * and contains all the methods for controlling the desired parameters, such as
@ -94,7 +63,17 @@ enum TIME_MEASURE {
*/ */
template <class T, APPROACH approach = BTCS_APPROACH, template <class T, APPROACH approach = BTCS_APPROACH,
SOLVER solver = THOMAS_ALGORITHM_SOLVER> SOLVER solver = THOMAS_ALGORITHM_SOLVER>
class Simulation { class Diffusion : public BaseSimulation {
private:
T timestep{-1};
int innerIterations{1};
int numThreads{omp_get_num_procs()};
UnfiormGrid<T> &grid;
Boundary<T> &bc;
const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
public: public:
/** /**
* @brief Set up a simulation environment. The timestep and number of * @brief Set up a simulation environment. The timestep and number of
@ -108,68 +87,7 @@ public:
* @param bc Valid boundary condition object * @param bc Valid boundary condition object
* @param approach Approach to solving the problem. Either FTCS or BTCS. * @param approach Approach to solving the problem. Either FTCS or BTCS.
*/ */
Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){}; Diffusion(UnfiormGrid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
/**
* @brief Set the option to output the results to a CSV file. Off by default.
*
*
* @param csv_output Valid output option. The following options can be set
* here:
* - CSV_OUTPUT_OFF: do not produce csv output
* - CSV_OUTPUT_ON: produce csv output with last
* concentration matrix
* - CSV_OUTPUT_VERBOSE: produce csv output with all
* concentration matrices
* - CSV_OUTPUT_XTREME: produce csv output with all
* concentration matrices and simulation environment
*/
void setOutputCSV(CSV_OUTPUT csv_output) {
if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
throw std::invalid_argument("Invalid CSV output option given!");
}
this->csv_output = csv_output;
}
/**
* @brief Set the options for outputting information to the console. Off by
* default.
*
* @param console_output Valid output option. The following options can be set
* here:
* - CONSOLE_OUTPUT_OFF: do not print any output to
* console
* - CONSOLE_OUTPUT_ON: print before and after
* concentrations to console
* - CONSOLE_OUTPUT_VERBOSE: print all concentration
* matrices to console
*/
void setOutputConsole(CONSOLE_OUTPUT console_output) {
if (console_output < CONSOLE_OUTPUT_OFF &&
console_output > CONSOLE_OUTPUT_VERBOSE) {
throw std::invalid_argument("Invalid console output option given!");
}
this->console_output = console_output;
}
/**
* @brief Set the Time Measure option. Off by default.
*
* @param time_measure The following options are allowed:
* - TIME_MEASURE_OFF: Time of simulation is not printed
* to console
* - TIME_MEASURE_ON: Time of simulation run is printed to
* console
*/
void setTimeMeasure(TIME_MEASURE time_measure) {
if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
throw std::invalid_argument("Invalid time measure option given!");
}
this->time_measure = time_measure;
}
/** /**
* @brief Setting the time step for each iteration step. Time step must be * @brief Setting the time step for each iteration step. Time step must be
@ -178,17 +96,15 @@ public:
* @param timestep Valid timestep greater than zero. * @param timestep Valid timestep greater than zero.
*/ */
void setTimestep(T timestep) { void setTimestep(T timestep) {
if (timestep <= 0) { tug_assert(timestep > 0, "Timestep has to be greater than zero.");
throw_invalid_argument("Timestep has to be greater than zero.");
}
if constexpr (approach == FTCS_APPROACH || if constexpr (approach == FTCS_APPROACH ||
approach == CRANK_NICOLSON_APPROACH) { approach == CRANK_NICOLSON_APPROACH) {
T cfl; T cfl;
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
const T deltaSquare = grid.getDelta(); const T deltaSquare = grid.getDeltaCol();
const T maxAlpha = grid.getAlpha().maxCoeff(); const T maxAlpha = grid.getAlphaX().maxCoeff();
// Courant-Friedrichs-Lewy condition // Courant-Friedrichs-Lewy condition
cfl = deltaSquare / (4 * maxAlpha); cfl = deltaSquare / (4 * maxAlpha);
@ -244,20 +160,6 @@ public:
*/ */
T getTimestep() const { return this->timestep; } T getTimestep() const { return this->timestep; }
/**
* @brief Set the desired iterations to be calculated. A value greater
* than zero must be specified here. Setting iterations is required.
*
* @param iterations Number of iterations to be simulated.
*/
void setIterations(int iterations) {
if (iterations <= 0) {
throw std::invalid_argument(
"Number of iterations must be greater than zero.");
}
this->iterations = iterations;
}
/** /**
* @brief Set the number of desired openMP Threads. * @brief Set the number of desired openMP Threads.
* *
@ -277,13 +179,6 @@ public:
} }
} }
/**
* @brief Return the currently set iterations to be calculated.
*
* @return int Number of iterations.
*/
int getIterations() const { return this->iterations; }
/** /**
* @brief Outputs the current concentrations of the grid on the console. * @brief Outputs the current concentrations of the grid on the console.
* *
@ -375,12 +270,8 @@ public:
* parameters. * parameters.
*/ */
void run() { void run() {
if (this->timestep == -1) { tug_assert(this->timestep > 0, "Timestep is not set!");
throw_invalid_argument("Timestep is not set!"); tug_assert(this->iterations > 0, "Number of iterations are not set!");
}
if (this->iterations == -1) {
throw_invalid_argument("Number of iterations are not set!");
}
std::string filename; std::string filename;
if (this->console_output > CONSOLE_OUTPUT_OFF) { if (this->console_output > CONSOLE_OUTPUT_OFF) {
@ -393,7 +284,6 @@ public:
auto begin = std::chrono::high_resolution_clock::now(); auto begin = std::chrono::high_resolution_clock::now();
if constexpr (approach == FTCS_APPROACH) { // FTCS case if constexpr (approach == FTCS_APPROACH) { // FTCS case
for (int i = 0; i < iterations * innerIterations; i++) { for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole(); printConcentrationsConsole();
@ -485,20 +375,5 @@ public:
<< milliseconds.count() << "ms" << std::endl; << milliseconds.count() << "ms" << std::endl;
} }
} }
private:
T 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<T> &grid;
Boundary<T> &bc;
const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
}; };
} // namespace tug } // namespace tug
#endif // SIMULATION_H_

392
include/tug/Grid.hpp
View File

@ -1,392 +0,0 @@
#ifndef GRID_H_
#define GRID_H_
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "Core/Matrix.hpp"
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <stdexcept>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class Grid {
public:
/**
* @brief Constructs a new 1D-Grid object of a given length, which holds a
* matrix with concentrations and a respective matrix of alpha coefficients.
* The domain length is per default the same as the length. The
* concentrations are all 20 by default and the alpha coefficients are 1.
*
* @param length Length of the 1D-Grid. Must be greater than 3.
*/
Grid(int length) : col(length), domainCol(length) {
if (length <= 3) {
throw std::invalid_argument(
"Given grid length too small. Must be greater than 3.");
}
this->dim = 1;
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
this->concentrations = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
}
/**
* @brief Constructs a new 2D-Grid object of given dimensions, which holds a
* matrix with concentrations and the respective matrices of alpha coefficient
* for each direction. The domain in x- and y-direction is per default equal
* to the col length and row length, respectively. The concentrations are all
* 20 by default across the entire grid and the alpha coefficients 1 in both
* directions.
*
* @param row Length of the 2D-Grid in y-direction. Must be greater than 3.
* @param col Length of the 2D-Grid in x-direction. Must be greater than 3.
*/
Grid(int _row, int _col)
: row(_row), col(_col), domainRow(_row), domainCol(_col) {
if (row <= 1 || col <= 1) {
throw std::invalid_argument(
"At least one dimension is 1. Use 1D grid for better results.");
}
this->dim = 2;
this->deltaRow =
static_cast<T>(this->domainRow) / static_cast<T>(this->row); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
this->concentrations = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaY = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
}
/**
* @brief Sets the concentrations matrix for a 1D or 2D-Grid.
*
* @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
* Matrix must have correct dimensions as defined in row and col. (Or length,
* in 1D case).
*/
void setConcentrations(const RowMajMat<T> &concentrations) {
if (concentrations.rows() != this->row ||
concentrations.cols() != this->col) {
throw std::invalid_argument(
"Given matrix of concentrations mismatch with Grid dimensions!");
}
this->concentrations = concentrations;
}
/**
* @brief Sets the concentrations matrix for a 1D or 2D-Grid.
*
* @param concentrations A pointer to an array holding the concentrations.
* Array must have correct dimensions as defined in row and col. (Or length,
* in 1D case). There is no check for correct dimensions, so be careful!
*/
void setConcentrations(T *concentrations) {
tug::RowMajMatMap<T> map(concentrations, this->row, this->col);
this->concentrations = map;
}
/**
* @brief Gets the concentrations matrix for a Grid.
*
* @return An Eigen3 matrix holding the concentrations and having
* the same dimensions as the grid.
*/
auto &getConcentrations() { return this->concentrations; }
/**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
* coefficients. Matrix columns must have same size as length of grid.
*/
void setAlpha(const RowMajMat<T> &alpha) {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probably "
"use 2D setter function!");
}
if (alpha.rows() != 1 || alpha.cols() != this->col) {
throw std::invalid_argument(
"Given matrix of alpha coefficients mismatch with Grid dimensions!");
}
this->alphaX = alpha;
}
/**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @param alpha A pointer to an array holding the alpha coefficients. Array
* must have correct dimensions as defined in length. There is no check for
* correct dimensions, so be careful!
*/
void setAlpha(T *alpha) {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probably "
"use 2D setter function!");
}
RowMajMatMap<T> map(alpha, 1, this->col);
this->alphaX = map;
}
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* dimensional.
*
* @param alphaX An Eigen3 MatrixX<T> holding the alpha coefficients in
* x-direction. Matrix must be of same size as the grid.
* @param alphaY An Eigen3 MatrixX<T> holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid.
*/
void setAlpha(const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY) {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, you should probably "
"use 1D setter function!");
}
if (alphaX.rows() != this->row || alphaX.cols() != this->col) {
throw std::invalid_argument(
"Given matrix of alpha coefficients in x-direction "
"mismatch with GRid dimensions!");
}
if (alphaY.rows() != this->row || alphaY.cols() != this->col) {
throw std::invalid_argument(
"Given matrix of alpha coefficients in y-direction "
"mismatch with GRid dimensions!");
}
this->alphaX = alphaX;
this->alphaY = alphaY;
}
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* dimensional.
*
* @param alphaX A pointer to an array holding the alpha coefficients in
* x-direction. Array must have correct dimensions as defined in row and col.
* There is no check for correct dimensions, so be careful!
* @param alphaY A pointer to an array holding the alpha coefficients in
* y-direction. Array must have correct dimensions as defined in row and col.
* There is no check for correct dimensions, so be careful!
*/
void setAlpha(T *alphaX, T *alphaY) {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, you should probably "
"use 1D setter function!");
}
RowMajMatMap<T> mapX(alphaX, this->row, this->col);
RowMajMatMap<T> mapY(alphaY, this->row, this->col);
this->alphaX = mapX;
this->alphaY = mapY;
}
/**
* @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @return A matrix with 1 row holding the alpha coefficients.
*/
const auto &getAlpha() const {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probably "
"use either getAlphaX() or getAlphaY()!");
}
return this->alphaX;
}
/**
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in x-direction.
*/
const auto &getAlphaX() const {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!");
}
return this->alphaX;
}
/**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in y-direction.
*/
const auto &getAlphaY() const {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!");
}
return this->alphaY;
}
/**
* @brief Gets the dimensions of the grid.
*
* @return Dimensions, either 1 or 2.
*/
int getDim() const { return this->dim; }
/**
* @brief Gets length of 1D grid. Must be one dimensional grid.
*
* @return Length of 1D grid.
*/
int getLength() const {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probably "
"use getRow() or getCol()!");
}
return col;
}
/**
* @brief Gets the number of rows of the grid.
*
* @return Number of rows.
*/
int getRow() const { return this->row; }
/**
* @brief Gets the number of columns of the grid.
*
* @return Number of columns.
*/
int getCol() const { return this->col; }
/**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
*
* @param domainLength A double value of the domain length. Must be positive.
*/
void setDomain(double domainLength) {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probaly "
"use the 2D domain setter!");
}
if (domainLength <= 0) {
throw std::invalid_argument("Given domain length is not positive!");
}
this->domainCol = domainLength;
this->deltaCol = double(this->domainCol) / double(this->col);
}
/**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
*
* @param domainRow A double value of the domain size in y-direction. Must
* be positive.
* @param domainCol A double value of the domain size in x-direction. Must
* be positive.
*/
void setDomain(double domainRow, double domainCol) {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, you should probably "
"use the 1D domain setter!");
}
if (domainRow <= 0 || domainCol <= 0) {
throw std::invalid_argument("Given domain size is not positive!");
}
this->domainRow = domainRow;
this->domainCol = domainCol;
this->deltaRow = double(this->domainRow) / double(this->row);
this->deltaCol = double(this->domainCol) / double(this->col);
}
/**
* @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
*
* @return Delta value.
*/
T getDelta() const {
if (dim != 1) {
throw std::invalid_argument(
"Grid is not one dimensional, you should probably "
"use the 2D delta getters");
}
return this->deltaCol;
}
/**
* @brief Gets the delta value in x-direction.
*
* @return Delta value in x-direction.
*/
T getDeltaCol() const { return this->deltaCol; }
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return Delta value in y-direction.
*/
T getDeltaRow() const {
if (dim != 2) {
throw std::invalid_argument(
"Grid is not two dimensional, meaning there is no "
"delta in y-direction!");
}
return this->deltaRow;
}
private:
int col; // number of grid columns
int row{1}; // number of grid rows
int dim; // 1D or 2D
T domainCol; // number of domain columns
T domainRow{0}; // number of domain rows
T deltaCol; // delta in x-direction (between columns)
T deltaRow{0}; // delta in y-direction (between rows)
RowMajMat<T> concentrations; // Matrix holding grid concentrations
RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
RowMajMat<T> alphaY; // Matrix holding alpha coefficients in y-direction
static constexpr T MAT_INIT_VAL = 0;
};
using Grid64 = Grid<double>;
using Grid32 = Grid<float>;
} // namespace tug
#endif // GRID_H_

30
include/tug/UniformAlpha.hpp Normal file
View File

@ -0,0 +1,30 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "tug/UniformGrid.hpp"
#include <Eigen/Core>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class UniformAlpha : public RowMajMat<T> {
public:
UniformAlpha(const UniformGrid<T> &grid, T init_value = 0)
: RowMajMat<T>::Constant(grid.rows(), grid.cols(), init_value) {}
UniformAlpha(const UniformGrid<T> &grid, T *alpha)
: tug::RowMajMatMap<T>(alpha, grid.rows(), grid.cols()) {}
};
} // namespace tug

117
include/tug/UniformGrid.hpp Normal file
View File

@ -0,0 +1,117 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "Core/Matrix.hpp"
#include "tug/Core/TugUtils.hpp"
#include <Eigen/Core>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class UniformGrid : public RowMajMat<T> {
public:
UniformGrid(Eigen::Index n_cells, T spatial_size, T initial_value)
: RowMajMat<T>(n_cells, initial_value), domainCol(spatial_size) {
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(n_cells);
}
UniformGrid(Eigen::Index row, Eigen::Index col, T spatial_row_size,
T spatial_col_size, T initial_value)
: RowMajMat<T>(row, col, initial_value), domainRow(spatial_row_size),
domainCol(spatial_col_size) {
this->deltaRow = static_cast<T>(this->domainRow) / static_cast<T>(row);
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(col);
}
UniformGrid(T *concentrations, Eigen::Index n_cells, T spatial_size)
: RowMajMatMap<T>(concentrations, 1, n_cells), domainCol(spatial_size) {
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(n_cells);
}
UniformGrid(T *concentrations, Eigen::Index row, Eigen::Index col,
T spatial_row_size, T spatial_col_size)
: RowMajMatMap<T>(concentrations, row, col), domainRow(spatial_row_size),
domainCol(spatial_col_size) {
this->deltaRow = static_cast<T>(this->domainRow) / static_cast<T>(row);
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(col);
}
/**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
*
* @param domainLength A double value of the domain length. Must be positive.
*/
void setDomain(double domainLength) {
tug_assert(this->dim == 1,
"Grid is not one dimensional, use 2D domain setter!");
tug_assert(domainLength > 0, "Given domain length is not positive!");
this->domainCol = domainLength;
this->deltaCol = double(this->domainCol) / double(this->cols());
}
/**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
*
* @param domainRow A double value of the domain size in y-direction. Must
* be positive.
* @param domainCol A double value of the domain size in x-direction. Must
* be positive.
*/
void setDomain(double domainRow, double domainCol) {
tug_assert(this->dim == 2,
"Grid is not two dimensional, use 1D domain setter!");
tug_assert(domainCol > 0,
"Given domain size in x-direction is not positive!");
tug_assert(domainRow > 0,
"Given domain size in y-direction is not positive!");
this->domainRow = domainRow;
this->domainCol = domainCol;
this->deltaRow =
double(this->domainRow) / double(this->concentrations.rows());
this->deltaCol =
double(this->domainCol) / double(this->concentrations.cols());
}
/**
* @brief Gets the delta value in x-direction.
*
* @return Delta value in x-direction.
*/
T getDeltaCol() const { return this->deltaCol; }
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return Delta value in y-direction.
*/
T getDeltaRow() const {
tug_assert(this->dim == 2,
"Grid is not two dimensional, there is no delta in "
"y-direction!");
return this->deltaRow;
}
private:
T domainCol; // number of domain columns
T domainRow{0}; // number of domain rows
T deltaCol; // delta in x-direction (between columns)
T deltaRow; // delta in y-direction (between rows)
};
using UniformGrid64 = UniformGrid<double>;
using UniformGrid32 = UniformGrid<float>;
} // namespace tug

10
naaice/BTCS_2D_NAAICE.cpp
View File

@ -5,8 +5,7 @@
#include <ostream> #include <ostream>
#include <stdexcept> #include <stdexcept>
#include <string> #include <string>
#include <string_view> #include <tug/Diffusion.hpp>
#include <tug/Simulation.hpp>
#include <vector> #include <vector>
#include "files.hpp" #include "files.hpp"
@ -105,15 +104,14 @@ int main(int argc, char *argv[]) {
// create a grid with a 5 x 10 field // create a grid with a 5 x 10 field
constexpr int row = 5; constexpr int row = 5;
constexpr int col = 10; constexpr int col = 10;
Grid64 grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
grid.setDomain(0.005, 0.01);
const auto init_values_vec = CSVToVector<double>(INPUT_CONC_FILE); const auto init_values_vec = CSVToVector<double>(INPUT_CONC_FILE);
Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col);
grid.setConcentrations(concentrations); UniformGrid64 grid(concentrations);
grid.setDomain(0.005, 0.01);
const double sum_init = concentrations.sum(); const double sum_init = concentrations.sum();
// // (optional) set alphas of the grid, e.g.: // // (optional) set alphas of the grid, e.g.:
@ -141,7 +139,7 @@ int main(int argc, char *argv[]) {
// // ************************ // // ************************
// set up a simulation environment // set up a simulation environment
Simulation simulation = Simulation(grid, bc); // grid,boundary Diffusion simulation(grid, bc); // grid,boundary
// set the timestep of the simulation // set the timestep of the simulation
simulation.setTimestep(360); // timestep simulation.setTimestep(360); // timestep

13
naaice/NAAICE_dble_vs_float.cpp
View File

@ -1,3 +1,4 @@
#include "tug/UniformGrid.hpp"
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <chrono> #include <chrono>
#include <cstdint> #include <cstdint>
@ -7,11 +8,11 @@
#include <stdexcept> #include <stdexcept>
#include <string> #include <string>
#include <string_view> #include <string_view>
#include <tug/Simulation.hpp>
#include <type_traits> #include <type_traits>
#include <vector> #include <vector>
#include "files.hpp" #include <files.hpp>
#include <tug/Diffusion.hpp>
using namespace tug; using namespace tug;
@ -114,15 +115,14 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
std::cout << name << " grid: " << ngrid << std::endl; std::cout << name << " grid: " << ngrid << std::endl;
Grid<T> grid(ngrid, ngrid);
// Grid64 grid(ngrid, ngrid); // Grid64 grid(ngrid, ngrid);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
grid.setDomain(0.1, 0.1);
Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::Constant(ngrid, ngrid, 0); Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::Constant(ngrid, ngrid, 0);
initConc64(50, 50) = 1E-6; initConc64(50, 50) = 1E-6;
grid.setConcentrations(initConc64);
UniformGrid<T> grid(initConc64.data(), ngrid, ngrid);
const T sum_init64 = initConc64.sum(); const T sum_init64 = initConc64.sum();
@ -142,8 +142,7 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
bc.setBoundarySideClosed(BC_SIDE_BOTTOM); bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
// set up a simulation environment // set up a simulation environment
Simulation Sim = Diffusion Sim(grid, bc); // grid_64,boundary,simulation-approach
Simulation<T, app>(grid, bc); // grid_64,boundary,simulation-approach
// Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER);

3
test/CMakeLists.txt
View File

@ -1,6 +1,7 @@
find_package(doctest REQUIRED) find_package(doctest REQUIRED)
add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp) # add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp)
add_executable(testTug setup.cpp testGrid.cpp)
target_link_libraries(testTug doctest::doctest tug) target_link_libraries(testTug doctest::doctest tug)
# get relative path of the CSV file # get relative path of the CSV file

17
test/testBoundary.cpp
View File

@ -1,3 +1,5 @@
#include <Eigen/Core>
#include <Eigen/src/Core/Matrix.h>
#include <doctest/doctest.h> #include <doctest/doctest.h>
#include <iostream> #include <iostream>
#include <stdio.h> #include <stdio.h>
@ -31,8 +33,10 @@ TEST_CASE("BoundaryElement") {
} }
TEST_CASE("Boundary Class") { TEST_CASE("Boundary Class") {
Grid grid1D = Grid64(10); Eigen::VectorXd conc(10);
Grid grid2D = Grid64(10, 12); UnfiormGrid grid1D = UniformGrid64(conc);
Eigen::MatrixXd conc2D(10, 12);
UnfiormGrid grid2D = UniformGrid64(conc2D);
Boundary boundary1D = Boundary(grid1D); Boundary boundary1D = Boundary(grid1D);
Boundary boundary2D = Boundary(grid2D); Boundary boundary2D = Boundary(grid2D);
vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0)); vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
@ -53,26 +57,21 @@ TEST_CASE("Boundary Class") {
CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1);
CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
BC_TYPE_CLOSED); BC_TYPE_CLOSED);
CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP));
CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM));
CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT));
CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP));
CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 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), CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
BC_TYPE_CONSTANT); BC_TYPE_CONSTANT);
CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
boundary1DVector[0].getType()); boundary1DVector[0].getType());
CHECK_NOTHROW(boundary1D.setInnerBoundary(0, inner_condition_value)); CHECK_NOTHROW(boundary1D.setInnerBoundary(0, inner_condition_value));
CHECK_THROWS(boundary1D.setInnerBoundary(0, 0, inner_condition_value));
CHECK_EQ(boundary1D.getInnerBoundary(0), innerBoundary); CHECK_EQ(boundary1D.getInnerBoundary(0), innerBoundary);
CHECK_EQ(boundary1D.getInnerBoundary(1).first, false); CHECK_EQ(boundary1D.getInnerBoundary(1).first, false);
} }
SUBCASE("Boundaries 2D case") { SUBCASE("Boundaries 2D case") {
CHECK_NOTHROW(Boundary boundary(grid1D)); CHECK_NOTHROW(Boundary boundary(grid2D));
CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10);
CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).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_TOP).size(), 12);
@ -85,13 +84,11 @@ TEST_CASE("Boundary Class") {
CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP));
CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 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), CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
BC_TYPE_CONSTANT); BC_TYPE_CONSTANT);
CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
boundary1DVector[0].getType()); boundary1DVector[0].getType());
CHECK_THROWS(boundary2D.setInnerBoundary(0, inner_condition_value));
CHECK_NOTHROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value)); CHECK_NOTHROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value));
CHECK_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary); CHECK_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary);
CHECK_EQ(boundary2D.getInnerBoundary(0, 2).first, false); CHECK_EQ(boundary2D.getInnerBoundary(0, 2).first, false);

50
test/testSimulation.cpp → test/testDiffusion.cpp
View File

@ -1,9 +1,8 @@
#include "TestUtils.hpp" #include "TestUtils.hpp"
#include <tug/Simulation.hpp> #include <tug/Diffusion.hpp>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <doctest/doctest.h> #include <doctest/doctest.h>
#include <stdio.h>
#include <string> #include <string>
// include the configured header file // include the configured header file
@ -13,19 +12,18 @@ using namespace Eigen;
using namespace std; using namespace std;
using namespace tug; using namespace tug;
Grid64 setupSimulation(double timestep, int iterations) { UniformGrid64 setupSimulation(double timestep, int iterations) {
int row = 11; int row = 11;
int col = 11; int col = 11;
int domain_row = 10; int domain_row = 10;
int domain_col = 10; int domain_col = 10;
// Grid // Grid
Grid grid = Grid64(row, col);
grid.setDomain(domain_row, domain_col);
MatrixXd concentrations = MatrixXd::Constant(row, col, 0); MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(5, 5) = 1; concentrations(5, 5) = 1;
grid.setConcentrations(concentrations);
UnfiormGrid grid = UniformGrid64(concentrations);
grid.setDomain(domain_row, domain_col);
MatrixXd alpha = MatrixXd::Constant(row, col, 1); MatrixXd alpha = MatrixXd::Constant(row, col, 1);
for (int i = 0; i < 5; i++) { for (int i = 0; i < 5; i++) {
@ -57,12 +55,12 @@ TEST_CASE("equality to reference matrix with FTCS") {
MatrixXd reference = CSV2Eigen(test_path); MatrixXd reference = CSV2Eigen(test_path);
cout << "FTCS Test: " << endl; cout << "FTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
// Simulation // Simulation
Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc); Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON); // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep); sim.setTimestep(timestep);
sim.setIterations(iterations); sim.setIterations(iterations);
@ -78,11 +76,11 @@ TEST_CASE("equality to reference matrix with BTCS") {
MatrixXd reference = CSV2Eigen(test_path); MatrixXd reference = CSV2Eigen(test_path);
cout << "BTCS Test: " << endl; cout << "BTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
// Simulation // Simulation
Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc); Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON); // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep); sim.setTimestep(timestep);
sim.setIterations(iterations); sim.setIterations(iterations);
@ -94,30 +92,31 @@ TEST_CASE("equality to reference matrix with BTCS") {
TEST_CASE("Initialize environment") { TEST_CASE("Initialize environment") {
int rc = 12; int rc = 12;
Grid64 grid(rc, rc); Eigen::MatrixXd concentrations(rc, rc);
UniformGrid64 grid(concentrations);
Boundary boundary(grid); Boundary boundary(grid);
CHECK_NOTHROW(Simulation sim(grid, boundary)); CHECK_NOTHROW(Diffusion sim(grid, boundary));
} }
TEST_CASE("Simulation environment") { TEST_CASE("Simulation environment") {
int rc = 12; int rc = 12;
Grid64 grid(rc, rc); Eigen::MatrixXd concentrations(rc, rc);
UniformGrid64 grid(concentrations);
grid.initAlpha();
Boundary boundary(grid); Boundary boundary(grid);
Simulation<double, tug::FTCS_APPROACH> sim(grid, boundary); Diffusion<double, tug::FTCS_APPROACH> sim(grid, boundary);
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_NOTHROW(sim.setIterations(2000));
CHECK_EQ(sim.getIterations(), 2000); CHECK_EQ(sim.getIterations(), 2000);
CHECK_THROWS(sim.setIterations(-300));
} }
SUBCASE("set timestep") { SUBCASE("set timestep") {
CHECK_NOTHROW(sim.setTimestep(0.1)); CHECK_NOTHROW(sim.setTimestep(0.1));
CHECK_EQ(sim.getTimestep(), 0.1); CHECK_EQ(sim.getTimestep(), 0.1);
CHECK_THROWS(sim.setTimestep(-0.3));
} }
} }
@ -126,13 +125,12 @@ TEST_CASE("Closed Boundaries - no change expected") {
constexpr std::uint32_t nrows = 5; constexpr std::uint32_t nrows = 5;
constexpr std::uint32_t ncols = 5; constexpr std::uint32_t ncols = 5;
tug::Grid64 grid(nrows, ncols);
auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
grid.setConcentrations(concentrations); tug::UniformGrid64 grid(concentrations);
grid.setAlpha(alphax, alphay); grid.setAlpha(alphax, alphay);
tug::Boundary bc(grid); tug::Boundary bc(grid);
@ -141,7 +139,7 @@ TEST_CASE("Closed Boundaries - no change expected") {
bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0); bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0);
bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0); bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0);
tug::Simulation<double> sim(grid, bc); tug::Diffusion<double> sim(grid, bc);
sim.setTimestep(1); sim.setTimestep(1);
sim.setIterations(1); sim.setIterations(1);
@ -156,20 +154,18 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") {
constexpr std::uint32_t nrows = 5; constexpr std::uint32_t nrows = 5;
constexpr std::uint32_t ncols = 5; constexpr std::uint32_t ncols = 5;
tug::Grid64 grid(nrows, ncols);
auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
grid.setConcentrations(concentrations); tug::UniformGrid64 grid(concentrations);
grid.setAlpha(alphax, alphay); grid.setAlpha(alphax, alphay);
tug::Boundary bc(grid); tug::Boundary bc(grid);
// inner // inner
bc.setInnerBoundary(2, 2, 0); bc.setInnerBoundary(2, 2, 0);
tug::Simulation<double> sim(grid, bc); tug::Diffusion<double> sim(grid, bc);
sim.setTimestep(1); sim.setTimestep(1);
sim.setIterations(1); sim.setIterations(1);
@ -184,4 +180,4 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") {
const bool less_zero = (grid.getConcentrations().array() < 0.0).any(); const bool less_zero = (grid.getConcentrations().array() < 0.0).any();
CHECK(less_zero == false); CHECK(less_zero == false);
} }

348
test/testGrid.cpp
View File

@ -1,268 +1,198 @@
#include "tug/UniformGrid.hpp"
#include <Eigen/Core> #include <Eigen/Core>
#include <Eigen/src/Core/Matrix.h>
#include <doctest/doctest.h> #include <doctest/doctest.h>
#include <tug/Grid.hpp> #include <tug/UniformGrid.hpp>
#include <vector> #include <vector>
using namespace Eigen; using namespace Eigen;
using namespace std; using namespace std;
using namespace tug; using namespace tug;
TEST_CASE("1D Grid, too small length") {
int l = 2;
CHECK_THROWS(Grid64(l));
}
TEST_CASE("2D Grid64, too small side") {
int r = 1;
int c = 4;
CHECK_THROWS(Grid64(r, c));
r = 4;
c = 1;
CHECK_THROWS(Grid64(r, c));
}
TEST_CASE("1D Grid64") { TEST_CASE("1D Grid64") {
int l = 12; int l = 12;
Grid64 grid(l); Eigen::VectorXd conc(l);
UniformGrid64 grid(l, l, 1);
SUBCASE("correct construction") { SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 1); CHECK_EQ(grid.getDim(), 1);
CHECK_EQ(grid.getLength(), l); CHECK_EQ(grid.getCol(), l);
CHECK_EQ(grid.getCol(), l); CHECK_EQ(grid.getCol(), l);
CHECK_EQ(grid.getRow(), 1); CHECK_EQ(grid.getRow(), 1);
CHECK_EQ(grid.getConcentrations().rows(), 1);
CHECK_EQ(grid.getConcentrations().cols(), l);
CHECK_EQ(grid.getAlpha().rows(), 1);
CHECK_EQ(grid.getAlpha().cols(), l);
CHECK_EQ(grid.getDeltaCol(), 1); CHECK_EQ(grid.getDeltaCol(), 1);
CHECK_THROWS(grid.getAlphaX()); // CHECK_EQ(grid.getConcentrations().rows(), 1);
CHECK_THROWS(grid.getAlphaY()); // CHECK_EQ(grid.getConcentrations().cols(), l);
CHECK_THROWS(grid.getDeltaRow()); // CHECK_EQ(grid.getAlphaX().rows(), 1);
// CHECK_EQ(grid.getAlphaX().cols(), l);
} }
SUBCASE("setting concentrations") { // SUBCASE("setting alpha") {
// correct concentrations matrix // // correct alpha matrix
MatrixXd concentrations = MatrixXd::Constant(1, l, 3); // MatrixXd alpha = MatrixXd::Constant(1, l, 3);
CHECK_NOTHROW(grid.setConcentrations(concentrations)); // CHECK_NOTHROW(grid.setAlpha(alpha));
// false concentrations matrix // grid.setAlpha(alpha);
MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); // CHECK_EQ(grid.getAlphaX(), alpha);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); // }
}
SUBCASE("setting alpha") {
// correct alpha matrix
MatrixXd alpha = MatrixXd::Constant(1, l, 3);
CHECK_NOTHROW(grid.setAlpha(alpha));
CHECK_THROWS(grid.setAlpha(alpha, alpha));
grid.setAlpha(alpha);
CHECK_EQ(grid.getAlpha(), alpha);
CHECK_THROWS(grid.getAlphaX());
CHECK_THROWS(grid.getAlphaY());
// false alpha matrix
MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
CHECK_THROWS(grid.setAlpha(wAlpha));
}
SUBCASE("setting domain") { SUBCASE("setting domain") {
int d = 8; int d = 8;
// set 1D domain // set 1D domain
CHECK_NOTHROW(grid.setDomain(d)); CHECK_NOTHROW(grid.setDomain(d));
// set 2D domain
CHECK_THROWS(grid.setDomain(d, d));
grid.setDomain(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
d = 0;
CHECK_THROWS(grid.setDomain(d));
d = -2;
CHECK_THROWS(grid.setDomain(d));
} }
} }
TEST_CASE("2D Grid64 quadratic") { // TEST_CASE("2D Grid64 quadratic") {
int rc = 12; // int rc = 12;
Grid64 grid(rc, rc); // Eigen::MatrixXd conc(rc, rc);
// Grid64 grid(conc);
// grid.initAlpha();
SUBCASE("correct construction") { // SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2); // CHECK_EQ(grid.getDim(), 2);
CHECK_THROWS(grid.getLength()); // CHECK_EQ(grid.getCol(), rc);
CHECK_EQ(grid.getCol(), rc); // CHECK_EQ(grid.getRow(), rc);
CHECK_EQ(grid.getRow(), rc);
CHECK_EQ(grid.getConcentrations().rows(), rc); // CHECK_EQ(grid.getConcentrations().rows(), rc);
CHECK_EQ(grid.getConcentrations().cols(), rc); // CHECK_EQ(grid.getConcentrations().cols(), rc);
CHECK_THROWS(grid.getAlpha());
CHECK_EQ(grid.getAlphaX().rows(), rc); // CHECK_EQ(grid.getAlphaX().rows(), rc);
CHECK_EQ(grid.getAlphaX().cols(), rc); // CHECK_EQ(grid.getAlphaX().cols(), rc);
CHECK_EQ(grid.getAlphaY().rows(), rc); // CHECK_EQ(grid.getAlphaY().rows(), rc);
CHECK_EQ(grid.getAlphaY().cols(), rc); // CHECK_EQ(grid.getAlphaY().cols(), rc);
CHECK_EQ(grid.getDeltaRow(), 1); // CHECK_EQ(grid.getDeltaRow(), 1);
CHECK_EQ(grid.getDeltaCol(), 1); // CHECK_EQ(grid.getDeltaCol(), 1);
} // }
SUBCASE("setting concentrations") { // SUBCASE("setting alphas") {
// correct concentrations matrix // // correct alpha matrices
MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); // MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
CHECK_NOTHROW(grid.setConcentrations(concentrations)); // MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// false concentrations matrix // grid.setAlpha(alphax, alphay);
MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); // CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); // CHECK_EQ(grid.getAlphaY(), alphay);
wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); // }
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
SUBCASE("setting alphas") { // SUBCASE("setting domain") {
// correct alpha matrices // int dr = 8;
MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); // int dc = 9;
MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.setAlpha(alphax)); // // set 2D domain
// CHECK_NOTHROW(grid.setDomain(dr, dc));
grid.setAlpha(alphax, alphay); // grid.setDomain(dr, dc);
CHECK_EQ(grid.getAlphaX(), alphax); // CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc));
CHECK_EQ(grid.getAlphaY(), alphay); // CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
CHECK_THROWS(grid.getAlpha()); // }
// }
// false alpha matrices // TEST_CASE("2D Grid64 non-quadratic") {
alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); // int r = 12;
CHECK_THROWS(grid.setAlpha(alphax, alphay)); // int c = 15;
alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); // Eigen::MatrixXd conc(r, c);
CHECK_THROWS(grid.setAlpha(alphax, alphay)); // Grid64 grid(conc);
} // grid.initAlpha();
SUBCASE("setting domain") { // SUBCASE("correct construction") {
int dr = 8; // CHECK_EQ(grid.getDim(), 2);
int dc = 9; // CHECK_EQ(grid.getCol(), c);
// CHECK_EQ(grid.getRow(), r);
// set 1D domain // CHECK_EQ(grid.getConcentrations().rows(), r);
CHECK_THROWS(grid.setDomain(dr)); // CHECK_EQ(grid.getConcentrations().cols(), c);
// set 2D domain // CHECK_EQ(grid.getAlphaX().rows(), r);
CHECK_NOTHROW(grid.setDomain(dr, dc)); // CHECK_EQ(grid.getAlphaX().cols(), c);
// CHECK_EQ(grid.getAlphaY().rows(), r);
// CHECK_EQ(grid.getAlphaY().cols(), c);
// CHECK_EQ(grid.getDeltaRow(), 1);
// CHECK_EQ(grid.getDeltaCol(), 1);
// }
grid.setDomain(dr, dc); // SUBCASE("setting alphas") {
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); // // correct alpha matrices
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); // MatrixXd alphax = MatrixXd::Constant(r, c, 2);
// MatrixXd alphay = MatrixXd::Constant(r, c, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// set too small domain // grid.setAlpha(alphax, alphay);
dr = 0; // CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_THROWS(grid.setDomain(dr, dc)); // CHECK_EQ(grid.getAlphaY(), alphay);
dr = 8; // }
dc = 0;
CHECK_THROWS(grid.setDomain(dr, dc));
dr = -2;
CHECK_THROWS(grid.setDomain(dr, dc));
}
}
TEST_CASE("2D Grid64 non-quadratic") { // SUBCASE("setting domain") {
int r = 12; // int dr = 8;
int c = 15; // int dc = 9;
Grid64 grid(r, c);
SUBCASE("correct construction") { // // set 2D domain
CHECK_EQ(grid.getDim(), 2); // CHECK_NOTHROW(grid.setDomain(dr, dc));
CHECK_THROWS(grid.getLength());
CHECK_EQ(grid.getCol(), c);
CHECK_EQ(grid.getRow(), r);
CHECK_EQ(grid.getConcentrations().rows(), r); // grid.setDomain(dr, dc);
CHECK_EQ(grid.getConcentrations().cols(), c); // CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
CHECK_THROWS(grid.getAlpha()); // CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
// }
// }
CHECK_EQ(grid.getAlphaX().rows(), r); // TEST_CASE("2D Grid64 non-quadratic from pointer") {
CHECK_EQ(grid.getAlphaX().cols(), c); // int r = 4;
CHECK_EQ(grid.getAlphaY().rows(), r); // int c = 5;
CHECK_EQ(grid.getAlphaY().cols(), c); // std::vector<double> concentrations(r * c);
CHECK_EQ(grid.getDeltaRow(), 1);
CHECK_EQ(grid.getDeltaCol(), 1);
}
SUBCASE("setting concentrations") { // for (int i = 0; i < r * c; i++) {
// correct concentrations matrix // concentrations[i] = i;
MatrixXd concentrations = MatrixXd::Constant(r, c, 2); // }
CHECK_NOTHROW(grid.setConcentrations(concentrations)); // Grid64 grid(concentrations.data(), r, c);
// grid.initAlpha();
// false concentrations matrix // SUBCASE("correct construction") {
MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); // CHECK_EQ(grid.getDim(), 2);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); // CHECK_EQ(grid.getCol(), c);
wConcentrations = MatrixXd::Constant(r + 3, c, 3); // CHECK_EQ(grid.getRow(), r);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
SUBCASE("setting alphas") { // CHECK_EQ(grid.getConcentrations().rows(), r);
// correct alpha matrices // CHECK_EQ(grid.getConcentrations().cols(), c);
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.setAlpha(alphax)); // CHECK_EQ(grid.getAlphaX().rows(), r);
// CHECK_EQ(grid.getAlphaX().cols(), c);
// CHECK_EQ(grid.getAlphaY().rows(), r);
// CHECK_EQ(grid.getAlphaY().cols(), c);
// CHECK_EQ(grid.getDeltaRow(), 1);
// CHECK_EQ(grid.getDeltaCol(), 1);
// }
grid.setAlpha(alphax, alphay); // SUBCASE("setting alphas") {
CHECK_EQ(grid.getAlphaX(), alphax); // // correct alpha matrices
CHECK_EQ(grid.getAlphaY(), alphay); // MatrixXd alphax = MatrixXd::Constant(r, c, 2);
CHECK_THROWS(grid.getAlpha()); // MatrixXd alphay = MatrixXd::Constant(r, c, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// false alpha matrices // grid.setAlpha(alphax, alphay);
alphax = MatrixXd::Constant(r + 3, c + 1, 3); // CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_THROWS(grid.setAlpha(alphax, alphay)); // CHECK_EQ(grid.getAlphaY(), alphay);
alphay = MatrixXd::Constant(r + 2, c + 1, 5); // }
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
SUBCASE("setting domain") { // SUBCASE("setting domain") {
int dr = 8; // int dr = 8;
int dc = 9; // int dc = 9;
// set 1D domain // // set 2D domain
CHECK_THROWS(grid.setDomain(dr)); // CHECK_NOTHROW(grid.setDomain(dr, dc));
// set 2D domain // grid.setDomain(dr, dc);
CHECK_NOTHROW(grid.setDomain(dr, dc)); // CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
// }
grid.setDomain(dr, dc); // SUBCASE("correct values") {
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); // CHECK_EQ(grid.getConcentrations()(0, 0), 0);
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); // CHECK_EQ(grid.getConcentrations()(0, 1), 1);
// CHECK_EQ(grid.getConcentrations()(1, 0), c);
// set too small domain // CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1);
dr = 0; // }
CHECK_THROWS(grid.setDomain(dr, dc)); // }
dr = 8;
dc = -1;
CHECK_THROWS(grid.setDomain(dr, dc));
dr = -2;
CHECK_THROWS(grid.setDomain(dr, dc));
}
SUBCASE("set concentration from pointer") {
std::vector<double> concentrations(r * c);
for (int i = 0; i < r * c; i++) {
concentrations[i] = i;
}
grid.setConcentrations(concentrations.data());
CHECK_EQ(grid.getConcentrations()(0, 0), 0);
CHECK_EQ(grid.getConcentrations()(0, 1), 1);
CHECK_EQ(grid.getConcentrations()(1, 0), c);
}
}