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,