max/TugJulia
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added new function simulate_base.

Browse Source 
- 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.
This commit is contained in:
Max Luebke 2022-02-28 15:09:46 +01:00
parent 6408fd89fe
commit 6f9d344cee
2 changed files with 113 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
src/BTCSDiffusion.cpp
View File

@ -1,4 +1,5 @@
#include "BTCSDiffusion.hpp" #include "BTCSDiffusion.hpp"
#include "BoundaryCondition.hpp"
#include <Eigen/SparseLU> #include <Eigen/SparseLU>
@ -74,71 +75,73 @@ void Diffusion::BTCSDiffusion::updateInternals() {
} }
} }
void Diffusion::BTCSDiffusion::simulate1D( void Diffusion::BTCSDiffusion::simulate_base(
Eigen::Map<DVectorRowMajor> &c, Diffusion::boundary_condition left, DVectorRowMajor &c, Eigen::Map<const BCVectorRowMajor> &bc,
Diffusion::boundary_condition right, Eigen::Map<const BCVectorRowMajor> &bc, Eigen::Map<const DVectorRowMajor> &alpha, double dx, double time_step,
Eigen::Map<const DVectorRowMajor> &alpha, double dx, int size) { int size, DVectorRowMajor &t0_c) {
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 // 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. // 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 // Then we will need ghost nodes. So this variable will give the count of
// ghost nodes. // 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 // set sizes of private and yet allocated vectors
b_vector.resize(size + bc_offset); // b_vector.resize(size + bc_offset);
x_vector.resize(size + bc_offset); // x_vector.resize(size + bc_offset);
/* // /*
* Begin to solve the equation system using LU solver of Eigen. // * Begin to solve the equation system using LU solver of Eigen.
* // *
* But first fill the A matrix and b vector. // * But first fill the A matrix and b vector.
*/ // */
// Set boundary condition for ghost nodes (for closed or flux system) or outer // // Set boundary condition for ghost nodes (for closed or flux system) or
// inlet nodes (constant boundary condition) // outer
A_matrix.resize(size + bc_offset, size + bc_offset); // // inlet nodes (constant boundary condition)
A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3)); // 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; // A_matrix.insert(0, 0) = 1;
b_vector[0] = // b_vector[0] =
(left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); // (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1; // A_matrix.insert(size + 1, size + 1) = 1;
b_vector[size + bc_offset - 1] = // b_vector[size + 1] =
(right_is_constant ? right.value // (right_is_constant ? right.value
: getBCFromFlux(right, c[size - 1], alpha[size - 1])); // : getBCFromFlux(right, c[size - 1], alpha[size -
// 1]));
// Start filling the A matrix // Start filling the A matrix
// =i= is used for equation system matrix and vector indexing // =i= is used for equation system matrix and vector indexing
// and =j= for indexing of c,alpha and bc // and =j= for indexing of c,alpha and bc
for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant; // for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant;
i++, j++) { // i++, j++) {
// if current grid cell is considered as constant boundary conditon // // if current grid cell is considered as constant boundary conditon
if (bc[j].type == Diffusion::BC_CONSTANT) { // if (bc[j].type == Diffusion::BC_CONSTANT) {
A_matrix.insert(i, i) = 1; // A_matrix.insert(i, i) = 1;
b_vector[i] = bc[j].value; // b_vector[i] = bc[j].value;
continue; // 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. - 2. * sx;
A_matrix.insert(i, i - 1) = sx; // A_matrix.insert(i, i - 1) = 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(); solveLES();
// write back result to input/output vector // 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( void Diffusion::BTCSDiffusion::simulate2D(
@ -165,7 +168,8 @@ void Diffusion::BTCSDiffusion::simulate2D(
fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant, fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant,
deltas[0], this->time_step / 2, bc.row(i)); 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(); solveLES();
@ -197,7 +201,8 @@ void Diffusion::BTCSDiffusion::simulate2D(
fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant, fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant,
deltas[1], this->time_step / 2, bc.col(i)); 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(); solveLES();
@ -211,81 +216,85 @@ void Diffusion::BTCSDiffusion::simulate2D(
c.transposeInPlace(); c.transposeInPlace();
} }
void Diffusion::BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, inline void Diffusion::BTCSDiffusion::fillMatrixFromRow(
int n_cols, int row, const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size,
bool left_constant, double dx, double time_step) {
bool right_constant,
double delta, double time_step,
const BCVectorRowMajor &bc) {
n_cols += 2; Diffusion::boundary_condition left = bc[0];
int offset = n_cols * row; 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) 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) 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++) { 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) { if (bc[k].type == Diffusion::BC_CONSTANT) {
A_matrix.insert(offset + j, offset + j) = 1; A_matrix.insert(j, j) = 1;
continue; continue;
} }
A_matrix.insert(offset + j, offset + j) = -1. - 2. * sx; A_matrix.insert(j, j) = -1. - 2. * sx;
A_matrix.insert(offset + j, offset + (j - 1)) = sx; A_matrix.insert(j, (j - 1)) = sx;
A_matrix.insert(offset + j, offset + (j + 1)) = sx; A_matrix.insert(j, (j + 1)) = sx;
} }
} }
void Diffusion::BTCSDiffusion::fillVectorFromRowADI( inline void Diffusion::BTCSDiffusion::fillVectorFromRowADI(
Eigen::Map<DMatrixRowMajor> &c, const Eigen::VectorXd alpha, int row, DVectorRowMajor &c, const Eigen::VectorXd alpha, const BCVectorRowMajor &bc,
double delta, boundary_condition left, boundary_condition right, DVectorRowMajor &t0_c, int size, double dx, double time_step) {
double time_step, const BCVectorRowMajor &bc) {
int ncol = c.cols(); Diffusion::boundary_condition left = bc[0];
int nrow = c.rows(); Diffusion::boundary_condition right = bc[size - 1];
int offset = ncol + 2;
if (left.type != Diffusion::BC_CONSTANT) { bool left_constant = (left.type == Diffusion::BC_CONSTANT);
// this is not correct currently.We will fix this when we are able to define bool right_constant = (right.type == Diffusion::BC_CONSTANT);
// FLUX boundary conditions
b_vector[offset * row] = getBCFromFlux(left, c(row, 0), alpha[0]);
}
if (right.type != Diffusion::BC_CONSTANT) { int b_size = b_vector.size();
b_vector[offset * row + (offset - 1)] =
getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]);
}
for (int j = 0; j < ncol; j++) { for (int j = 0; j < size; j++) {
boundary_condition tmp_bc = bc[j]; boundary_condition tmp_bc = bc[j];
if (tmp_bc.type == Diffusion::BC_CONSTANT) { if (tmp_bc.type == Diffusion::BC_CONSTANT) {
b_vector[offset * row + (j + 1)] = tmp_bc.value; b_vector[j + 1] = tmp_bc.value;
continue; continue;
} }
double y_values[3]; // double y_values[3];
y_values[0] = // y_values[0] =
(row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j), alpha[j])); // (row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j),
y_values[1] = c(row, j); // alpha[j]));
y_values[2] = // y_values[1] = c(row, j);
(row != nrow - 1 ? c(row + 1, j) // y_values[2] =
: getBCFromFlux(tmp_bc, c(row, j), alpha[j])); // (row != nrow - 1 ? c(row + 1, j)
// : getBCFromFlux(tmp_bc, c(row, j), alpha[j]));
double t0_c = double t0_c_j = time_step * alpha[j] * (t0_c[j] / (dx * dx));
time_step * alpha[j] * b_vector[j + 1] = -c[j] - t0_c_j;
((y_values[0] - 2 * y_values[1] + y_values[2]) / (delta * delta)); }
b_vector[offset * row + (j + 1)] = -c(row, j) - (t0_c);
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]);
} }
} }

23
src/BTCSDiffusion.hpp
View File

@ -7,6 +7,7 @@
#include <Eigen/src/Core/Map.h> #include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h> #include <Eigen/src/Core/util/Constants.h>
#include <cstddef>
#include <tuple> #include <tuple>
#include <type_traits> #include <type_traits>
#include <vector> #include <vector>
@ -112,24 +113,22 @@ private:
Eigen::RowMajor> Eigen::RowMajor>
BCVectorRowMajor; BCVectorRowMajor;
void simulate1D(Eigen::Map<DVectorRowMajor> &c, void simulate_base(DVectorRowMajor &c,
Diffusion::boundary_condition left,
Diffusion::boundary_condition right,
Eigen::Map<const BCVectorRowMajor> &bc, Eigen::Map<const BCVectorRowMajor> &bc,
Eigen::Map<const DVectorRowMajor> &alpha, double dx, Eigen::Map<const DVectorRowMajor> &alpha, double dx,
int size); double time_step, int size, DVectorRowMajor &t0_c);
void simulate2D(Eigen::Map<DMatrixRowMajor> &c, void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
Eigen::Map<const DMatrixRowMajor> &alpha, Eigen::Map<const DMatrixRowMajor> &alpha,
Eigen::Map<const BCMatrixRowMajor> &bc); Eigen::Map<const BCMatrixRowMajor> &bc);
void fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, int row, inline void fillMatrixFromRow(const DVectorRowMajor &alpha,
bool left_constant, bool right_constant, double delta, const BCVectorRowMajor &bc, int size, double dx,
double time_step, const BCVectorRowMajor &bc); double time_step);
void fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c, inline void fillVectorFromRowADI(DVectorRowMajor &c,
const Eigen::VectorXd alpha, int row, double delta, const Eigen::VectorXd alpha,
Diffusion::boundary_condition left, const BCVectorRowMajor &bc,
Diffusion::boundary_condition right, DVectorRowMajor &t0_c, int size, double dx,
double time_step, const BCVectorRowMajor &bc); double time_step);
void simulate3D(std::vector<double> &c); void simulate3D(std::vector<double> &c);
inline double getBCFromFlux(Diffusion::boundary_condition bc, inline double getBCFromFlux(Diffusion::boundary_condition bc,
double nearest_value, double neighbor_alpha); double nearest_value, double neighbor_alpha);