proposal implementation for MDL merge request

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]

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;

src/FTCS.cpp

@@ -274,38 +274,40 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     double deltaCol = grid.getDeltaCol();
 
     // MDL: here we have to compute the max time step
-    double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
+    // double deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
-    double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+    // double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
     
-    double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare;
+    // double minDelta2 = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare;
-    double maxAlphaX = grid.getAlphaX().maxCoeff();
+    // double maxAlphaX = grid.getAlphaX().maxCoeff();
-    double maxAlphaY = grid.getAlphaY().maxCoeff();
+    // double maxAlphaY = grid.getAlphaY().maxCoeff();
-    double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
+    // double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
     
-    double CFL_MDL = minDelta2 / (4*maxAlpha); // Formula from Marco --> seems to be unstable
+    // 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_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 MDL: " << CFL_MDL << endl;
-    cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl;
+    // cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl;
-    double required_dt = timestep;
+    // double required_dt = timestep;
-    cout << "FTCS_2D :: required dt=" << required_dt <<  endl;
+    // cout << "FTCS_2D :: required dt=" << required_dt <<  endl;
 
-    int inner_iterations = 1;
+    // int inner_iterations = 1;
-    double allowed_dt = timestep;
+    // double timestep = timestep;
-    if (required_dt > CFL_MDL) {
+    // if (required_dt > CFL_MDL) {
 
-	inner_iterations = (int) ceil(required_dt/CFL_MDL);
+    //   inner_iterations = (int)ceil(required_dt / CFL_MDL);
-	allowed_dt = required_dt/(double)inner_iterations;
+    //   timestep = required_dt / (double)inner_iterations;
-	    
+
-	cout << "FTCS_2D :: Required " << inner_iterations << " inner iterations with dt=" << allowed_dt <<  endl;
+    //   cout << "FTCS_2D :: Required " << inner_iterations
-    } else {
+    //        << " inner iterations with dt=" << timestep << endl;
-	cout << "FTCS_2D :: No inner iterations required, dt=" << required_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){
+    // 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 &timestep) {
 	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 &timestep) {
 #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 &timestep) {
 #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 &timestep) {
 #pragma omp parallel for
 	for (int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-            + allowed_dt / (deltaRow*deltaRow) 
+            + timestep / (deltaRow*deltaRow) 
 	    * (
 	       calcVerticalChangeTopBoundary(grid, bc, row, col)
 	       )
-            + allowed_dt / (deltaCol*deltaCol) 
+            + timestep / (deltaCol*deltaCol) 
 	    * (
 	       calcHorizontalChange(grid, row, col)
 	       )
@@ -387,11 +389,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 #pragma omp parallel for
 	for(int col=1; col<colMax-1;col++){
 	    concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-		+ allowed_dt / (deltaRow*deltaRow) 
+		+ timestep / (deltaRow*deltaRow) 
                 * (
 		   calcVerticalChangeBottomBoundary(grid, bc, row, col)
 		   )
-		+ allowed_dt / (deltaCol*deltaCol) 
+		+ timestep / (deltaCol*deltaCol) 
                 * (
 		   calcHorizontalChange(grid, row, col)
 		   )
@@ -403,11 +405,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 	row = 0;
 	col = 0;
 	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)
             * (
 	       calcVerticalChangeTopBoundary(grid, bc, row, col)
 	       )
@@ -418,11 +420,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 	row = 0;
 	col = colMax-1;
 	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)
             * (
 	       calcVerticalChangeTopBoundary(grid, bc, row, col)
 	       )
@@ -433,11 +435,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 	row = rowMax-1;
 	col = 0;
 	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)
             * (
 	       calcVerticalChangeBottomBoundary(grid, bc, row, col)
 	       )
@@ -448,11 +450,11 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 	row = rowMax-1;
 	col = colMax-1;
 	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)
             * (
 	       calcVerticalChangeBottomBoundary(grid, bc, row, col)
 	       )
@@ -460,7 +462,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 
 	// overwrite obsolete concentrations 
 	grid.setConcentrations(concentrations_t1);
-    }
+    // }
 }
 
 

src/Simulation.cpp

@@ -1,3 +1,5 @@
+#include <cmath>
+#include <cstddef>
 #include <filesystem>
 #include <stdexcept>
 #include <string>
@@ -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();
             }

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