diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 56bc933..d3fc8a4 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -42,4 +42,4 @@ lint: - make only: refs: - - master + - main diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 274dd05..4d0d195 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -1,4 +1,5 @@ #include "BTCSDiffusion.hpp" +#include "BoundaryCondition.hpp" #include @@ -6,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -14,6 +16,8 @@ #include +#define BTCS_MAX_DEP_PER_CELL 3 + Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { assert(dim <= 3); @@ -73,219 +77,216 @@ void Diffusion::BTCSDiffusion::updateInternals() { deltas[i] = (double)domain_size[i] / grid_cells[i]; } } +void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c, + const BCVectorRowMajor &bc, + const DVectorRowMajor &alpha, + double dx, double time_step, + int size, + const DVectorRowMajor &t0_c) { -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) { + reserveMemory(size, BTCS_MAX_DEP_PER_CELL); - 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. - // Then we will need ghost nodes. So this variable will give the count of - // ghost nodes. - int bc_offset = !left_is_constant + !right_is_constant; - ; - - // set sizes of private and yet allocated vectors - b_vector.resize(size + bc_offset); - x_vector.resize(size + bc_offset); - - /* - * Begin to solve the equation system using LU solver of Eigen. - * - * But first fill the A matrix and b vector. - */ - - // Set boundary condition for ghost nodes (for closed or flux system) or outer - // inlet nodes (constant boundary condition) - A_matrix.resize(size + bc_offset, size + bc_offset); - A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3)); - - A_matrix.insert(0, 0) = 1; - b_vector[0] = - (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - - A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1; - b_vector[size + bc_offset - 1] = - (right_is_constant ? right.value - : getBCFromFlux(right, c[size - 1], alpha[size - 1])); - - // Start filling the A matrix - // =i= is used for equation system matrix and vector indexing - // and =j= for indexing of c,alpha and bc - for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant; - i++, j++) { - - // if current grid cell is considered as constant boundary conditon - if (bc[j].type == Diffusion::BC_CONSTANT) { - A_matrix.insert(i, i) = 1; - b_vector[i] = bc[j].value; - continue; - } - - double sx = (alpha[j] * time_step) / (dx * dx); - - A_matrix.insert(i, i) = -1. - 2. * sx; - A_matrix.insert(i, i - 1) = sx; - A_matrix.insert(i, i + 1) = sx; - - b_vector[i] = -c[j]; - } + fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step); + fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, + dx, time_step); solveLES(); - // write back result to input/output vector - c = x_vector.segment(!left_is_constant, c.size()); + c = x_vector.segment(1, 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; + + DVectorRowMajor input_field = c.row(0); + + simulate_base(input_field, bc, alpha, dx, time_step, size, + Eigen::VectorXd::Constant(size, 0)); + + c.row(0) << input_field; } void Diffusion::BTCSDiffusion::simulate2D( Eigen::Map &c, Eigen::Map &alpha, Eigen::Map &bc) { + int n_rows, n_cols; + double dx; + DMatrixRowMajor t0_c; + double local_dt = this->time_step / 2.; + dx = this->deltas[0]; DMatrixRowMajor tmp_vector; - int n_cols = c.cols(); - unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]); + n_rows = this->grid_cells[1]; + n_cols = this->grid_cells[0]; - A_matrix.resize(size, size); - A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); + t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); - 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; - - 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)); + for (int i = 0; i < n_rows; i++) { + DVectorRowMajor input_field = c.row(i); + simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols, + t0_c.row(i)); + c.row(i) << input_field; } - solveLES(); + dx = this->deltas[1]; - tmp_vector = x_vector; - tmp_vector.transposeInPlace(); - tmp_vector.conservativeResize(c.rows(), c.cols() + 2); + t0_c = + calc_t0_c(c.transpose(), alpha.transpose(), bc.transpose(), local_dt, dx); - Eigen::Map tmp(tmp_vector.data(), c.rows(), c.cols() + 2); - - c = tmp_vector.block(0, 1, c.rows(), c.cols()); - c.transposeInPlace(); - - size = (this->grid_cells[0] * (this->grid_cells[1] + 2)); - - A_matrix.resize(size, size); - A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); - - b_vector.resize(size); - x_vector.resize(size); - - 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; - - 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)); + for (int i = 0; i < n_cols; i++) { + DVectorRowMajor input_field = c.col(i); + simulate_base(input_field, bc.col(i), alpha.col(i), dx, local_dt, n_rows, + t0_c.row(i)); + c.col(i) << input_field.transpose(); } - - solveLES(); - - 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.transposeInPlace(); } -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) { +Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c( + const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, + const BCMatrixRowMajor &bc, double time_step, double dx) { - n_cols += 2; - int offset = n_cols * row; + int n_rows = this->grid_cells[1]; + int n_cols = this->grid_cells[0]; - A_matrix.insert(offset, offset) = 1; + DMatrixRowMajor t0_c(n_rows, n_cols); + + double y_values[3]; + + // first, iterate over first row + for (int j = 0; j < n_cols; j++) { + y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j)); + y_values[1] = c(0, j); + y_values[2] = c(1, j); + + t0_c(0, j) = time_step * alpha(0, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + + // then iterate over inlet + for (int i = 1; i < n_rows - 1; i++) { + for (int j = 0; j < n_cols; j++) { + + y_values[0] = c(i - 1, j); + y_values[1] = c(i, j); + y_values[2] = c(i + 1, j); + + t0_c(i, j) = time_step * alpha(i, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + } + + int end = n_rows - 1; + + // and finally over last row + for (int j = 0; j < n_cols; j++) { + y_values[0] = c(end - 1, j); + y_values[1] = c(end, j); + y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j)); + + t0_c(end, j) = time_step * alpha(end, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + + return t0_c; +} + +inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( + const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, + double dx, double time_step) { + + Diffusion::boundary_condition left = bc[0]; + Diffusion::boundary_condition right = bc[size - 1]; + + bool left_constant = (left.type == Diffusion::BC_CONSTANT); + bool right_constant = (right.type == Diffusion::BC_CONSTANT); + + int A_size = A_matrix.cols(); + + A_matrix.insert(0, 0) = 1; if (left_constant) - A_matrix.insert(offset + 1, offset + 1) = 1; + A_matrix.insert(1, 1) = 1; - A_matrix.insert(offset + (n_cols - 1), offset + (n_cols - 1)) = 1; + A_matrix.insert(A_size - 1, A_size - 1) = 1; if (right_constant) - A_matrix.insert(offset + (n_cols - 2), offset + (n_cols - 2)) = 1; + A_matrix.insert(A_size - 2, A_size - 2) = 1; - for (int j = 1 + left_constant, k = j - 1; j < n_cols - (1 - right_constant); + for (int j = 1 + left_constant, k = j - 1; k < size - right_constant; j++, k++) { - double sx = (alpha[j - 1] * time_step) / (delta * delta); + double sx = (alpha[k] * time_step) / (dx * dx); if (bc[k].type == Diffusion::BC_CONSTANT) { - A_matrix.insert(offset + j, offset + j) = 1; + A_matrix.insert(j, j) = 1; continue; } - A_matrix.insert(offset + j, offset + j) = -1. - 2. * sx; - A_matrix.insert(offset + j, offset + (j - 1)) = sx; - A_matrix.insert(offset + j, offset + (j + 1)) = sx; + A_matrix.insert(j, j) = -1. - 2. * sx; + A_matrix.insert(j, (j - 1)) = sx; + A_matrix.insert(j, (j + 1)) = sx; } } -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) { +inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI( + const DVectorRowMajor &c, const DVectorRowMajor alpha, + const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, + double dx, double time_step) { - int ncol = c.cols(); - int nrow = c.rows(); - int offset = ncol + 2; + Diffusion::boundary_condition left = bc[0]; + Diffusion::boundary_condition right = bc[size - 1]; - 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]); - } + bool left_constant = (left.type == Diffusion::BC_CONSTANT); + bool right_constant = (right.type == Diffusion::BC_CONSTANT); - if (right.type != Diffusion::BC_CONSTANT) { - b_vector[offset * row + (offset - 1)] = - getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]); - } + int b_size = b_vector.size(); - for (int j = 0; j < ncol; j++) { + for (int j = 0; j < size; j++) { boundary_condition tmp_bc = bc[j]; if (tmp_bc.type == Diffusion::BC_CONSTANT) { - b_vector[offset * row + (j + 1)] = tmp_bc.value; + b_vector[j + 1] = tmp_bc.value; continue; } - double y_values[3]; - y_values[0] = - (row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j), alpha[j])); - y_values[1] = c(row, j); - y_values[2] = - (row != nrow - 1 ? c(row + 1, j) - : getBCFromFlux(tmp_bc, c(row, j), alpha[j])); + // double y_values[3]; + // y_values[0] = + // (row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j), + // alpha[j])); + // y_values[1] = c(row, j); + // y_values[2] = + // (row != nrow - 1 ? c(row + 1, j) + // : getBCFromFlux(tmp_bc, c(row, j), alpha[j])); - double t0_c = - time_step * alpha[j] * - ((y_values[0] - 2 * y_values[1] + y_values[2]) / (delta * delta)); - b_vector[offset * row + (j + 1)] = -c(row, j) - (t0_c); + double t0_c_j = time_step * alpha[j] * (t0_c[j] / (dx * dx)); + b_vector[j + 1] = -c[j] - t0_c_j; + } + + if (!left_constant) { + // this is not correct currently.We will fix this when we are able to define + // FLUX boundary conditions + b_vector[0] = getBCFromFlux(left, c[0], alpha[0]); + } + + if (!right_constant) { + b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]); } } @@ -300,8 +301,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 d8c7cf7..1c9556a 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -7,10 +7,13 @@ #include #include #include +#include #include #include #include +#define BTCS_MAX_DEP_PER_CELL 3 + namespace Diffusion { /*! * Class implementing a solution for a 1/2/3D diffusion equation using backward @@ -112,25 +115,33 @@ private: Eigen::RowMajor> BCVectorRowMajor; + void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc, + const DVectorRowMajor &alpha, double dx, double time_step, + int size, const DVectorRowMajor &t0_c); + void simulate1D(Eigen::Map &c, - Diffusion::boundary_condition left, - Diffusion::boundary_condition right, - Eigen::Map &bc, - Eigen::Map &alpha, double dx, - int size); + Eigen::Map &alpha, + Eigen::Map &bc); + void simulate2D(Eigen::Map &c, 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, const BCVectorRowMajor &bc); - void fillVectorFromRowADI(Eigen::Map &c, - const Eigen::VectorXd alpha, int row, double delta, - Diffusion::boundary_condition left, - Diffusion::boundary_condition right, - double time_step, const BCVectorRowMajor &bc); + DMatrixRowMajor calc_t0_c(const DMatrixRowMajor &c, + const DMatrixRowMajor &alpha, + const BCMatrixRowMajor &bc, double time_step, double dx); + + inline void fillMatrixFromRow(const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, + double time_step); + inline void fillVectorFromRowADI(const DVectorRowMajor &c, + const DVectorRowMajor alpha, + const BCVectorRowMajor &bc, + 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();