From 9a760bd9d98eccd38c82a5fa3363c1fdece42698 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 28 Feb 2022 20:23:50 +0100 Subject: [PATCH] Rewrite simulate1D function --- src/BTCSDiffusion.cpp | 143 ++++++++++++++++++++++++++---------------- src/BTCSDiffusion.hpp | 25 +++++--- 2 files changed, 107 insertions(+), 61 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 019c471..d615d75 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -15,6 +15,8 @@ #include +#define BTCS_MAX_DEP_PER_CELL 3 + Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { assert(dim <= 3); @@ -144,81 +146,116 @@ void Diffusion::BTCSDiffusion::simulate_base( // c = x_vector.segment(!left_is_constant, c.size()); } +inline void Diffusion::BTCSDiffusion::reserveMemory(int size, + int max_count_per_line) { + size += 2; + + A_matrix.resize(size, size); + A_matrix.reserve(Eigen::VectorXi::Constant(size, max_count_per_line)); + + b_vector.resize(size); + x_vector.resize(size); +} + +void Diffusion::BTCSDiffusion::simulate1D( + Eigen::Map &c, Eigen::Map &alpha, + Eigen::Map &bc) { + + int size = this->grid_cells[0]; + double dx = this->deltas[0]; + double time_step = this->time_step; + + reserveMemory(size, BTCS_MAX_DEP_PER_CELL); + + fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step); + fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, + dx, time_step); + + std::cout << A_matrix << std::endl; + + solveLES(); + + c = x_vector.segment(1, size); +} + void Diffusion::BTCSDiffusion::simulate2D( Eigen::Map &c, Eigen::Map &alpha, Eigen::Map &bc) { - double local_dt = this->time_step / 2.; - DMatrixRowMajor tmp_vector; + // double local_dt = this->time_step / 2.; + // DMatrixRowMajor tmp_vector; - int n_cols = c.cols(); - unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]); + // int n_cols = c.cols(); + // unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]); - A_matrix.resize(size, size); - A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); + // A_matrix.resize(size, size); + // A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); - b_vector.resize(size); - x_vector.resize(size); + // b_vector.resize(size); + // x_vector.resize(size); - for (int i = 0; i < c.rows(); i++) { - boundary_condition left = bc(i, 0); - bool left_constant = left.type == Diffusion::BC_CONSTANT; - boundary_condition right = bc(i, n_cols - 1); - bool right_constant = right.type == Diffusion::BC_CONSTANT; + // for (int i = 0; i < c.rows(); i++) { + // boundary_condition left = bc(i, 0); + // bool left_constant = left.type == Diffusion::BC_CONSTANT; + // boundary_condition right = bc(i, n_cols - 1); + // 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, bc.row(i)); - fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt, - bc.row(i)); - } + // fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant, + // 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(); + // solveLES(); - tmp_vector = x_vector; - tmp_vector.transposeInPlace(); - tmp_vector.conservativeResize(c.rows(), c.cols() + 2); + // tmp_vector = x_vector; + // tmp_vector.transposeInPlace(); + // tmp_vector.conservativeResize(c.rows(), c.cols() + 2); - Eigen::Map tmp(tmp_vector.data(), c.rows(), c.cols() + 2); + // Eigen::Map tmp(tmp_vector.data(), c.rows(), c.cols() + 2); - c = tmp_vector.block(0, 1, c.rows(), c.cols()); - c.transposeInPlace(); + // c = tmp_vector.block(0, 1, c.rows(), c.cols()); + // c.transposeInPlace(); - size = (this->grid_cells[0] * (this->grid_cells[1] + 2)); + // size = (this->grid_cells[0] * (this->grid_cells[1] + 2)); - A_matrix.resize(size, size); - A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); + // A_matrix.resize(size, size); + // A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); - b_vector.resize(size); - x_vector.resize(size); + // b_vector.resize(size); + // x_vector.resize(size); - n_cols = c.cols(); + // n_cols = c.cols(); - for (int i = 0; i < c.rows(); i++) { - 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; + // for (int i = 0; i < c.rows(); i++) { + // 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, bc.col(i)); - fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt, - bc.col(i)); - } + // fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant, + // 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(); + // solveLES(); - tmp_vector = x_vector; - tmp_vector.transposeInPlace(); - tmp_vector.conservativeResize(c.rows(), c.cols() + 2); + // tmp_vector = x_vector; + // tmp_vector.transposeInPlace(); + // tmp_vector.conservativeResize(c.rows(), c.cols() + 2); - c = tmp_vector.block(0, 1, c.rows(), c.cols()); + // c = tmp_vector.block(0, 1, c.rows(), c.cols()); - c.transposeInPlace(); + // c.transposeInPlace(); } inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( - const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, - double dx, double time_step) { + const Eigen::VectorXd &alpha, + const Eigen::Vector &bc, + int size, double dx, double time_step) { Diffusion::boundary_condition left = bc[0]; Diffusion::boundary_condition right = bc[size - 1]; @@ -254,8 +291,9 @@ inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( } inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI( - DVectorRowMajor &c, const Eigen::VectorXd alpha, const BCVectorRowMajor &bc, - DVectorRowMajor &t0_c, int size, double dx, double time_step) { + const DVectorRowMajor &c, const Eigen::VectorXd alpha, + const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, + double dx, double time_step) { Diffusion::boundary_condition left = bc[0]; Diffusion::boundary_condition right = bc[size - 1]; @@ -309,8 +347,7 @@ void Diffusion::BTCSDiffusion::simulate(double *c, double *alpha, 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]); + simulate1D(c_in, alpha_in, bc_in); } if (this->grid_dim == 2) { Eigen::Map c_in(c, this->grid_cells[1], diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 6316f2d..72a54b9 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -12,6 +12,8 @@ #include #include +#define BTCS_MAX_DEP_PER_CELL 3 + namespace Diffusion { /*! * Class implementing a solution for a 1/2/3D diffusion equation using backward @@ -113,23 +115,30 @@ private: Eigen::RowMajor> BCVectorRowMajor; - void simulate_base(DVectorRowMajor &c, - Eigen::Map &bc, + void simulate_base(DVectorRowMajor &c, Eigen::Map &bc, Eigen::Map &alpha, double dx, double time_step, int size, DVectorRowMajor &t0_c); + + void simulate1D(Eigen::Map &c, + Eigen::Map &alpha, + Eigen::Map &bc); + void simulate2D(Eigen::Map &c, Eigen::Map &alpha, Eigen::Map &bc); - inline void fillMatrixFromRow(const DVectorRowMajor &alpha, - const BCVectorRowMajor &bc, int size, double dx, - double time_step); - inline void fillVectorFromRowADI(DVectorRowMajor &c, + inline void fillMatrixFromRow( + const Eigen::VectorXd &alpha, + const Eigen::Vector &bc, + int size, double dx, double time_step); + inline void fillVectorFromRowADI(const DVectorRowMajor &c, const Eigen::VectorXd alpha, const BCVectorRowMajor &bc, - DVectorRowMajor &t0_c, int size, double dx, - double time_step); + const DVectorRowMajor &t0_c, int size, + double dx, double time_step); void simulate3D(std::vector &c); + + inline void reserveMemory(int size, int max_count_per_line); inline double getBCFromFlux(Diffusion::boundary_condition bc, double nearest_value, double neighbor_alpha); void solveLES();