From 6f9d344cee6979a0031bf569c1a7dbe8738af1bf Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 28 Feb 2022 15:09:46 +0100 Subject: [PATCH] Added new function simulate_base. - With this new function we abstract the actual filling of the A_Matrix and b_vector into processes which are indepent of the dimension. - This code will not run and so the pipeline will fail. --- src/BTCSDiffusion.cpp | 191 ++++++++++++++++++++++-------------------- src/BTCSDiffusion.hpp | 27 +++--- 2 files changed, 113 insertions(+), 105 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 274dd05..019c471 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -1,4 +1,5 @@ #include "BTCSDiffusion.hpp" +#include "BoundaryCondition.hpp" #include @@ -74,71 +75,73 @@ void Diffusion::BTCSDiffusion::updateInternals() { } } -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 == Diffusion::BC_CONSTANT); - bool right_is_constant = (right.type == Diffusion::BC_CONSTANT); +void Diffusion::BTCSDiffusion::simulate_base( + DVectorRowMajor &c, Eigen::Map &bc, + Eigen::Map &alpha, double dx, double time_step, + int size, DVectorRowMajor &t0_c) { // 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; - ; + // 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); + // 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. - */ + // /* + // * 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)); + // // 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(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])); + // A_matrix.insert(size + 1, size + 1) = 1; + // b_vector[size + 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++) { + // 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; - } + // // 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); + // 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; + // 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]; - } + // b_vector[i] = -c[j]; + // } + + fillMatrixFromRow(alpha, bc, size, dx, time_step); + fillVectorFromRowADI(c, alpha, bc, t0_c, 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(!left_is_constant, c.size()); } void Diffusion::BTCSDiffusion::simulate2D( @@ -165,7 +168,8 @@ void Diffusion::BTCSDiffusion::simulate2D( 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)); + fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt, + bc.row(i)); } solveLES(); @@ -197,7 +201,8 @@ void Diffusion::BTCSDiffusion::simulate2D( 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)); + fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt, + bc.col(i)); } solveLES(); @@ -211,81 +216,85 @@ void Diffusion::BTCSDiffusion::simulate2D( 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) { +inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( + const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, + double dx, double time_step) { - n_cols += 2; - int offset = n_cols * row; + Diffusion::boundary_condition left = bc[0]; + Diffusion::boundary_condition right = bc[size - 1]; - A_matrix.insert(offset, offset) = 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; j < size - (1 - 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( + DVectorRowMajor &c, const Eigen::VectorXd alpha, const BCVectorRowMajor &bc, + 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, b_vector[1], alpha[0]); + } + + if (!right_constant) { + b_vector[b_size - 1] = + getBCFromFlux(right, b_vector[size - 2], alpha[size - 1]); } } diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index d8c7cf7..6316f2d 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -112,24 +113,22 @@ private: Eigen::RowMajor> BCVectorRowMajor; - void simulate1D(Eigen::Map &c, - Diffusion::boundary_condition left, - Diffusion::boundary_condition right, - Eigen::Map &bc, - Eigen::Map &alpha, double dx, - int size); + void simulate_base(DVectorRowMajor &c, + Eigen::Map &bc, + Eigen::Map &alpha, double dx, + double time_step, int size, DVectorRowMajor &t0_c); 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); + inline void fillMatrixFromRow(const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, + double time_step); + inline void fillVectorFromRowADI(DVectorRowMajor &c, + const Eigen::VectorXd alpha, + const BCVectorRowMajor &bc, + DVectorRowMajor &t0_c, int size, double dx, + double time_step); void simulate3D(std::vector &c); inline double getBCFromFlux(Diffusion::boundary_condition bc, double nearest_value, double neighbor_alpha);