From 6f9d344cee6979a0031bf569c1a7dbe8738af1bf Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 28 Feb 2022 15:09:46 +0100 Subject: [PATCH 01/12] 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); From 9a760bd9d98eccd38c82a5fa3363c1fdece42698 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 28 Feb 2022 20:23:50 +0100 Subject: [PATCH 02/12] 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(); From b7b37e92318ffcf31a71dbb0f850e6a73dab3b54 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:08:24 +0100 Subject: [PATCH 03/12] Update indexing + Bug fix - Wrong stopping criteria @ filling of matrix - Fill left and right side of b_vector with values from c instead of b_vector --- src/BTCSDiffusion.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index d615d75..814a6a0 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -171,8 +171,6 @@ void Diffusion::BTCSDiffusion::simulate1D( 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); @@ -275,7 +273,7 @@ inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( if (right_constant) A_matrix.insert(A_size - 2, A_size - 2) = 1; - for (int j = 1 + left_constant, k = j - 1; j < size - (1 - right_constant); + for (int j = 1 + left_constant, k = j - 1; k < size - right_constant; j++, k++) { double sx = (alpha[k] * time_step) / (dx * dx); @@ -327,12 +325,12 @@ inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI( 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]); + b_vector[0] = getBCFromFlux(left, c[0], alpha[0]); } if (!right_constant) { b_vector[b_size - 1] = - getBCFromFlux(right, b_vector[size - 2], alpha[size - 1]); + getBCFromFlux(right, c[size - 1], alpha[size - 1]); } } From d65fcd453b25e0bafef4067ea04df75866477719 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:10:08 +0100 Subject: [PATCH 04/12] Upodate .gitlab-ci to run lint only @ main --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From d0b75496c7c6d37f71c79ccb22a36aa491c3677c Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:18:55 +0100 Subject: [PATCH 05/12] Remove simulate_base function. --- src/BTCSDiffusion.cpp | 108 +++++++++++++++++++++--------------------- src/BTCSDiffusion.hpp | 7 +-- 2 files changed, 58 insertions(+), 57 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 814a6a0..71ff04b 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -77,74 +77,74 @@ void Diffusion::BTCSDiffusion::updateInternals() { } } -void Diffusion::BTCSDiffusion::simulate_base( - DVectorRowMajor &c, Eigen::Map &bc, - Eigen::Map &alpha, double dx, double time_step, - int size, DVectorRowMajor &t0_c) { +// 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; - // ; +// // 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); +// // 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. - // */ +// // /* +// // * 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 + 1, size + 1) = 1; - // b_vector[size + 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++) { +// // 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; - // } +// // // 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); +// fillMatrixFromRow(alpha, bc, size, dx, time_step); +// fillVectorFromRowADI(c, alpha, bc, t0_c, size, dx, time_step); - solveLES(); +// solveLES(); - // write back result to input/output vector - // c = x_vector.segment(!left_is_constant, c.size()); -} +// // write back result to input/output vector +// // c = x_vector.segment(!left_is_constant, c.size()); +// } inline void Diffusion::BTCSDiffusion::reserveMemory(int size, int max_count_per_line) { diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 72a54b9..ff64189 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -115,9 +115,10 @@ private: Eigen::RowMajor> BCVectorRowMajor; - void simulate_base(DVectorRowMajor &c, Eigen::Map &bc, - Eigen::Map &alpha, double dx, - double time_step, int size, DVectorRowMajor &t0_c); + // 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, From 9ec382877ec1879a9fdb24bd61946e6a840eaf50 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:19:06 +0100 Subject: [PATCH 06/12] Fix function parameters. - Use private data types instead of plain Eigen types --- src/BTCSDiffusion.cpp | 6 +++--- src/BTCSDiffusion.hpp | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 71ff04b..f18916b 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -251,8 +251,8 @@ void Diffusion::BTCSDiffusion::simulate2D( } inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( - const Eigen::VectorXd &alpha, - const Eigen::Vector &bc, + const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, double time_step) { Diffusion::boundary_condition left = bc[0]; @@ -289,7 +289,7 @@ inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( } inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI( - const DVectorRowMajor &c, const Eigen::VectorXd alpha, + const DVectorRowMajor &c, const DVectorRowMajor alpha, const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, double dx, double time_step) { diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index ff64189..fcbab40 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -128,12 +128,12 @@ private: Eigen::Map &alpha, Eigen::Map &bc); - inline void fillMatrixFromRow( - const Eigen::VectorXd &alpha, - const Eigen::Vector &bc, - int size, double dx, double time_step); + inline void fillMatrixFromRow(const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, + double time_step); + inline void fillVectorFromRowADI(const DVectorRowMajor &c, - const Eigen::VectorXd alpha, + const DVectorRowMajor alpha, const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, double dx, double time_step); From a5a66f1403ba1039befce5b933d841c1414efdbc Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:25:13 +0100 Subject: [PATCH 07/12] Update: pass additional information as parameter to simulate1D --- src/BTCSDiffusion.cpp | 10 ++++------ src/BTCSDiffusion.hpp | 3 ++- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index f18916b..555b213 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -159,11 +159,8 @@ inline void Diffusion::BTCSDiffusion::reserveMemory(int 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; + Eigen::Map &bc, int size, double dx, + double time_step) { reserveMemory(size, BTCS_MAX_DEP_PER_CELL); @@ -345,7 +342,8 @@ 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, alpha_in, bc_in); + simulate1D(c_in, alpha_in, bc_in, this->grid_cells[0], this->deltas[0], + this->time_step); } 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 fcbab40..f6ce0f2 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -122,7 +122,8 @@ private: void simulate1D(Eigen::Map &c, Eigen::Map &alpha, - Eigen::Map &bc); + Eigen::Map &bc, int size, double dx, + double time_step); void simulate2D(Eigen::Map &c, Eigen::Map &alpha, From fb5ee6431e15e8403c0d455397419a3dece423f6 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 11:38:08 +0100 Subject: [PATCH 08/12] Update: also pass t0_c to simulate_1D --- src/BTCSDiffusion.cpp | 12 ++++++------ src/BTCSDiffusion.hpp | 3 ++- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 555b213..d8c875f 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -159,14 +159,13 @@ inline void Diffusion::BTCSDiffusion::reserveMemory(int size, void Diffusion::BTCSDiffusion::simulate1D( Eigen::Map &c, Eigen::Map &alpha, - Eigen::Map &bc, int size, double dx, - double time_step) { + Eigen::Map &bc, const DVectorRowMajor &t0_c, + int size, double dx, double 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); + fillVectorFromRowADI(c, alpha, bc, t0_c, size, dx, time_step); solveLES(); @@ -342,8 +341,9 @@ 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, alpha_in, bc_in, this->grid_cells[0], this->deltas[0], - this->time_step); + simulate1D(c_in, alpha_in, bc_in, + Eigen::VectorXd::Constant(this->grid_cells[0], 0), + this->grid_cells[0], this->deltas[0], this->time_step); } 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 f6ce0f2..e17137b 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -122,7 +122,8 @@ private: void simulate1D(Eigen::Map &c, Eigen::Map &alpha, - Eigen::Map &bc, int size, double dx, + Eigen::Map &bc, + const DVectorRowMajor &t0_c, int size, double dx, double time_step); void simulate2D(Eigen::Map &c, From d0072f9f32d3ec2bf96617637408c779ee969b21 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 13:01:03 +0100 Subject: [PATCH 09/12] Revert to commit d65fcd4. --- src/BTCSDiffusion.cpp | 130 +++++++++++++++++++++--------------------- src/BTCSDiffusion.hpp | 12 ++-- 2 files changed, 69 insertions(+), 73 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index d8c875f..5da512c 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -77,74 +77,74 @@ void Diffusion::BTCSDiffusion::updateInternals() { } } -// void Diffusion::BTCSDiffusion::simulate_base( -// DVectorRowMajor &c, Eigen::Map &bc, -// Eigen::Map &alpha, double dx, double time_step, -// int size, DVectorRowMajor &t0_c) { +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; -// // ; + // 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); + // 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. -// // */ + // /* + // * 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 + 1, size + 1) = 1; -// // b_vector[size + 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++) { + // 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; -// // } + // // 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); + fillMatrixFromRow(alpha, bc, size, dx, time_step); + fillVectorFromRowADI(c, alpha, bc, t0_c, size, dx, time_step); -// solveLES(); + solveLES(); -// // write back result to input/output vector -// // c = x_vector.segment(!left_is_constant, c.size()); -// } + // write back result to input/output vector + // c = x_vector.segment(!left_is_constant, c.size()); +} inline void Diffusion::BTCSDiffusion::reserveMemory(int size, int max_count_per_line) { @@ -159,13 +159,17 @@ inline void Diffusion::BTCSDiffusion::reserveMemory(int size, void Diffusion::BTCSDiffusion::simulate1D( Eigen::Map &c, Eigen::Map &alpha, - Eigen::Map &bc, const DVectorRowMajor &t0_c, - int size, double dx, double time_step) { + 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, t0_c, size, dx, time_step); + fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, + dx, time_step); solveLES(); @@ -247,9 +251,8 @@ void Diffusion::BTCSDiffusion::simulate2D( } inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( - const DVectorRowMajor &alpha, - const BCVectorRowMajor &bc, - int size, double dx, double time_step) { + 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]; @@ -325,8 +328,7 @@ inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI( } if (!right_constant) { - b_vector[b_size - 1] = - getBCFromFlux(right, c[size - 1], alpha[size - 1]); + b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]); } } @@ -341,9 +343,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, alpha_in, bc_in, - Eigen::VectorXd::Constant(this->grid_cells[0], 0), - this->grid_cells[0], this->deltas[0], this->time_step); + 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 e17137b..7097fa2 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -115,16 +115,13 @@ private: Eigen::RowMajor> BCVectorRowMajor; - // void simulate_base(DVectorRowMajor &c, Eigen::Map - // &bc, - // Eigen::Map &alpha, double dx, - // double time_step, int size, DVectorRowMajor &t0_c); + 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, - const DVectorRowMajor &t0_c, int size, double dx, - double time_step); + Eigen::Map &bc); void simulate2D(Eigen::Map &c, Eigen::Map &alpha, @@ -133,7 +130,6 @@ private: 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, From 9d3ee1f913029f7570510f25c35b63ccf2ec1612 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 13:13:33 +0100 Subject: [PATCH 10/12] Use simulate_base for actual solving of les. --- src/BTCSDiffusion.cpp | 85 ++++++++----------------------------------- src/BTCSDiffusion.hpp | 6 +-- 2 files changed, 18 insertions(+), 73 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 5da512c..2f6b477 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -76,74 +76,22 @@ 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::simulate_base( - DVectorRowMajor &c, Eigen::Map &bc, - Eigen::Map &alpha, double dx, double time_step, - int size, DVectorRowMajor &t0_c) { + reserveMemory(size, BTCS_MAX_DEP_PER_CELL); - // 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 + 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++) { - - // // 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, bc, size, dx, time_step); - fillVectorFromRowADI(c, alpha, bc, t0_c, size, dx, time_step); + 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, @@ -165,15 +113,12 @@ void Diffusion::BTCSDiffusion::simulate1D( double dx = this->deltas[0]; double time_step = this->time_step; - reserveMemory(size, BTCS_MAX_DEP_PER_CELL); + DVectorRowMajor input_field = c.row(0); - fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step); - fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, - dx, time_step); + simulate_base(input_field, bc, alpha, dx, time_step, size, + Eigen::VectorXd::Constant(size, 0)); - solveLES(); - - c = x_vector.segment(1, size); + c.row(0) << input_field; } void Diffusion::BTCSDiffusion::simulate2D( diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 7097fa2..6657f97 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -115,9 +115,9 @@ private: Eigen::RowMajor> BCVectorRowMajor; - void simulate_base(DVectorRowMajor &c, Eigen::Map &bc, - Eigen::Map &alpha, double dx, - double time_step, int size, DVectorRowMajor &t0_c); + 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, Eigen::Map &alpha, From 9c1afe8e2d46313d22f8085e0c2155253322d89e Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 13:56:32 +0100 Subject: [PATCH 11/12] Implement 2D-row-wise in one direction --- src/BTCSDiffusion.cpp | 86 +++++++++++++++++++++++++++++++++---------- src/BTCSDiffusion.hpp | 8 +++- 2 files changed, 72 insertions(+), 22 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 2f6b477..a015eb2 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -125,32 +126,28 @@ void Diffusion::BTCSDiffusion::simulate2D( Eigen::Map &c, Eigen::Map &alpha, Eigen::Map &bc) { - // double local_dt = this->time_step / 2.; - // DMatrixRowMajor tmp_vector; + int n_rows, n_cols; + + double local_dt = this->time_step / 2.; + double dx = this->deltas[0]; + DMatrixRowMajor tmp_vector; + + n_rows = this->grid_cells[1]; + n_cols = this->grid_cells[0]; - // 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)); + DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); - // b_vector.resize(size); - // x_vector.resize(size); + std::cout << t0_c.row(0) << std::endl; - // 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 < 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; + } - // 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(); @@ -195,6 +192,55 @@ void Diffusion::BTCSDiffusion::simulate2D( // c.transposeInPlace(); } +Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c( + const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, + const BCMatrixRowMajor &bc, double time_step, double dx) { + + int n_rows = this->grid_cells[1]; + int n_cols = this->grid_cells[0]; + + 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) { diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 6657f97..1c9556a 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -116,8 +116,8 @@ private: BCVectorRowMajor; void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc, - const DVectorRowMajor &alpha, double dx, - double time_step, int size, const DVectorRowMajor &t0_c); + const DVectorRowMajor &alpha, double dx, double time_step, + int size, const DVectorRowMajor &t0_c); void simulate1D(Eigen::Map &c, Eigen::Map &alpha, @@ -127,6 +127,10 @@ private: Eigen::Map &alpha, Eigen::Map &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); From d4a87261518f8b141afedb31faea8e2e497502bd Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 1 Mar 2022 14:05:23 +0100 Subject: [PATCH 12/12] Implement 2D-row-wise in both directions --- src/BTCSDiffusion.cpp | 71 ++++++++++++------------------------------- 1 file changed, 19 insertions(+), 52 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index a015eb2..4d0d195 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -127,69 +127,36 @@ void Diffusion::BTCSDiffusion::simulate2D( Eigen::Map &bc) { int n_rows, n_cols; + double dx; + DMatrixRowMajor t0_c; double local_dt = this->time_step / 2.; - double dx = this->deltas[0]; + dx = this->deltas[0]; DMatrixRowMajor tmp_vector; n_rows = this->grid_cells[1]; n_cols = this->grid_cells[0]; - // unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]); - - DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); - - std::cout << t0_c.row(0) << std::endl; + t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); 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)); + 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)); - // } - - // 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(); + 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(); + } } Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c( @@ -221,8 +188,8 @@ Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c( 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); + t0_c(i, j) = time_step * alpha(i, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); } } @@ -230,12 +197,12 @@ Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c( // and finally over last row for (int j = 0; j < n_cols; j++) { - y_values[0] = c(end-1,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); + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); } return t0_c;