From 25f8c3fe6e08704c5ee12b4e1f33a55abe2dcfb4 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 7 Aug 2023 12:35:11 +0200 Subject: [PATCH] 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) {