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

Merge branch 'hannes-philipp' into 11-comment-and-error-handling

Browse Source
This commit is contained in:
Hannes Signer 2023-08-05 16:21:12 +02:00
commit eb339667c6
5 changed files with 233 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
examples/reference-FTCS_2D_closed.cpp
View File

@ -5,8 +5,8 @@
using namespace std;
int main(int argc, char *argv[]) {
int row = 11;
int col = 11;
int row = 50;
int col = 50;
int domain_row = 10;
int domain_col = 10;
@ -45,14 +45,11 @@ int main(int argc, char *argv[]) {
// Simulation
Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
sim.setTimestep(0.001);
sim.setIterations(7000);
sim.setOutputCSV(CSV_OUTPUT_ON);
sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setIterations(100);
sim.setOutputCSV(CSV_OUTPUT_OFF);
sim.setOutputConsole(CONSOLE_OUTPUT_OFF);
// RUN
sim.run();
cout << grid.getConcentrations() << endl;
}

7
src/FTCS.cpp
View File

@ -9,6 +9,7 @@
#include <cstddef>
#include <tug/Boundary.hpp>
#include <iostream>
#include <omp.h>
using namespace std;
@ -276,6 +277,8 @@ void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// inner cells
// these are independent of the boundary condition type
omp_set_num_threads(10);
#pragma omp parallel for
for (int row = 1; row < rowMax-1; row++) {
for (int col = 1; col < colMax-1; col++) {
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
@ -295,6 +298,7 @@ void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// left without corners / looping over rows
// hold column constant at index 0
int col = 0;
#pragma omp parallel for
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
@ -311,6 +315,7 @@ void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// right without corners / looping over rows
// hold column constant at max index
col = colMax-1;
#pragma omp parallel for
for (int row = 1; row < rowMax-1; row++) {
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep / (deltaCol*deltaCol)
@ -328,6 +333,7 @@ void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// top without corners / looping over columns
// hold row constant at index 0
int row = 0;
#pragma omp parallel for
for (int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)
@ -344,6 +350,7 @@ void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
// bottom without corners / looping over columns
// hold row constant at max index
row = rowMax-1;
#pragma omp parallel for
for(int col=1; col<colMax-1;col++){
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
+ timestep / (deltaRow*deltaRow)

6
src/Grid.cpp
View File

@ -123,19 +123,19 @@ void Grid::setDomain(int domainLength) {
if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probaly use the 2D domain setter!");
}
if (domainLength < 1) {
if (domainLength <= 0) {
throw_invalid_argument("Given domain length is not positive!");
}
this->domainCol = domainLength;
this->deltaCol = double(this->domainCol)/this->col;
this->deltaCol = double(this->domainCol)/double(this->col);
}
void Grid::setDomain(int domainRow, int domainCol) {
if (dim != 2) {
throw_invalid_argument("Grid is not two dimensional, you should probably use the 1D domain setter!");
}
if (domainRow < 1 || domainCol < 1) {
if (domainRow <= 0 || domainCol <= 0) {
throw_invalid_argument("Given domain size is not positive!");
}

5
src/Simulation.cpp
View File

@ -154,7 +154,7 @@ void Simulation::run() {
}
if (approach == FTCS_APPROACH) {
auto begin = std::chrono::high_resolution_clock::now();
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
@ -165,6 +165,9 @@ void Simulation::run() {
FTCS(grid, bc, timestep);
}
auto end = std::chrono::high_resolution_clock::now();
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
std::cout << milliseconds.count() << endl;
} else if (approach == BTCS_APPROACH) {

219
test/testGrid.cpp
View File

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