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

Refactor simulate function signature

Browse Source
This commit is contained in:
Max Lübke 2022-02-28 14:10:53 +01:00
parent 89e6345e7e
commit d7e240c6a8
5 changed files with 164 additions and 140 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
src/BTCSDiffusion.cpp
View File

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

80
src/BTCSDiffusion.hpp
View File

@ -1,6 +1,8 @@
#ifndef BTCSDIFFUSION_H_ #ifndef BTCSDIFFUSION_H_
#define BTCSDIFFUSION_H_ #define BTCSDIFFUSION_H_
#include "BoundaryCondition.hpp"
#include <Eigen/Sparse> #include <Eigen/Sparse>
#include <Eigen/src/Core/Map.h> #include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
@ -9,11 +11,7 @@
#include <type_traits> #include <type_traits>
#include <vector> #include <vector>
/*! namespace Diffusion {
* Defines both types of boundary condition as a datatype.
*/
typedef int bctype;
/*! /*!
* Class implementing a solution for a 1/2/3D diffusion equation using backward * Class implementing a solution for a 1/2/3D diffusion equation using backward
* euler. * euler.
@ -21,21 +19,6 @@ typedef int bctype;
class BTCSDiffusion { class BTCSDiffusion {
public: 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. * Creates a diffusion module.
* *
@ -108,7 +91,7 @@ public:
* continious memory (row major). * continious memory (row major).
* @param alpha Vector of diffusion coefficients for each grid element. * @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 * Set the timestep of the simulation
@ -117,53 +100,46 @@ public:
*/ */
void setTimestep(double time_step); 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: 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> typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
DMatrixRowMajor; DMatrixRowMajor;
typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor> typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
DVectorRowMajor; 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, void simulate1D(Eigen::Map<DVectorRowMajor> &c,
boundary_condition right, const std::vector<double> &alpha, Diffusion::boundary_condition left,
double dx, int size); Diffusion::boundary_condition right,
Eigen::Map<const BCVectorRowMajor> &bc,
Eigen::Map<const DVectorRowMajor> &alpha, double dx,
int size);
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);
void fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, int row, void fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols, int row,
bool left_constant, bool right_constant, double delta, bool left_constant, bool right_constant, double delta,
double time_step); double time_step, const BCVectorRowMajor &bc);
void fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c, void fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
const Eigen::VectorXd alpha, int row, double delta, const Eigen::VectorXd alpha, int row, double delta,
boundary_condition left, boundary_condition right, Diffusion::boundary_condition left,
double time_step); Diffusion::boundary_condition right,
double time_step, const BCVectorRowMajor &bc);
void simulate3D(std::vector<double> &c); void simulate3D(std::vector<double> &c);
inline double getBCFromFlux(boundary_condition bc, double nearest_value, inline double getBCFromFlux(Diffusion::boundary_condition bc,
double neighbor_alpha); double nearest_value, double neighbor_alpha);
void solveLES(); void solveLES();
void updateInternals(); void updateInternals();
// std::vector<boundary_condition> bc; // 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::SparseMatrix<double> A_matrix;
Eigen::VectorXd b_vector; Eigen::VectorXd b_vector;
@ -176,5 +152,5 @@ private:
std::vector<double> domain_size; std::vector<double> domain_size;
std::vector<double> deltas; std::vector<double> deltas;
}; };
} // namespace Diffusion
#endif // BTCSDIFFUSION_H_ #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 "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max #include <algorithm> // for copy, max
#include <cmath>
#include <iomanip> #include <iomanip>
#include <iostream> // for std #include <iostream> // for std
#include <vector> // for vector #include <vector> // for vector
using namespace std; using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
// dimension of grid // dimension of grid
@ -15,6 +19,7 @@ int main(int argc, char *argv[]) {
// create input + diffusion coefficients for each grid cell // create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n, 1 * pow(10, -1)); std::vector<double> alpha(n, 1 * pow(10, -1));
std::vector<double> field(n, 1 * std::pow(10, -6)); std::vector<double> field(n, 1 * std::pow(10, -6));
std::vector<boundary_condition> bc(n, {0,0});
// create instance of diffusion module // create instance of diffusion module
BTCSDiffusion diffu(dim); BTCSDiffusion diffu(dim);
@ -22,8 +27,9 @@ int main(int argc, char *argv[]) {
diffu.setXDimensions(1, n); diffu.setXDimensions(1, n);
// set the boundary condition for the left ghost cell to dirichlet // set the boundary condition for the left ghost cell to dirichlet
diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, bc[0] = {Diffusion::BC_CONSTANT, 5*std::pow(10,-6)};
5. * std::pow(10, -6)); // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT,
// 5. * std::pow(10, -6));
// set timestep for simulation to 1 second // set timestep for simulation to 1 second
diffu.setTimestep(1.); diffu.setTimestep(1.);
@ -33,7 +39,7 @@ int main(int argc, char *argv[]) {
// loop 100 times // loop 100 times
// output is currently generated by the method itself // output is currently generated by the method itself
for (int i = 0; i < 100; i++) { for (int i = 0; i < 100; i++) {
diffu.simulate(field, alpha); diffu.simulate(field.data(), alpha.data(), bc.data());
cout << "Iteration: " << i << "\n\n"; 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 "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max #include <algorithm> // for copy, max
#include <cmath>
#include <iomanip> #include <iomanip>
#include <iostream> // for std #include <iostream> // for std
#include <vector> // for vector #include <vector> // for vector
using namespace std; using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
// dimension of grid // dimension of grid
@ -16,6 +20,7 @@ int main(int argc, char *argv[]) {
// create input + diffusion coefficients for each grid cell // create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n * m, 1 * pow(10, -1)); std::vector<double> alpha(n * m, 1 * pow(10, -1));
std::vector<double> field(n * m, 1 * std::pow(10, -6)); 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 // create instance of diffusion module
BTCSDiffusion diffu(dim); BTCSDiffusion diffu(dim);
@ -23,9 +28,8 @@ int main(int argc, char *argv[]) {
diffu.setXDimensions(1, n); diffu.setXDimensions(1, n);
diffu.setYDimensions(1, m); diffu.setYDimensions(1, m);
// set the boundary condition for the left ghost cell to dirichlet //set inlet to higher concentration without setting bc
diffu.setBoundaryCondition(2, 2, BTCSDiffusion::BC_CONSTANT, field[12] = 5*std::pow(10,-3);
5. * std::pow(10, -6));
// set timestep for simulation to 1 second // set timestep for simulation to 1 second
diffu.setTimestep(1.); diffu.setTimestep(1.);
@ -33,7 +37,7 @@ int main(int argc, char *argv[]) {
cout << setprecision(12); cout << setprecision(12);
for (int t = 0; t < 10; t++) { for (int t = 0; t < 10; t++) {
diffu.simulate(field, alpha); diffu.simulate(field.data(), alpha.data(), bc.data());
cout << "Iteration: " << t << "\n\n"; cout << "Iteration: " << t << "\n\n";