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