From d88d7956a557b9a832181be049f7c2ea676aab48 Mon Sep 17 00:00:00 2001 From: philippun Date: Thu, 20 Jul 2023 11:47:24 +0200 Subject: [PATCH] calculated several examples and implemented csv out --- examples/FTCS_2D_proto_example.cpp | 16 +++++++++-- include/tug/Boundary.hpp | 2 +- include/tug/Simulation.hpp | 6 ++++ proto/FTCS.ipynb | 18 ++++++++---- src/Boundary.cpp | 2 +- src/Simulation.cpp | 45 ++++++++++++++++++++++-------- 6 files changed, 68 insertions(+), 21 deletions(-) diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index c599448..dec4e34 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -25,6 +25,9 @@ int main(int argc, char *argv[]) { // (optional) set the concentrations, e.g.: // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value // grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(20,20,0); + // concentrations(0,0) = 2000; + grid.setConcentrations(concentrations); // (optional) set alphas of the grid, e.g.: // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value @@ -42,6 +45,12 @@ int main(int argc, char *argv[]) { // (optional) set boundary condition values for one side, e.g.: // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values + VectorXd bc_zero_values = VectorXd::Constant(20,0); + bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + VectorXd bc_front_values = VectorXd::Constant(20,2000); + bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); // ************************ @@ -52,10 +61,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(0.01); // timestep + simulation.setTimestep(0.1); // timestep // (optional) set the number of iterations - simulation.setIterations(20); + simulation.setIterations(1000); + + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); // **** RUN SIMULATION **** diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index fc7543a..d633503 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -42,7 +42,7 @@ class Boundary { * @param side * @param values */ - void setBoundaryConditionValue(BC_SIDE side, VectorXd &values); + void setBoundaryConditionValue(BC_SIDE side, VectorXd values); /** * @brief Get the Boundary Condition Value object diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 8e27308..f0a9ebf 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -1,5 +1,7 @@ #include "Boundary.hpp" +using namespace std; + enum APPROACH { FTCS_APPROACH, BTCS_APPROACH @@ -57,6 +59,10 @@ class Simulation { */ auto getIterations(); + void printConcentrationsConsole(); + + void printConcentrationsCSV(string ident); + /** * @brief * diff --git a/proto/FTCS.ipynb b/proto/FTCS.ipynb index f3efad0..b64c934 100644 --- a/proto/FTCS.ipynb +++ b/proto/FTCS.ipynb @@ -113,10 +113,12 @@ " for j in range(1, grid_size['x']-1): #columns\n", " C_t1[i,j] = C_t[i,j] \\\n", " + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n", - " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) \n", + " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n", " + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n", - " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) \n", + " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n", " \n", " # boundary conditions\n", @@ -128,7 +130,8 @@ " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n", " + 2 * alpha_x[i,j] * bc_left) \\\n", " + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n", - " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) \n", + " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \n", "\n", " # right without corners / looping over rows\n", @@ -140,7 +143,8 @@ " - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n", " + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n", " + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n", - " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) \n", + " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n", "\n", " # top without corners / looping over columns\n", @@ -151,7 +155,8 @@ " - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n", " + 2 * alpha_y[0,j] * bc_top) \\\n", " + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n", - " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) \n", + " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n", "\n", " # bottom without corners / looping over columns\n", @@ -163,7 +168,8 @@ " - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n", " + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j]) \\\n", " + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n", - " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", + " - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) \n", + " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n", "\n", " # corner top left\n", diff --git a/src/Boundary.cpp b/src/Boundary.cpp index b464f0d..516c422 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -28,7 +28,7 @@ BC_TYPE Boundary::getBoundaryConditionType() { return this->type; } -void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd &values) { +void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd values) { if (type != BC_TYPE_CONSTANT) { // TODO check if correct way for handling warning cerr << "Values will not be used, wrong BC_TYPE!"; diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 44d1736..53ce95b 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -1,14 +1,12 @@ #include #include +#include + #include "FTCS.cpp" using namespace std; -// auto FTCS(Grid &grid, Boundary &bc, double timestep) { - -// } - Simulation::Simulation(Grid grid, Boundary bc, APPROACH approach) : grid(grid), bc(bc) { //probably to DEBUG assignment of grid and bc this->grid = grid; @@ -45,19 +43,44 @@ auto Simulation::getIterations() { return this->iterations; } +void Simulation::printConcentrationsConsole() { + cout << "Concentrations:" << endl; + cout << grid.getConcentrations() << endl; + cout << endl; +} + +void Simulation::printConcentrationsCSV(string ident) { + ofstream file; + + string filename = "output-" + ident + ".csv"; + // string directory = "output/"; + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } + + IOFormat do_not_align(StreamPrecision, DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << endl; + file << endl << endl; + file.close(); +} void Simulation::run() { if (approach == FTCS_APPROACH) { - cout << grid.getConcentrations() << endl; - cout << endl; + printConcentrationsConsole(); for (int i = 0; i < iterations; i++) { - grid.setConcentrations(FTCS(grid, bc, timestep)); - if (i == 10) { - cout << grid.getConcentrations() << endl; - cout << endl; + if (csv_output == CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV("test"); } + grid.setConcentrations(FTCS(grid, bc, timestep)); + // if (i != 0 && i % 200 == 0) { + // printConcentrationsConsole(); + // } + } + printConcentrationsConsole(); + if (csv_output >= CSV_OUTPUT_ON) { + printConcentrationsCSV("test"); } - cout << grid.getConcentrations() << endl; } else if (approach == BTCS_APPROACH) { for (int i = 0; i < iterations; i++) { //TODO