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

implemented BTCS 2D constant case

Browse Source
This commit is contained in:
philippun 2023-08-10 22:21:15 +02:00
parent 435314ba61
commit fdb5c436ea
10 changed files with 268 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
examples/BTCS_2D_proto_example.cpp Normal file
View File

@ -0,0 +1,80 @@
#include <tug/Simulation.hpp>
int main(int argc, char *argv[]) {
// EASY_PROFILER_ENABLE;
// profiler::startListen();
// **************
// **** GRID ****
// **************
// profiler::startListen();
// create a grid with a 20 x 20 field
int row = 20;
int col = 20;
Grid grid = Grid(row,col);
// (optional) set the domain, e.g.:
// grid.setDomain(20, 20);
// (optional) set the concentrations, e.g.:
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
// grid.setConcentrations(concentrations);
MatrixXd concentrations = MatrixXd::Constant(row,col,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
// MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value
// grid.setAlpha(alphax, alphay);
// ******************
// **** BOUNDARY ****
// ******************
// create a boundary with constant values
Boundary bc = Boundary(grid);
bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
// (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);
// ************************
// **** SIMULATION ENV ****
// ************************
// set up a simulation environment
Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
// set the timestep of the simulation
simulation.setTimestep(0.1); // timestep
// set the number of iterations
simulation.setIterations(100);
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
// **** RUN SIMULATION ****
// run the simulation
// EASY_BLOCK("SIMULATION")
simulation.run();
// EASY_END_BLOCK;
// profiler::dumpBlocksToFile("test_profile.prof");
// profiler::stopListen();
}

2
examples/CMakeLists.txt
View File

@ -2,6 +2,7 @@ add_executable(first_example first_example.cpp)
add_executable(second_example second_example.cpp)
add_executable(boundary_example1D boundary_example1D.cpp)
add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
@ -9,6 +10,7 @@ target_link_libraries(first_example tug)
target_link_libraries(second_example tug)
target_link_libraries(boundary_example1D tug)
target_link_libraries(FTCS_2D_proto_example tug)
target_link_libraries(BTCS_2D_proto_example tug)
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)

5
include/tug/Boundary.hpp
View File

@ -32,7 +32,7 @@ enum BC_SIDE {
*/
class BoundaryElement {
public:
// bc type closed
/**
* @brief Construct a new Boundary Element object for the closed case.
* The boundary type is here automatically set to the type
@ -154,6 +154,9 @@ class Boundary {
*/
vector<BoundaryElement> getBoundarySide(BC_SIDE side);
// TODO write documentation and tests for this method
VectorXd getBoundarySideValues(BC_SIDE side);
/**
* @brief Returns the boundary condition of a specified element on a given side.
*

28
proto/BTCS.ipynb
View File

@ -98,9 +98,9 @@
"metadata": {},
"outputs": [],
"source": [
"##############\n",
"### CLOSED ###\n",
"##############\n",
"################\n",
"### CONSTANT ###\n",
"################\n",
"\n",
"# creates the coeffiecient matrix for a given row\n",
"def create_coeff_matrix(alpha_x, row_index, s_x):\n",
@ -153,8 +153,6 @@
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
"\n",
" # first row\n",
" # TODO check formula for alpha boundary cases\n",
" # added factor of 2 in front of alpha_y!!!\n",
" if row_index == 0:\n",
" for i in range(0, sv.shape[0]):\n",
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
@ -168,9 +166,7 @@
"\n",
"\n",
" # last row\n",
" # TODO check formula for alpha boundary cases\n",
" # added factor of 2 in front of alpha_y!!!\n",
" if row_index == sv.shape[0]-1:\n",
" if row_index == concentrations.shape[0]-1:\n",
" for i in range(0, sv.shape[0]):\n",
" i = i\n",
" sv[i] = s_y * alpha_y[row_index,i]*bc_bottom \\\n",
@ -184,12 +180,10 @@
"\n",
"\n",
" # first column (additional fixed concentration change from perpendicular dimension)\n",
" # TODO check correctness of additional summand\n",
" i = 0\n",
" sv[0] += 2 * s_x * alpha_x[row_index,0] * bc_left # alpha_x[row_index,0]?\n",
"\n",
" # last column (additional fixed concentration change from perpendicular direction)\n",
" # TODO check correctness of additional summand\n",
" n = sv.shape[0]-1\n",
" sv[n] += 2 * s_x * alpha_x[row_index,n] * bc_right\n",
"\n",
@ -233,9 +227,9 @@
"metadata": {},
"outputs": [],
"source": [
"################\n",
"### CONSTANT ###\n",
"################\n",
"##############\n",
"### CLOSED ###\n",
"##############\n",
"\n",
"# creates the coeffiecient matrix for a given row\n",
"def create_coeff_matrix(alpha_x, row_index, s_x):\n",
@ -246,7 +240,7 @@
" cm = np.zeros((alpha_x.shape[1],alpha_x.shape[1]))\n",
"\n",
" # top left\n",
" cm[0,0] = 1 + s_x * (alpha_interblock(alpha_x[r,0], alpha_x[r,1])) # alpha_x[0,0]\n",
" cm[0,0] = 1 + s_x * alpha_interblock(alpha_x[r,0], alpha_x[r,1]) # alpha_x[0,0]\n",
" cm[0,1] = -s_x * alpha_interblock(alpha_x[r, 0], alpha_x[r, 1])\n",
"\n",
" # inner cells\n",
@ -288,8 +282,6 @@
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
"\n",
" # first row\n",
" # TODO check formula for alpha boundary cases\n",
" # added factor of 2 in front of alpha_y!!!\n",
" if row_index == 0:\n",
" for i in range(0, sv.shape[0]):\n",
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
@ -301,8 +293,6 @@
"\n",
"\n",
" # last row\n",
" # TODO check formula for alpha boundary cases\n",
" # added factor of 2 in front of alpha_y!!!\n",
" if row_index == sv.shape[0]-1:\n",
" for i in range(0, sv.shape[0]):\n",
" i = i\n",
@ -601,7 +591,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
"version": "3.11.4"
},
"orig_nbformat": 4
},

161
src/BTCSv2.cpp
View File

@ -1,41 +1,176 @@
#include "TugUtils.hpp"
#include "FTCS.cpp"
#include <arm_neon.h>
#include <cstddef>
#include <tug/Boundary.hpp>
using namespace Eigen;
// calculates arithmetic or harmonic mean of alpha between two cells
static double calcAlphaIntercell(double &alpha1, double &alpha2, bool useHarmonic = true) {
if (useHarmonic) {
return double(2) / ((double(1)/alpha1) + (double(1)/alpha2));
} else {
return 0.5 * (alpha1 + alpha2);
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int dim, int rowIndex, double sx) {
// square matrix of column^2 dimension for the coefficients
SparseMatrix<double> cm(dim, dim);
cm.reserve(VectorXi::Constant(dim, 3));
// left column
cm.insert(0,0) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
+ 2 * alpha(rowIndex,0));
cm.insert(0,1) = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
// inner columns
int n = dim-1;
for (int i = 1; i < n; i++) {
cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i));
cm.insert(i,i) = 1 + sx * (
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1))
+ calcAlphaIntercell(alpha(rowIndex,i-1), alpha(i,1))
);
cm.insert(i,i+1) = -sx * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1));
}
// right column
cm.insert(n,n-1) = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
cm.insert(n,n) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n))
+ 2 * alpha(rowIndex,n));
cm.makeCompressed();
return cm;
}
static MatrixXd createCoeffMatrix(Grid &grid, int rowIndex, double sx) {
static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
VectorXd &bcLeft, VectorXd &bcRight, VectorXd &bcTop,
VectorXd &bcBottom, int length, int rowIndex, double sx, double sy) {
// square matrix of column^2 dimension for the coefficients
int dim = grid.getCol();
SparseMatrix<double, RowMajor> cm(dim, dim);
VectorXd sv(length);
int numRows = concentrations.rows();
// top left
cm.coeffRef(0,0) = 1 + sx * (calcAlphaIntercell(grid.getAlphaX()(rowIndex,0), grid.getAlphaX()(rowIndex,1)));
// inner rows
if (rowIndex > 0 && rowIndex < numRows-1) {
for (int i = 0; i < length; i++) {
sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
* concentrations(rowIndex+1,i)
+ (
1 - sy * (
calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
+ calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
)
) * concentrations(rowIndex,i)
+ sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
* concentrations(rowIndex-1,i)
;
}
}
// first row
if (rowIndex == 0) {
for (int i = 0; i < length; i++) {
sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
* concentrations(rowIndex,i)
+ (
1 - sy * (
calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
+ 2 * alphaY(rowIndex,i)
)
) * concentrations(rowIndex,i)
+ sy * alphaY(rowIndex,i) * bcTop(i)
;
}
}
// last row
if (rowIndex == numRows-1) {
for (int i = 0; i < length; i++) {
sv(i) = sy * alphaY(rowIndex,i) * bcBottom(i)
+ (
1 - sy * (
2 * alphaY(rowIndex,i)
+ calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
)
) * concentrations(rowIndex,i)
+ sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
* concentrations(rowIndex-1,i)
;
}
}
// first column -> additional fixed concentration change from perpendicular dimension
sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(0);
// last column -> additional fixed concentration change from perpendicular dimension
sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(length-1);
return sv;
}
// BTCS solution for 1D grid
static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
// TODO
}
static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
SparseLU<SparseMatrix<double>> solver;
solver.analyzePattern(A);
solver.factorize(A);
return solver.solve(b);
}
// BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
VectorXd row_t1(colMax);
SparseMatrix<double> A;
VectorXd b;
MatrixXd alphaX = grid.getAlphaX();
MatrixXd alphaY = grid.getAlphaY();
VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
VectorXd bcTop = bc.getBoundarySideValues(BC_SIDE_TOP);
VectorXd bcBottom = bc.getBoundarySideValues(BC_SIDE_BOTTOM);
MatrixXd concentrations = grid.getConcentrations();
for (int i = 0; i < rowMax; i++) {
A = createCoeffMatrix(alphaX, colMax, i, sx);
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy);
row_t1 = solve(A, b);
for (int j = 0; j < colMax; j++) {
concentrations_t1(i,j) = row_t1(j);
}
}
concentrations_t1.transposeInPlace();
concentrations.transposeInPlace();
alphaX.transposeInPlace();
alphaY.transposeInPlace();
for (int i = 0; i < colMax; i++) {
A = createCoeffMatrix(alphaY, rowMax, i, sy);
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
bcLeft, bcRight, rowMax, i, sy, sx);
row_t1 = solve(A, b);
for (int j = 0; j < rowMax; j++) {
concentrations(i,j) = row_t1(j);
}
}
concentrations.transposeInPlace();
grid.setConcentrations(concentrations);
}

11
src/Boundary.cpp
View File

@ -110,6 +110,17 @@ vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
return this->boundaries[side];
}
VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
int length = boundaries[side].size();
VectorXd values(length);
for (int i = 0; i < length; i++) {
values(i) = getBoundaryElementValue(side, i);
}
return values;
}
BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) {
if((boundaries[side].size() < index) || index < 0){
throw_invalid_argument("Index is selected either too large or too small.");

2
src/CMakeLists.txt
View File

@ -1,4 +1,4 @@
add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp)
add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp)
target_link_libraries(tug Eigen3::Eigen)

10
src/FTCS.cpp
View File

@ -14,16 +14,6 @@
using namespace std;
// calculates arithmetic or harmonic mean of alpha between two cells
static double calcAlphaIntercell(double &alpha1, double &alpha2, bool useHarmonic = true) {
if (useHarmonic) {
return double(2) / ((double(1)/alpha1) + (double(1)/alpha2));
} else {
return 0.5 * (alpha1 + alpha2);
}
}
// calculates horizontal change on one cell independent of boundary type
static double calcHorizontalChange(Grid &grid, int &row, int &col) {

9
src/Simulation.cpp
View File

@ -7,8 +7,7 @@
#include <fstream>
#include "FTCS.cpp"
#include "TugUtils.hpp"
#include "BTCSv2.cpp"
#include <tug/progressbar.hpp>
@ -234,7 +233,7 @@ void Simulation::run() {
} else if (approach == BTCS_APPROACH) {
for (int i = 0; i < iterations; i++) {
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
@ -242,8 +241,8 @@ void Simulation::run() {
printConcentrationsCSV(filename);
}
//TODO
break;
BTCS(this->grid, this->bc, this->timestep);
}
}

9
src/TugUtils.hpp
View File

@ -24,3 +24,12 @@
duration.count(); \
})
#endif // BTCSUTILS_H_
// calculates arithmetic or harmonic mean of alpha between two cells
static double calcAlphaIntercell(double &alpha1, double &alpha2, bool useHarmonic = true) {
if (useHarmonic) {
return double(2) / ((double(1)/alpha1) + (double(1)/alpha2));
} else {
return 0.5 * (alpha1 + alpha2);
}
}