From e1a08ea555554a53deba6d041f23475577fe8423 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 14 Feb 2022 16:46:49 +0100 Subject: [PATCH] Use Eigen::Matrix for internal BC representation --- src/BTCSDiffusion.cpp | 39 ++++++++++++++++----------------------- src/BTCSDiffusion.hpp | 5 +++-- src/main.cpp | 2 +- src/main_2D.cpp | 10 +++++----- 4 files changed, 25 insertions(+), 31 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index dc68093..161df32 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -14,9 +14,9 @@ #include -const int BTCSDiffusion::BC_CONSTANT = 0; -const int BTCSDiffusion::BC_CLOSED = 1; -const int BTCSDiffusion::BC_FLUX = 2; +const int BTCSDiffusion::BC_CLOSED = 0; +const int BTCSDiffusion::BC_FLUX = 1; +const int BTCSDiffusion::BC_CONSTANT = 2; BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { assert(dim <= 3); @@ -65,13 +65,7 @@ void BTCSDiffusion::updateInternals() { deltas[i] = (double)domain_size[i] / grid_cells[i]; } - int cells = 1; - - for (int i = 0; i < grid_dim; i++) { - cells *= (grid_cells[i] + 2); - } - - bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0}); + bc.resize((grid_dim > 1 ? grid_cells[1] : 1), grid_cells[0]); } void BTCSDiffusion::simulate1D(Eigen::Map &c, @@ -121,9 +115,9 @@ void BTCSDiffusion::simulate1D(Eigen::Map &c, i++, j++) { // if current grid cell is considered as constant boundary conditon - if (bc[j].type == BTCSDiffusion::BC_CONSTANT) { + if (bc(1,j).type == BTCSDiffusion::BC_CONSTANT) { A_matrix.insert(i, i) = 1; - b_vector[i] = bc[j].value; + b_vector[i] = bc(1,j).value; continue; } @@ -157,9 +151,9 @@ void BTCSDiffusion::simulate2D(Eigen::Map &c, x_vector.resize(size); for (int i = 0; i < c.rows(); i++) { - boundary_condition left = bc[i * n_cols]; + boundary_condition left = bc(i, 0); bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT; - boundary_condition right = bc[((i + 1) * n_cols) - 1]; + boundary_condition right = bc(i, n_cols - 1); bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT; fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant, @@ -186,13 +180,12 @@ void BTCSDiffusion::simulate2D(Eigen::Map &c, b_vector.resize(size); x_vector.resize(size); - int bottom_offset = bc.size() - (this->grid_cells[0]); n_cols = c.cols(); for (int i = 0; i < c.rows(); i++) { - boundary_condition left = bc[i]; + boundary_condition left = bc(0,i); bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT; - boundary_condition right = bc[bottom_offset + i]; + boundary_condition right = bc(n_cols-1,i); bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT; fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant, @@ -233,7 +226,7 @@ void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, j++, k++) { double sx = (alpha[j - 1] * time_step) / (delta * delta); - if (this->bc[row * (n_cols - 2) + k].type == BTCSDiffusion::BC_CONSTANT) { + if (this->bc(row, k).type == BTCSDiffusion::BC_CONSTANT) { A_matrix.insert(offset + j, offset + j) = 1; continue; } @@ -265,7 +258,7 @@ void BTCSDiffusion::fillVectorFromRow2D(Eigen::Map &c, } for (int j = 1; j < offset - 1; j++) { - boundary_condition tmp_bc = this->bc[ncol * row + (j - 1)]; + boundary_condition tmp_bc = this->bc(row, j-1); if (tmp_bc.type == BTCSDiffusion::BC_CONSTANT) { b_vector[offset * row + j] = tmp_bc.value; @@ -297,7 +290,7 @@ void BTCSDiffusion::simulate(std::vector &c, if (this->grid_dim == 1) { assert(c.size() == grid_cells[0]); Eigen::Map c_in(c.data(), this->grid_cells[0]); - simulate1D(c_in, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0], + simulate1D(c_in, bc(1,0), bc(1,bc.cols()-1), alpha, this->deltas[0], this->grid_cells[0]); } if (this->grid_dim == 2) { @@ -330,10 +323,10 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, return val; } -void BTCSDiffusion::setBoundaryCondition(int index, bctype type, double value) { +void BTCSDiffusion::setBoundaryCondition(int i, int j, bctype type, double value) { - bc[index].type = type; - bc[index].value = value; + bc(i,j).type = type; + bc(i,j).value = value; } inline void BTCSDiffusion::solveLES() { diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 70fe052..2e9b8fd 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -129,7 +129,7 @@ public: * during solving. For flux value refers to a gradient of change for this grid * cell. For closed this value has no effect since a gradient of 0 is used. */ - void setBoundaryCondition(int index, bctype type, double value); + void setBoundaryCondition(int row, int column, bctype type, double value); private: typedef struct boundary_condition { @@ -161,7 +161,8 @@ private: void solveLES(); void updateInternals(); - std::vector bc; + // std::vector bc; + Eigen::Matrix bc; Eigen::SparseMatrix A_matrix; Eigen::VectorXd b_vector; diff --git a/src/main.cpp b/src/main.cpp index 522c5c6..6f03403 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -22,7 +22,7 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(0, BTCSDiffusion::BC_CONSTANT, + diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, 5. * std::pow(10, -6)); // set timestep for simulation to 1 second diff --git a/src/main_2D.cpp b/src/main_2D.cpp index d278425..49878ab 100644 --- a/src/main_2D.cpp +++ b/src/main_2D.cpp @@ -1,5 +1,5 @@ #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET -#include // for copy, max +#include // for copy, max #include #include // for std #include // for vector @@ -14,8 +14,8 @@ int main(int argc, char *argv[]) { int m = 5; // create input + diffusion coefficients for each grid cell - std::vector alpha(n*m, 1 * pow(10, -1)); - std::vector field(n*m, 1 * std::pow(10, -6)); + std::vector alpha(n * m, 1 * pow(10, -1)); + std::vector field(n * m, 1 * std::pow(10, -6)); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -24,7 +24,7 @@ int main(int argc, char *argv[]) { diffu.setYDimensions(1, m); // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(0, BTCSDiffusion::BC_CONSTANT, + diffu.setBoundaryCondition(2, 2, BTCSDiffusion::BC_CONSTANT, 5. * std::pow(10, -6)); // set timestep for simulation to 1 second @@ -43,7 +43,7 @@ int main(int argc, char *argv[]) { for (int i = 0; i < m; i++) { // iterate through columns for (int j = 0; j < n; j++) { - cout << left << std::setw(20) << field[i*n + j]; + cout << left << std::setw(20) << field[i * n + j]; } cout << "\n"; }