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

refactor: format all source files to LLVM standard

Browse Source
This commit is contained in:
Max Lübke 2023-09-14 10:24:03 +02:00
parent 81774e72c1
commit 0d34752837
24 changed files with 2413 additions and 2416 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
examples/BTCS_1D_proto_example.cpp
View File

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

130
examples/BTCS_2D_proto_example.cpp
View File

@ -1,84 +1,82 @@
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
// EASY_PROFILER_ENABLE; // EASY_PROFILER_ENABLE;
// profiler::startListen(); // profiler::startListen();
// ************** // **************
// **** GRID **** // **** GRID ****
// ************** // **************
// profiler::startListen(); // profiler::startListen();
// 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;
Grid grid = Grid(row,col); Grid grid = Grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
// grid.setDomain(20, 20); // grid.setDomain(20, 20);
// (optional) set the concentrations, e.g.: // (optional) set the concentrations, e.g.:
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); //
// 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); grid.setConcentrations(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
// MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value
// grid.setAlpha(alphax, alphay); // grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
// ****************** // create a boundary with constant values
// **** BOUNDARY **** Boundary bc = Boundary(grid);
// ****************** bc.setBoundarySideClosed(BC_SIDE_LEFT);
bc.setBoundarySideClosed(BC_SIDE_RIGHT);
bc.setBoundarySideClosed(BC_SIDE_TOP);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
// bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
// bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
// bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
// bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
// create a boundary with constant values // (optional) set boundary condition values for one side, e.g.:
Boundary bc = Boundary(grid); // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
bc.setBoundarySideClosed(BC_SIDE_LEFT); // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values
bc.setBoundarySideClosed(BC_SIDE_RIGHT); // VectorXd bc_zero_values = VectorXd::Constant(20,0);
bc.setBoundarySideClosed(BC_SIDE_TOP); // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM); // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
// bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // VectorXd bc_front_values = VectorXd::Constant(20,2000);
// bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
// bc.setBoundarySideConstant(BC_SIDE_TOP, 0); // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
// bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
// ************************
// **** SIMULATION ENV ****
// ************************
// (optional) set boundary condition values for one side, e.g.: // set up a simulation environment
// VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value Simulation simulation =
// bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
// VectorXd bc_zero_values = VectorXd::Constant(20,0);
// bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
// bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
// VectorXd bc_front_values = VectorXd::Constant(20,2000);
// bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
// bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
// ************************ // set the number of iterations
// **** SIMULATION ENV **** simulation.setIterations(300);
// ************************
// set up a simulation environment // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach // CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_XTREME);
// set the timestep of the simulation // **** RUN SIMULATION ****
simulation.setTimestep(0.1); // timestep
// set the number of iterations // run the simulation
simulation.setIterations(300);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] // EASY_BLOCK("SIMULATION")
simulation.setOutputCSV(CSV_OUTPUT_XTREME); simulation.run();
// EASY_END_BLOCK;
// **** RUN SIMULATION **** // profiler::dumpBlocksToFile("test_profile.prof");
// profiler::stopListen();
// run the simulation
// EASY_BLOCK("SIMULATION")
simulation.run();
// EASY_END_BLOCK;
// profiler::dumpBlocksToFile("test_profile.prof");
// profiler::stopListen();
} }

32
examples/CRNI_2D_proto_example.cpp
View File

@ -1,24 +1,24 @@
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int row = 20; int row = 20;
int col = 20; int col = 20;
Grid grid(row, col); Grid grid(row, col);
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); grid.setConcentrations(concentrations);
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
bc.setBoundarySideClosed(BC_SIDE_LEFT); bc.setBoundarySideClosed(BC_SIDE_LEFT);
bc.setBoundarySideClosed(BC_SIDE_RIGHT); bc.setBoundarySideClosed(BC_SIDE_RIGHT);
bc.setBoundarySideClosed(BC_SIDE_TOP); bc.setBoundarySideClosed(BC_SIDE_TOP);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM); bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH);
simulation.setTimestep(0.1); simulation.setTimestep(0.1);
simulation.setIterations(50); simulation.setIterations(50);
simulation.setOutputCSV(CSV_OUTPUT_XTREME); simulation.setOutputCSV(CSV_OUTPUT_XTREME);
simulation.run(); simulation.run();
} }

62
examples/FTCS_1D_proto_example.cpp
View File

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

114
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -1,8 +1,9 @@
/** /**
* @file FTCS_2D_proto_closed_mdl.cpp * @file FTCS_2D_proto_closed_mdl.cpp
* @author Hannes Signer, Philipp Ungrund, MDL * @author Hannes Signer, Philipp Ungrund, MDL
* @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary condition; optional command line argument: number of cols and rows * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary
* * condition; optional command line argument: number of cols and rows
*
*/ */
#include <cstdlib> #include <cstdlib>
@ -11,75 +12,76 @@
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int row = 64; int row = 64;
if (argc == 2) { if (argc == 2) {
// no cmd line argument, take col=row=64 // no cmd line argument, take col=row=64
row = atoi(argv[1]); row = atoi(argv[1]);
} }
int col=row; int col = row;
std::cout << "Nrow =" << row << std::endl; std::cout << "Nrow =" << row << std::endl;
// ************** // **************
// **** GRID **** // **** GRID ****
// ************** // **************
// 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;
Grid grid = Grid(row,col); Grid grid = Grid(row, col);
// (optional) set the domain, e.g.: // (optional) set the domain, e.g.:
// grid.setDomain(20, 20); // grid.setDomain(20, 20);
// (optional) set the concentrations, e.g.: // (optional) set the concentrations, e.g.:
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); //
MatrixXd concentrations = MatrixXd::Constant(row,col,0); // #row,#col,value
concentrations(n2,n2) = 1; MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(n2,n2+1) = 1; concentrations(n2, n2) = 1;
concentrations(n2+1,n2) = 1; concentrations(n2, n2 + 1) = 1;
concentrations(n2+1,n2+1) = 1; concentrations(n2 + 1, n2) = 1;
grid.setConcentrations(concentrations); concentrations(n2 + 1, n2 + 1) = 1;
grid.setConcentrations(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
MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value
grid.setAlpha(alphax, alphay); grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
// ****************** // create a boundary with constant values
// **** BOUNDARY **** Boundary bc = Boundary(grid);
// ******************
// create a boundary with constant values // (optional) set boundary condition values for one side, e.g.:
Boundary bc = Boundary(grid); bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values
bc.setBoundarySideClosed(BC_SIDE_RIGHT);
bc.setBoundarySideClosed(BC_SIDE_TOP);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
// (optional) set boundary condition values for one side, e.g.: // ************************
bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values // **** SIMULATION ENV ****
bc.setBoundarySideClosed(BC_SIDE_RIGHT); // ************************
bc.setBoundarySideClosed(BC_SIDE_TOP);
bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
// set up a simulation environment
Simulation simulation =
Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// ************************ // set the timestep of the simulation
// **** SIMULATION ENV **** simulation.setTimestep(10000); // timestep
// ************************
// set up a simulation environment // set the number of iterations
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach simulation.setIterations(100);
// set the timestep of the simulation // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
simulation.setTimestep(10000); // timestep // CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// set the number of iterations // **** RUN SIMULATION ****
simulation.setIterations(100);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] // run the simulation
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); simulation.run();
// **** RUN SIMULATION ****
// run the simulation return 0;
simulation.run();
return 0;
} }

130
examples/FTCS_2D_proto_example.cpp
View File

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

105
examples/FTCS_2D_proto_example_mdl.cpp
View File

@ -1,77 +1,78 @@
/** /**
* @file FTCS_2D_proto_example.cpp * @file FTCS_2D_proto_example.cpp
* @author Hannes Signer, Philipp Ungrund * @author Hannes Signer, Philipp Ungrund
* @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach * @brief Creates a prototypical standard TUG simulation in 2D with FTCS
* and constant boundary condition * approach and constant boundary condition
* *
*/ */
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
// **************
// **** GRID ****
// **************
// create a grid with a 20 x 20 field // **************
int row = 64; // **** GRID ****
int col = 64; // **************
int n2 = row/2-1;
Grid grid = Grid(row,col);
// (optional) set the domain, e.g.: // create a grid with a 20 x 20 field
// grid.setDomain(20, 20); int row = 64;
int col = 64;
int n2 = row / 2 - 1;
Grid grid = Grid(row, col);
// (optional) set the concentrations, e.g.: // (optional) set the domain, e.g.:
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // grid.setDomain(20, 20);
MatrixXd concentrations = MatrixXd::Constant(row,col,0);
concentrations(n2,n2) = 1;
concentrations(n2,n2+1) = 1;
concentrations(n2+1,n2) = 1;
concentrations(n2+1,n2+1) = 1;
grid.setConcentrations(concentrations);
// (optional) set alphas of the grid, e.g.: // (optional) set the concentrations, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); //
MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value // #row,#col,value
grid.setAlpha(alphax, alphay); MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(n2, n2) = 1;
concentrations(n2, n2 + 1) = 1;
concentrations(n2 + 1, n2) = 1;
concentrations(n2 + 1, n2 + 1) = 1;
grid.setConcentrations(concentrations);
// (optional) set alphas of the grid, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value
MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value
grid.setAlpha(alphax, alphay);
// ****************** // ******************
// **** BOUNDARY **** // **** BOUNDARY ****
// ****************** // ******************
// create a boundary with constant values // create a boundary with constant values
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
// (optional) set boundary condition values for one side, e.g.: // (optional) set boundary condition values for one side, e.g.:
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
bc.setBoundarySideConstant(BC_SIDE_TOP, 0); bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
// ************************
// **** SIMULATION ENV ****
// ************************
// ************************ // set up a simulation environment
// **** SIMULATION ENV **** Simulation simulation =
// ************************ Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
// set up a simulation environment // (optional) set the timestep of the simulation
Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach simulation.setTimestep(1000); // timestep
// (optional) set the timestep of the simulation // (optional) set the number of iterations
simulation.setTimestep(1000); // timestep simulation.setIterations(5);
// (optional) set the number of iterations // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
simulation.setIterations(5); // CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_OFF);
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] // **** RUN SIMULATION ****
simulation.setOutputCSV(CSV_OUTPUT_OFF);
// **** RUN SIMULATION ****
// run the simulation // run the simulation
simulation.run(); simulation.run();
return 0; return 0;
} }

93
examples/profiling_openmp.cpp
View File

@ -1,67 +1,66 @@
#include <chrono>
#include <fstream>
#include <iostream>
#include <string> #include <string>
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
#include <chrono>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int n[] = {2000}; int n[] = {2000};
int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
int iterations[1] = {1}; int iterations[1] = {1};
int repetition = 10; int repetition = 10;
for (int l = 0; l < size(threads); l++) {
for(int l=0; l<size(threads); l++){ // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
//string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
ofstream myfile; ofstream myfile;
myfile.open("speedup_1000.csv", std::ios::app); myfile.open("speedup_1000.csv", std::ios::app);
myfile << "Number threads: " << threads[l] << endl; myfile << "Number threads: " << threads[l] << endl;
for (int i = 0; i < size(n); i++){ for (int i = 0; i < size(n); i++) {
cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
//myfile << "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++){ for (int j = 0; j < size(iterations); j++) {
cout << "Iterations: " << iterations[j] << endl; cout << "Iterations: " << iterations[j] << endl;
//myfile << "Iterations: " << iterations[j] << endl; // myfile << "Iterations: " << iterations[j] << endl;
for (int k = 0; k < repetition; k++){ for (int k = 0; k < repetition; k++) {
cout << "Wiederholung: " << k << endl; cout << "Wiederholung: " << k << endl;
Grid grid = Grid(n[i], n[i]); Grid grid = Grid(n[i], n[i]);
grid.setDomain(1, 1); grid.setDomain(1, 1);
MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
concentrations(n[i]/2,n[i]/2) = 1; concentrations(n[i] / 2, n[i] / 2) = 1;
grid.setConcentrations(concentrations); grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
Simulation sim = Simulation(grid, bc, BTCS_APPROACH); Simulation sim = Simulation(grid, bc, BTCS_APPROACH);
sim.setSolver(THOMAS_ALGORITHM_SOLVER); sim.setSolver(THOMAS_ALGORITHM_SOLVER);
if(argc == 2){ if (argc == 2) {
int numThreads = atoi(argv[1]); int numThreads = atoi(argv[1]);
sim.setNumberThreads(numThreads); sim.setNumberThreads(numThreads);
} } else {
else{ sim.setNumberThreads(threads[l]);
sim.setNumberThreads(threads[l]); }
}
sim.setTimestep(0.01); sim.setTimestep(0.01);
sim.setIterations(iterations[j]); sim.setIterations(iterations[j]);
sim.setOutputCSV(CSV_OUTPUT_OFF); sim.setOutputCSV(CSV_OUTPUT_OFF);
auto begin = std::chrono::high_resolution_clock::now(); auto begin = std::chrono::high_resolution_clock::now();
sim.run(); sim.run();
auto end = std::chrono::high_resolution_clock::now(); auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin); auto milliseconds =
myfile << milliseconds.count() << endl; std::chrono::duration_cast<std::chrono::milliseconds>(end -
} begin);
myfile << milliseconds.count() << endl;
} }
cout << endl; }
myfile << endl; cout << endl;
myfile << endl;
} }
myfile.close(); myfile.close();
} }
} }

93
examples/profiling_speedup.cpp
View File

@ -1,67 +1,66 @@
#include <chrono>
#include <fstream>
#include <iostream>
#include <string> #include <string>
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
#include <chrono>
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int n[] = {2000}; int n[] = {2000};
int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
int iterations[1] = {1}; int iterations[1] = {1};
int repetition = 10; int repetition = 10;
for (int l = 0; l < size(threads); l++) {
for(int l=0; l<size(threads); l++){ // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
//string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
ofstream myfile; ofstream myfile;
myfile.open("speedup_1000.csv", std::ios::app); myfile.open("speedup_1000.csv", std::ios::app);
myfile << "Number threads: " << threads[l] << endl; myfile << "Number threads: " << threads[l] << endl;
for (int i = 0; i < size(n); i++){ for (int i = 0; i < size(n); i++) {
cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
//myfile << "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++){ for (int j = 0; j < size(iterations); j++) {
cout << "Iterations: " << iterations[j] << endl; cout << "Iterations: " << iterations[j] << endl;
//myfile << "Iterations: " << iterations[j] << endl; // myfile << "Iterations: " << iterations[j] << endl;
for (int k = 0; k < repetition; k++){ for (int k = 0; k < repetition; k++) {
cout << "Wiederholung: " << k << endl; cout << "Wiederholung: " << k << endl;
Grid grid = Grid(n[i], n[i]); Grid grid = Grid(n[i], n[i]);
grid.setDomain(1, 1); grid.setDomain(1, 1);
MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
concentrations(n[i]/2,n[i]/2) = 1; concentrations(n[i] / 2, n[i] / 2) = 1;
grid.setConcentrations(concentrations); grid.setConcentrations(concentrations);
MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
Boundary bc = Boundary(grid); Boundary bc = Boundary(grid);
Simulation sim = Simulation(grid, bc, BTCS_APPROACH); Simulation sim = Simulation(grid, bc, BTCS_APPROACH);
sim.setSolver(THOMAS_ALGORITHM_SOLVER); sim.setSolver(THOMAS_ALGORITHM_SOLVER);
if(argc == 2){ if (argc == 2) {
int numThreads = atoi(argv[1]); int numThreads = atoi(argv[1]);
sim.setNumberThreads(numThreads); sim.setNumberThreads(numThreads);
} } else {
else{ sim.setNumberThreads(threads[l]);
sim.setNumberThreads(threads[l]); }
}
sim.setTimestep(0.01); sim.setTimestep(0.01);
sim.setIterations(iterations[j]); sim.setIterations(iterations[j]);
sim.setOutputCSV(CSV_OUTPUT_OFF); sim.setOutputCSV(CSV_OUTPUT_OFF);
auto begin = std::chrono::high_resolution_clock::now(); auto begin = std::chrono::high_resolution_clock::now();
sim.run(); sim.run();
auto end = std::chrono::high_resolution_clock::now(); auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin); auto milliseconds =
myfile << milliseconds.count() << endl; std::chrono::duration_cast<std::chrono::milliseconds>(end -
} begin);
myfile << milliseconds.count() << endl;
} }
cout << endl; }
myfile << endl; cout << endl;
myfile << endl;
} }
myfile.close(); myfile.close();
} }
} }

74
examples/reference-FTCS_2D_closed.cpp
View File

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

333
include/tug/Boundary.hpp
View File

@ -1,37 +1,29 @@
/** /**
* @file Boundary.hpp * @file Boundary.hpp
* @brief API of Boundary class, that holds all information for each boundary condition * @brief API of Boundary class, that holds all information for each boundary
* at the edges of the diffusion grid. * condition at the edges of the diffusion grid.
* *
*/ */
#ifndef BOUNDARY_H_ #ifndef BOUNDARY_H_
#define BOUNDARY_H_ #define BOUNDARY_H_
#include <cstddef>
#include "Grid.hpp" #include "Grid.hpp"
#include <cstddef>
using namespace std; using namespace std;
using namespace Eigen; using namespace Eigen;
/** /**
* @brief Enum defining the two implemented boundary conditions. * @brief Enum defining the two implemented boundary conditions.
* *
*/ */
enum BC_TYPE { enum BC_TYPE { BC_TYPE_CLOSED, BC_TYPE_CONSTANT };
BC_TYPE_CLOSED,
BC_TYPE_CONSTANT
};
/** /**
* @brief Enum defining all 4 possible sides to a 1D and 2D grid. * @brief Enum defining all 4 possible sides to a 1D and 2D grid.
* *
*/ */
enum BC_SIDE { enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM };
BC_SIDE_LEFT,
BC_SIDE_RIGHT,
BC_SIDE_TOP,
BC_SIDE_BOTTOM
};
/** /**
* This class defines the boundary conditions of individual boundary elements. * This class defines the boundary conditions of individual boundary elements.
@ -39,177 +31,184 @@ enum BC_SIDE {
* The class serves as an auxiliary class for structuring the Boundary class. * The class serves as an auxiliary class for structuring the Boundary class.
*/ */
class BoundaryElement { class BoundaryElement {
public: public:
/**
/** * @brief Construct a new Boundary Element object for the closed case.
* @brief Construct a new Boundary Element object for the closed case. * The boundary type is here automatically set to the type
* The boundary type is here automatically set to the type * BC_TYPE_CLOSED, where the value takes -1 and does not hold any
* BC_TYPE_CLOSED, where the value takes -1 and does not hold any * physical meaning.
* physical meaning. */
*/ BoundaryElement();
BoundaryElement();
/** /**
* @brief Construct a new Boundary Element object for the constant case. * @brief Construct a new Boundary Element object for the constant case.
* The boundary type is automatically set to the type * The boundary type is automatically set to the type
* BC_TYPE_CONSTANT. * BC_TYPE_CONSTANT.
* *
* @param value Value of the constant concentration to be assumed at the * @param value Value of the constant concentration to be assumed at the
* corresponding boundary element. * corresponding boundary element.
*/ */
BoundaryElement(double value); BoundaryElement(double value);
/** /**
* @brief Allows changing the boundary type of a corresponding * @brief Allows changing the boundary type of a corresponding
* BoundaryElement object. * BoundaryElement object.
* *
* @param type Type of boundary condition. Either BC_TYPE_CONSTANT or * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or
BC_TYPE_CLOSED. BC_TYPE_CLOSED.
*/ */
void setType(BC_TYPE type); void setType(BC_TYPE type);
/**
* @brief Sets the value of a boundary condition for the constant case.
*
* @param value Concentration to be considered constant for the
* corresponding boundary element.
*/
void setValue(double value);
/** /**
* @brief Return the type of the boundary condition, i.e. whether the * @brief Sets the value of a boundary condition for the constant case.
* boundary is considered closed or constant. *
* * @param value Concentration to be considered constant for the
* @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or * corresponding boundary element.
BC_TYPE_CONSTANT. */
*/ void setValue(double value);
BC_TYPE getType();
/** /**
* @brief Return the concentration value for the constant boundary condition. * @brief Return the type of the boundary condition, i.e. whether the
* * boundary is considered closed or constant.
* @return double Value of the concentration. *
*/ * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or
double getValue(); BC_TYPE_CONSTANT.
*/
BC_TYPE getType();
private: /**
BC_TYPE type; * @brief Return the concentration value for the constant boundary condition.
double value; *
* @return double Value of the concentration.
*/
double getValue();
private:
BC_TYPE type;
double value;
}; };
/** /**
* This class implements the functionality and management of the boundary * This class implements the functionality and management of the boundary
* conditions in the grid to be simulated. * conditions in the grid to be simulated.
*/ */
class Boundary { class Boundary {
public: public:
/** /**
* @brief Creates a boundary object based on the passed grid object and * @brief Creates a boundary object based on the passed grid object and
* initializes the boundaries as closed. * initializes the boundaries as closed.
* *
* @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(Grid grid); Boundary(Grid grid);
/** /**
* @brief Sets all elements of the specified boundary side to the boundary * @brief Sets all elements of the specified boundary side to the boundary
* condition closed. * condition closed.
* *
* @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);
/** /**
* @brief Sets all elements of the specified boundary side to the boundary * @brief Sets all elements of the specified boundary side to the boundary
* condition constant. Thereby the concentration values of the * condition constant. Thereby the concentration values of the
* boundaries are set to the passed value. * boundaries are set to the passed value.
* *
* @param side Side to be set to constant, e.g. BC_SIDE_LEFT. * @param side Side to be set to constant, e.g. BC_SIDE_LEFT.
* @param value Concentration to be set for all elements of the specified page. * @param value Concentration to be set for all elements of the specified
*/ * page.
void setBoundarySideConstant(BC_SIDE side, double value); */
void setBoundarySideConstant(BC_SIDE side, double value);
/** /**
* @brief Specifically sets the boundary element of the specified side * @brief Specifically sets the boundary element of the specified side
* defined by the index to the boundary condition closed. * defined by the index to the boundary condition closed.
* *
* @param side Side in which an element is to be defined as closed. * @param side Side in which an element is to be defined as closed.
* @param index Index of the boundary element on the corresponding * @param index Index of the boundary element on the corresponding
* boundary side. Must index an element of the corresponding side. * boundary side. Must index an element of the corresponding
*/ * side.
void setBoundaryElementClosed(BC_SIDE side, int index); */
void setBoundaryElementClosed(BC_SIDE side, int index);
/** /**
* @brief Specifically sets the boundary element of the specified side * @brief Specifically sets the boundary element of the specified side
* defined by the index to the boundary condition constant with the * defined by the index to the boundary condition constant with the
given concentration value. given concentration value.
* *
* @param side Side in which an element is to be defined as constant. * @param side Side in which an element is to be defined as constant.
* @param index Index of the boundary element on the corresponding * @param index Index of the boundary element on the corresponding
* boundary side. Must index an element of the corresponding side. * boundary side. Must index an element of the corresponding
* @param value Concentration value to which the boundary element should be set. side.
*/ * @param value Concentration value to which the boundary element should be
void setBoundaryElementConstant(BC_SIDE side, int index, double value); set.
*/
void setBoundaryElementConstant(BC_SIDE side, int index, double value);
/** /**
* @brief Returns the boundary condition of a specified side as a vector * @brief Returns the boundary condition of a specified side as a vector
* of BoundarsElement objects. * of BoundarsElement objects.
* *
* @param side Boundary side from which the boundary conditions are to be returned. * @param side Boundary side from which the boundary conditions are to be
* @return vector<BoundaryElement> Contains the boundary conditions as BoundaryElement objects. * returned.
*/ * @return vector<BoundaryElement> Contains the boundary conditions as
const vector<BoundaryElement> getBoundarySide(BC_SIDE side); * BoundaryElement objects.
*/
const vector<BoundaryElement> getBoundarySide(BC_SIDE side);
/** /**
* @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some
boundary is closed. specific boundary is closed.
* *
* @param side Boundary side for which the values are to be returned. * @param side Boundary side for which the values are to be returned.
* @return VectorXd Vector with values as doubles. * @return VectorXd Vector with values as doubles.
*/ */
VectorXd getBoundarySideValues(BC_SIDE side); VectorXd getBoundarySideValues(BC_SIDE side);
/** /**
* @brief Returns the boundary condition of a specified element on a given side. * @brief Returns the boundary condition of a specified element on a given
* * side.
* @param side Boundary side in which the boundary condition is located. *
* @param index Index of the boundary element on the corresponding * @param side Boundary side in which the boundary condition is located.
* boundary side. Must index an element of the corresponding side. * @param index Index of the boundary element on the corresponding
* @return BoundaryElement Boundary condition as a BoundaryElement object. * boundary side. Must index an element of the corresponding
*/ * side.
BoundaryElement getBoundaryElement(BC_SIDE side, int index); * @return BoundaryElement Boundary condition as a BoundaryElement object.
*/
BoundaryElement getBoundaryElement(BC_SIDE side, int index);
/** /**
* @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED or * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED
BC_TYPE_CONSTANT. or BC_TYPE_CONSTANT.
* *
* @param side Boundary side in which the boundary condition type is located. * @param side Boundary side in which the boundary condition type is located.
* @param index Index of the boundary element on the corresponding * @param index Index of the boundary element on the corresponding
* boundary side. Must index an element of the corresponding side. * boundary side. Must index an element of the corresponding
* @return BC_TYPE Boundary Type of the corresponding boundary condition. side.
*/ * @return BC_TYPE Boundary Type of the corresponding boundary condition.
BC_TYPE getBoundaryElementType(BC_SIDE side, int index); */
BC_TYPE getBoundaryElementType(BC_SIDE side, int index);
/** /**
* @brief Returns the concentration value of a corresponding * @brief Returns the concentration value of a corresponding
* BoundaryElement object if it is a constant boundary condition. * BoundaryElement object if it is a constant boundary condition.
* *
* @param side Boundary side in which the boundary condition value is * @param side Boundary side in which the boundary condition value is
* located. * located.
* @param index Index of the boundary element on the corresponding * @param index Index of the boundary element on the corresponding
* boundary side. Must index an element of the corresponding * boundary side. Must index an element of the corresponding
* side. * side.
* @return double Concentration of the corresponding BoundaryElement object. * @return double Concentration of the corresponding BoundaryElement object.
*/ */
double getBoundaryElementValue(BC_SIDE side, int index); double getBoundaryElementValue(BC_SIDE side, int index);
private: private:
Grid grid; // Boundary is directly dependent on the dimensions of a predefined Grid grid; // Boundary is directly dependent on the dimensions of a predefined
vector<vector<BoundaryElement>> boundaries; // Vector with Boundary Element information vector<vector<BoundaryElement>>
boundaries; // Vector with Boundary Element information
}; };
#endif #endif

293
include/tug/Grid.hpp
View File

@ -1,8 +1,8 @@
/** /**
* @file Grid.hpp * @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a * @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients. * respective matrix/matrices of alpha coefficients.
* *
*/ */
#include <Eigen/Core> #include <Eigen/Core>
@ -11,165 +11,166 @@
using namespace Eigen; using namespace Eigen;
class Grid { class Grid {
public: public:
/**
* @brief Constructs a new 1D-Grid object of a given length, which holds a
* matrix with concentrations and a respective matrix of alpha coefficients.
* The domain length is per default the same as the length. The
* concentrations are all 20 by default and the alpha coefficients are 1.
*
* @param length Length of the 1D-Grid. Must be greater than 3.
*/
Grid(int length);
/** /**
* @brief Constructs a new 1D-Grid object of a given length, which holds a matrix * @brief Constructs a new 2D-Grid object of given dimensions, which holds a
* with concentrations and a respective matrix of alpha coefficients. * matrix with concentrations and the respective matrices of alpha coefficient
* The domain length is per default the same as the length. The concentrations * for each direction. The domain in x- and y-direction is per default equal
* are all 20 by default and the alpha coefficients are 1. * 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
* @param length Length of the 1D-Grid. Must be greater than 3. * directions.
*/ *
Grid(int length); * @param row Length of the 2D-Grid in y-direction. Must be greater than 3.
* @param col Length of the 2D-Grid in x-direction. Must be greater than 3.
*/
Grid(int row, int col);
/** /**
* @brief Constructs a new 2D-Grid object of given dimensions, which holds a matrix * @brief Sets the concentrations matrix for a 1D or 2D-Grid.
* with concentrations and the respective matrices of alpha coefficient for *
* each direction. The domain in x- and y-direction is per default equal to * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix
* the col length and row length, respectively. * must have correct dimensions as defined in row and col. (Or length, in 1D
* The concentrations are all 20 by default across the entire grid and the * case).
* alpha coefficients 1 in both directions. */
* void setConcentrations(MatrixXd concentrations);
* @param row Length of the 2D-Grid in y-direction. Must be greater than 3.
* @param col Length of the 2D-Grid in x-direction. Must be greater than 3.
*/
Grid(int row, int col);
/** /**
* @brief Sets the concentrations matrix for a 1D or 2D-Grid. * @brief Gets the concentrations matrix for a Grid.
* *
* @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix must * @return MatrixXd An Eigen3 matrix holding the concentrations and having the
* have correct dimensions as defined in row and col. (Or length, * same dimensions as the grid.
* in 1D case). */
*/ const MatrixXd getConcentrations();
void setConcentrations(MatrixXd concentrations);
/** /**
* @brief Gets the concentrations matrix for a Grid. * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
* * dimensional.
* @return MatrixXd An Eigen3 matrix holding the concentrations and having the *
* same dimensions as the grid. * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
*/ * Matrix columns must have same size as length of grid.
const MatrixXd getConcentrations(); */
void setAlpha(MatrixXd alpha);
/** /**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one dimensional. * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* * dimensional.
* @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. Matrix *
* columns must have same size as length of grid. * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in
*/ * x-direction. Matrix must be of same size as the grid.
void setAlpha(MatrixXd alpha); * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid.
*/
void setAlpha(MatrixXd alphaX, MatrixXd alphaY);
/** /**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two dimensional. * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
* * dimensional.
* @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in x-direction. *
* Matrix must be of same size as the grid. * @return MatrixXd A matrix with 1 row holding the alpha coefficients.
* @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in y-direction. */
* Matrix must be of same size as the grid. const MatrixXd getAlpha();
*/
void setAlpha(MatrixXd alphaX, MatrixXd alphaY);
/** /**
* @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one dimensional. * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
* * Grid must be two dimensional.
* @return MatrixXd A matrix with 1 row holding the alpha coefficients. *
*/ * @return MatrixXd A matrix holding the alpha coefficients in x-direction.
const MatrixXd getAlpha(); */
const MatrixXd getAlphaX();
/** /**
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. Grid must be * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
* two dimensional. * Grid must be two dimensional.
* *
* @return MatrixXd A matrix holding the alpha coefficients in x-direction. * @return MatrixXd A matrix holding the alpha coefficients in y-direction.
*/ */
const MatrixXd getAlphaX(); const MatrixXd getAlphaY();
/** /**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. Grid must be * @brief Gets the dimensions of the grid.
* two dimensional. *
* * @return int Dimensions, either 1 or 2.
* @return MatrixXd A matrix holding the alpha coefficients in y-direction. */
*/ int getDim();
const MatrixXd getAlphaY();
/** /**
* @brief Gets the dimensions of the grid. * @brief Gets length of 1D grid. Must be one dimensional grid.
* *
* @return int Dimensions, either 1 or 2. * @return int Length of 1D grid.
*/ */
int getDim(); int getLength();
/** /**
* @brief Gets length of 1D grid. Must be one dimensional grid. * @brief Gets the number of rows of the grid.
* *
* @return int Length of 1D grid. * @return int Number of rows.
*/ */
int getLength(); int getRow();
/** /**
* @brief Gets the number of rows of the grid. * @brief Gets the number of columns of the grid.
* *
* @return int Number of rows. * @return int Number of columns.
*/ */
int getRow(); int getCol();
/** /**
* @brief Gets the number of columns of the grid. * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
* *
* @return int Number of columns. * @param domainLength A double value of the domain length. Must be positive.
*/ */
int getCol(); void setDomain(double domainLength);
/** /**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
* *
* @param domainLength A double value of the domain length. Must be positive. * @param domainRow A double value of the domain size in y-direction. Must be
*/ * positive.
void setDomain(double domainLength); * @param domainCol A double value of the domain size in x-direction. Must be
* positive.
*/
void setDomain(double domainRow, double domainCol);
/** /**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
* *
* @param domainRow A double value of the domain size in y-direction. Must be positive. * @return double Delta value.
* @param domainCol A double value of the domain size in x-direction. Must be positive. */
*/ double getDelta();
void setDomain(double domainRow,double domainCol);
/** /**
* @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * @brief Gets the delta value in x-direction.
* *
* @return double Delta value. * @return double Delta value in x-direction.
*/ */
double getDelta(); double getDeltaCol();
/** /**
* @brief Gets the delta value in x-direction. * @brief Gets the delta value in y-direction. Must be two dimensional grid.
* *
* @return double Delta value in x-direction. * @return double Delta value in y-direction.
*/ */
double getDeltaCol(); double getDeltaRow();
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return double Delta value in y-direction.
*/
double getDeltaRow();
private:
int col; // number of grid columns
int row; // number of grid rows
int dim; // 1D or 2D
double domainCol; // number of domain columns
double domainRow; // number of domain rows
double deltaCol; // delta in x-direction (between columns)
double deltaRow; // delta in y-direction (between rows)
MatrixXd concentrations; // Matrix holding grid concentrations
MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
private:
int col; // number of grid columns
int row; // number of grid rows
int dim; // 1D or 2D
double domainCol; // number of domain columns
double domainRow; // number of domain rows
double deltaCol; // delta in x-direction (between columns)
double deltaRow; // delta in y-direction (between rows)
MatrixXd concentrations; // Matrix holding grid concentrations
MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
}; };

315
include/tug/Simulation.hpp
View File

@ -1,8 +1,8 @@
/** /**
* @file Simulation.hpp * @file Simulation.hpp
* @brief API of Simulation class, that holds all information regarding a specific simulation * @brief API of Simulation class, that holds all information regarding a
* run like its timestep, number of iterations and output options. Simulation object * specific simulation run like its timestep, number of iterations and output
* also holds a predefined Grid and Boundary object. * options. Simulation object also holds a predefined Grid and Boundary object.
* *
*/ */
#include "Boundary.hpp" #include "Boundary.hpp"
@ -11,52 +11,54 @@
using namespace std; using namespace std;
/** /**
* @brief Enum defining the two implemented solution approaches. * @brief Enum defining the two implemented solution approaches.
* *
*/ */
enum APPROACH { enum APPROACH {
FTCS_APPROACH, // Forward Time-Centered Space FTCS_APPROACH, // Forward Time-Centered Space
BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver
CRANK_NICOLSON_APPROACH CRANK_NICOLSON_APPROACH
}; };
/** /**
* @brief Enum defining the Linear Equation solvers * @brief Enum defining the Linear Equation solvers
* *
*/ */
enum SOLVER { enum SOLVER {
EIGEN_LU_SOLVER, // EigenLU solver EIGEN_LU_SOLVER, // EigenLU solver
THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for tridiagonal matrices THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for
// tridiagonal matrices
}; };
/** /**
* @brief Enum holding different options for .csv output. * @brief Enum holding different options for .csv output.
* *
*/ */
enum CSV_OUTPUT { enum CSV_OUTPUT {
CSV_OUTPUT_OFF, // do not produce csv output CSV_OUTPUT_OFF, // do not produce csv output
CSV_OUTPUT_ON, // produce csv output with last concentration matrix CSV_OUTPUT_ON, // produce csv output with last concentration matrix
CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices
CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary conditions at beginning CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary
// conditions at beginning
}; };
/** /**
* @brief Enum holding different options for console output. * @brief Enum holding different options for console output.
* *
*/ */
enum CONSOLE_OUTPUT { enum CONSOLE_OUTPUT {
CONSOLE_OUTPUT_OFF, // do not print any output to console CONSOLE_OUTPUT_OFF, // do not print any output to console
CONSOLE_OUTPUT_ON, // print before and after concentrations to console CONSOLE_OUTPUT_ON, // print before and after concentrations to console
CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console
}; };
/** /**
* @brief Enum holding different options for time measurement. * @brief Enum holding different options for time measurement.
* *
*/ */
enum TIME_MEASURE { enum TIME_MEASURE {
TIME_MEASURE_OFF, // do not print any time measures TIME_MEASURE_OFF, // do not print any time measures
TIME_MEASURE_ON // print time measure after last iteration TIME_MEASURE_ON // print time measure after last iteration
}; };
/** /**
@ -66,148 +68,155 @@ enum TIME_MEASURE {
* *
*/ */
class Simulation { class Simulation {
public: public:
/** /**
* @brief Set up a simulation environment. The timestep and number of iterations * @brief Set up a simulation environment. The timestep and number of
* must be set. For the BTCS approach, the Thomas algorithm is used as * iterations must be set. For the BTCS approach, the Thomas algorithm is used
* the default linear equation solver as this is faster for tridiagonal * as the default linear equation solver as this is faster for tridiagonal
* matrices. CSV output, console output and time measure are off by default. * matrices. CSV output, console output and time measure are off by
* Also, the number of cores is set to the maximum number of cores -1 by default. * default. Also, the number of cores is set to the maximum number of cores -1
* * by default.
* @param grid Valid grid object *
* @param bc Valid boundary condition object * @param grid Valid grid object
* @param approach Approach to solving the problem. Either FTCS or BTCS. * @param bc Valid boundary condition object
*/ * @param approach Approach to solving the problem. Either FTCS or BTCS.
Simulation(Grid &grid, Boundary &bc, APPROACH approach); */
Simulation(Grid &grid, Boundary &bc, APPROACH approach);
/** /**
* @brief Set the option to output the results to a CSV file. Off by default. * @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 * @param csv_output Valid output option. The following options can be set
* here: * here:
* - CSV_OUTPUT_OFF: do not produce csv output * - CSV_OUTPUT_OFF: do not produce csv output
* - CSV_OUTPUT_ON: produce csv output with last * - CSV_OUTPUT_ON: produce csv output with last
* concentration matrix * concentration matrix
* - CSV_OUTPUT_VERBOSE: produce csv output with all * - CSV_OUTPUT_VERBOSE: produce csv output with all
* concentration matrices * concentration matrices
* - CSV_OUTPUT_XTREME: produce csv output with all * - CSV_OUTPUT_XTREME: produce csv output with all
* concentration matrices and simulation environment * concentration matrices and simulation environment
*/ */
void setOutputCSV(CSV_OUTPUT csv_output); void setOutputCSV(CSV_OUTPUT csv_output);
/** /**
* @brief Set the options for outputting information to the console. Off by default. * @brief Set the options for outputting information to the console. Off by
* * default.
* @param console_output Valid output option. The following options can be set *
* here: * @param console_output Valid output option. The following options can be set
* - CONSOLE_OUTPUT_OFF: do not print any output to console * here:
* - CONSOLE_OUTPUT_ON: print before and after concentrations to console * - CONSOLE_OUTPUT_OFF: do not print any output to
* - CONSOLE_OUTPUT_VERBOSE: print all concentration matrices to console * console
*/ * - CONSOLE_OUTPUT_ON: print before and after
void setOutputConsole(CONSOLE_OUTPUT console_output); * concentrations to console
* - CONSOLE_OUTPUT_VERBOSE: print all concentration
* matrices to console
*/
void setOutputConsole(CONSOLE_OUTPUT console_output);
/** /**
* @brief Set the Time Measure option. Off by default. * @brief Set the Time Measure option. Off by default.
* *
* @param time_measure The following options are allowed: * @param time_measure The following options are allowed:
* - TIME_MEASURE_OFF: Time of simulation is not printed to console * - TIME_MEASURE_OFF: Time of simulation is not printed
* - TIME_MEASURE_ON: Time of simulation run is printed to console * to console
*/ * - TIME_MEASURE_ON: Time of simulation run is printed to
void setTimeMeasure(TIME_MEASURE time_measure); * console
*/
void setTimeMeasure(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
* greater than zero. Setting the timestep is required. * greater than zero. Setting the timestep is required.
* *
* @param timestep Valid timestep greater than zero. * @param timestep Valid timestep greater than zero.
*/ */
void setTimestep(double timestep); void setTimestep(double timestep);
/** /**
* @brief Currently set time step is returned. * @brief Currently set time step is returned.
* *
* @return double timestep * @return double timestep
*/ */
double getTimestep(); double getTimestep();
/** /**
* @brief Set the desired iterations to be calculated. A value greater * @brief Set the desired iterations to be calculated. A value greater
* than zero must be specified here. Setting iterations is required. * than zero must be specified here. Setting iterations is required.
* *
* @param iterations Number of iterations to be simulated. * @param iterations Number of iterations to be simulated.
*/ */
void setIterations(int iterations); void setIterations(int iterations);
/** /**
* @brief Set the desired linear equation solver to be used for BTCS approach. Without effect * @brief Set the desired linear equation solver to be used for BTCS approach.
* in case of FTCS approach. * Without effect in case of FTCS approach.
* *
* @param solver Solver to be used. Default is Thomas Algorithm as it is more efficient for * @param solver Solver to be used. Default is Thomas Algorithm as it is more
* tridiagonal Matrices. * efficient for tridiagonal Matrices.
*/ */
void setSolver(SOLVER solver); void setSolver(SOLVER solver);
/** /**
* @brief Set the number of desired openMP Threads. * @brief Set the number of desired openMP Threads.
* *
* @param num_threads Number of desired threads. Must have a value between * @param num_threads Number of desired threads. Must have a value between
* 1 and the maximum available number of processors. The maximum number of * 1 and the maximum available number of processors. The
* processors is set as the default case during Simulation construction. * maximum number of processors is set as the default case during Simulation
*/ * construction.
void setNumberThreads(int num_threads); */
void setNumberThreads(int num_threads);
/** /**
* @brief Return the currently set iterations to be calculated. * @brief Return the currently set iterations to be calculated.
* *
* @return int Number of iterations. * @return int Number of iterations.
*/ */
int getIterations(); int getIterations();
/** /**
* @brief Outputs the current concentrations of the grid on the console. * @brief Outputs the current concentrations of the grid on the console.
* *
*/ */
void printConcentrationsConsole(); void printConcentrationsConsole();
/** /**
* @brief Creates a CSV file with a name containing the current simulation * @brief Creates a CSV file with a name containing the current simulation
* parameters. If the data name already exists, an additional counter is * parameters. If the data name already exists, an additional counter
* appended to the name. The name of the file is built up as follows: * is appended to the name. The name of the file is built up as follows:
* <approach> + <number rows> + <number columns> + <number of iterations>+<counter>.csv * <approach> + <number rows> + <number columns> + <number of
* * iterations>+<counter>.csv
* @return string Filename with configured simulation parameters. *
*/ * @return string Filename with configured simulation parameters.
string createCSVfile(); */
string createCSVfile();
/** /**
* @brief Writes the currently calculated concentration values of the grid * @brief Writes the currently calculated concentration values of the grid
* into the CSV file with the passed filename. * into the CSV file with the passed filename.
* *
* @param filename Name of the file to which the concentration values are * @param filename Name of the file to which the concentration values are
* to be written. * to be written.
*/ */
void printConcentrationsCSV(string filename); void printConcentrationsCSV(string filename);
/** /**
* @brief Method starts the simulation process with the previously set * @brief Method starts the simulation process with the previously set
* parameters. * parameters.
*/ */
void run(); void run();
private: private:
double timestep;
double timestep; int iterations;
int iterations; int innerIterations;
int innerIterations; int numThreads;
int numThreads; CSV_OUTPUT csv_output;
CSV_OUTPUT csv_output; CONSOLE_OUTPUT console_output;
CONSOLE_OUTPUT console_output; TIME_MEASURE time_measure;
TIME_MEASURE time_measure;
Grid &grid;
Boundary &bc;
APPROACH approach;
SOLVER solver;
Grid &grid;
Boundary &bc;
APPROACH approach;
SOLVER solver;
}; };

658
src/BTCSv2.cpp
View File

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

210
src/Boundary.cpp
View File

@ -1,162 +1,162 @@
#include "TugUtils.cpp" #include "TugUtils.cpp"
#include <iostream> #include <iostream>
#include <omp.h> #include <omp.h>
#include <tug/Boundary.hpp>
#include <stdexcept> #include <stdexcept>
#include <tug/Boundary.hpp>
using namespace std; using namespace std;
BoundaryElement::BoundaryElement() { BoundaryElement::BoundaryElement() {
this->type = BC_TYPE_CLOSED; this->type = BC_TYPE_CLOSED;
this->value = -1; // without meaning in closed case this->value = -1; // without meaning in closed case
} }
BoundaryElement::BoundaryElement(double value) { BoundaryElement::BoundaryElement(double value) {
this->type = BC_TYPE_CONSTANT; this->type = BC_TYPE_CONSTANT;
this->value = value; this->value = value;
} }
void BoundaryElement::setType(BC_TYPE type) { void BoundaryElement::setType(BC_TYPE type) { this->type = type; }
this->type = type;
}
void BoundaryElement::setValue(double value) { void BoundaryElement::setValue(double value) {
if(value < 0){ if (value < 0) {
throw_invalid_argument("No negative concentration allowed."); throw_invalid_argument("No negative concentration allowed.");
} }
if(type == BC_TYPE_CLOSED){ if (type == BC_TYPE_CLOSED) {
throw_invalid_argument( throw_invalid_argument(
"No constant boundary concentrations can be set for closed " "No constant boundary concentrations can be set for closed "
"boundaries. Please change type first."); "boundaries. Please change type first.");
} }
this->value = value; this->value = value;
} }
BC_TYPE BoundaryElement::getType() { BC_TYPE BoundaryElement::getType() { return this->type; }
return this->type;
}
double BoundaryElement::getValue() { double BoundaryElement::getValue() { return this->value; }
return this->value;
}
Boundary::Boundary(Grid grid) : grid(grid) { Boundary::Boundary(Grid grid) : grid(grid) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
this->boundaries = vector<vector<BoundaryElement>>(2); // in 1D only left and right boundary this->boundaries = vector<vector<BoundaryElement>>(
2); // in 1D only left and right boundary
this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {
this->boundaries = vector<vector<BoundaryElement>>(4); this->boundaries = vector<vector<BoundaryElement>>(4);
this->boundaries[BC_SIDE_LEFT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement()); this->boundaries[BC_SIDE_LEFT] =
this->boundaries[BC_SIDE_RIGHT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement()); vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_TOP] = vector<BoundaryElement>(grid.getCol(), BoundaryElement()); this->boundaries[BC_SIDE_RIGHT] =
this->boundaries[BC_SIDE_BOTTOM] = vector<BoundaryElement>(grid.getCol(), BoundaryElement()); vector<BoundaryElement>(grid.getRow(), BoundaryElement());
} this->boundaries[BC_SIDE_TOP] =
vector<BoundaryElement>(grid.getCol(), BoundaryElement());
this->boundaries[BC_SIDE_BOTTOM] =
vector<BoundaryElement>(grid.getCol(), BoundaryElement());
}
} }
void Boundary::setBoundarySideClosed(BC_SIDE side) { void Boundary::setBoundarySideClosed(BC_SIDE side) {
if(grid.getDim() == 1){ if (grid.getDim() == 1) {
if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument( throw_invalid_argument(
"For the one-dimensional case, only the BC_SIDE_LEFT and " "For the one-dimensional case, only the BC_SIDE_LEFT and "
"BC_SIDE_RIGHT borders exist."); "BC_SIDE_RIGHT borders exist.");
}
} }
}
int n; int n;
if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) {
n = grid.getRow(); n = grid.getRow();
} else { } else {
n = grid.getCol(); n = grid.getCol();
} }
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement()); this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement());
} }
void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
if(grid.getDim() == 1){ if (grid.getDim() == 1) {
if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument( throw_invalid_argument(
"For the one-dimensional case, only the BC_SIDE_LEFT and " "For the one-dimensional case, only the BC_SIDE_LEFT and "
"BC_SIDE_RIGHT borders exist."); "BC_SIDE_RIGHT borders exist.");
}
} }
}
int n; int n;
if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) {
n = grid.getRow(); n = grid.getRow();
} else { } else {
n = grid.getCol(); n = grid.getCol();
} }
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement(value)); this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement(value));
} }
void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { void Boundary::setBoundaryElementClosed(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){ if ((boundaries[side].size() < index) || index < 0) {
throw_invalid_argument("Index is selected either too large or too small."); throw_invalid_argument("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);
} }
void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) { void Boundary::setBoundaryElementConstant(BC_SIDE side, int index,
// tests whether the index really points to an element of the boundary side. double value) {
if((boundaries[side].size() < index) || index < 0){ // tests whether the index really points to an element of the boundary side.
throw_invalid_argument("Index is selected either too large or too small."); if ((boundaries[side].size() < index) || index < 0) {
} throw_invalid_argument("Index is selected either too large or too small.");
this->boundaries[side][index].setType(BC_TYPE_CONSTANT); }
this->boundaries[side][index].setValue(value); this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
this->boundaries[side][index].setValue(value);
} }
const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) { const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
if(grid.getDim() == 1){ if (grid.getDim() == 1) {
if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument( throw_invalid_argument(
"For the one-dimensional trap, only the BC_SIDE_LEFT and " "For the one-dimensional trap, only the BC_SIDE_LEFT and "
"BC_SIDE_RIGHT borders exist."); "BC_SIDE_RIGHT borders exist.");
}
} }
return this->boundaries[side]; }
return this->boundaries[side];
} }
VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
int length = boundaries[side].size(); int length = boundaries[side].size();
VectorXd values(length); VectorXd values(length);
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {
values(i) = -1; values(i) = -1;
continue; continue;
}
values(i) = getBoundaryElementValue(side, i);
} }
values(i) = getBoundaryElementValue(side, i);
}
return values; return values;
} }
BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) {
if((boundaries[side].size() < index) || index < 0){ if ((boundaries[side].size() < index) || index < 0) {
throw_invalid_argument("Index is selected either too large or too small."); throw_invalid_argument("Index is selected either too large or too small.");
} }
return this->boundaries[side][index]; return this->boundaries[side][index];
} }
BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) {
if((boundaries[side].size() < index) || index < 0){ if ((boundaries[side].size() < index) || index < 0) {
throw_invalid_argument("Index is selected either too large or too small."); throw_invalid_argument("Index is selected either too large or too small.");
} }
return this->boundaries[side][index].getType(); return this->boundaries[side][index].getType();
} }
double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { double Boundary::getBoundaryElementValue(BC_SIDE side, int index) {
if((boundaries[side].size() < index) || index < 0){ if ((boundaries[side].size() < index) || index < 0) {
throw_invalid_argument("Index is selected either too large or too small."); throw_invalid_argument("Index is selected either too large or too small.");
} }
if(boundaries[side][index].getType() != BC_TYPE_CONSTANT){ if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) {
throw_invalid_argument("A value can only be output if it is a constant boundary condition."); throw_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();
} }

646
src/FTCS.cpp
View File

@ -1,433 +1,389 @@
/** /**
* @file FTCS.cpp * @file FTCS.cpp
* @brief Implementation of heterogenous FTCS (forward time-centered space) solution * @brief Implementation of heterogenous FTCS (forward time-centered space)
* of diffusion equation in 1D and 2D space. * solution of diffusion equation in 1D and 2D space.
* *
*/ */
#include "TugUtils.cpp" #include "TugUtils.cpp"
#include <cstddef> #include <cstddef>
#include <tug/Boundary.hpp>
#include <iostream> #include <iostream>
#include <omp.h> #include <omp.h>
#include <tug/Boundary.hpp>
using namespace std; using namespace std;
// calculates horizontal change on one cell independent of boundary type // calculates horizontal change on one cell independent of boundary type
static double calcHorizontalChange(Grid &grid, int &row, int &col) { static double calcHorizontalChange(Grid &grid, int &row, int &col) {
double result = double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) *
* grid.getConcentrations()(row,col+1) grid.getConcentrations()(row, col + 1) -
- ( (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) +
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
) grid.getAlphaX()(row, col))) *
* grid.getConcentrations()(row,col) grid.getConcentrations()(row, col) +
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
* grid.getConcentrations()(row,col-1); grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1);
return result; return result;
} }
// calculates vertical change on one cell independent of boundary type // calculates vertical change on one cell independent of boundary type
static double calcVerticalChange(Grid &grid, int &row, int &col) { static double calcVerticalChange(Grid &grid, int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row+1,col)
- (
calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+ calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
)
* grid.getConcentrations()(row,col)
+ calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
* grid.getConcentrations()(row-1,col);
return result; double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col))) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row - 1, col);
return result;
} }
// calculates horizontal change on one cell with a constant left boundary // calculates horizontal change on one cell with a constant left boundary
static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc,
int &row, int &col) {
double result = double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) *
* grid.getConcentrations()(row,col+1) grid.getConcentrations()(row, col + 1) -
- ( (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) +
+ 2 * grid.getAlphaX()(row,col) 2 * grid.getAlphaX()(row, col)) *
) grid.getConcentrations()(row, col) +
* grid.getConcentrations()(row,col) 2 * grid.getAlphaX()(row, col) *
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row); bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
return result; return result;
} }
// calculates horizontal change on one cell with a closed left boundary // calculates horizontal change on one cell with a closed left boundary
static double calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { static double calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row,
int &col) {
double result =
calcAlphaIntercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col))
* (grid.getConcentrations()(row, col+1) - grid.getConcentrations()(row, col));
return result;
}
double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col + 1) -
grid.getConcentrations()(row, col));
return result;
}
// 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
static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, int &row, int &col) { static double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc,
if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { int &row, int &col) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) { return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
return calcHorizontalChangeLeftBoundaryClosed(grid, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) {
} else { return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
throw_invalid_argument("Undefined Boundary Condition Type!"); } else {
} throw_invalid_argument("Undefined Boundary Condition Type!");
}
} }
// calculates horizontal change on one cell with a constant right boundary // calculates horizontal change on one cell with a constant right boundary
static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { static double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
Boundary &bc, int &row,
int &col) {
double result = double result = 2 * grid.getAlphaX()(row, col) *
2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
- ( (calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) +
+ 2 * grid.getAlphaX()(row,col) 2 * grid.getAlphaX()(row, col)) *
) grid.getConcentrations()(row, col) +
* grid.getConcentrations()(row,col) calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
+ calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) grid.getAlphaX()(row, col)) *
* grid.getConcentrations()(row,col-1); grid.getConcentrations()(row, col - 1);
return result; return result;
} }
// calculates horizontal change on one cell with a closed right boundary // calculates horizontal change on one cell with a closed right boundary
static double calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { static double calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row,
int &col) {
double result =
- (calcAlphaIntercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row, col-1)));
return result;
}
double result = -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col) -
grid.getConcentrations()(row, col - 1)));
return result;
}
// 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
static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, int &row, int &col) { static double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc,
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { int &row, int &col) {
return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
} else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) { return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
return calcHorizontalChangeRightBoundaryClosed(grid, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) {
} else { return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
throw_invalid_argument("Undefined Boundary Condition Type!"); } else {
} throw_invalid_argument("Undefined Boundary Condition Type!");
}
} }
// calculates vertical change on one cell with a constant top boundary // calculates vertical change on one cell with a constant top boundary
static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc,
int &row, int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col))
* grid.getConcentrations()(row+1,col)
- (
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col))
+ 2 * grid.getAlphaY()(row, col)
)
* grid.getConcentrations()(row, col)
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
return result; double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
2 * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row, col) +
2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_TOP, col);
return result;
} }
// calculates vertical change on one cell with a closed top boundary // calculates vertical change on one cell with a closed top boundary
static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, int &col) { static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
int &col) {
double result =
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getConcentrations()(row, col))
* (grid.getConcentrations()(row+1, col) - grid.getConcentrations()(row, col));
return result; double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getConcentrations()(row, col)) *
(grid.getConcentrations()(row + 1, col) -
grid.getConcentrations()(row, col));
return result;
} }
// 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
static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row, int &col) { static double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, int &row,
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { int &col) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col); if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
return calcVerticalChangeTopBoundaryClosed(grid, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
} else { return calcVerticalChangeTopBoundaryClosed(grid, row, col);
throw_invalid_argument("Undefined Boundary Condition Type!"); } else {
} throw_invalid_argument("Undefined Boundary Condition Type!");
}
} }
// calculates vertical change on one cell with a constant bottom boundary // calculates vertical change on one cell with a constant bottom boundary
static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc,
int &row, int &col) {
double result = double result = 2 * grid.getAlphaY()(row, col) *
2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
- ( (calcAlphaIntercell(grid.getAlphaY()(row, col),
calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) grid.getAlphaY()(row - 1, col)) +
+ 2 * grid.getAlphaY()(row, col) 2 * grid.getAlphaY()(row, col)) *
) grid.getConcentrations()(row, col) +
* grid.getConcentrations()(row, col) calcAlphaIntercell(grid.getAlphaY()(row, col),
+ calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) grid.getAlphaY()(row - 1, col)) *
* grid.getConcentrations()(row-1,col); grid.getConcentrations()(row - 1, col);
return result; return result;
} }
// calculates vertical change on one cell with a closed bottom boundary // calculates vertical change on one cell with a closed bottom boundary
static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row,
int &col) {
double result = double result = -(calcAlphaIntercell(grid.getAlphaY()(row, col),
- (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) grid.getAlphaY()(row - 1, col)) *
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col))); (grid.getConcentrations()(row, col) -
grid.getConcentrations()(row - 1, col)));
return result; return result;
} }
// 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
static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, int &row, int &col) { static double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc,
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { int &row, int &col) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col); if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) { return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
return calcVerticalChangeBottomBoundaryClosed(grid, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
} else { return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
throw_invalid_argument("Undefined Boundary Condition Type!"); } else {
} throw_invalid_argument("Undefined Boundary Condition Type!");
}
} }
// FTCS solution for 1D grid // FTCS solution for 1D grid
static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) { static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
int colMax = grid.getCol(); int colMax = grid.getCol();
double deltaCol = grid.getDeltaCol(); double deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0);
// only one row in 1D case -> row constant at index 0 // only one row in 1D case -> row constant at index 0
int row = 0; int row = 0;
// 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) = grid.getConcentrations()(row, col) +
+ timestep / (deltaCol*deltaCol) timestep / (deltaCol * deltaCol) *
* ( (calcHorizontalChange(grid, row, col));
calcHorizontalChange(grid, row, col) }
)
;
}
// 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) = grid.getConcentrations()(row,col) concentrations_t1(row, col) =
+ timestep / (deltaCol*deltaCol) grid.getConcentrations()(row, col) +
* ( timestep / (deltaCol * deltaCol) *
calcHorizontalChangeLeftBoundary(grid, bc, row, col) (calcHorizontalChangeLeftBoundary(grid, bc, row, col));
)
;
// right boundary; hold column constant at max index
col = colMax - 1;
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col));
// right boundary; hold column constant at max index // overwrite obsolete concentrations
col = colMax-1; grid.setConcentrations(concentrations_t1);
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
;
// overwrite obsolete concentrations
grid.setConcentrations(concentrations_t1);
} }
// FTCS solution for 2D grid // FTCS solution for 2D grid
static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep, int numThreads) { static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
int rowMax = grid.getRow(); int numThreads) {
int colMax = grid.getCol(); int rowMax = grid.getRow();
double deltaRow = grid.getDeltaRow(); int colMax = grid.getCol();
double deltaCol = grid.getDeltaCol(); double deltaRow = grid.getDeltaRow();
double deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
// inner cells
// these are independent of the boundary condition type
// omp_set_num_threads(10);
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax-1; row++) {
for (int col = 1; col < colMax-1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChange(grid, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChange(grid, row, col)
)
;
}
}
// boundary conditions
// left without corners / looping over rows
// hold column constant at index 0
int col = 0;
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChange(grid, row, col)
)
;
}
// right without corners / looping over rows
// hold column constant at max index
col = colMax-1;
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChange(grid, row, col)
)
;
}
// top without corners / looping over columns
// hold row constant at index 0
int row = 0;
#pragma omp parallel for num_threads(numThreads)
for (int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChange(grid, row, col)
)
;
}
// bottom without corners / looping over columns
// hold row constant at max index
row = rowMax-1;
#pragma omp parallel for num_threads(numThreads)
for(int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChange(grid, row, col)
)
;
}
// corner top left
// hold row and column constant at 0
row = 0;
col = 0;
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
// corner top right
// hold row constant at 0 and column constant at max index
row = 0;
col = colMax-1;
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
// corner bottom left // inner cells
// hold row constant at max index and column constant at 0 // these are independent of the boundary condition type
row = rowMax-1; // omp_set_num_threads(10);
col = 0; #pragma omp parallel for num_threads(numThreads)
concentrations_t1(row,col) = grid.getConcentrations()(row,col) for (int row = 1; row < rowMax - 1; row++) {
+ timestep/(deltaCol*deltaCol) for (int col = 1; col < colMax - 1; col++) {
* ( concentrations_t1(row, col) = grid.getConcentrations()(row, col) +
calcHorizontalChangeLeftBoundary(grid, bc, row, col) timestep / (deltaRow * deltaRow) *
) (calcVerticalChange(grid, row, col)) +
+ timestep/(deltaRow*deltaRow) timestep / (deltaCol * deltaCol) *
* ( (calcHorizontalChange(grid, row, col));
calcVerticalChangeBottomBoundary(grid, bc, row, col) }
) }
;
// corner bottom right // boundary conditions
// hold row and column constant at max index // left without corners / looping over rows
row = rowMax-1; // hold column constant at index 0
col = colMax-1; int col = 0;
concentrations_t1(row,col) = grid.getConcentrations()(row,col) #pragma omp parallel for num_threads(numThreads)
+ timestep/(deltaCol*deltaCol) for (int row = 1; row < rowMax - 1; row++) {
* ( concentrations_t1(row, col) =
calcHorizontalChangeRightBoundary(grid, bc, row, col) grid.getConcentrations()(row, col) +
) timestep / (deltaCol * deltaCol) *
+ timestep/(deltaRow*deltaRow) (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
* ( timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
calcVerticalChangeBottomBoundary(grid, bc, row, col) }
)
;
// overwrite obsolete concentrations // right without corners / looping over rows
grid.setConcentrations(concentrations_t1); // hold column constant at max index
// } col = colMax - 1;
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) {
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
}
// top without corners / looping over columns
// hold row constant at index 0
int row = 0;
#pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
}
// bottom without corners / looping over columns
// hold row constant at max index
row = rowMax - 1;
#pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
}
// corner top left
// hold row and column constant at 0
row = 0;
col = 0;
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col));
// corner top right
// hold row constant at 0 and column constant at max index
row = 0;
col = colMax - 1;
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col));
// corner bottom left
// hold row constant at max index and column constant at 0
row = rowMax - 1;
col = 0;
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col));
// corner bottom right
// hold row and column constant at max index
row = rowMax - 1;
col = colMax - 1;
concentrations_t1(row, col) =
grid.getConcentrations()(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col));
// overwrite obsolete concentrations
grid.setConcentrations(concentrations_t1);
// }
} }
// entry point; differentiate between 1D and 2D grid // entry point; differentiate between 1D and 2D grid
static void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) { static void FTCS(Grid &grid, Boundary &bc, double &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) {
FTCS_2D(grid, bc, timestep, numThreads); FTCS_2D(grid, bc, timestep, numThreads);
} else { } else {
throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); throw_invalid_argument(
} "Error: Only 1- and 2-dimensional grids are defined!");
}
} }

213
src/Grid.cpp
View File

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

511
src/Simulation.cpp
View File

@ -1,10 +1,10 @@
#include <cmath> #include <cmath>
#include <cstddef> #include <cstddef>
#include <filesystem> #include <filesystem>
#include <omp.h>
#include <stdexcept> #include <stdexcept>
#include <string> #include <string>
#include <tug/Simulation.hpp> #include <tug/Simulation.hpp>
#include <omp.h>
#include <fstream> #include <fstream>
@ -13,319 +13,330 @@
#include "BTCSv2.cpp" #include "BTCSv2.cpp"
using namespace std; using namespace std;
Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach)
: grid(grid), bc(bc) {
this->approach = approach; this->approach = approach;
this->solver = THOMAS_ALGORITHM_SOLVER; this->solver = THOMAS_ALGORITHM_SOLVER;
this->timestep = -1; // error per default; needs to be set this->timestep = -1; // error per default; needs to be set
this->iterations = -1; this->iterations = -1;
this->innerIterations = 1; this->innerIterations = 1;
this->numThreads = omp_get_num_procs()-1; this->numThreads = omp_get_num_procs() - 1;
this->csv_output = CSV_OUTPUT_OFF; this->csv_output = CSV_OUTPUT_OFF;
this->console_output = CONSOLE_OUTPUT_OFF; this->console_output = CONSOLE_OUTPUT_OFF;
this->time_measure = TIME_MEASURE_OFF; this->time_measure = TIME_MEASURE_OFF;
} }
void Simulation::setOutputCSV(CSV_OUTPUT csv_output) { void Simulation::setOutputCSV(CSV_OUTPUT csv_output) {
if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
throw_invalid_argument("Invalid CSV output option given!"); throw_invalid_argument("Invalid CSV output option given!");
} }
this->csv_output = csv_output; this->csv_output = csv_output;
} }
void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) {
if (console_output < CONSOLE_OUTPUT_OFF && console_output > CONSOLE_OUTPUT_VERBOSE) { if (console_output < CONSOLE_OUTPUT_OFF &&
throw_invalid_argument("Invalid console output option given!"); console_output > CONSOLE_OUTPUT_VERBOSE) {
} throw_invalid_argument("Invalid console output option given!");
}
this->console_output = console_output; this->console_output = console_output;
} }
void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { void Simulation::setTimeMeasure(TIME_MEASURE time_measure) {
if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
throw_invalid_argument("Invalid time measure option given!"); throw_invalid_argument("Invalid time measure option given!");
} }
this->time_measure = time_measure; this->time_measure = time_measure;
} }
void Simulation::setTimestep(double timestep) { void Simulation::setTimestep(double timestep) {
if(timestep <= 0){ if (timestep <= 0) {
throw_invalid_argument("Timestep has to be greater than zero."); throw_invalid_argument("Timestep has to be greater than zero.");
} }
if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) {
double deltaRowSquare; double deltaRowSquare;
double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
double minDeltaSquare; double minDeltaSquare;
double maxAlphaX, maxAlphaY, maxAlpha; double maxAlphaX, maxAlphaY, maxAlpha;
string dim; string dim;
if (grid.getDim() == 2) { if (grid.getDim() == 2) {
dim = "2D"; dim = "2D";
deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; minDeltaSquare =
maxAlphaX = grid.getAlphaX().maxCoeff(); (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare;
maxAlphaY = grid.getAlphaY().maxCoeff(); maxAlphaX = grid.getAlphaX().maxCoeff();
maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; maxAlphaY = grid.getAlphaY().maxCoeff();
maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
} else if (grid.getDim() == 1) { } else if (grid.getDim() == 1) {
dim = "1D"; dim = "1D";
minDeltaSquare = deltaColSquare; minDeltaSquare = deltaColSquare;
maxAlpha = grid.getAlpha().maxCoeff(); maxAlpha = grid.getAlpha().maxCoeff();
} else {
throw_invalid_argument("Critical error: Undefined number of dimensions!");
}
// Courant-Friedrichs-Lewy condition
double cfl = minDeltaSquare / (4*maxAlpha);
// stability equation from Wikipedia; might be useful if applied cfl does not work in some cases
// double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare)));
string approachPrefix = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl << endl;
cout << approachPrefix << "_" << dim << " :: required dt=" << timestep << endl;
if (timestep > cfl) {
this->innerIterations = (int)ceil(timestep / cfl);
this->timestep = timestep / (double)innerIterations;
cerr << "Warning :: Timestep was adjusted, because of stability "
"conditions. Time duration was approximately preserved by "
"adjusting internal number of iterations."
<< endl;
cout << approachPrefix << "_" << dim << " :: Required " << this->innerIterations
<< " inner iterations with dt=" << this->timestep << endl;
} else {
this->timestep = timestep;
cout << approachPrefix << "_" << dim << " :: No inner iterations required, dt=" << timestep << endl;
}
} else { } else {
this->timestep = timestep; throw_invalid_argument("Critical error: Undefined number of dimensions!");
} }
// Courant-Friedrichs-Lewy condition
double cfl = minDeltaSquare / (4 * maxAlpha);
// stability equation from Wikipedia; might be useful if applied cfl does
// not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha *
// ((1/deltaRowSquare) + (1/deltaColSquare)));
string approachPrefix =
(approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
<< endl;
cout << approachPrefix << "_" << dim << " :: required dt=" << timestep
<< endl;
if (timestep > cfl) {
this->innerIterations = (int)ceil(timestep / cfl);
this->timestep = timestep / (double)innerIterations;
cerr << "Warning :: Timestep was adjusted, because of stability "
"conditions. Time duration was approximately preserved by "
"adjusting internal number of iterations."
<< endl;
cout << approachPrefix << "_" << dim << " :: Required "
<< this->innerIterations
<< " inner iterations with dt=" << this->timestep << endl;
} else {
this->timestep = timestep;
cout << approachPrefix << "_" << dim
<< " :: No inner iterations required, dt=" << timestep << endl;
}
} else {
this->timestep = timestep;
}
} }
double Simulation::getTimestep() { double Simulation::getTimestep() { return this->timestep; }
return this->timestep;
}
void Simulation::setIterations(int iterations) { void Simulation::setIterations(int iterations) {
if(iterations <= 0){ if (iterations <= 0) {
throw_invalid_argument("Number of iterations must be greater than zero."); throw_invalid_argument("Number of iterations must be greater than zero.");
} }
this->iterations = iterations; this->iterations = iterations;
} }
void Simulation::setSolver(SOLVER solver) { void Simulation::setSolver(SOLVER solver) {
if (this->approach == FTCS_APPROACH) { if (this->approach == FTCS_APPROACH) {
cerr << "Warning: Solver was set, but FTCS approach initialized. Setting the solver " cerr << "Warning: Solver was set, but FTCS approach initialized. Setting "
"is thus without effect." "the solver "
<< endl; "is thus without effect."
} << endl;
}
this->solver = solver; this->solver = solver;
} }
void Simulation::setNumberThreads(int numThreads){ void Simulation::setNumberThreads(int numThreads) {
if(numThreads>0 && numThreads<=omp_get_num_procs()){ if (numThreads > 0 && numThreads <= omp_get_num_procs()) {
this->numThreads=numThreads; this->numThreads = numThreads;
} } else {
else{ int maxThreadNumber = omp_get_num_procs();
int maxThreadNumber = omp_get_num_procs(); string outputMessage =
string outputMessage = "Number of threads exceeds the number of processor cores (" "Number of threads exceeds the number of processor cores (" +
+ to_string(maxThreadNumber) + ") or is less than 1."; to_string(maxThreadNumber) + ") or is less than 1.";
throw_invalid_argument(outputMessage); throw_invalid_argument(outputMessage);
} }
} }
int Simulation::getIterations() { int Simulation::getIterations() { return this->iterations; }
return this->iterations;
}
void Simulation::printConcentrationsConsole() { void Simulation::printConcentrationsConsole() {
cout << grid.getConcentrations() << endl; cout << grid.getConcentrations() << endl;
cout << endl; cout << endl;
} }
string Simulation::createCSVfile() { string Simulation::createCSVfile() {
ofstream file; ofstream file;
int appendIdent = 0; int appendIdent = 0;
string appendIdentString; string appendIdentString;
// string approachString = (approach == 0) ? "FTCS" : "BTCS"; // string approachString = (approach == 0) ? "FTCS" : "BTCS";
string approachString = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); string approachString =
string row = to_string(grid.getRow()); (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
string col = to_string(grid.getCol()); string row = to_string(grid.getRow());
string numIterations = to_string(iterations); string col = to_string(grid.getCol());
string numIterations = to_string(iterations);
string filename = approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; string filename =
approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
while (filesystem::exists(filename)) { while (filesystem::exists(filename)) {
appendIdent += 1; appendIdent += 1;
appendIdentString = to_string(appendIdent); appendIdentString = to_string(appendIdent);
filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv"; filename = approachString + "_" + row + "_" + col + "_" + numIterations +
} "-" + appendIdentString + ".csv";
}
file.open(filename); file.open(filename);
if (!file) { if (!file) {
exit(1); exit(1);
} }
// adds lines at the beginning of verbose output csv that represent the boundary conditions and their values // adds lines at the beginning of verbose output csv that represent the
// -1 in case of closed boundary // boundary conditions and their values -1 in case of closed boundary
if (csv_output == CSV_OUTPUT_XTREME) { if (csv_output == CSV_OUTPUT_XTREME) {
IOFormat one_row(StreamPrecision, DontAlignCols, "", " "); IOFormat one_row(StreamPrecision, DontAlignCols, "", " ");
file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) << endl; // boundary left file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) << endl; // boundary right << endl; // boundary left
file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) << endl; // boundary top file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row)
file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) << endl; // boundary bottom << endl; // boundary right
file << endl << endl; file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row)
} << endl; // boundary top
file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row)
<< endl; // boundary bottom
file << endl << endl;
}
file.close(); file.close();
return filename; return filename;
} }
void Simulation::printConcentrationsCSV(string filename) { void Simulation::printConcentrationsCSV(string filename) {
ofstream file; ofstream file;
file.open(filename, std::ios_base::app); file.open(filename, std::ios_base::app);
if (!file) { if (!file) {
exit(1); exit(1);
} }
IOFormat do_not_align(StreamPrecision, DontAlignCols); IOFormat do_not_align(StreamPrecision, DontAlignCols);
file << grid.getConcentrations().format(do_not_align) << endl; file << grid.getConcentrations().format(do_not_align) << endl;
file << endl << endl; file << endl << endl;
file.close(); file.close();
} }
void Simulation::run() { void Simulation::run() {
if (this->timestep == -1) { if (this->timestep == -1) {
throw_invalid_argument("Timestep is not set!"); throw_invalid_argument("Timestep is not set!");
} }
if (this->iterations == -1) { if (this->iterations == -1) {
throw_invalid_argument("Number of iterations are not set!"); throw_invalid_argument("Number of iterations are not set!");
} }
string filename; string filename;
if (this->console_output > CONSOLE_OUTPUT_OFF) { if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole();
}
if (this->csv_output > CSV_OUTPUT_OFF) {
filename = createCSVfile();
}
auto begin = std::chrono::high_resolution_clock::now();
if (approach == FTCS_APPROACH) { // FTCS case
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole(); printConcentrationsConsole();
} }
if (this->csv_output > CSV_OUTPUT_OFF) { if (csv_output >= CSV_OUTPUT_VERBOSE) {
filename = createCSVfile();
}
auto begin = std::chrono::high_resolution_clock::now();
if (approach == FTCS_APPROACH) { // FTCS case
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
// if (i % (iterations * innerIterations / 100) == 0) {
// double percentage = (double)i / ((double)iterations * (double)innerIterations) * 100;
// if ((int)percentage % 10 == 0) {
// cout << "Progress: " << percentage << "%" << endl;
// }
// }
}
} else if (approach == BTCS_APPROACH) { // BTCS case
if (solver == EIGEN_LU_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads);
}
} else if (solver == THOMAS_ALGORITHM_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
}
}
} else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
double beta = 0.5;
// TODO this implementation is very inefficient!
// a separate implementation that sets up a specific tridiagonal matrix for Crank-Nicolson would be better
MatrixXd concentrations;
MatrixXd concentrationsFTCS;
MatrixXd concentrationsResult;
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
concentrations = grid.getConcentrations();
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsFTCS = grid.getConcentrations();
grid.setConcentrations(concentrations);
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsResult = beta * concentrationsFTCS + (1-beta) * grid.getConcentrations();
grid.setConcentrations(concentrationsResult);
}
}
auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole();
}
if (this->csv_output > CSV_OUTPUT_OFF) {
printConcentrationsCSV(filename); printConcentrationsCSV(filename);
}
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
// if (i % (iterations * innerIterations / 100) == 0) {
// double percentage = (double)i / ((double)iterations *
// (double)innerIterations) * 100; if ((int)percentage % 10 == 0) {
// cout << "Progress: " << percentage << "%" << endl;
// }
// }
} }
if (this->time_measure > TIME_MEASURE_OFF) {
string approachString = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); } else if (approach == BTCS_APPROACH) { // BTCS case
string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << endl; if (solver == EIGEN_LU_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads);
}
} else if (solver == THOMAS_ALGORITHM_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
}
} }
} else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
double beta = 0.5;
// TODO this implementation is very inefficient!
// a separate implementation that sets up a specific tridiagonal matrix for
// Crank-Nicolson would be better
MatrixXd concentrations;
MatrixXd concentrationsFTCS;
MatrixXd concentrationsResult;
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
concentrations = grid.getConcentrations();
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsFTCS = grid.getConcentrations();
grid.setConcentrations(concentrations);
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsResult =
beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations();
grid.setConcentrations(concentrationsResult);
}
}
auto end = std::chrono::high_resolution_clock::now();
auto milliseconds =
std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole();
}
if (this->csv_output > CSV_OUTPUT_OFF) {
printConcentrationsCSV(filename);
}
if (this->time_measure > TIME_MEASURE_OFF) {
string approachString =
(approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
cout << approachString << dimString << ":: run() finished in "
<< milliseconds.count() << "ms" << endl;
}
} }
#endif #endif

15
src/TugUtils.cpp
View File

@ -1,7 +1,7 @@
#include <chrono> #include <chrono>
#include <fstream>
#include <stdexcept> #include <stdexcept>
#include <string> #include <string>
#include <fstream>
using namespace std; using namespace std;
@ -29,10 +29,11 @@ using namespace std;
}) })
// calculates arithmetic or harmonic mean of alpha between two cells // calculates arithmetic or harmonic mean of alpha between two cells
static double calcAlphaIntercell(const double &alpha1, const double &alpha2, bool useHarmonic = true) { static double calcAlphaIntercell(const double &alpha1, const double &alpha2,
if (useHarmonic) { bool useHarmonic = true) {
return double(2) / ((double(1)/alpha1) + (double(1)/alpha2)); if (useHarmonic) {
} else { return double(2) / ((double(1) / alpha1) + (double(1) / alpha2));
return 0.5 * (alpha1 + alpha2); } else {
} return 0.5 * (alpha1 + alpha2);
}
} }

52
test/TestUtils.cpp
View File

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

122
test/testBoundary.cpp
View File

@ -1,68 +1,74 @@
#include <stdio.h>
#include <doctest/doctest.h> #include <doctest/doctest.h>
#include <tug/Boundary.hpp>
#include <string>
#include <typeinfo>
#include <iostream> #include <iostream>
#include <stdio.h>
#include <string>
#include <tug/Boundary.hpp>
#include <typeinfo>
TEST_CASE("BoundaryElement"){ TEST_CASE("BoundaryElement") {
SUBCASE("Closed case"){ SUBCASE("Closed case") {
BoundaryElement boundaryElementClosed = BoundaryElement(); BoundaryElement boundaryElementClosed = BoundaryElement();
CHECK_NOTHROW(BoundaryElement()); CHECK_NOTHROW(BoundaryElement());
CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN)); CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN));
CHECK_THROWS(boundaryElementClosed.setValue(0.2)); CHECK_THROWS(boundaryElementClosed.setValue(0.2));
} }
SUBCASE("Constant case"){ SUBCASE("Constant case") {
BoundaryElement boundaryElementConstant = BoundaryElement(0.1); BoundaryElement boundaryElementConstant = BoundaryElement(0.1);
CHECK_NOTHROW(BoundaryElement(0.1)); CHECK_NOTHROW(BoundaryElement(0.1));
CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT);
CHECK_EQ(boundaryElementConstant.getValue(), 0.1); CHECK_EQ(boundaryElementConstant.getValue(), 0.1);
CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); CHECK_NOTHROW(boundaryElementConstant.setValue(0.2));
CHECK_EQ(boundaryElementConstant.getValue(), 0.2); CHECK_EQ(boundaryElementConstant.getValue(), 0.2);
} }
} }
TEST_CASE("Boundary Class"){ TEST_CASE("Boundary Class") {
Grid grid1D = Grid(10); Grid grid1D = Grid(10);
Grid grid2D = Grid(10, 12); Grid grid2D = Grid(10, 12);
Boundary boundary1D = Boundary(grid1D); Boundary boundary1D = Boundary(grid1D);
Boundary boundary2D = Boundary(grid2D); Boundary boundary2D = Boundary(grid2D);
vector<BoundaryElement> boundary1DVector(1, BoundaryElement(1.0)); vector<BoundaryElement> boundary1DVector(1, BoundaryElement(1.0));
SUBCASE("Boundaries 1D case"){ SUBCASE("Boundaries 1D case") {
CHECK_NOTHROW(Boundary boundary(grid1D)); CHECK_NOTHROW(Boundary boundary(grid1D));
CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1);
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), BC_TYPE_CLOSED); CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); BC_TYPE_CLOSED);
CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP));
CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM));
CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT));
CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP));
CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2));
CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
} BC_TYPE_CONSTANT);
CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
boundary1DVector[0].getType());
}
SUBCASE("Boundaries 2D case"){ SUBCASE("Boundaries 2D case") {
CHECK_NOTHROW(Boundary boundary(grid1D)); CHECK_NOTHROW(Boundary boundary(grid1D));
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);
CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12);
CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); BC_TYPE_CLOSED);
CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP));
CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM));
CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT));
CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP));
CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12));
CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
} BC_TYPE_CONSTANT);
CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
boundary1DVector[0].getType());
}
} }

26
test/testFTCS.cpp
View File

@ -1,19 +1,19 @@
#include <doctest/doctest.h>
#include <../src/FTCS.cpp> #include <../src/FTCS.cpp>
#include <doctest/doctest.h>
#include <limits> #include <limits>
TEST_CASE("Maths") { TEST_CASE("Maths") {
SUBCASE("mean between two alphas") { SUBCASE("mean between two alphas") {
double alpha1 = 10; double alpha1 = 10;
double alpha2 = 20; double alpha2 = 20;
double average = 15; double average = 15;
double harmonicMean = double(2) / ((double(1)/alpha1)+(double(1)/alpha2)); double harmonicMean =
double(2) / ((double(1) / alpha1) + (double(1) / alpha2));
// double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) - harmonicMean);
// CHECK(difference < std::numeric_limits<double>::epsilon());
CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean);
CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average);
}
// double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) -
// harmonicMean); CHECK(difference <
// std::numeric_limits<double>::epsilon());
CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean);
CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average);
}
} }

382
test/testGrid.cpp
View File

@ -1,251 +1,249 @@
#include <Eigen/Core>
#include <doctest/doctest.h> #include <doctest/doctest.h>
#include <tug/Grid.hpp> #include <tug/Grid.hpp>
#include <Eigen/Core>
TEST_CASE("1D Grid, too small length") { TEST_CASE("1D Grid, too small length") {
int l = 2; int l = 2;
CHECK_THROWS(Grid(l)); CHECK_THROWS(Grid(l));
} }
TEST_CASE("2D Grid, too small side") { TEST_CASE("2D Grid, too small side") {
int r = 2; int r = 2;
int c = 4; int c = 4;
CHECK_THROWS(Grid(r, c)); CHECK_THROWS(Grid(r, c));
r = 4; r = 4;
c = 2; c = 2;
CHECK_THROWS(Grid(r, c)); CHECK_THROWS(Grid(r, c));
} }
TEST_CASE("1D Grid") { TEST_CASE("1D Grid") {
int l = 12; int l = 12;
Grid grid(l); Grid grid(l);
SUBCASE("correct construction") { SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 1); CHECK_EQ(grid.getDim(), 1);
CHECK_EQ(grid.getLength(), l); CHECK_EQ(grid.getLength(), 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().rows(), 1);
CHECK_EQ(grid.getConcentrations().cols(), l); CHECK_EQ(grid.getConcentrations().cols(), l);
CHECK_EQ(grid.getAlpha().rows(), 1); CHECK_EQ(grid.getAlpha().rows(), 1);
CHECK_EQ(grid.getAlpha().cols(), l); CHECK_EQ(grid.getAlpha().cols(), l);
CHECK_EQ(grid.getDeltaCol(), 1); CHECK_EQ(grid.getDeltaCol(), 1);
CHECK_THROWS(grid.getAlphaX()); CHECK_THROWS(grid.getAlphaX());
CHECK_THROWS(grid.getAlphaY()); CHECK_THROWS(grid.getAlphaY());
CHECK_THROWS(grid.getDeltaRow()); CHECK_THROWS(grid.getDeltaRow());
} }
SUBCASE("setting concentrations") { SUBCASE("setting concentrations") {
// correct concentrations matrix // correct concentrations matrix
MatrixXd concentrations = MatrixXd::Constant(1, l, 3); MatrixXd concentrations = MatrixXd::Constant(1, l, 3);
CHECK_NOTHROW(grid.setConcentrations(concentrations)); CHECK_NOTHROW(grid.setConcentrations(concentrations));
// false concentrations matrix // false concentrations matrix
MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); CHECK_THROWS(grid.setConcentrations(wConcentrations));
} }
SUBCASE("setting alpha") { SUBCASE("setting alpha") {
// correct alpha matrix // correct alpha matrix
MatrixXd alpha = MatrixXd::Constant(1, l, 3); MatrixXd alpha = MatrixXd::Constant(1, l, 3);
CHECK_NOTHROW(grid.setAlpha(alpha)); CHECK_NOTHROW(grid.setAlpha(alpha));
CHECK_THROWS(grid.setAlpha(alpha, alpha)); CHECK_THROWS(grid.setAlpha(alpha, alpha));
grid.setAlpha(alpha); grid.setAlpha(alpha);
CHECK_EQ(grid.getAlpha(), alpha); CHECK_EQ(grid.getAlpha(), alpha);
CHECK_THROWS(grid.getAlphaX()); CHECK_THROWS(grid.getAlphaX());
CHECK_THROWS(grid.getAlphaY()); CHECK_THROWS(grid.getAlphaY());
// false alpha matrix // false alpha matrix
MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
CHECK_THROWS(grid.setAlpha(wAlpha)); 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 // set 2D domain
CHECK_THROWS(grid.setDomain(d, d)); 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()); CHECK_THROWS(grid.getDeltaRow());
// set too small domain // set too small domain
d = 0; d = 0;
CHECK_THROWS(grid.setDomain(d)); CHECK_THROWS(grid.setDomain(d));
d = -2; d = -2;
CHECK_THROWS(grid.setDomain(d)); CHECK_THROWS(grid.setDomain(d));
} }
} }
TEST_CASE("2D Grid quadratic") { TEST_CASE("2D Grid quadratic") {
int rc = 12; int rc = 12;
Grid grid(rc, rc); Grid grid(rc, rc);
SUBCASE("correct construction") { SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2); CHECK_EQ(grid.getDim(), 2);
CHECK_THROWS(grid.getLength()); 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_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 concentrations") {
// correct concentrations matrix // correct concentrations matrix
MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2);
CHECK_NOTHROW(grid.setConcentrations(concentrations)); CHECK_NOTHROW(grid.setConcentrations(concentrations));
// false concentrations matrix
MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(rc + 3, rc, 4);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
// false concentrations matrix SUBCASE("setting alphas") {
MatrixXd wConcentrations = MatrixXd::Constant(rc, rc+3, 1); // correct alpha matrices
CHECK_THROWS(grid.setConcentrations(wConcentrations)); MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
wConcentrations = MatrixXd::Constant(rc+3, rc, 4); MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
wConcentrations = MatrixXd::Constant(rc+2, rc+4, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
SUBCASE("setting alphas") { CHECK_THROWS(grid.setAlpha(alphax));
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.setAlpha(alphax)); grid.setAlpha(alphax, alphay);
CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_EQ(grid.getAlphaY(), alphay);
CHECK_THROWS(grid.getAlpha());
grid.setAlpha(alphax, alphay); // false alpha matrices
CHECK_EQ(grid.getAlphaX(), alphax); alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
CHECK_EQ(grid.getAlphaY(), alphay); CHECK_THROWS(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.getAlpha()); alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
// false alpha matrices SUBCASE("setting domain") {
alphax = MatrixXd::Constant(rc+3, rc+1, 3); int dr = 8;
CHECK_THROWS(grid.setAlpha(alphax, alphay)); int dc = 9;
alphay = MatrixXd::Constant(rc+2, rc+1, 3);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
SUBCASE("setting domain") { // set 1D domain
int dr = 8; CHECK_THROWS(grid.setDomain(dr));
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(rc));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
grid.setDomain(dr, dc); // set too small domain
CHECK_EQ(grid.getDeltaCol(), double(dc)/double(rc)); dr = 0;
CHECK_EQ(grid.getDeltaRow(), double(dr)/double(rc)); CHECK_THROWS(grid.setDomain(dr, dc));
dr = 8;
// set too small domain dc = 0;
dr = 0; CHECK_THROWS(grid.setDomain(dr, dc));
CHECK_THROWS(grid.setDomain(dr, dc)); dr = -2;
dr = 8; CHECK_THROWS(grid.setDomain(dr, dc));
dc = 0; }
CHECK_THROWS(grid.setDomain(dr, dc));
dr = -2;
CHECK_THROWS(grid.setDomain(dr, dc));
}
} }
TEST_CASE("2D Grid non-quadratic") { TEST_CASE("2D Grid non-quadratic") {
int r = 12; int r = 12;
int c = 15; int c = 15;
Grid grid(r, c); Grid grid(r, c);
SUBCASE("correct construction") { SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2); CHECK_EQ(grid.getDim(), 2);
CHECK_THROWS(grid.getLength()); CHECK_THROWS(grid.getLength());
CHECK_EQ(grid.getCol(), c); CHECK_EQ(grid.getCol(), c);
CHECK_EQ(grid.getRow(), r); CHECK_EQ(grid.getRow(), r);
CHECK_EQ(grid.getConcentrations().rows(), r); CHECK_EQ(grid.getConcentrations().rows(), r);
CHECK_EQ(grid.getConcentrations().cols(), c); CHECK_EQ(grid.getConcentrations().cols(), c);
CHECK_THROWS(grid.getAlpha()); CHECK_THROWS(grid.getAlpha());
CHECK_EQ(grid.getAlphaX().rows(), r); CHECK_EQ(grid.getAlphaX().rows(), r);
CHECK_EQ(grid.getAlphaX().cols(), c); CHECK_EQ(grid.getAlphaX().cols(), c);
CHECK_EQ(grid.getAlphaY().rows(), r); CHECK_EQ(grid.getAlphaY().rows(), r);
CHECK_EQ(grid.getAlphaY().cols(), c); CHECK_EQ(grid.getAlphaY().cols(), c);
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 concentrations") {
// correct concentrations matrix // correct concentrations matrix
MatrixXd concentrations = MatrixXd::Constant(r, c, 2); MatrixXd concentrations = MatrixXd::Constant(r, c, 2);
CHECK_NOTHROW(grid.setConcentrations(concentrations)); CHECK_NOTHROW(grid.setConcentrations(concentrations));
// false concentrations matrix
MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(r + 3, c, 3);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
// false concentrations matrix SUBCASE("setting alphas") {
MatrixXd wConcentrations = MatrixXd::Constant(r, c+3, 6); // correct alpha matrices
CHECK_THROWS(grid.setConcentrations(wConcentrations)); MatrixXd alphax = MatrixXd::Constant(r, c, 2);
wConcentrations = MatrixXd::Constant(r+3, c, 3); MatrixXd alphay = MatrixXd::Constant(r, c, 4);
CHECK_THROWS(grid.setConcentrations(wConcentrations)); CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
wConcentrations = MatrixXd::Constant(r+2, c+4, 2);
CHECK_THROWS(grid.setConcentrations(wConcentrations));
}
SUBCASE("setting alphas") { CHECK_THROWS(grid.setAlpha(alphax));
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.setAlpha(alphax)); grid.setAlpha(alphax, alphay);
CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_EQ(grid.getAlphaY(), alphay);
CHECK_THROWS(grid.getAlpha());
grid.setAlpha(alphax, alphay); // false alpha matrices
CHECK_EQ(grid.getAlphaX(), alphax); alphax = MatrixXd::Constant(r + 3, c + 1, 3);
CHECK_EQ(grid.getAlphaY(), alphay); CHECK_THROWS(grid.setAlpha(alphax, alphay));
CHECK_THROWS(grid.getAlpha()); alphay = MatrixXd::Constant(r + 2, c + 1, 5);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
// false alpha matrices SUBCASE("setting domain") {
alphax = MatrixXd::Constant(r+3, c+1, 3); int dr = 8;
CHECK_THROWS(grid.setAlpha(alphax, alphay)); int dc = 9;
alphay = MatrixXd::Constant(r+2, c+1, 5);
CHECK_THROWS(grid.setAlpha(alphax, alphay));
}
SUBCASE("setting domain") { // set 1D domain
int dr = 8; CHECK_THROWS(grid.setDomain(dr));
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); // set too small domain
CHECK_EQ(grid.getDeltaCol(), double(dc)/double(c)); dr = 0;
CHECK_EQ(grid.getDeltaRow(), double(dr)/double(r)); CHECK_THROWS(grid.setDomain(dr, dc));
dr = 8;
// set too small domain dc = -1;
dr = 0; CHECK_THROWS(grid.setDomain(dr, dc));
CHECK_THROWS(grid.setDomain(dr, dc)); dr = -2;
dr = 8; CHECK_THROWS(grid.setDomain(dr, dc));
dc = -1; }
CHECK_THROWS(grid.setDomain(dr, dc));
dr = -2;
CHECK_THROWS(grid.setDomain(dr, dc));
}
} }

154
test/testSimulation.cpp
View File

@ -1,110 +1,104 @@
#include <stdio.h>
#include <doctest/doctest.h>
#include <tug/Simulation.hpp>
#include "TestUtils.cpp" #include "TestUtils.cpp"
#include <doctest/doctest.h>
#include <stdio.h>
#include <string> #include <string>
#include <tug/Simulation.hpp>
// include the configured header file // include the configured header file
#include <testSimulation.hpp> #include <testSimulation.hpp>
static Grid setupSimulation(APPROACH approach, double timestep, int iterations) { static Grid setupSimulation(APPROACH approach, double timestep,
int row = 11; int iterations) {
int col = 11; int row = 11;
int domain_row = 10; int col = 11;
int domain_col = 10; int domain_row = 10;
int domain_col = 10;
// Grid
Grid grid = Grid(row, col);
grid.setDomain(domain_row, domain_col);
// Grid MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
Grid grid = Grid(row, col); concentrations(5, 5) = 1;
grid.setDomain(domain_row, domain_col); grid.setConcentrations(concentrations);
MatrixXd concentrations = MatrixXd::Constant(row, col, 0); MatrixXd alpha = MatrixXd::Constant(row, col, 1);
concentrations(5,5) = 1; for (int i = 0; i < 5; i++) {
grid.setConcentrations(concentrations); for (int j = 0; j < 6; j++) {
alpha(i, j) = 0.01;
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
for (int i = 0; i < 5; i++) {
for (int j = 0; j < 6; j++) {
alpha(i, j) = 0.01;
}
} }
for (int i = 0; i < 5; i++) { }
for (int j = 6; j < 11; j++) { for (int i = 0; i < 5; i++) {
alpha(i, j) = 0.001; 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++) { for (int i = 5; i < 11; i++) {
alpha(i, j) = 0.1; for (int j = 6; j < 11; j++) {
} alpha(i, j) = 0.1;
} }
grid.setAlpha(alpha, alpha); }
grid.setAlpha(alpha, alpha);
// Boundary
Boundary bc = Boundary(grid);
// Boundary // Simulation
Boundary bc = Boundary(grid); Simulation sim = Simulation(grid, bc, approach);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep);
sim.setIterations(iterations);
sim.run();
// RUN
// Simulation return grid;
Simulation sim = Simulation(grid, bc, approach);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep);
sim.setIterations(iterations);
sim.run();
// RUN
return grid;
} }
TEST_CASE("equality to reference matrix with FTCS") { TEST_CASE("equality to reference matrix with FTCS") {
// set string from the header file // set string from the header file
string test_path = testSimulationCSVDir; string test_path = testSimulationCSVDir;
MatrixXd reference = CSV2Eigen(test_path); MatrixXd reference = CSV2Eigen(test_path);
cout << "FTCS Test: " << endl; cout << "FTCS Test: " << endl;
Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000);
cout << endl; cout << endl;
CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true);
} }
TEST_CASE("equality to reference matrix with BTCS") { TEST_CASE("equality to reference matrix with BTCS") {
// set string from the header file // set string from the header file
string test_path = testSimulationCSVDir; string test_path = testSimulationCSVDir;
MatrixXd reference = CSV2Eigen(test_path); MatrixXd reference = CSV2Eigen(test_path);
cout << "BTCS Test: " << endl; cout << "BTCS Test: " << endl;
Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); Grid grid = setupSimulation(BTCS_APPROACH, 1, 7);
cout << endl; cout << endl;
CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true);
} }
TEST_CASE("Initialize environment"){ TEST_CASE("Initialize environment") {
int rc = 12; int rc = 12;
Grid grid(rc, rc); Grid grid(rc, rc);
Boundary boundary(grid); Boundary boundary(grid);
CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH));
} }
TEST_CASE("Simulation environment"){ TEST_CASE("Simulation environment") {
int rc = 12; int rc = 12;
Grid grid(rc, rc); Grid grid(rc, rc);
Boundary boundary(grid); Boundary boundary(grid);
Simulation sim(grid, boundary, FTCS_APPROACH); Simulation sim(grid, boundary, FTCS_APPROACH);
SUBCASE("default paremeters") { SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); }
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)); 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)); CHECK_THROWS(sim.setTimestep(-0.3));
} }
} }