mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 17:38:23 +01:00
Merge branch 'ml' into 'dev'
Implemented 1D Diffusion See merge request mluebke/diffusion!1
This commit is contained in:
Max Lübke
commit
ef831a2cd4
129
src/BTCSDiffusion.cpp
Normal file
129
src/BTCSDiffusion.cpp
Normal file
@ -0,0 +1,129 @@
|
||||
#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>
|
||||
|
||||
const BCSide BTCSDiffusion::LEFT = 0;
|
||||
const BCSide BTCSDiffusion::RIGHT = 1;
|
||||
|
||||
BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
|
||||
this->grid_dim = 1;
|
||||
|
||||
// per default use Neumann condition with gradient of 0 at the end of the grid
|
||||
this->bc.resize(2, -1);
|
||||
}
|
||||
BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
|
||||
|
||||
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);
|
||||
}
|
||||
BTCSDiffusion::BTCSDiffusion(int x, int y, int z)
|
||||
: dim_x(x), dim_y(y), dim_z(z) {
|
||||
|
||||
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);
|
||||
|
||||
// calculate size needed for A matrix and b,x vectors
|
||||
int size = this->dim_x + 2;
|
||||
|
||||
Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
|
||||
Eigen::VectorXd x_out(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
|
||||
*
|
||||
* 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());
|
||||
|
||||
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
|
||||
solver;
|
||||
solver.analyzePattern(A);
|
||||
|
||||
solver.factorize(A);
|
||||
|
||||
std::cout << solver.lastErrorMessage() << std::endl;
|
||||
|
||||
x_out = solver.solve(b);
|
||||
|
||||
std::cout << std::setprecision(10) << x_out << std::endl << std::endl;
|
||||
|
||||
for (int i = 0; i < c.size(); i++) {
|
||||
c[i] = x_out[i + 1];
|
||||
}
|
||||
}
|
||||
89
src/BTCSDiffusion.hpp
Normal file
89
src/BTCSDiffusion.hpp
Normal file
@ -0,0 +1,89 @@
|
||||
#ifndef BTCSDIFFUSION_H_
|
||||
#define BTCSDIFFUSION_H_
|
||||
|
||||
#include <Eigen/Sparse>
|
||||
#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;
|
||||
|
||||
/*!
|
||||
* Class implementing a solution for a 1/2/3D diffusion equation using backward
|
||||
* euler.
|
||||
*/
|
||||
class BTCSDiffusion {
|
||||
|
||||
public:
|
||||
/*!
|
||||
* Set left boundary condition.
|
||||
*/
|
||||
static const BCSide LEFT;
|
||||
|
||||
/*!
|
||||
* Set right boundary condition.
|
||||
*/
|
||||
static const BCSide RIGHT;
|
||||
|
||||
/*!
|
||||
* Create 1D-diffusion module.
|
||||
*
|
||||
* @param x Count of cells in x direction.
|
||||
*/
|
||||
BTCSDiffusion(int x);
|
||||
|
||||
/*!
|
||||
* Currently not implemented: Create 2D-diffusion module.
|
||||
*
|
||||
* @param x Count of cells in x direction.
|
||||
* @param y Count of cells in y direction.
|
||||
*/
|
||||
BTCSDiffusion(int x, int y);
|
||||
|
||||
/*!
|
||||
* Currently not implemented: Create 3D-diffusion module.
|
||||
*
|
||||
* @param x Count of cells in x direction.
|
||||
* @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);
|
||||
|
||||
/*!
|
||||
* With given ghost zones simulate diffusion. Only 1D allowed at this moment.
|
||||
*
|
||||
* @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);
|
||||
|
||||
private:
|
||||
std::vector<double> bc;
|
||||
int grid_dim;
|
||||
int dim_x;
|
||||
int dim_y;
|
||||
int dim_z;
|
||||
};
|
||||
|
||||
#endif // BTCSDIFFUSION_H_
|
||||
@ -1,4 +1,4 @@
|
||||
add_library(diffusion OBJECT diffusion.cpp diffusion.hpp)
|
||||
add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp)
|
||||
target_link_libraries(diffusion Eigen3::Eigen)
|
||||
|
||||
add_executable(test main.cpp)
|
||||
|
||||
@ -1,122 +0,0 @@
|
||||
#include "diffusion.hpp"
|
||||
|
||||
#include <Eigen/SparseLU>
|
||||
#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/SparseCholesky>
|
||||
#include<Eigen/SparseQR>
|
||||
#include <Eigen/src/SparseQR/SparseQR.h>
|
||||
#include <iostream>
|
||||
|
||||
void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
|
||||
double timestep) {
|
||||
|
||||
double dx = 1. / x;
|
||||
double dy = 1. / y;
|
||||
|
||||
int size = (x * y) - 4;
|
||||
|
||||
int local_x = x - 2;
|
||||
|
||||
Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
|
||||
Eigen::VectorXd x_out(size);
|
||||
std::vector<T> tripletList;
|
||||
tripletList.reserve(size * 5);
|
||||
|
||||
for (int i = x - 1 ; i < 2*x - 3 ; i++) {
|
||||
double sx = (alpha[i + 2] * timestep) / (dx * dx);
|
||||
double sy = (alpha[i + 2] * timestep) / (dy * dy);
|
||||
|
||||
|
||||
tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
|
||||
tripletList.push_back(T(i, i - (x - 1), sy));
|
||||
tripletList.push_back(T(i, i + x, sy));
|
||||
tripletList.push_back(T(i, i + 1, sx));
|
||||
tripletList.push_back(T(i, i - 1, sx));
|
||||
|
||||
b[i] = -c[i+2];
|
||||
}
|
||||
|
||||
for (int i = 2*x - 1; i < (y-2)*x-3; i++) {
|
||||
double sx = (alpha[i + 2] * timestep) / (dx * dx);
|
||||
double sy = (alpha[i + 2] * timestep) / (dy * dy);
|
||||
|
||||
tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
|
||||
tripletList.push_back(T(i, i - x, sy));
|
||||
tripletList.push_back(T(i, i + x, sy));
|
||||
tripletList.push_back(T(i, i + 1, sx));
|
||||
tripletList.push_back(T(i, i - 1, sx));
|
||||
|
||||
b[i] = -c[i+2];
|
||||
}
|
||||
|
||||
for (int i = (y-2)*x-1; i < (y-1)*x-3; i++) {
|
||||
double sx = (alpha[i + 2] * timestep) / (dx * dx);
|
||||
double sy = (alpha[i + 2] * timestep) / (dy * dy);
|
||||
|
||||
tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
|
||||
tripletList.push_back(T(i, i - x, sy));
|
||||
tripletList.push_back(T(i, i + (x-1), sy));
|
||||
tripletList.push_back(T(i, i + 1, sx));
|
||||
tripletList.push_back(T(i, i - 1, sx));
|
||||
|
||||
b[i] = -c[i+2];
|
||||
}
|
||||
|
||||
// for (int i = 0; i < (size-local_x); i++) {
|
||||
// int current = local_x + i;
|
||||
// double sx = (alpha[current] * timestep) / (dx * dx);
|
||||
// double sy = (alpha[current] * timestep) / (dy * dy);
|
||||
|
||||
// tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy)));
|
||||
// tripletList.push_back(T(i, current + local_x, sy));
|
||||
// tripletList.push_back(T(i, current - local_x, sy));
|
||||
// tripletList.push_back(T(i, current + 1, sx));
|
||||
// tripletList.push_back(T(i, current - 1, sx));
|
||||
|
||||
// std::cout << current << std::endl;
|
||||
|
||||
// b[i] = -c[x+i];
|
||||
// }
|
||||
|
||||
// for (int i = 0; i < y; i++) {
|
||||
// for (int j = 0; j < x; j++) {
|
||||
// double sx = (alpha[i * y + j] * timestep) / (dx * dx);
|
||||
// double sy = (alpha[i * y + j] * timestep) / (dy * dy);
|
||||
|
||||
// tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 *
|
||||
// sy))); tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy));
|
||||
// tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
|
||||
// tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
|
||||
// tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
|
||||
|
||||
// b[(i * x) + j] = -c[(i * x) + j];
|
||||
// }
|
||||
// }
|
||||
|
||||
std::cout << b << std::endl;
|
||||
|
||||
Eigen::SparseMatrix<double> A(size, (x * y) - 4);
|
||||
A.setFromTriplets(tripletList.begin(), tripletList.end());
|
||||
|
||||
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
|
||||
|
||||
// Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
|
||||
// solver;
|
||||
solver.analyzePattern(A);
|
||||
|
||||
solver.factorize(A);
|
||||
|
||||
std::cout << A << std::endl;
|
||||
std::cout << solver.lastErrorMessage() << std::endl;
|
||||
|
||||
x_out = solver.solve(b);
|
||||
|
||||
std::cout << x_out << std::endl;
|
||||
}
|
||||
@ -1,11 +0,0 @@
|
||||
#ifndef DIFFUSION_H_
|
||||
#define DIFFUSION_H_
|
||||
|
||||
#include <Eigen/SparseCore>
|
||||
#include <vector>
|
||||
|
||||
typedef Eigen::Triplet<double> T;
|
||||
|
||||
extern void BTCS2D(int x, int y, std::vector<double> &c,
|
||||
std::vector<double> &alpha, double timestep);
|
||||
#endif // DIFFUSION_H_
|
||||
27
src/main.cpp
27
src/main.cpp
@ -1,4 +1,4 @@
|
||||
#include "diffusion.hpp"
|
||||
#include "BTCSDiffusion.hpp"
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
@ -7,16 +7,25 @@ using namespace std;
|
||||
|
||||
int main(int argc, char *argv[]) {
|
||||
|
||||
int x = 5;
|
||||
int y = 5;
|
||||
int x = 20;
|
||||
|
||||
std::vector<double> alpha(x * y, -1.5 * pow(10, -2));
|
||||
std::vector<double> input(x * y, 0);
|
||||
input[x + 1] = 5.55 * std::pow(10, -6);
|
||||
input[x + 2] = 5.5556554 * std::pow(10, -6);
|
||||
input[x + 3] = 5.234564213 * std::pow(10, -6);
|
||||
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;
|
||||
|
||||
BTCS2D(x, y, input, alpha, 10.);
|
||||
bc_left.push_back(5. * std::pow(10, -6));
|
||||
bc_right.push_back(-1);
|
||||
|
||||
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);
|
||||
|
||||
for (int i = 0; i < 100; i++) {
|
||||
diffu.simulate(input, alpha, 1.);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user