From f4924ac8b2965e291ba205d5832e437d628248e4 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 1 Aug 2023 16:56:25 +0200 Subject: [PATCH] add: Implementation of max time step --- examples/reference-FTCS_2D_closed.cpp | 4 ++-- src/Simulation.cpp | 18 +++++++++++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index ad0c466..89c3b29 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -40,8 +40,8 @@ int main(int argc, char *argv[]) { // Simulation Simulation sim = Simulation(grid, bc, FTCS_APPROACH); - sim.setTimestep(0.001); - sim.setIterations(7000); + //sim.setTimestep(0.001); + sim.setIterations(2); sim.setOutputCSV(CSV_OUTPUT_VERBOSE); diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 6d7bc6f..e6aea04 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -13,7 +13,23 @@ Simulation::Simulation(Grid grid, Boundary bc, APPROACH approach) : grid(grid), this->approach = approach; //TODO calculate max time step - this->timestep = 0.01; + + 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 maxStableTimestep = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable + double maxStableTimestep = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + + cout << maxStableTimestep << endl; + + this->timestep = maxStableTimestep; + + this->iterations = 1000; this->csv_output = CSV_OUTPUT_OFF; }