max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 09:28:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'fix_confusion_inner_outer' into 'hannes-philipp'

Browse Source 
fix confusion between CFL/inner and outer iterations

See merge request naaice/tug!11
This commit is contained in:
Philipp Ungrund 2023-08-07 17:23:04 +02:00
commit 5b01149642
6 changed files with 272 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
examples/CMakeLists.txt
View File

@ -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)

32
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -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);
@ -59,14 +67,14 @@ 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
simulation.setTimestep(1000); // timestep
// set the timestep of the simulation
simulation.setTimestep(10000); // timestep
// (optional) set the number of iterations
simulation.setIterations(5);
// set the number of iterations
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 ****

1
include/tug/Simulation.hpp
View File

@ -148,6 +148,7 @@ class Simulation {
double timestep;
int iterations;
int innerIterations;
CSV_OUTPUT csv_output;
CONSOLE_OUTPUT console_output;
TIME_MEASURE time_measure;

304
src/FTCS.cpp
View File

@ -273,160 +273,196 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
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 timestep = timestep;
// if (required_dt > CFL_MDL) {
// 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){
// 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)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChange(grid, row, col)
)
+ timestep / (deltaCol*deltaCol)
calcVerticalChange(grid, 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)
+ 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)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep / (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)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep / (deltaRow*deltaRow)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep / (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)
)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
* (
calcHorizontalChange(grid, row, col)
)
* (
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)
+ timestep / (deltaRow*deltaRow)
* (
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
+ timestep / (deltaCol*deltaCol)
calcVerticalChangeBottomBoundary(grid, bc, row, col)
)
+ timestep / (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)
+ 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 top left
// hold row and column constant at 0
row = 0;
col = 0;
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+ timestep/(deltaCol*deltaCol)
// 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)
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
calcVerticalChangeBottomBoundary(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)
// 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)
calcHorizontalChangeRightBoundary(grid, bc, row, col)
)
+ timestep/(deltaRow*deltaRow)
* (
calcVerticalChangeTopBoundary(grid, bc, row, col)
)
;
calcVerticalChangeBottomBoundary(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);
// overwrite obsolete concentrations
grid.setConcentrations(concentrations_t1);
// }
}

101
src/Simulation.cpp
View File

@ -1,3 +1,5 @@
#include <cmath>
#include <cstddef>
#include <filesystem>
#include <stdexcept>
#include <string>
@ -11,28 +13,20 @@
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;
//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;
this->csv_output = CSV_OUTPUT_OFF;
this->console_output = CONSOLE_OUTPUT_OFF;
this->time_measure = TIME_MEASURE_OFF;
@ -64,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() {
@ -151,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();
@ -161,7 +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;
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
@ -169,11 +220,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) {

8
test/testSimulation.cpp
View File

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