max/TugJulia
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implemented 1D diffusion with new data structure

Browse Source
This commit is contained in:
Max Lübke 2021-12-13 19:39:22 +01:00
parent b41ad0659c
commit 816a1cc67a
4 changed files with 168 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -7,3 +7,4 @@ build/
compile_commands.json
.cache/
.ccls-cache/
/iwyu/

174
src/BTCSDiffusion.cpp
View File

@ -1,129 +1,127 @@
#include "BTCSDiffusion.hpp"
#include <Eigen/SparseCholesky>
#include <Eigen/SparseLU>
#include <Eigen/SparseQR>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <Eigen/src/OrderingMethods/Ordering.h>
#include <Eigen/src/SparseCholesky/SimplicialCholesky.h>
#include <Eigen/src/SparseCore/SparseMap.h>
#include <Eigen/src/SparseCore/SparseMatrix.h>
#include <Eigen/src/SparseCore/SparseMatrixBase.h>
#include <Eigen/src/SparseLU/SparseLU.h>
#include <Eigen/src/SparseQR/SparseQR.h>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <tuple>
#include <vector>
const BCSide BTCSDiffusion::LEFT = 0;
const BCSide BTCSDiffusion::RIGHT = 1;
const int BTCSDiffusion::BC_NEUMANN = 0;
const int BTCSDiffusion::BC_DIRICHLET = 1;
BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
this->grid_dim = 1;
this->dx = 1. / (x - 1);
// per default use Neumann condition with gradient of 0 at the end of the grid
this->bc.resize(2, -1);
this->bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_NEUMANN, 0.));
}
BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
this->grid_dim = 2;
// this->grid_dim = 2;
this->bc.reserve(x * 2 + y * 2);
// per default use Neumann condition with gradient of 0 at the end of the grid
std::fill(this->bc.begin(), this->bc.end(), -1);
// this->bc.reserve(x * 2 + y * 2);
// // per default use Neumann condition with gradient of 0 at the end of the
// grid std::fill(this->bc.begin(), this->bc.end(), -1);
}
BTCSDiffusion::BTCSDiffusion(int x, int y, int z)
: dim_x(x), dim_y(y), dim_z(z) {
BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) {
this->grid_dim = 3;
// this->grid_dim = 3;
// TODO: reserve memory for boundary conditions
}
void BTCSDiffusion::setBoundaryCondition(std::vector<double> input,
BCSide side) {
if (this->grid_dim == 1) {
bc[side] = input[0];
}
}
void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
double timestep) {
// calculate dx
double dx = 1. / (this->dim_x - 1);
void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
double bc_right,
const std::vector<double> &alpha, double dx,
int size) {
// calculate size needed for A matrix and b,x vectors
int size = this->dim_x + 2;
// we need 2 more grid cells for ghost cells
size = size + 2;
Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
Eigen::VectorXd x_out(size);
// set sizes of private and yet allocated vectors
b_vector.resize(size);
x_vector.resize(size);
/*
* Initalization of matrix A
* This is done by triplets. See:
* https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
*/
std::vector<T> tripletList;
tripletList.reserve(c.size() * 3 + bc.size());
int A_line = 0;
// For all concentrations create one row in matrix A
for (int i = 1; i < this->dim_x + 1; i++) {
double sx = (alpha[i - 1] * timestep) / (dx * dx);
tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
tripletList.push_back(T(A_line, i - 1, sx));
tripletList.push_back(T(A_line, i + 1, sx));
b[A_line] = -c[i - 1];
A_line++;
}
// append left and right boundary conditions/ghost zones
tripletList.push_back(T(A_line, 0, 1));
// if value is -1 apply Neumann condition with given gradient
// TODO: set specific gradient
if (bc[0] == -1)
b[A_line] = c[0];
// else apply given Dirichlet condition
else
b[A_line] = this->bc[0];
A_line++;
tripletList.push_back(T(A_line, size - 1, 1));
// b[A_line] = bc[1];
if (bc[1] == -1)
b[A_line] = c[c.size() - 1];
else
b[A_line] = this->bc[1];
/*
* Begin to solve the equation system
* Begin to solve the equation system using LU solver of Eigen.
*
* But first fill the A matrix and b vector.
*
* At this point there is some debugging output in the code.
* TODO: remove output
*/
Eigen::SparseMatrix<double> A(size, size);
A.setFromTriplets(tripletList.begin(), tripletList.end());
A_matrix.resize(size, size);
A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
A_matrix.insert(0, 0) = 1;
A_matrix.insert(size - 1, size - 1) = 1;
b_vector[0] = bc_left;
b_vector[size - 1] = bc_right;
for (int i = 1; i < this->n_x + 1; i++) {
double sx = (alpha[i - 1] * 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];
}
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;
solver.analyzePattern(A);
solver.analyzePattern(A_matrix);
solver.factorize(A);
solver.factorize(A_matrix);
std::cout << solver.lastErrorMessage() << std::endl;
x_out = solver.solve(b);
x_vector = solver.solve(b_vector);
std::cout << std::setprecision(10) << x_out << std::endl << std::endl;
std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
for (int i = 0; i < c.size(); i++) {
c[i] = x_out[i + 1];
c[i] = x_vector[i + 1];
}
}
void BTCSDiffusion::setTimestep(double time_step) {
this->time_step = 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]);
simulate1D(c, bc_left, bc_right, alpha, this->dx, this->n_x);
}
}
double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
double neighbor_alpha) {
double val = -1;
int type = std::get<0>(bc[index]);
if (type == BTCSDiffusion::BC_NEUMANN) {
val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha *
std::get<1>(bc[index]);
} else if (type == BTCSDiffusion::BC_DIRICHLET) {
val = std::get<1>(bc[index]);
} else {
// TODO: implement error handling here. Type was set to wrong value.
}
return val;
}
void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) {
std::get<0>(bc[index]) = type;
std::get<1>(bc[index]) = val;
}

96
src/BTCSDiffusion.hpp
View File

@ -2,19 +2,31 @@
#define BTCSDIFFUSION_H_
#include <Eigen/Sparse>
#include <tuple>
#include <vector>
/*!
* Type defining the side of given boundary condition.
*/
typedef int BCSide;
/*!
* Datatype to fill the sparse matrix which is used to solve the equation
* system.
*/
typedef Eigen::Triplet<double> T;
/*!
* Defines both types of boundary condition as a datatype.
*/
typedef int bctype;
/*!
* A boundary condition consists of two features. A type and the according
* value. Here we can differentiate between:
*
* - Neumann boundary conditon: type BC_NEUMANN with the value defining the
* gradient
* - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
* the boundary condition
*/
typedef std::vector<std::tuple<bctype, double>> boundary_condition;
/*!
* Class implementing a solution for a 1/2/3D diffusion equation using backward
* euler.
@ -23,21 +35,20 @@ class BTCSDiffusion {
public:
/*!
* Set left boundary condition.
* Defines a Neumann boundary condition.
*/
static const BCSide LEFT;
static const int BC_NEUMANN;
/*!
* Set right boundary condition.
* Defines a Dirichlet boundary condition.
*/
static const BCSide RIGHT;
static const int BC_DIRICHLET;
/*!
* Create 1D-diffusion module.
*
* @param x Count of cells in x direction.
*/
BTCSDiffusion(int x);
explicit BTCSDiffusion(int x);
/*!
* Currently not implemented: Create 2D-diffusion module.
@ -45,7 +56,7 @@ public:
* @param x Count of cells in x direction.
* @param y Count of cells in y direction.
*/
BTCSDiffusion(int x, int y);
explicit BTCSDiffusion(int x, int y);
/*!
* Currently not implemented: Create 3D-diffusion module.
@ -54,18 +65,7 @@ public:
* @param y Count of cells in y direction.
* @param z Count of cells in z direction.
*/
BTCSDiffusion(int x, int y, int z);
/*!
* Sets internal boundary condition at the end of the grid/ghost zones.
* Currently only implemented for 1D diffusion.
*
* @param input Vector containing all the values to initialize the ghost
* zones.
* @param side Sets the side of the boundary condition. See BCSide for more
* information.
*/
void setBoundaryCondition(std::vector<double> input, BCSide side);
explicit BTCSDiffusion(int x, int y, int z);
/*!
* With given ghost zones simulate diffusion. Only 1D allowed at this moment.
@ -73,17 +73,51 @@ public:
* @param c Vector describing the concentration of one solution of the grid as
* continious memory (Row-wise).
* @param alpha Vector of diffusioncoefficients for each grid element.
* @param timestep Time (in seconds ?) to simulate.
*/
void simulate(std::vector<double> &c, std::vector<double> &alpha,
double timestep);
void simulate(std::vector<double> &c, const std::vector<double> &alpha);
/*!
* Set the timestep of the simulation
*
* @param time_step Time step (in seconds ???)
*/
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 boundary condition vector.
* @param val Value of the boundary condition (gradient for Neumann, exact
* value for Dirichlet).
* @param Type of the grid cell.
*/
void setBoundaryCondition(int index, double val, bctype type);
private:
std::vector<double> bc;
void simulate1D(std::vector<double> &c, double bc_left, double bc_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);
boundary_condition bc;
Eigen::SparseMatrix<double> A_matrix;
Eigen::VectorXd b_vector;
Eigen::VectorXd x_vector;
double time_step;
int grid_dim;
int dim_x;
int dim_y;
int dim_z;
int n_x;
double dx;
int n_y;
double dy;
int n_z;
double dz;
};
#endif // BTCSDIFFUSION_H_

30
src/main.cpp
View File

@ -1,30 +1,32 @@
#include "BTCSDiffusion.hpp"
#include <cmath>
#include <iostream>
#include <vector>
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include <algorithm> // for copy, max
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
int main(int argc, char *argv[]) {
// count of grid cells
int x = 20;
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(x, 1 * pow(10, -1));
std::vector<double> input(x, 1 * std::pow(10, -6));
std::vector<double> bc_left, bc_right;
bc_left.push_back(5. * std::pow(10, -6));
bc_right.push_back(-1);
// create instance of diffusion module
BTCSDiffusion diffu(x);
diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);
// we don't need this since Neumann condition with gradient of 0 is set per
// default
// diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT);
// set the boundary condition for the left ghost cell to dirichlet
diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
BTCSDiffusion::BC_DIRICHLET);
// set timestep for simulation to 1 second
diffu.setTimestep(1.);
// loop 100 times
// output is currently generated by the method itself
for (int i = 0; i < 100; i++) {
diffu.simulate(input, alpha, 1.);
diffu.simulate(input, alpha);
}
return 0;