From 2c2851a0377308586ecb637ca96911be62782aad Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Fri, 4 Aug 2023 00:23:40 +0200 Subject: [PATCH 1/7] comment and input validation of Boundary class --- include/tug/Boundary.hpp | 187 +++++++++++++++++++++++++++++++-------- src/Boundary.cpp | 22 +++++ 2 files changed, 174 insertions(+), 35 deletions(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index b38024a..f865cef 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -20,71 +20,188 @@ enum BC_SIDE { BC_SIDE_BOTTOM }; +/** + * This class defines the boundary conditions of individual boundary elements. + * These can be flexibly used and combined later in other classes. + * The class serves as an auxiliary class for structuring the Boundary class. + */ class BoundaryElement { public: // bc type closed - BoundaryElement(); + /** + * @brief Construct a new Boundary Element object for the closed case. + * The boundary type is here automatically set to the type + * BC_TYPE_CLOSED, where the value takes NaN. + */ + BoundaryElement(); - // bc type constant - BoundaryElement(double value); + /** + * @brief Construct a new Boundary Element object for the constant case. + * The boundary type is automatically set to the type + * BC_TYPE_CONSTANT. + * + * @param value Value of the constant concentration to be assumed at the + * corresponding boundary element. + */ + BoundaryElement(double value); // TODO negative concentration allowed? - void setType(BC_TYPE type); + /** + * @brief Allows changing the boundary type of a corresponding + * BoundaryElement object. + * + * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or + BC_TYPE_CLOSED. + */ + void setType(BC_TYPE type); + + /** + * @brief Sets the value of a boundary condition for the constant case. + * + * @param value Concentration to be considered constant for the + * corresponding boundary element. + */ + void setValue(double value); - void setValue(double value); + /** + * @brief Return the type of the boundary condition, i.e. whether the + * boundary is considered closed or constant. + * + * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or + BC_TYPE_CONSTANT. + */ + BC_TYPE getType(); - BC_TYPE getType(); - - double getValue(); + /** + * @brief Return the concentration value for the constant boundary condition. + * + * @return double Value of the concentration. + */ + double getValue(); private: BC_TYPE type; double value; }; -class BoundaryWall { - public: - BoundaryWall(int length); - void setWall(BC_TYPE type, double value = NAN); +// TODO can be deleted? +// class BoundaryWall { +// public: +// BoundaryWall(int length); - vector getWall(); +// void setWall(BC_TYPE type, double value = NAN); - void setBoundaryElement(int index, BC_TYPE type, double value = NAN); +// vector getWall(); - BoundaryElement getBoundaryElement(); +// void setBoundaryElement(int index, BC_TYPE type, double value = NAN); - private: - BC_SIDE side; - int length; - vector wall; +// BoundaryElement getBoundaryElement(); -}; +// private: +// BC_SIDE side; +// int length; +// vector wall; +// }; + +/** + * This class implements the functionality and management of the boundary + * conditions in the grid to be simulated. + * This class implements the functionality and management of the boundary + * conditions in the grid to be simulated. + */ class Boundary { public: + /** + * @brief Creates a boundary object based on the passed grid object and + * initializes the boundaries as closed. + * + * @param grid Grid object on the basis of which the simulation is to take place. + */ + Boundary(Grid grid); - /** - * @brief Construct a new Boundary object - * - * @param grid - */ - Boundary(Grid grid); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition closed. + * + * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. + */ + void setBoundarySideClosed(BC_SIDE side); - void setBoundarySideClosed(BC_SIDE side); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition constant. Thereby the concentration values of the + * boundaries are set to the passed value. + * + * @param side Side to be set to constant, e.g. BC_SIDE_LEFT. + * @param value Concentration to be set for all elements of the specified page. + */ + void setBoundarySideConstant(BC_SIDE side, double value); - void setBoundarySideConstant(BC_SIDE side, double value); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition closed. + * + * @param side Side in which an element is to be defined as closed. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding side. + */ + void setBoundaryElementClosed(BC_SIDE side, int index); - void setBoundaryElementClosed(BC_SIDE side, int index); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition constant with the + given concentration value. + * + * @param side Side in which an element is to be defined as constant. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding side. + * @param value Concentration value to which the boundary element should be set. + */ + void setBoundaryElementConstant(BC_SIDE side, int index, double value); - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + /** + * @brief Returns the boundary condition of a specified side as a vector + * of BoundarsElement objects. + * + * @param side Boundary side from which the boundaryconditions are to be returned. + * @return vector Contains the boundary conditions as BoundaryElement objects. + */ + vector getBoundarySide(BC_SIDE side); - vector getBoundarySide(BC_SIDE side); + /** + * @brief Returns the boundary condition of a specified element on a given side. + * + * @param side Boundary side in which the boundary condition is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding side. + * @return BoundaryElement Boundary condition as a BoundaryElement object. + */ + BoundaryElement getBoundaryElement(BC_SIDE side, int index); - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + /** + * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED or + BC_TYPE_CONSTANT. + * + * @param side Boundary side in which the boundary condition type is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding side. + * @return BC_TYPE Boundary Type of the corresponding boundary condition. + */ + BC_TYPE getBoundaryElementType(BC_SIDE side, int index); - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); - - double getBoundaryElementValue(BC_SIDE side, int index); + /** + * @brief Returns the concentration value of a corresponding + * BoundaryElement object if it is a constant boundary condition. + * + * @param side Boundary side in which the boundary condition value is + * located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + * @return double Concentration of the corresponding BoundaryElement object. + */ + double getBoundaryElementValue(BC_SIDE side, int index); private: Grid grid; diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 87752ea..2abc938 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,3 +1,4 @@ +#include "TugUtils.hpp" #include "tug/BoundaryCondition.hpp" #include #include @@ -7,6 +8,7 @@ using namespace std; BoundaryElement::BoundaryElement() { + this->type = BC_TYPE_CLOSED; this->value = NAN; } @@ -59,10 +61,18 @@ void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { } void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { + // tests whether the index really points to an element of the boundary side. + if((boundaries[side].size() < index) || index < 0){ + throw_invalid_argument("Index is selected either too large or too small."); + } this->boundaries[side][index].setType(BC_TYPE_CLOSED); } void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) { + // tests whether the index really points to an element of the boundary side. + if((boundaries[side].size() < index) || index < 0){ + throw_invalid_argument("Index is selected either too large or too small."); + } this->boundaries[side][index].setType(BC_TYPE_CONSTANT); this->boundaries[side][index].setValue(value); } @@ -72,14 +82,26 @@ vector Boundary::getBoundarySide(BC_SIDE side) { } BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { + if((boundaries[side].size() < index) || index < 0){ + throw_invalid_argument("Index is selected either too large or too small."); + } return this->boundaries[side][index]; } BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { + if((boundaries[side].size() < index) || index < 0){ + throw_invalid_argument("Index is selected either too large or too small."); + } return this->boundaries[side][index].getType(); } double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { + if((boundaries[side].size() < index) || index < 0){ + throw_invalid_argument("Index is selected either too large or too small."); + } + if(boundaries[side][index].getType() != BC_TYPE_CONSTANT){ + throw_invalid_argument("A value can only be output if it is a constant boundary condition."); + } return this->boundaries[side][index].getValue(); } From 154091e405f40ddd67280b082b4658d985b9b0c5 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Fri, 4 Aug 2023 15:39:02 +0200 Subject: [PATCH 2/7] add user input validation --- include/tug/Boundary.hpp | 28 +++++++--------------------- src/Boundary.cpp | 3 +++ 2 files changed, 10 insertions(+), 21 deletions(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index f865cef..98d0181 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -1,3 +1,9 @@ +/** + * @file Boundary.hpp + * @brief + * + * + */ #ifndef BOUNDARY_H_ #define BOUNDARY_H_ @@ -43,7 +49,7 @@ class BoundaryElement { * @param value Value of the constant concentration to be assumed at the * corresponding boundary element. */ - BoundaryElement(double value); // TODO negative concentration allowed? + BoundaryElement(double value); /** * @brief Allows changing the boundary type of a corresponding @@ -84,26 +90,6 @@ class BoundaryElement { }; -// TODO can be deleted? -// class BoundaryWall { -// public: -// BoundaryWall(int length); - -// void setWall(BC_TYPE type, double value = NAN); - -// vector getWall(); - -// void setBoundaryElement(int index, BC_TYPE type, double value = NAN); - -// BoundaryElement getBoundaryElement(); - -// private: -// BC_SIDE side; -// int length; -// vector wall; - -// }; - /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 2abc938..18f3daf 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -23,6 +23,9 @@ void BoundaryElement::setType(BC_TYPE type) { } void BoundaryElement::setValue(double value) { + if(value < 0){ + throw_invalid_argument("No negative concentration allowed."); + } this->value = value; } From aa4ce6a086ea342693a57f7932c88b82171dab10 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Sat, 5 Aug 2023 21:28:44 +0200 Subject: [PATCH 3/7] comments for Simulation files --- include/tug/Simulation.hpp | 171 ++++++++++++++++++++++--------------- src/Simulation.cpp | 6 ++ 2 files changed, 109 insertions(+), 68 deletions(-) diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 3d123b7..9d334eb 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -1,3 +1,7 @@ +/** + * @file Simulation.hpp + * @brief + */ #include "Boundary.hpp" #include @@ -27,87 +31,118 @@ enum TIME_MEASURE { TIME_MEASURE_VERBOSE // print time measures after each iteration }; +/** + * @brief The class forms the interface for performing the diffusion simulations + * and contains all the methods for controlling the desired parameters, such as + * time step, number of simulations, etc. + * + */ class Simulation { public: + /** + * @brief Set up a runnable simulation environment with the largest stable + * time step and 1000 iterations by passing the required parameters. + * + * @param grid Valid grid object + * @param bc Valid boundary condition object + * @param approach Approach to solving the problem. Either FTCS or BTCS. + */ + Simulation(Grid &grid, Boundary &bc, APPROACH approach); - /** - * @brief Construct a new Simulation object - * - * @param grid - * @param bc - * @param aproach - */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); + /** + * @brief Set the option to output the results to a CSV file. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - CSV_OUTPUT_OFF: do not produce csv output + * - CSV_OUTPUT_ON: produce csv output with last + * concentration matrix + * - CSV_OUTPUT_VERBOSE: produce csv output with all + * concentration matrices + * - CSV_OUTPUT_XTREME: produce csv output with all + * concentration matrices and simulation environment + */ + void setOutputCSV(CSV_OUTPUT csv_output); - /** - * @brief - * - * @param csv_output - */ - void setOutputCSV(CSV_OUTPUT csv_output); + /** + * @brief Set the options for outputting information to the console. + * + * @param console_output Valid output option. The following options can be set + * here: + * - CONSOLE_OUTPUT_OFF: do not print any output to console + * - CONSOLE_OUTPUT_ON: print before and after concentrations to console + * - CONSOLE_OUTPUT_VERBOSE: print all concentration matrices to console + */ + void setOutputConsole(CONSOLE_OUTPUT console_output); - /** - * @brief Set the Output Console object - * - * @param console_output - */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + /** + * @brief Set the Time Measure object + * + * @param time_measure + */ + void setTimeMeasure(TIME_MEASURE time_measure); - /** - * @brief Set the Time Measure object - * - * @param time_measure - */ - void setTimeMeasure(TIME_MEASURE time_measure); + /** + * @brief Setting the time step for each iteration step. Time step must be + * greater than zero. + * + * @param timestep Valid timestep greater than zero. + */ + void setTimestep(double timestep); - /** - * @brief Set the Timestep object - * - * @param timetstep - */ - void setTimestep(double timetstep); + /** + * @brief Currently set time step is returned. + * + * @return double timestep + */ + double getTimestep(); - /** - * @brief Get the Timestep object - * - */ - double getTimestep(); + /** + * @brief Set the desired iterations to be calculated. A value greater + * than zero must be specified here. + * + * @param iterations Number of iterations to be simulated. + */ + void setIterations(int iterations); - /** - * @brief Set the Iterations object - * - * @param iterations - */ - void setIterations(int iterations); + /** + * @brief Return the currently set iterations to be calculated. + * + * @return int Number of iterations. + */ + int getIterations(); - /** - * @brief Get the Iterations object - * - * @return auto - */ - int getIterations(); + /** + * @brief Outputs the current concentrations of the grid on the console. + * + */ + void printConcentrationsConsole(); - /** - * @brief Print the current concentrations of the grid to standard out. - * - */ - void printConcentrationsConsole(); + /** + * @brief Creates a CSV file with a name containing the current simulation + * parameters. If the data name already exists, an additional counter is + * appended to the name. The name of the file is built up as follows: + * + + + -.csv + * + * @return string Filename with given simulation parameter. + */ + string createCSVfile(); - /** - * @brief - * - * @return string - */ - string createCSVfile(); + /** + * @brief Writes the currently calculated concentration values of the grid + * into the CSV file with the passed filename. + * + * @param filename Name of the file to which the concentration values are + * to be written. + */ + void printConcentrationsCSV(string filename); - void printConcentrationsCSV(string filename); - - /** - * @brief - * - * @return Grid - */ - void run(); + /** + * @brief Method starts the simulation process with the previously set + * parameters. + */ + void run(); private: diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 5fd552c..91deee8 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -65,6 +65,9 @@ void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { void Simulation::setTimestep(double timestep) { //TODO check timestep in FTCS for max value + if(timestep <= 0){ + throw_invalid_argument("Timestep has to be greater than zero."); + } this->timestep = timestep; } @@ -73,6 +76,9 @@ double Simulation::getTimestep() { } void Simulation::setIterations(int iterations) { + if(iterations <= 0){ + throw_invalid_argument("Number of iterations must be greater than zero."); + } this->iterations = iterations; } From 30bc6766046595f648c141c303f4ec3c0e052516 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Sun, 6 Aug 2023 19:24:17 +0200 Subject: [PATCH 4/7] add tests for Simulation class --- include/tug/Simulation.hpp | 2 +- test/testSimulation.cpp | 38 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 9d334eb..6ec69c4 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -123,7 +123,7 @@ class Simulation { * @brief Creates a CSV file with a name containing the current simulation * parameters. If the data name already exists, an additional counter is * appended to the name. The name of the file is built up as follows: - * + + + -.csv + * + + + +.csv * * @return string Filename with given simulation parameter. */ diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index aaac3fa..70ad0cf 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -63,3 +63,41 @@ TEST_CASE("equality to reference matrix") { Grid grid = setupSimulation(); CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); } + +TEST_CASE("Initialize environment"){ + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); + + CHECK_NOTHROW(Simulation sim(grid, boundary, FTCS_APPROACH)); + } + +TEST_CASE("Simulation environment"){ + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); + Simulation sim(grid, boundary, FTCS_APPROACH); + + SUBCASE("default paremeters"){ + CHECK_EQ(sim.getIterations(), 1000); + } + + SUBCASE("set iterations"){ + CHECK_NOTHROW(sim.setIterations(2000)); + CHECK_EQ(sim.getIterations(), 2000); + CHECK_THROWS(sim.setIterations(-300)); + } + + SUBCASE("set timestep"){ + CHECK_NOTHROW(sim.setTimestep(0.1)); + CHECK_EQ(sim.getTimestep(), 0.1); + CHECK_THROWS(sim.setTimestep(-0.3)); + } + + SUBCASE("filename"){ + string s1 = sim.createCSVfile(); + string s2 = "FTCS_12_12_1000"; + CHECK_EQ(s1.find(s2) != std::string::npos, true); + } +} + From 0b19b7c1973adab55e0014c063bef82e132835e1 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Sun, 6 Aug 2023 22:14:03 +0200 Subject: [PATCH 5/7] add test cases for Boundary class and additional input validation --- src/Boundary.cpp | 26 +++++++++++++++++ test/testBoundary.cpp | 68 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+) create mode 100644 test/testBoundary.cpp diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 18f3daf..da3a3e4 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -26,6 +26,11 @@ void BoundaryElement::setValue(double value) { if(value < 0){ throw_invalid_argument("No negative concentration allowed."); } + if(type == BC_TYPE_CLOSED){ + throw_invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } this->value = value; } @@ -56,10 +61,24 @@ Boundary::Boundary(Grid grid) : grid(grid) { } void Boundary::setBoundarySideClosed(BC_SIDE side) { + if(grid.getDim() == 1){ + if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ + throw_invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } this->boundaries[side] = vector(grid.getRow(), BoundaryElement()); } void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { + if(grid.getDim() == 1){ + if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ + throw_invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } this->boundaries[side] = vector(grid.getRow(), BoundaryElement(value)); } @@ -81,6 +100,13 @@ void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) } vector Boundary::getBoundarySide(BC_SIDE side) { + if(grid.getDim() == 1){ + if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ + throw_invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } return this->boundaries[side]; } diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp new file mode 100644 index 0000000..47ca1fa --- /dev/null +++ b/test/testBoundary.cpp @@ -0,0 +1,68 @@ +#include +#include +#include +#include +#include +#include + +TEST_CASE("BoundaryElement"){ + + SUBCASE("Closed case"){ + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); + CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); + CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN)); + CHECK_THROWS(boundaryElementClosed.setValue(0.2)); + } + + SUBCASE("Constant case"){ + BoundaryElement boundaryElementConstant = BoundaryElement(0.1); + CHECK_NOTHROW(BoundaryElement(0.1)); + CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); + CHECK_EQ(boundaryElementConstant.getValue(), 0.1); + CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); + CHECK_EQ(boundaryElementConstant.getValue(), 0.2); + } +} + +TEST_CASE("Boundary Class"){ + Grid grid1D = Grid(10); + Grid grid2D = Grid(10, 12); + Boundary boundary1D = Boundary(grid1D); + Boundary boundary2D = Boundary(grid2D); + vector boundary1DVector(1, BoundaryElement(1.0)); + + SUBCASE("Boundaries 1D case"){ + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); + CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); + } + + SUBCASE("Boundaries 2D case"){ + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); + CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); + } +} \ No newline at end of file From 25f8c3fe6e08704c5ee12b4e1f33a55abe2dcfb4 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 7 Aug 2023 12:35:11 +0200 Subject: [PATCH 6/7] MDL: distinguish between "inner" (= due to CFL) and "outer" iterations --- examples/CMakeLists.txt | 3 + examples/FTCS_2D_proto_closed_mdl.cpp | 28 ++- src/FTCS.cpp | 318 ++++++++++++++------------ src/Simulation.cpp | 63 +++-- 4 files changed, 241 insertions(+), 171 deletions(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index eb2dda4..f59d901 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -13,3 +13,6 @@ target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(reference-FTCS_2D_closed tug) # target_link_libraries(FTCS_2D_proto_example easy_profiler) + +add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) +target_link_libraries(FTCS_2D_proto_closed_mdl tug) diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 442b64a..1c6ffea 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -1,22 +1,30 @@ /** - * @file FTCS_2D_proto_example.cpp - * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach - * and constant boundary condition + * @file FTCS_2D_proto_closed_mdl.cpp + * @author Hannes Signer, Philipp Ungrund, MDL + * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary condition; optional command line argument: number of cols and rows * */ +#include +#include #include int main(int argc, char *argv[]) { - + + int row = 64; + + if (argc == 2) { + // no cmd line argument, take col=row=64 + row = atoi(argv[1]); + } + int col=row; + + std::cout << "Nrow =" << row << std::endl; // ************** // **** GRID **** // ************** // create a grid with a 20 x 20 field - int row = 64; - int col = 64; int n2 = row/2-1; Grid grid = Grid(row,col); @@ -60,13 +68,13 @@ int main(int argc, char *argv[]) { Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach // (optional) set the timestep of the simulation - simulation.setTimestep(1000); // timestep + simulation.setTimestep(10000); // timestep // (optional) set the number of iterations - simulation.setIterations(5); + simulation.setIterations(100); // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); // **** RUN SIMULATION **** diff --git a/src/FTCS.cpp b/src/FTCS.cpp index fa9f546..11e8652 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -273,160 +273,194 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { double deltaRow = grid.getDeltaRow(); double deltaCol = grid.getDeltaCol(); - // matrix for concentrations at time t+1 - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); + // MDL: here we have to compute the max time step + double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + + double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + double maxAlphaX = grid.getAlphaX().maxCoeff(); + double maxAlphaY = grid.getAlphaY().maxCoeff(); + double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + + double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable + double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + + cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl; + cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl; + double required_dt = timestep; + cout << "FTCS_2D :: required dt=" << required_dt << endl; - // 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) - + timestep / (deltaRow*deltaRow) + int inner_iterations = 1; + double allowed_dt = timestep; + if (required_dt > CFL_MDL) { + + inner_iterations = (int) ceil(required_dt/CFL_MDL); + allowed_dt = required_dt/(double)inner_iterations; + + cout << "FTCS_2D :: Required " << inner_iterations << " inner iterations with dt=" << allowed_dt << endl; + } else { + cout << "FTCS_2D :: No inner iterations required, dt=" << required_dt << endl; + } + + // we loop for inner iterations + for (int it =0; it < inner_iterations; ++it){ + + cout << "FTCS_2D :: iteration " << it+1 << "/" << inner_iterations << endl; + // matrix for concentrations at time t+1 + 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 + for (int row = 1; row < rowMax-1; row++) { + for (int col = 1; col < colMax-1; col++) { + concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + allowed_dt / (deltaRow*deltaRow) * ( - calcVerticalChange(grid, row, col) - ) - + timestep / (deltaCol*deltaCol) + calcVerticalChange(grid, row, col) + ) + + allowed_dt / (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 - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row, col) = grid.getConcentrations()(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 + for (int row = 1; row < rowMax-1; row++) { + concentrations_t1(row, col) = grid.getConcentrations()(row,col) + + allowed_dt / (deltaCol*deltaCol) * ( - calcHorizontalChangeLeftBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) + calcHorizontalChangeLeftBoundary(grid, bc, row, col) + ) + + allowed_dt / (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 - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) + calcVerticalChange(grid, row, col) + ) + ; + } + + // 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) + + allowed_dt / (deltaCol*deltaCol) * ( - calcHorizontalChangeRightBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) + calcHorizontalChangeRightBoundary(grid, bc, row, col) + ) + + allowed_dt / (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 - for (int col=1; colapproach = approach; - //TODO calculate max time step - double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - - double minDelta = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - double maxAlphaX = grid.getAlphaX().maxCoeff(); - double maxAlphaY = grid.getAlphaY().maxCoeff(); - double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - - double maxStableTimestepMdl = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable - double maxStableTimestep = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia - - // cout << "Max stable time step MDL: " << maxStableTimestepMdl << endl; - // cout << "Max stable time step: " << maxStableTimestep << endl; - - this->timestep = maxStableTimestep; + // MDL no: we need to distinguish between "required dt" and + // "number of (outer) iterations" at which the user needs an + // output and the actual CFL-allowed timestep and consequently the + // number of "inner" iterations which the explicit FTCS needs to + // reach them. The following, at least at the moment, cannot be + // computed here since "timestep" is not yet set when this + // function is called. I brought everything into "FTCS_2D"! - this->iterations = 1000; + // TODO calculate max time step + + // double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + // double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + + // double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + // double maxAlphaX = grid.getAlphaX().maxCoeff(); + // double maxAlphaY = grid.getAlphaY().maxCoeff(); + // double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + + // double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable + // double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + + // cout << "Sim :: CFL condition MDL: " << CFL_MDL << endl; + // double required_dt = this->timestep; + // cout << "Sim :: required dt=" << required_dt << endl; + // cout << "Sim :: CFL condition Wiki: " << CFL_Wiki << endl; + + // if (required_dt > CFL_MDL) { + + // int inner_iterations = (int) ceil(required_dt/CFL_MDL); + // double allowed_dt = required_dt/(double)inner_iterations; + + // cout << "Sim :: Required " << inner_iterations << " inner iterations with dt=" << allowed_dt << endl; + // this->timestep = allowed_dt; + // this->iterations = inner_iterations; + // } else { + // cout << "Sim :: No inner iterations required, dt=" << required_dt << endl; + // } + this->csv_output = CSV_OUTPUT_OFF; this->console_output = CONSOLE_OUTPUT_OFF; this->time_measure = TIME_MEASURE_OFF; @@ -162,6 +183,8 @@ void Simulation::run() { if (approach == FTCS_APPROACH) { auto begin = std::chrono::high_resolution_clock::now(); for (int i = 0; i < iterations; i++) { + // MDL: distinguish between "outer" and "inner" iterations + std::cout << ":: run(): Outer iteration " << i+1 << "/" << iterations << endl; if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); } @@ -169,11 +192,13 @@ void Simulation::run() { printConcentrationsCSV(filename); } - FTCS(grid, bc, timestep); + FTCS(this->grid, this->bc, this->timestep); } auto end = std::chrono::high_resolution_clock::now(); auto milliseconds = std::chrono::duration_cast(end - begin); - std::cout << milliseconds.count() << endl; + + // MDL: meaningful stdout messages + std::cout << ":: run() finished in " << milliseconds.count() << "ms" << endl; } else if (approach == BTCS_APPROACH) { From e1b703849072c7439f2e39f20ccbc8b80fe21316 Mon Sep 17 00:00:00 2001 From: philippun Date: Mon, 7 Aug 2023 16:51:44 +0200 Subject: [PATCH 7/7] proposal implementation for MDL merge request --- examples/FTCS_2D_proto_closed_mdl.cpp | 4 +- include/tug/Simulation.hpp | 1 + src/FTCS.cpp | 88 ++++++++++++----------- src/Simulation.cpp | 100 ++++++++++++++++---------- test/testSimulation.cpp | 8 +-- 5 files changed, 113 insertions(+), 88 deletions(-) diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 1c6ffea..4ed9599 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -67,10 +67,10 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // (optional) set the timestep of the simulation + // set the timestep of the simulation simulation.setTimestep(10000); // timestep - // (optional) set the number of iterations + // set the number of iterations simulation.setIterations(100); // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 6ec69c4..772c0e6 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -148,6 +148,7 @@ class Simulation { double timestep; int iterations; + int innerIterations; CSV_OUTPUT csv_output; CONSOLE_OUTPUT console_output; TIME_MEASURE time_measure; diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 11e8652..8d946e9 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -274,38 +274,40 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { double deltaCol = grid.getDeltaCol(); // MDL: here we have to compute the max time step - double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + // double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + // double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - double maxAlphaX = grid.getAlphaX().maxCoeff(); - double maxAlphaY = grid.getAlphaY().maxCoeff(); - double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + // double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + // double maxAlphaX = grid.getAlphaX().maxCoeff(); + // double maxAlphaY = grid.getAlphaY().maxCoeff(); + // double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable - double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + // double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable + // double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia - cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl; - cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl; - double required_dt = timestep; - cout << "FTCS_2D :: required dt=" << required_dt << endl; + // cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl; + // cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl; + // double required_dt = timestep; + // cout << "FTCS_2D :: required dt=" << required_dt << endl; - int inner_iterations = 1; - double allowed_dt = timestep; - if (required_dt > CFL_MDL) { + // int inner_iterations = 1; + // double timestep = timestep; + // if (required_dt > CFL_MDL) { - inner_iterations = (int) ceil(required_dt/CFL_MDL); - allowed_dt = required_dt/(double)inner_iterations; - - cout << "FTCS_2D :: Required " << inner_iterations << " inner iterations with dt=" << allowed_dt << endl; - } else { - cout << "FTCS_2D :: No inner iterations required, dt=" << required_dt << endl; - } + // inner_iterations = (int)ceil(required_dt / CFL_MDL); + // timestep = required_dt / (double)inner_iterations; + + // cout << "FTCS_2D :: Required " << inner_iterations + // << " inner iterations with dt=" << timestep << endl; + // } else { + // cout << "FTCS_2D :: No inner iterations required, dt=" << required_dt + // << endl; + // } // we loop for inner iterations - for (int it =0; it < inner_iterations; ++it){ + // for (int it =0; it < inner_iterations; ++it){ - cout << "FTCS_2D :: iteration " << it+1 << "/" << inner_iterations << endl; + // cout << "FTCS_2D :: iteration " << it+1 << "/" << inner_iterations << endl; // matrix for concentrations at time t+1 MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); @@ -316,11 +318,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { for (int row = 1; row < rowMax-1; row++) { for (int col = 1; col < colMax-1; col++) { concentrations_t1(row, col) = grid.getConcentrations()(row, col) - + allowed_dt / (deltaRow*deltaRow) + + timestep / (deltaRow*deltaRow) * ( calcVerticalChange(grid, row, col) ) - + allowed_dt / (deltaCol*deltaCol) + + timestep / (deltaCol*deltaCol) * ( calcHorizontalChange(grid, row, col) ) @@ -335,11 +337,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { #pragma omp parallel for for (int row = 1; row < rowMax-1; row++) { concentrations_t1(row, col) = grid.getConcentrations()(row,col) - + allowed_dt / (deltaCol*deltaCol) + + timestep / (deltaCol*deltaCol) * ( calcHorizontalChangeLeftBoundary(grid, bc, row, col) ) - + allowed_dt / (deltaRow*deltaRow) + + timestep / (deltaRow*deltaRow) * ( calcVerticalChange(grid, row, col) ) @@ -352,11 +354,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { #pragma omp parallel for for (int row = 1; row < rowMax-1; row++) { concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + allowed_dt / (deltaCol*deltaCol) + + timestep / (deltaCol*deltaCol) * ( calcHorizontalChangeRightBoundary(grid, bc, row, col) ) - + allowed_dt / (deltaRow*deltaRow) + + timestep / (deltaRow*deltaRow) * ( calcVerticalChange(grid, row, col) ) @@ -370,11 +372,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep) { #pragma omp parallel for for (int col=1; col +#include #include #include #include @@ -13,7 +15,9 @@ using namespace std; Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { this->approach = approach; - + this->timestep = -1; // error per default + this->iterations = -1; + this->innerIterations = 1; // MDL no: we need to distinguish between "required dt" and // "number of (outer) iterations" at which the user needs an @@ -22,37 +26,6 @@ Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid) // reach them. The following, at least at the moment, cannot be // computed here since "timestep" is not yet set when this // function is called. I brought everything into "FTCS_2D"! - - - // TODO calculate max time step - - // double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - // double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - - // double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - // double maxAlphaX = grid.getAlphaX().maxCoeff(); - // double maxAlphaY = grid.getAlphaY().maxCoeff(); - // double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - - // double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable - // double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia - - // cout << "Sim :: CFL condition MDL: " << CFL_MDL << endl; - // double required_dt = this->timestep; - // cout << "Sim :: required dt=" << required_dt << endl; - // cout << "Sim :: CFL condition Wiki: " << CFL_Wiki << endl; - - // if (required_dt > CFL_MDL) { - - // int inner_iterations = (int) ceil(required_dt/CFL_MDL); - // double allowed_dt = required_dt/(double)inner_iterations; - - // cout << "Sim :: Required " << inner_iterations << " inner iterations with dt=" << allowed_dt << endl; - // this->timestep = allowed_dt; - // this->iterations = inner_iterations; - // } else { - // cout << "Sim :: No inner iterations required, dt=" << required_dt << endl; - // } this->csv_output = CSV_OUTPUT_OFF; this->console_output = CONSOLE_OUTPUT_OFF; @@ -85,11 +58,59 @@ void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { } void Simulation::setTimestep(double timestep) { - //TODO check timestep in FTCS for max value if(timestep <= 0){ throw_invalid_argument("Timestep has to be greater than zero."); } - this->timestep = timestep; + + double deltaRowSquare; + double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + double minDeltaSquare; + double maxAlphaX, maxAlphaY, maxAlpha; + if (grid.getDim() == 2) { + + deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + + minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + maxAlphaX = grid.getAlphaX().maxCoeff(); + maxAlphaY = grid.getAlphaY().maxCoeff(); + maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + + } else if (grid.getDim() == 1) { + minDeltaSquare = deltaColSquare; + maxAlpha = grid.getAlpha().maxCoeff(); + + + } else { + throw_invalid_argument("Critical error: Undefined number of dimensions!"); + } + + // TODO check formula 1D case + double CFL_MDL = minDeltaSquare / (4*maxAlpha); // Formula from Marco --> seems to be unstable + double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + + cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl; + cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl; + cout << "FTCS_2D :: required dt=" << timestep << endl; + + if (timestep > CFL_MDL) { + + this->innerIterations = (int)ceil(timestep / CFL_MDL); + 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 << "FTCS_2D :: Required " << this->innerIterations + << " inner iterations with dt=" << this->timestep << endl; + + } else { + + this->timestep = timestep; + cout << "FTCS_2D :: No inner iterations required, dt=" << timestep << endl; + + } + } double Simulation::getTimestep() { @@ -172,6 +193,13 @@ void Simulation::printConcentrationsCSV(string filename) { } void Simulation::run() { + if (this->timestep == -1) { + throw_invalid_argument("Timestep is not set!"); + } + if (this->iterations == -1) { + throw_invalid_argument("Number of iterations are not set!"); + } + string filename; if (this->console_output > CONSOLE_OUTPUT_OFF) { printConcentrationsConsole(); @@ -182,9 +210,9 @@ void Simulation::run() { if (approach == FTCS_APPROACH) { auto begin = std::chrono::high_resolution_clock::now(); - for (int i = 0; i < iterations; i++) { + for (int i = 0; i < iterations * innerIterations; i++) { // MDL: distinguish between "outer" and "inner" iterations - std::cout << ":: run(): Outer iteration " << i+1 << "/" << iterations << endl; + // std::cout << ":: run(): Outer iteration " << i+1 << "/" << iterations << endl; if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); } diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index 70ad0cf..44d6f98 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -79,7 +79,7 @@ TEST_CASE("Simulation environment"){ Simulation sim(grid, boundary, FTCS_APPROACH); SUBCASE("default paremeters"){ - CHECK_EQ(sim.getIterations(), 1000); + CHECK_EQ(sim.getIterations(), -1); } SUBCASE("set iterations"){ @@ -93,11 +93,5 @@ TEST_CASE("Simulation environment"){ CHECK_EQ(sim.getTimestep(), 0.1); CHECK_THROWS(sim.setTimestep(-0.3)); } - - SUBCASE("filename"){ - string s1 = sim.createCSVfile(); - string s2 = "FTCS_12_12_1000"; - CHECK_EQ(s1.find(s2) != std::string::npos, true); - } }