mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
MDL: distinguish between "inner" (= due to CFL) and "outer" iterations
This commit is contained in:
Marco De Lucia
parent
ff611d7a97
commit
25f8c3fe6e
@ -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)
|
||||
|
||||
@ -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 <cstdlib>
|
||||
#include <iostream>
|
||||
#include <tug/Simulation.hpp>
|
||||
|
||||
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 ****
|
||||
|
||||
|
||||
318
src/FTCS.cpp
318
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; col<colMax-1;col++){
|
||||
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; col<colMax-1;col++){
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChange(grid, row, col)
|
||||
)
|
||||
+ allowed_dt / (deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChange(grid, row, col)
|
||||
)
|
||||
;
|
||||
}
|
||||
|
||||
// bottom without corners / looping over columns
|
||||
// hold row constant at max index
|
||||
row = rowMax-1;
|
||||
#pragma omp parallel for
|
||||
for(int col=1; col<colMax-1;col++){
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
}
|
||||
|
||||
// bottom without corners / looping over columns
|
||||
// hold row constant at max index
|
||||
row = rowMax-1;
|
||||
#pragma omp parallel for
|
||||
for(int col=1; col<colMax-1;col++){
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
|
||||
+ allowed_dt / (deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChange(grid, row, col)
|
||||
)
|
||||
;
|
||||
calcHorizontalChange(grid, row, col)
|
||||
)
|
||||
;
|
||||
}
|
||||
|
||||
// corner top left
|
||||
// hold row and column constant at 0
|
||||
row = 0;
|
||||
col = 0;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ allowed_dt/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner top right
|
||||
// hold row constant at 0 and column constant at max index
|
||||
row = 0;
|
||||
col = colMax-1;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ allowed_dt/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner bottom left
|
||||
// hold row constant at max index and column constant at 0
|
||||
row = rowMax-1;
|
||||
col = 0;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ allowed_dt/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner bottom right
|
||||
// hold row and column constant at max index
|
||||
row = rowMax-1;
|
||||
col = colMax-1;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ allowed_dt/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ allowed_dt/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// overwrite obsolete concentrations
|
||||
grid.setConcentrations(concentrations_t1);
|
||||
}
|
||||
|
||||
// corner top left
|
||||
// hold row and column constant at 0
|
||||
row = 0;
|
||||
col = 0;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner top right
|
||||
// hold row constant at 0 and column constant at max index
|
||||
row = 0;
|
||||
col = colMax-1;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner bottom left
|
||||
// hold row constant at max index and column constant at 0
|
||||
row = rowMax-1;
|
||||
col = 0;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// corner bottom right
|
||||
// hold row and column constant at max index
|
||||
row = rowMax-1;
|
||||
col = colMax-1;
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
// overwrite obsolete concentrations
|
||||
grid.setConcentrations(concentrations_t1);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -11,28 +11,49 @@
|
||||
using namespace std;
|
||||
|
||||
Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) {
|
||||
|
||||
this->approach = 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<std::chrono::milliseconds>(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) {
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user