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

calculated several examples and implemented csv out

Browse Source
This commit is contained in:
philippun 2023-07-20 11:47:24 +02:00
parent 67f289c1f8
commit d88d7956a5
6 changed files with 68 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
examples/FTCS_2D_proto_example.cpp
View File

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

2
include/tug/Boundary.hpp
View File

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

6
include/tug/Simulation.hpp
View File

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

18
proto/FTCS.ipynb
View File

@ -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",

2
src/Boundary.cpp
View File

@ -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!";

45
src/Simulation.cpp
View File

@ -1,14 +1,12 @@
#include <stdexcept>
#include <tug/Simulation.hpp>
#include <fstream>
#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