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

Merge branch 'func_signature' into '2D'

Browse Source 
Refactor simulate function signature

See merge request mluebke/diffusion!7
This commit is contained in:
Max Lübke 2022-02-28 14:10:53 +01:00
commit 0ebdbe638d
5 changed files with 164 additions and 140 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

167
src/BTCSDiffusion.cpp
View File

@ -14,11 +14,7 @@
#include <iostream>
const int BTCSDiffusion::BC_CLOSED = 0;
const int BTCSDiffusion::BC_FLUX = 1;
const int BTCSDiffusion::BC_CONSTANT = 2;
BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
assert(dim <= 3);
grid_cells.resize(dim, 1);
@ -26,8 +22,8 @@ BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
deltas.resize(dim, 1);
}
void BTCSDiffusion::setXDimensions(double domain_size,
unsigned int n_grid_cells) {
void Diffusion::BTCSDiffusion::setXDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 0);
this->domain_size[0] = domain_size;
this->grid_cells[0] = n_grid_cells;
@ -35,8 +31,8 @@ void BTCSDiffusion::setXDimensions(double domain_size,
updateInternals();
}
void BTCSDiffusion::setYDimensions(double domain_size,
unsigned int n_grid_cells) {
void Diffusion::BTCSDiffusion::setYDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 1);
this->domain_size[1] = domain_size;
this->grid_cells[1] = n_grid_cells;
@ -44,8 +40,8 @@ void BTCSDiffusion::setYDimensions(double domain_size,
updateInternals();
}
void BTCSDiffusion::setZDimensions(double domain_size,
unsigned int n_grid_cells) {
void Diffusion::BTCSDiffusion::setZDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 2);
this->domain_size[2] = domain_size;
this->grid_cells[2] = n_grid_cells;
@ -53,29 +49,38 @@ void BTCSDiffusion::setZDimensions(double domain_size,
updateInternals();
}
unsigned int BTCSDiffusion::getXGridCellsN() { return this->grid_cells[0]; }
unsigned int BTCSDiffusion::getYGridCellsN() { return this->grid_cells[1]; }
unsigned int BTCSDiffusion::getZGridCellsN() { return this->grid_cells[2]; }
unsigned int BTCSDiffusion::getXDomainSize() { return this->domain_size[0]; }
unsigned int BTCSDiffusion::getYDomainSize() { return this->domain_size[1]; }
unsigned int BTCSDiffusion::getZDomainSize() { return this->domain_size[2]; }
unsigned int Diffusion::BTCSDiffusion::getXGridCellsN() {
return this->grid_cells[0];
}
unsigned int Diffusion::BTCSDiffusion::getYGridCellsN() {
return this->grid_cells[1];
}
unsigned int Diffusion::BTCSDiffusion::getZGridCellsN() {
return this->grid_cells[2];
}
unsigned int Diffusion::BTCSDiffusion::getXDomainSize() {
return this->domain_size[0];
}
unsigned int Diffusion::BTCSDiffusion::getYDomainSize() {
return this->domain_size[1];
}
unsigned int Diffusion::BTCSDiffusion::getZDomainSize() {
return this->domain_size[2];
}
void BTCSDiffusion::updateInternals() {
void Diffusion::BTCSDiffusion::updateInternals() {
for (int i = 0; i < grid_dim; i++) {
deltas[i] = (double)domain_size[i] / grid_cells[i];
}
bc.resize((grid_dim > 1 ? grid_cells[1] : 1), grid_cells[0]);
}
void BTCSDiffusion::simulate1D(Eigen::Map<DVectorRowMajor> &c,
boundary_condition left,
boundary_condition right,
const std::vector<double> &alpha, double dx,
int size) {
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) {
bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
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.
@ -115,9 +120,9 @@ void BTCSDiffusion::simulate1D(Eigen::Map<DVectorRowMajor> &c,
i++, j++) {
// if current grid cell is considered as constant boundary conditon
if (bc(1,j).type == BTCSDiffusion::BC_CONSTANT) {
if (bc[j].type == Diffusion::BC_CONSTANT) {
A_matrix.insert(i, i) = 1;
b_vector[i] = bc(1,j).value;
b_vector[i] = bc[j].value;
continue;
}
@ -136,10 +141,11 @@ void BTCSDiffusion::simulate1D(Eigen::Map<DVectorRowMajor> &c,
c = x_vector.segment(!left_is_constant, c.size());
}
void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
Eigen::Map<const DMatrixRowMajor> &alpha) {
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.;
double local_dt = this->time_step / 2.;
DMatrixRowMajor tmp_vector;
int n_cols = c.cols();
@ -153,13 +159,13 @@ void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
for (int i = 0; i < c.rows(); i++) {
boundary_condition left = bc(i, 0);
bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT;
bool left_constant = left.type == Diffusion::BC_CONSTANT;
boundary_condition right = bc(i, n_cols - 1);
bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT;
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);
fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt);
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();
@ -184,14 +190,14 @@ void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
n_cols = c.cols();
for (int i = 0; i < c.rows(); i++) {
boundary_condition left = bc(0,i);
bool left_constant = left.type == BTCSDiffusion::BC_CONSTANT;
boundary_condition right = bc(n_cols-1,i);
bool right_constant = right.type == BTCSDiffusion::BC_CONSTANT;
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);
fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt);
deltas[1], this->time_step / 2, bc.col(i));
fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt, bc.col(i));
}
solveLES();
@ -205,10 +211,12 @@ void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
c.transposeInPlace();
}
void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols,
int row, bool left_constant,
bool right_constant, double delta,
double time_step) {
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) {
n_cols += 2;
int offset = n_cols * row;
@ -227,7 +235,7 @@ void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols,
j++, k++) {
double sx = (alpha[j - 1] * time_step) / (delta * delta);
if (this->bc(row, k).type == BTCSDiffusion::BC_CONSTANT) {
if (bc[k].type == Diffusion::BC_CONSTANT) {
A_matrix.insert(offset + j, offset + j) = 1;
continue;
}
@ -238,32 +246,31 @@ void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols,
}
}
void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
const Eigen::VectorXd alpha, int row,
double delta, boundary_condition left,
boundary_condition right,
double time_step) {
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) {
int ncol = c.cols();
int nrow = c.rows();
int offset = ncol + 2;
if (left.type != BTCSDiffusion::BC_CONSTANT) {
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]);
}
if (right.type != BTCSDiffusion::BC_CONSTANT) {
if (right.type != Diffusion::BC_CONSTANT) {
b_vector[offset * row + (offset - 1)] =
getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]);
}
for (int j = 0; j < ncol; j++) {
boundary_condition tmp_bc = this->bc(row, j);
boundary_condition tmp_bc = bc[j];
if (tmp_bc.type == BTCSDiffusion::BC_CONSTANT) {
b_vector[offset * row + (j+1)] = tmp_bc.value;
if (tmp_bc.type == Diffusion::BC_CONSTANT) {
b_vector[offset * row + (j + 1)] = tmp_bc.value;
continue;
}
@ -282,39 +289,43 @@ void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
}
}
void BTCSDiffusion::setTimestep(double time_step) {
void Diffusion::BTCSDiffusion::setTimestep(double time_step) {
this->time_step = time_step;
}
void BTCSDiffusion::simulate(std::vector<double> &c,
const std::vector<double> &alpha) {
void Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
Diffusion::boundary_condition *bc) {
if (this->grid_dim == 1) {
assert(c.size() == grid_cells[0]);
Eigen::Map<DVectorRowMajor> c_in(c.data(), this->grid_cells[0]);
simulate1D(c_in, bc(1,0), bc(1,bc.cols()-1), alpha, this->deltas[0],
this->grid_cells[0]);
Eigen::Map<DVectorRowMajor> c_in(c, this->grid_cells[0]);
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]);
}
if (this->grid_dim == 2) {
assert(c.size() == grid_cells[0] * grid_cells[1]);
Eigen::Map<DMatrixRowMajor> c_in(c.data(), this->grid_cells[1],
Eigen::Map<DMatrixRowMajor> c_in(c, this->grid_cells[1],
this->grid_cells[0]);
Eigen::Map<const DMatrixRowMajor> alpha_in(
alpha.data(), this->grid_cells[1], this->grid_cells[0]);
Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, this->grid_cells[1],
this->grid_cells[0]);
simulate2D(c_in, alpha_in);
Eigen::Map<const BCMatrixRowMajor> bc_in(bc, this->grid_cells[1],
this->grid_cells[0]);
simulate2D(c_in, alpha_in, bc_in);
}
}
inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
double neighbor_c,
double neighbor_alpha) {
inline double Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc,
double neighbor_c,
double neighbor_alpha) {
double val;
if (bc.type == BTCSDiffusion::BC_CLOSED) {
if (bc.type == Diffusion::BC_CLOSED) {
val = neighbor_c;
} else if (bc.type == BTCSDiffusion::BC_FLUX) {
} else if (bc.type == Diffusion::BC_FLUX) {
// TODO
// val = bc[index].value;
} else {
@ -324,13 +335,7 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
return val;
}
void BTCSDiffusion::setBoundaryCondition(int i, int j, bctype type, double value) {
bc(i,j).type = type;
bc(i,j).value = value;
}
inline void BTCSDiffusion::solveLES() {
inline void Diffusion::BTCSDiffusion::solveLES() {
// start to solve
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;

80
src/BTCSDiffusion.hpp
View File

@ -1,6 +1,8 @@
#ifndef BTCSDIFFUSION_H_
#define BTCSDIFFUSION_H_
#include "BoundaryCondition.hpp"
#include <Eigen/Sparse>
#include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h>
@ -9,11 +11,7 @@
#include <type_traits>
#include <vector>
/*!
* Defines both types of boundary condition as a datatype.
*/
typedef int bctype;
namespace Diffusion {
/*!
* Class implementing a solution for a 1/2/3D diffusion equation using backward
* euler.
@ -21,21 +19,6 @@ typedef int bctype;
class BTCSDiffusion {
public:
/*!
* Defines a constant/Dirichlet boundary condition.
*/
static const int BC_CONSTANT;
/*!
* Defines a closed/Neumann boundary condition.
*/
static const int BC_CLOSED;
/*!
* Defines a flux/Cauchy boundary condition.
*/
static const int BC_FLUX;
/*!
* Creates a diffusion module.
*
@ -108,7 +91,7 @@ public:
* continious memory (row major).
* @param alpha Vector of diffusion coefficients for each grid element.
*/
void simulate(std::vector<double> &c, const std::vector<double> &alpha);
void simulate(double *c, double *alpha, Diffusion::boundary_condition *bc);
/*!
* Set the timestep of the simulation
@ -117,53 +100,46 @@ public:
*/
void setTimestep(double time_step);
/*!
* Set the boundary condition of the given grid. This is done by defining an
* index (exact order still to be determined), the type of the boundary
* condition and the according value.
*
* @param index Index of the grid cell the boundary condition is applied to.
* @param type Type of the boundary condition. Must be constant, closed or
* flux.
* @param value For constant boundary conditions this value is set
* during solving. For flux value refers to a gradient of change for this grid
* cell. For closed this value has no effect since a gradient of 0 is used.
*/
void setBoundaryCondition(int row, int column, bctype type, double value);
private:
typedef struct boundary_condition {
bctype type;
double value;
} boundary_condition;
typedef Eigen::Triplet<double> T;
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
DMatrixRowMajor;
typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
DVectorRowMajor;
typedef Eigen::Matrix<Diffusion::boundary_condition, Eigen::Dynamic,
Eigen::Dynamic, Eigen::RowMajor>
BCMatrixRowMajor;
typedef Eigen::Matrix<Diffusion::boundary_condition, 1, Eigen::Dynamic,
Eigen::RowMajor>
BCVectorRowMajor;
void simulate1D(Eigen::Map<DVectorRowMajor> &c, boundary_condition left,
boundary_condition right, const std::vector<double> &alpha,
double dx, int size);
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);
void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
Eigen::Map<const DMatrixRowMajor> &alpha);
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);
double time_step, const BCVectorRowMajor &bc);
void fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
const Eigen::VectorXd alpha, int row, double delta,
boundary_condition left, boundary_condition right,
double time_step);
Diffusion::boundary_condition left,
Diffusion::boundary_condition right,
double time_step, const BCVectorRowMajor &bc);
void simulate3D(std::vector<double> &c);
inline double getBCFromFlux(boundary_condition bc, double nearest_value,
double neighbor_alpha);
inline double getBCFromFlux(Diffusion::boundary_condition bc,
double nearest_value, double neighbor_alpha);
void solveLES();
void updateInternals();
// std::vector<boundary_condition> bc;
Eigen::Matrix<boundary_condition, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> bc;
// Eigen::Matrix<boundary_condition, Eigen::Dynamic, Eigen::Dynamic,
// Eigen::RowMajor>
// bc;
Eigen::SparseMatrix<double> A_matrix;
Eigen::VectorXd b_vector;
@ -176,5 +152,5 @@ private:
std::vector<double> domain_size;
std::vector<double> deltas;
};
} // namespace Diffusion
#endif // BTCSDIFFUSION_H_

33
src/BoundaryCondition.hpp Normal file
View File

@ -0,0 +1,33 @@
#ifndef BOUNDARYCONDITION_H_
#define BOUNDARYCONDITION_H_
namespace Diffusion {
/*!
* Defines both types of boundary condition as a datatype.
*/
typedef int bctype;
/*!
* Defines a closed/Neumann boundary condition.
*/
static const bctype BC_CLOSED = 0;
/*!
* Defines a flux/Cauchy boundary condition.
*/
static const bctype BC_FLUX = 1;
/*!
* Defines a constant/Dirichlet boundary condition.
*/
static const bctype BC_CONSTANT = 2;
typedef struct boundary_condition {
bctype type;
double value;
} boundary_condition;
} // namespace Diffusion
#endif // BOUNDARYCONDITION_H_

12
src/main.cpp
View File

@ -1,10 +1,14 @@
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
@ -15,6 +19,7 @@ int main(int argc, char *argv[]) {
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n, 1 * pow(10, -1));
std::vector<double> field(n, 1 * std::pow(10, -6));
std::vector<boundary_condition> bc(n, {0,0});
// create instance of diffusion module
BTCSDiffusion diffu(dim);
@ -22,8 +27,9 @@ int main(int argc, char *argv[]) {
diffu.setXDimensions(1, n);
// set the boundary condition for the left ghost cell to dirichlet
diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT,
5. * std::pow(10, -6));
bc[0] = {Diffusion::BC_CONSTANT, 5*std::pow(10,-6)};
// diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT,
// 5. * std::pow(10, -6));
// set timestep for simulation to 1 second
diffu.setTimestep(1.);
@ -33,7 +39,7 @@ int main(int argc, char *argv[]) {
// loop 100 times
// output is currently generated by the method itself
for (int i = 0; i < 100; i++) {
diffu.simulate(field, alpha);
diffu.simulate(field.data(), alpha.data(), bc.data());
cout << "Iteration: " << i << "\n\n";

12
src/main_2D.cpp
View File

@ -1,10 +1,14 @@
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
@ -16,6 +20,7 @@ int main(int argc, char *argv[]) {
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n * m, 1 * pow(10, -1));
std::vector<double> field(n * m, 1 * std::pow(10, -6));
std::vector<boundary_condition> bc(n*m, {0,0});
// create instance of diffusion module
BTCSDiffusion diffu(dim);
@ -23,9 +28,8 @@ int main(int argc, char *argv[]) {
diffu.setXDimensions(1, n);
diffu.setYDimensions(1, m);
// set the boundary condition for the left ghost cell to dirichlet
diffu.setBoundaryCondition(2, 2, BTCSDiffusion::BC_CONSTANT,
5. * std::pow(10, -6));
//set inlet to higher concentration without setting bc
field[12] = 5*std::pow(10,-3);
// set timestep for simulation to 1 second
diffu.setTimestep(1.);
@ -33,7 +37,7 @@ int main(int argc, char *argv[]) {
cout << setprecision(12);
for (int t = 0; t < 10; t++) {
diffu.simulate(field, alpha);
diffu.simulate(field.data(), alpha.data(), bc.data());
cout << "Iteration: " << t << "\n\n";