mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 17:38:23 +01:00
Merge branch 'row_wise_2D' into '2D'
Implemented diffusion row-wise. See merge request mluebke/diffusion!8
This commit is contained in:
Max Lübke
commit
3f1747fba9
@ -42,4 +42,4 @@ lint:
|
||||
- make
|
||||
only:
|
||||
refs:
|
||||
- master
|
||||
- main
|
||||
|
||||
@ -1,4 +1,5 @@
|
||||
#include "BTCSDiffusion.hpp"
|
||||
#include "BoundaryCondition.hpp"
|
||||
|
||||
#include <Eigen/SparseLU>
|
||||
|
||||
@ -6,6 +7,7 @@
|
||||
#include <Eigen/src/Core/Matrix.h>
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
#include <cstddef>
|
||||
#include <iomanip>
|
||||
#include <iterator>
|
||||
#include <ostream>
|
||||
@ -14,6 +16,8 @@
|
||||
|
||||
#include <iostream>
|
||||
|
||||
#define BTCS_MAX_DEP_PER_CELL 3
|
||||
|
||||
Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
|
||||
assert(dim <= 3);
|
||||
|
||||
@ -73,219 +77,216 @@ 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::simulate1D(
|
||||
Eigen::Map<DVectorRowMajor> &c, Diffusion::boundary_condition left,
|
||||
Diffusion::boundary_condition right, Eigen::Map<const BCVectorRowMajor> &bc,
|
||||
Eigen::Map<const DVectorRowMajor> &alpha, double dx, int size) {
|
||||
reserveMemory(size, BTCS_MAX_DEP_PER_CELL);
|
||||
|
||||
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
|
||||
// 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 + 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]));
|
||||
|
||||
// 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.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,
|
||||
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<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
|
||||
Eigen::Map<const BCVectorRowMajor> &bc) {
|
||||
|
||||
int size = this->grid_cells[0];
|
||||
double dx = this->deltas[0];
|
||||
double time_step = this->time_step;
|
||||
|
||||
DVectorRowMajor input_field = c.row(0);
|
||||
|
||||
simulate_base(input_field, bc, alpha, dx, time_step, size,
|
||||
Eigen::VectorXd::Constant(size, 0));
|
||||
|
||||
c.row(0) << input_field;
|
||||
}
|
||||
|
||||
void Diffusion::BTCSDiffusion::simulate2D(
|
||||
Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
|
||||
Eigen::Map<const BCMatrixRowMajor> &bc) {
|
||||
|
||||
int n_rows, n_cols;
|
||||
double dx;
|
||||
DMatrixRowMajor t0_c;
|
||||
|
||||
double local_dt = this->time_step / 2.;
|
||||
dx = this->deltas[0];
|
||||
DMatrixRowMajor tmp_vector;
|
||||
|
||||
int n_cols = c.cols();
|
||||
unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]);
|
||||
n_rows = this->grid_cells[1];
|
||||
n_cols = this->grid_cells[0];
|
||||
|
||||
A_matrix.resize(size, size);
|
||||
A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
|
||||
t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
|
||||
|
||||
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;
|
||||
|
||||
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));
|
||||
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;
|
||||
}
|
||||
|
||||
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<Eigen::MatrixXd> 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));
|
||||
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();
|
||||
}
|
||||
|
||||
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();
|
||||
}
|
||||
|
||||
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) {
|
||||
Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c(
|
||||
const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
|
||||
const BCMatrixRowMajor &bc, double time_step, double dx) {
|
||||
|
||||
n_cols += 2;
|
||||
int offset = n_cols * row;
|
||||
int n_rows = this->grid_cells[1];
|
||||
int n_cols = this->grid_cells[0];
|
||||
|
||||
A_matrix.insert(offset, offset) = 1;
|
||||
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) {
|
||||
|
||||
Diffusion::boundary_condition left = bc[0];
|
||||
Diffusion::boundary_condition right = bc[size - 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; k < size - 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<DMatrixRowMajor> &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(
|
||||
const DVectorRowMajor &c, const DVectorRowMajor alpha,
|
||||
const BCVectorRowMajor &bc, const 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, c[0], alpha[0]);
|
||||
}
|
||||
|
||||
if (!right_constant) {
|
||||
b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]);
|
||||
}
|
||||
}
|
||||
|
||||
@ -300,8 +301,7 @@ void Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
|
||||
Eigen::Map<const DVectorRowMajor> alpha_in(alpha, this->grid_cells[0]);
|
||||
Eigen::Map<const BCVectorRowMajor> 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<DMatrixRowMajor> c_in(c, this->grid_cells[1],
|
||||
|
||||
@ -7,10 +7,13 @@
|
||||
#include <Eigen/src/Core/Map.h>
|
||||
#include <Eigen/src/Core/Matrix.h>
|
||||
#include <Eigen/src/Core/util/Constants.h>
|
||||
#include <cstddef>
|
||||
#include <tuple>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
#define BTCS_MAX_DEP_PER_CELL 3
|
||||
|
||||
namespace Diffusion {
|
||||
/*!
|
||||
* Class implementing a solution for a 1/2/3D diffusion equation using backward
|
||||
@ -112,25 +115,33 @@ private:
|
||||
Eigen::RowMajor>
|
||||
BCVectorRowMajor;
|
||||
|
||||
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<DVectorRowMajor> &c,
|
||||
Diffusion::boundary_condition left,
|
||||
Diffusion::boundary_condition right,
|
||||
Eigen::Map<const BCVectorRowMajor> &bc,
|
||||
Eigen::Map<const DVectorRowMajor> &alpha, double dx,
|
||||
int size);
|
||||
Eigen::Map<const DVectorRowMajor> &alpha,
|
||||
Eigen::Map<const BCVectorRowMajor> &bc);
|
||||
|
||||
void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
|
||||
Eigen::Map<const DMatrixRowMajor> &alpha,
|
||||
Eigen::Map<const BCMatrixRowMajor> &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<DMatrixRowMajor> &c,
|
||||
const Eigen::VectorXd alpha, int row, double delta,
|
||||
Diffusion::boundary_condition left,
|
||||
Diffusion::boundary_condition right,
|
||||
double time_step, const BCVectorRowMajor &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);
|
||||
inline void fillVectorFromRowADI(const DVectorRowMajor &c,
|
||||
const DVectorRowMajor alpha,
|
||||
const BCVectorRowMajor &bc,
|
||||
const DVectorRowMajor &t0_c, int size,
|
||||
double dx, double time_step);
|
||||
void simulate3D(std::vector<double> &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();
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user