From d7e240c6a8ef25e0b4e724ccf08c7882391d1a12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Mon, 28 Feb 2022 14:10:53 +0100 Subject: [PATCH] Refactor simulate function signature --- src/BTCSDiffusion.cpp | 167 ++++++++++++++++++++------------------ src/BTCSDiffusion.hpp | 80 +++++++----------- src/BoundaryCondition.hpp | 33 ++++++++ src/main.cpp | 12 ++- src/main_2D.cpp | 12 ++- 5 files changed, 164 insertions(+), 140 deletions(-) create mode 100644 src/BoundaryCondition.hpp diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index a2b414c..274dd05 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -14,11 +14,7 @@ #include -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) { +Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { assert(dim <= 3); grid_cells.resize(dim, 1); @@ -26,8 +22,8 @@ BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { deltas.resize(dim, 1); } -void BTCSDiffusion::setXDimensions(double domain_size, - unsigned int n_grid_cells) { +void Diffusion::BTCSDiffusion::setXDimensions(double domain_size, + unsigned int n_grid_cells) { assert(this->grid_dim > 0); this->domain_size[0] = domain_size; this->grid_cells[0] = n_grid_cells; @@ -35,8 +31,8 @@ void BTCSDiffusion::setXDimensions(double domain_size, updateInternals(); } -void BTCSDiffusion::setYDimensions(double domain_size, - unsigned int n_grid_cells) { +void Diffusion::BTCSDiffusion::setYDimensions(double domain_size, + unsigned int n_grid_cells) { assert(this->grid_dim > 1); this->domain_size[1] = domain_size; this->grid_cells[1] = n_grid_cells; @@ -44,8 +40,8 @@ void BTCSDiffusion::setYDimensions(double domain_size, updateInternals(); } -void BTCSDiffusion::setZDimensions(double domain_size, - unsigned int n_grid_cells) { +void Diffusion::BTCSDiffusion::setZDimensions(double domain_size, + unsigned int n_grid_cells) { assert(this->grid_dim > 2); this->domain_size[2] = domain_size; this->grid_cells[2] = n_grid_cells; @@ -53,29 +49,38 @@ void BTCSDiffusion::setZDimensions(double domain_size, updateInternals(); } -unsigned int BTCSDiffusion::getXGridCellsN() { return this->grid_cells[0]; } -unsigned int BTCSDiffusion::getYGridCellsN() { return this->grid_cells[1]; } -unsigned int BTCSDiffusion::getZGridCellsN() { return this->grid_cells[2]; } -unsigned int BTCSDiffusion::getXDomainSize() { return this->domain_size[0]; } -unsigned int BTCSDiffusion::getYDomainSize() { return this->domain_size[1]; } -unsigned int BTCSDiffusion::getZDomainSize() { return this->domain_size[2]; } +unsigned int Diffusion::BTCSDiffusion::getXGridCellsN() { + return this->grid_cells[0]; +} +unsigned int Diffusion::BTCSDiffusion::getYGridCellsN() { + return this->grid_cells[1]; +} +unsigned int Diffusion::BTCSDiffusion::getZGridCellsN() { + return this->grid_cells[2]; +} +unsigned int Diffusion::BTCSDiffusion::getXDomainSize() { + return this->domain_size[0]; +} +unsigned int Diffusion::BTCSDiffusion::getYDomainSize() { + return this->domain_size[1]; +} +unsigned int Diffusion::BTCSDiffusion::getZDomainSize() { + return this->domain_size[2]; +} -void BTCSDiffusion::updateInternals() { +void Diffusion::BTCSDiffusion::updateInternals() { for (int i = 0; i < grid_dim; i++) { deltas[i] = (double)domain_size[i] / grid_cells[i]; } - - bc.resize((grid_dim > 1 ? grid_cells[1] : 1), grid_cells[0]); } -void BTCSDiffusion::simulate1D(Eigen::Map &c, - boundary_condition left, - boundary_condition right, - const std::vector &alpha, double dx, - int size) { +void Diffusion::BTCSDiffusion::simulate1D( + Eigen::Map &c, Diffusion::boundary_condition left, + Diffusion::boundary_condition right, Eigen::Map &bc, + Eigen::Map &alpha, double dx, int size) { - bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT); - bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT); + bool left_is_constant = (left.type == Diffusion::BC_CONSTANT); + bool right_is_constant = (right.type == Diffusion::BC_CONSTANT); // The sizes for matrix and vectors of the equation system is defined by the // actual size of the input vector and if the system is (partially) closed. @@ -115,9 +120,9 @@ void BTCSDiffusion::simulate1D(Eigen::Map &c, i++, j++) { // if current grid cell is considered as constant boundary conditon - if (bc(1,j).type == BTCSDiffusion::BC_CONSTANT) { + if (bc[j].type == Diffusion::BC_CONSTANT) { A_matrix.insert(i, i) = 1; - b_vector[i] = bc(1,j).value; + b_vector[i] = bc[j].value; continue; } @@ -136,10 +141,11 @@ void BTCSDiffusion::simulate1D(Eigen::Map &c, c = x_vector.segment(!left_is_constant, c.size()); } -void BTCSDiffusion::simulate2D(Eigen::Map &c, - Eigen::Map &alpha) { +void Diffusion::BTCSDiffusion::simulate2D( + Eigen::Map &c, Eigen::Map &alpha, + Eigen::Map &bc) { - double local_dt = this->time_step/ 2.; + double local_dt = this->time_step / 2.; DMatrixRowMajor tmp_vector; int n_cols = c.cols(); @@ -153,13 +159,13 @@ void BTCSDiffusion::simulate2D(Eigen::Map &c, for (int i = 0; i < c.rows(); i++) { boundary_condition left = bc(i, 0); - bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT; + bool left_constant = left.type == Diffusion::BC_CONSTANT; boundary_condition right = bc(i, n_cols - 1); - bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT; + bool right_constant = right.type == Diffusion::BC_CONSTANT; fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant, - deltas[0], this->time_step / 2); - fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt); + deltas[0], this->time_step / 2, bc.row(i)); + fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt, bc.row(i)); } solveLES(); @@ -184,14 +190,14 @@ void BTCSDiffusion::simulate2D(Eigen::Map &c, n_cols = c.cols(); for (int i = 0; i < c.rows(); i++) { - boundary_condition left = bc(0,i); - bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT; - boundary_condition right = bc(n_cols-1,i); - bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT; + boundary_condition left = bc(0, i); + bool left_constant = left.type == Diffusion::BC_CONSTANT; + boundary_condition right = bc(n_cols - 1, i); + bool right_constant = right.type == Diffusion::BC_CONSTANT; fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant, - deltas[1], this->time_step / 2); - fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt); + deltas[1], this->time_step / 2, bc.col(i)); + fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt, bc.col(i)); } solveLES(); @@ -205,10 +211,12 @@ void BTCSDiffusion::simulate2D(Eigen::Map &c, c.transposeInPlace(); } -void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, - int row, bool left_constant, - bool right_constant, double delta, - double time_step) { +void Diffusion::BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, + int n_cols, int row, + bool left_constant, + bool right_constant, + double delta, double time_step, + const BCVectorRowMajor &bc) { n_cols += 2; int offset = n_cols * row; @@ -227,7 +235,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, k).type == BTCSDiffusion::BC_CONSTANT) { + if (bc[k].type == Diffusion::BC_CONSTANT) { A_matrix.insert(offset + j, offset + j) = 1; continue; } @@ -238,32 +246,31 @@ void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, } } -void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map &c, - const Eigen::VectorXd alpha, int row, - double delta, boundary_condition left, - boundary_condition right, - double time_step) { +void Diffusion::BTCSDiffusion::fillVectorFromRowADI( + Eigen::Map &c, const Eigen::VectorXd alpha, int row, + double delta, boundary_condition left, boundary_condition right, + double time_step, const BCVectorRowMajor &bc) { int ncol = c.cols(); int nrow = c.rows(); int offset = ncol + 2; - if (left.type != BTCSDiffusion::BC_CONSTANT) { + if (left.type != Diffusion::BC_CONSTANT) { // this is not correct currently.We will fix this when we are able to define // FLUX boundary conditions b_vector[offset * row] = getBCFromFlux(left, c(row, 0), alpha[0]); } - if (right.type != BTCSDiffusion::BC_CONSTANT) { + if (right.type != Diffusion::BC_CONSTANT) { b_vector[offset * row + (offset - 1)] = getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]); } for (int j = 0; j < ncol; j++) { - boundary_condition tmp_bc = this->bc(row, j); + boundary_condition tmp_bc = bc[j]; - if (tmp_bc.type == BTCSDiffusion::BC_CONSTANT) { - b_vector[offset * row + (j+1)] = tmp_bc.value; + if (tmp_bc.type == Diffusion::BC_CONSTANT) { + b_vector[offset * row + (j + 1)] = tmp_bc.value; continue; } @@ -282,39 +289,43 @@ void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map &c, } } -void BTCSDiffusion::setTimestep(double time_step) { +void Diffusion::BTCSDiffusion::setTimestep(double time_step) { this->time_step = time_step; } -void BTCSDiffusion::simulate(std::vector &c, - const std::vector &alpha) { +void Diffusion::BTCSDiffusion::simulate(double *c, double *alpha, + Diffusion::boundary_condition *bc) { 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(1,0), bc(1,bc.cols()-1), alpha, this->deltas[0], - this->grid_cells[0]); + Eigen::Map c_in(c, this->grid_cells[0]); + Eigen::Map alpha_in(alpha, this->grid_cells[0]); + Eigen::Map bc_in(bc, this->grid_cells[0]); + + simulate1D(c_in, bc[0], bc[this->grid_cells[0] - 1], bc_in, alpha_in, + this->deltas[0], this->grid_cells[0]); } if (this->grid_dim == 2) { - assert(c.size() == grid_cells[0] * grid_cells[1]); - Eigen::Map c_in(c.data(), this->grid_cells[1], + Eigen::Map c_in(c, this->grid_cells[1], this->grid_cells[0]); - Eigen::Map alpha_in( - alpha.data(), this->grid_cells[1], this->grid_cells[0]); + Eigen::Map alpha_in(alpha, this->grid_cells[1], + this->grid_cells[0]); - simulate2D(c_in, alpha_in); + Eigen::Map bc_in(bc, this->grid_cells[1], + this->grid_cells[0]); + + simulate2D(c_in, alpha_in, bc_in); } } -inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, - double neighbor_c, - double neighbor_alpha) { +inline double Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc, + double neighbor_c, + double neighbor_alpha) { double val; - if (bc.type == BTCSDiffusion::BC_CLOSED) { + if (bc.type == Diffusion::BC_CLOSED) { val = neighbor_c; - } else if (bc.type == BTCSDiffusion::BC_FLUX) { + } else if (bc.type == Diffusion::BC_FLUX) { // TODO // val = bc[index].value; } else { @@ -324,13 +335,7 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, return val; } -void BTCSDiffusion::setBoundaryCondition(int i, int j, bctype type, double value) { - - bc(i,j).type = type; - bc(i,j).value = value; -} - -inline void BTCSDiffusion::solveLES() { +inline void Diffusion::BTCSDiffusion::solveLES() { // start to solve Eigen::SparseLU, Eigen::COLAMDOrdering> solver; diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 8d12660..d8c7cf7 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -1,6 +1,8 @@ #ifndef BTCSDIFFUSION_H_ #define BTCSDIFFUSION_H_ +#include "BoundaryCondition.hpp" + #include #include #include @@ -9,11 +11,7 @@ #include #include -/*! - * Defines both types of boundary condition as a datatype. - */ -typedef int bctype; - +namespace Diffusion { /*! * Class implementing a solution for a 1/2/3D diffusion equation using backward * euler. @@ -21,21 +19,6 @@ typedef int bctype; class BTCSDiffusion { public: - /*! - * Defines a constant/Dirichlet boundary condition. - */ - static const int BC_CONSTANT; - - /*! - * Defines a closed/Neumann boundary condition. - */ - static const int BC_CLOSED; - - /*! - * Defines a flux/Cauchy boundary condition. - */ - static const int BC_FLUX; - /*! * Creates a diffusion module. * @@ -108,7 +91,7 @@ public: * continious memory (row major). * @param alpha Vector of diffusion coefficients for each grid element. */ - void simulate(std::vector &c, const std::vector &alpha); + void simulate(double *c, double *alpha, Diffusion::boundary_condition *bc); /*! * Set the timestep of the simulation @@ -117,53 +100,46 @@ public: */ void setTimestep(double time_step); - /*! - * Set the boundary condition of the given grid. This is done by defining an - * index (exact order still to be determined), the type of the boundary - * condition and the according value. - * - * @param index Index of the grid cell the boundary condition is applied to. - * @param type Type of the boundary condition. Must be constant, closed or - * flux. - * @param value For constant boundary conditions this value is set - * 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 row, int column, bctype type, double value); - private: - typedef struct boundary_condition { - bctype type; - double value; - } boundary_condition; - typedef Eigen::Triplet T; - typedef Eigen::Matrix DMatrixRowMajor; typedef Eigen::Matrix DVectorRowMajor; + typedef Eigen::Matrix + BCMatrixRowMajor; + typedef Eigen::Matrix + BCVectorRowMajor; - void simulate1D(Eigen::Map &c, boundary_condition left, - boundary_condition right, const std::vector &alpha, - double dx, int size); + void simulate1D(Eigen::Map &c, + Diffusion::boundary_condition left, + Diffusion::boundary_condition right, + Eigen::Map &bc, + Eigen::Map &alpha, double dx, + int size); void simulate2D(Eigen::Map &c, - Eigen::Map &alpha); + Eigen::Map &alpha, + Eigen::Map &bc); void fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, int row, bool left_constant, bool right_constant, double delta, - double time_step); + double time_step, const BCVectorRowMajor &bc); void fillVectorFromRowADI(Eigen::Map &c, const Eigen::VectorXd alpha, int row, double delta, - boundary_condition left, boundary_condition right, - double time_step); + Diffusion::boundary_condition left, + Diffusion::boundary_condition right, + double time_step, const BCVectorRowMajor &bc); void simulate3D(std::vector &c); - inline double getBCFromFlux(boundary_condition bc, double nearest_value, - double neighbor_alpha); + inline double getBCFromFlux(Diffusion::boundary_condition bc, + double nearest_value, double neighbor_alpha); void solveLES(); void updateInternals(); // std::vector bc; - Eigen::Matrix bc; + // Eigen::Matrix + // bc; Eigen::SparseMatrix A_matrix; Eigen::VectorXd b_vector; @@ -176,5 +152,5 @@ private: std::vector domain_size; std::vector deltas; }; - +} // namespace Diffusion #endif // BTCSDIFFUSION_H_ diff --git a/src/BoundaryCondition.hpp b/src/BoundaryCondition.hpp new file mode 100644 index 0000000..e47482d --- /dev/null +++ b/src/BoundaryCondition.hpp @@ -0,0 +1,33 @@ +#ifndef BOUNDARYCONDITION_H_ +#define BOUNDARYCONDITION_H_ + +namespace Diffusion { + +/*! + * Defines both types of boundary condition as a datatype. + */ +typedef int bctype; + +/*! + * Defines a closed/Neumann boundary condition. + */ +static const bctype BC_CLOSED = 0; + +/*! + * Defines a flux/Cauchy boundary condition. + */ +static const bctype BC_FLUX = 1; + +/*! + * Defines a constant/Dirichlet boundary condition. + */ +static const bctype BC_CONSTANT = 2; + +typedef struct boundary_condition { + bctype type; + double value; +} boundary_condition; + +} // namespace Diffusion + +#endif // BOUNDARYCONDITION_H_ diff --git a/src/main.cpp b/src/main.cpp index 6f03403..6432527 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,10 +1,14 @@ #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "BoundaryCondition.hpp" #include // for copy, max +#include #include #include // for std #include // for vector using namespace std; +using namespace Diffusion; + int main(int argc, char *argv[]) { // dimension of grid @@ -15,6 +19,7 @@ int main(int argc, char *argv[]) { // create input + diffusion coefficients for each grid cell std::vector alpha(n, 1 * pow(10, -1)); std::vector field(n, 1 * std::pow(10, -6)); + std::vector bc(n, {0,0}); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -22,8 +27,9 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, - 5. * std::pow(10, -6)); + bc[0] = {Diffusion::BC_CONSTANT, 5*std::pow(10,-6)}; + // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, + // 5. * std::pow(10, -6)); // set timestep for simulation to 1 second diffu.setTimestep(1.); @@ -33,7 +39,7 @@ int main(int argc, char *argv[]) { // loop 100 times // output is currently generated by the method itself for (int i = 0; i < 100; i++) { - diffu.simulate(field, alpha); + diffu.simulate(field.data(), alpha.data(), bc.data()); cout << "Iteration: " << i << "\n\n"; diff --git a/src/main_2D.cpp b/src/main_2D.cpp index fbcfaf5..e67004c 100644 --- a/src/main_2D.cpp +++ b/src/main_2D.cpp @@ -1,10 +1,14 @@ #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "BoundaryCondition.hpp" #include // for copy, max +#include #include #include // for std #include // for vector using namespace std; +using namespace Diffusion; + int main(int argc, char *argv[]) { // dimension of grid @@ -16,6 +20,7 @@ int main(int argc, char *argv[]) { // 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 bc(n*m, {0,0}); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -23,9 +28,8 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); diffu.setYDimensions(1, m); - // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(2, 2, BTCSDiffusion::BC_CONSTANT, - 5. * std::pow(10, -6)); + //set inlet to higher concentration without setting bc + field[12] = 5*std::pow(10,-3); // set timestep for simulation to 1 second diffu.setTimestep(1.); @@ -33,7 +37,7 @@ int main(int argc, char *argv[]) { cout << setprecision(12); for (int t = 0; t < 10; t++) { - diffu.simulate(field, alpha); + diffu.simulate(field.data(), alpha.data(), bc.data()); cout << "Iteration: " << t << "\n\n";