max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 17:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Update simulation loop to new indexing

Browse Source
This commit is contained in:
Max Luebke 2022-01-21 14:09:30 +01:00
parent 5606b559c7
commit 2cd8a1e5d8
2 changed files with 49 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
src/BTCSDiffusion.cpp
View File

@ -79,17 +79,24 @@ void BTCSDiffusion::updateInternals() {
// // TODO: reserve memory for boundary conditions
// }
void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
double bc_right,
void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
boundary_condition right,
const std::vector<double> &alpha, double dx,
int size) {
bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
int loop_end = size + !right_is_constant;
// we need 2 more grid cells for ghost cells
// size = size + 2;
int bc_offset = !left_is_constant + !right_is_constant;
;
// set sizes of private and yet allocated vectors
b_vector.resize(size + 2);
x_vector.resize(size + 2);
b_vector.resize(size + bc_offset);
x_vector.resize(size + bc_offset);
/*
* Begin to solve the equation system using LU solver of Eigen.
@ -100,23 +107,27 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
* TODO: remove output
*/
A_matrix.resize(size + 2, size + 2);
A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, 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(size + 1, size + 1) = 1;
b_vector[0] = (left_is_constant ? c[0] : getBCFromFlux(left, c[0], alpha[0]));
b_vector[0] = bc_left;
b_vector[size + 1] = bc_right;
A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1;
b_vector[size + bc_offset - 1] =
(right_is_constant ? c[size - 1] : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
for (int i = 1; i < size + 1; i++) {
double sx = (alpha[i - 1] * time_step) / (dx * dx);
// A_matrix.insert(0, 0) = 1;
// A_matrix.insert(size + 1, size + 1) = 1;
for (int i = 1; i < size - right_is_constant; i++) {
double sx = (alpha[i - !(left_is_constant)] * 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[i - 1];
b_vector[i] = -c[i - !(left_is_constant)];
}
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@ -132,7 +143,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
for (int i = 0; i < c.size(); i++) {
c[i] = x_vector[i + 1];
c[i] = x_vector[i + !left_is_constant];
}
}
@ -143,27 +154,24 @@ void BTCSDiffusion::setTimestep(double time_step) {
void BTCSDiffusion::simulate(std::vector<double> &c,
const std::vector<double> &alpha) {
if (this->grid_dim == 1) {
double bc_left = getBCFromTuple(0, c[0], alpha[0]);
double bc_right =
getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
// double bc_left = getBCFromTuple(0, c[0], alpha[0]);
// double bc_right =
// getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
simulate1D(c, bc_left, bc_right, alpha, this->deltas[0],
this->grid_cells[0]);
simulate1D(c, bc[0], bc[1], alpha, this->deltas[0], this->grid_cells[0]);
}
}
double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
double neighbor_alpha) {
double val = -1;
int type = bc[index].type;
inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, double neighbor_c,
double neighbor_alpha) {
if (type == BTCSDiffusion::BC_CLOSED) {
double val;
if (bc.type == BTCSDiffusion::BC_CLOSED) {
val = neighbor_c;
// val = neighbor_c + (this->time_step / (this->deltas[0] *
// this->deltas[0])) *
// neighbor_alpha * std::get<1>(bc[index]);
} else if (type == BTCSDiffusion::BC_CONSTANT) {
val = bc[index].value;
} else if (bc.type == BTCSDiffusion::BC_FLUX) {
//TODO
// val = bc[index].value;
} else {
// TODO: implement error handling here. Type was set to wrong value.
}

26
src/BTCSDiffusion.hpp
View File

@ -5,7 +5,6 @@
#include <tuple>
#include <vector>
/*!
* Defines both types of boundary condition as a datatype.
*/
@ -42,10 +41,10 @@ public:
* - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
* the boundary condition
*/
typedef struct boundary_condition {
bctype type;
double value;
} boundary_condition;
typedef struct boundary_condition {
bctype type;
double value;
} boundary_condition;
/*!
* A boundary condition consists of two features. A type and the according
@ -58,11 +57,11 @@ public:
*/
// typedef std::vector<std::tuple<bctype, double>> boundary_condition;
/*!
* Datatype to fill the sparse matrix which is used to solve the equation
* system.
*/
typedef Eigen::Triplet<double> T;
/*!
* Datatype to fill the sparse matrix which is used to solve the equation
* system.
*/
typedef Eigen::Triplet<double> T;
/*!
* Create 1D-diffusion module.
@ -122,12 +121,13 @@ public:
void setBoundaryCondition(int index, double val, bctype type);
private:
void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
const std::vector<double> &alpha, double dx, int size);
void simulate1D(std::vector<double> &c, boundary_condition left,
boundary_condition right, const std::vector<double> &alpha,
double dx, int size);
void simulate2D(std::vector<double> &c);
void simulate3D(std::vector<double> &c);
double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
inline double getBCFromFlux(boundary_condition bc, double nearest_value, double neighbor_alpha);
void updateInternals();