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

Implement 2D-row-wise in one direction

Browse Source
This commit is contained in:
Max Luebke 2022-03-01 13:56:32 +01:00
parent 9d3ee1f913
commit 9c1afe8e2d
2 changed files with 72 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
src/BTCSDiffusion.cpp
View File

@ -7,6 +7,7 @@
#include <Eigen/src/Core/Matrix.h>
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <iomanip>
#include <iterator>
#include <ostream>
@ -125,32 +126,28 @@ void Diffusion::BTCSDiffusion::simulate2D(
Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
Eigen::Map<const BCMatrixRowMajor> &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) {

8
src/BTCSDiffusion.hpp
View File

@ -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<DVectorRowMajor> &c,
Eigen::Map<const DVectorRowMajor> &alpha,
@ -127,6 +127,10 @@ private:
Eigen::Map<const DMatrixRowMajor> &alpha,
Eigen::Map<const BCMatrixRowMajor> &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);