From 00dad33b0f54a1d123043e2be78f07d50b72b4ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 5 Nov 2021 10:07:54 +0100 Subject: [PATCH 01/12] use new indexing method for linear system --- src/diffusion.cpp | 92 ++++++++++------------------------------------- 1 file changed, 19 insertions(+), 73 deletions(-) diff --git a/src/diffusion.cpp b/src/diffusion.cpp index 7ecad7f..cda6f4b 100644 --- a/src/diffusion.cpp +++ b/src/diffusion.cpp @@ -1,6 +1,8 @@ #include "diffusion.hpp" +#include #include +#include #include #include #include @@ -9,8 +11,6 @@ #include #include #include -#include -#include #include #include @@ -20,92 +20,38 @@ void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, double dx = 1. / x; double dy = 1. / y; - int size = (x * y) - 4; - - int local_x = x - 2; + int size = (x - 2) * (y - 2); Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); - Eigen::VectorXd x_out(size); + Eigen::VectorXd x_out(x * y); std::vector 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); + int A_line = 0; + for (int i = 1; i < y - 1; i++) { + for (int j = 1; j < x - 1; j++) { + double sx = (alpha[i * x + j] * timestep) / (dx * dx); + double sy = (alpha[i * x + j] * 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)); + tripletList.push_back(T(A_line, i * x + j, (1 + 2 * sx + 2 * sy))); + tripletList.push_back(T(A_line, (i - 1) * x + j, sy)); + tripletList.push_back(T(A_line, (i + 1) * x + j, sy)); + tripletList.push_back(T(A_line, i * x + (j + 1), sx)); + tripletList.push_back(T(A_line, i * x + (j - 1), sx)); - b[i] = -c[i+2]; + b[A_line] = -c[i*x+j]; + A_line++; + } } - 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 A(size, (x * y) - 4); A.setFromTriplets(tripletList.begin(), tripletList.end()); -Eigen::SparseQR, Eigen::COLAMDOrdering> solver; + Eigen::SparseQR, Eigen::COLAMDOrdering> + solver; // Eigen::SparseLU, Eigen::COLAMDOrdering> // solver; From 8dff22c7bf55c83e1acdd9adef84d27d25c51b18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 5 Nov 2021 11:05:19 +0100 Subject: [PATCH 02/12] 1105 --- src/diffusion.cpp | 8 ++++++-- src/main.cpp | 16 ++++++++-------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/diffusion.cpp b/src/diffusion.cpp index cda6f4b..52c71e0 100644 --- a/src/diffusion.cpp +++ b/src/diffusion.cpp @@ -13,6 +13,7 @@ #include #include #include +#include void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, double timestep) { @@ -34,7 +35,10 @@ void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, double sx = (alpha[i * x + j] * timestep) / (dx * dx); double sy = (alpha[i * x + j] * timestep) / (dy * dy); - tripletList.push_back(T(A_line, i * x + j, (1 + 2 * sx + 2 * sy))); + tripletList.push_back(T(A_line, i * x + j, (1. + 2. * sx + 2. * sy))); + + std::cout << sx << std::endl; + tripletList.push_back(T(A_line, (i - 1) * x + j, sy)); tripletList.push_back(T(A_line, (i + 1) * x + j, sy)); tripletList.push_back(T(A_line, i * x + (j + 1), sx)); @@ -47,7 +51,7 @@ void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, std::cout << b << std::endl; - Eigen::SparseMatrix A(size, (x * y) - 4); + Eigen::SparseMatrix A(size, x*y); A.setFromTriplets(tripletList.begin(), tripletList.end()); Eigen::SparseQR, Eigen::COLAMDOrdering> diff --git a/src/main.cpp b/src/main.cpp index 8269d9b..be66bcc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -7,16 +7,16 @@ using namespace std; int main(int argc, char *argv[]) { - int x = 5; - int y = 5; + int x = 4; + int y = 4; - std::vector alpha(x * y, -1.5 * pow(10, -2)); - std::vector 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 alpha(x * y, 1 * pow(10, -9)); + std::vector input(x * y, 1 * std::pow(10,-6)); + input[x + 1] = 5 * std::pow(10, -6); + // input[x + 2] = 5.5556554 * std::pow(10, -6); + // input[x + 3] = 5.234564213 * std::pow(10, -6); - BTCS2D(x, y, input, alpha, 10.); + BTCS2D(x, y, input, alpha, 1.); return 0; } From d18857510569ddf75dd2bc08f275192b7a16a87f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 11 Nov 2021 11:53:14 +0100 Subject: [PATCH 03/12] implement 1D --- src/diffusion.cpp | 63 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 61 insertions(+), 2 deletions(-) diff --git a/src/diffusion.cpp b/src/diffusion.cpp index 52c71e0..726fdf2 100644 --- a/src/diffusion.cpp +++ b/src/diffusion.cpp @@ -13,8 +13,67 @@ #include #include #include +#include #include +void BTCS1D(int x, std::vector &c, std::vector &alpha, + double timestep, std::vector &bc) { + double dx = 1. / x; + + int size = x + 2; + + Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); + Eigen::VectorXd x_out(size); + std::vector tripletList; + tripletList.reserve(c.size() * 3 + bc.size()); + + int A_line = 0; + + for (int i = 1; i < 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++; + } + + tripletList.push_back(T(A_line, 0, 1)); + b[A_line] = bc[0]; + + A_line++; + tripletList.push_back(T(A_line, size-1, 1)); + b[A_line] = bc[1]; + + // std::cout << b << std::endl; + + Eigen::SparseMatrix A(size, size); + A.setFromTriplets(tripletList.begin(), tripletList.end()); + + // std::cout << A << std::endl; + Eigen::SparseQR, Eigen::COLAMDOrdering> + solver; + + // Eigen::SparseLU, Eigen::COLAMDOrdering> + // 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]; + } +} + void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, double timestep) { @@ -44,14 +103,14 @@ void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, tripletList.push_back(T(A_line, i * x + (j + 1), sx)); tripletList.push_back(T(A_line, i * x + (j - 1), sx)); - b[A_line] = -c[i*x+j]; + b[A_line] = -c[i * x + j]; A_line++; } } std::cout << b << std::endl; - Eigen::SparseMatrix A(size, x*y); + Eigen::SparseMatrix A(size, x * y); A.setFromTriplets(tripletList.begin(), tripletList.end()); Eigen::SparseQR, Eigen::COLAMDOrdering> From 3fa31a8f40a1e520d98b25a5974811593f8b5d86 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 23 Nov 2021 14:12:37 +0100 Subject: [PATCH 04/12] push old state --- src/diffusion.cpp | 3 ++- src/diffusion.hpp | 3 +++ src/main.cpp | 16 ++++++++++------ 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/diffusion.cpp b/src/diffusion.cpp index 726fdf2..5b398a2 100644 --- a/src/diffusion.cpp +++ b/src/diffusion.cpp @@ -46,7 +46,8 @@ void BTCS1D(int x, std::vector &c, std::vector &alpha, A_line++; tripletList.push_back(T(A_line, size-1, 1)); - b[A_line] = bc[1]; + // b[A_line] = bc[1]; + b[A_line] = c[c.size()-1]; // std::cout << b << std::endl; diff --git a/src/diffusion.hpp b/src/diffusion.hpp index 6645f8c..9d80e45 100644 --- a/src/diffusion.hpp +++ b/src/diffusion.hpp @@ -6,6 +6,9 @@ typedef Eigen::Triplet T; +extern void BTCS1D(int x, std::vector &c, std::vector &alpha, + double timestep, std::vector &bc); + extern void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, double timestep); #endif // DIFFUSION_H_ diff --git a/src/main.cpp b/src/main.cpp index be66bcc..3613d17 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -7,16 +7,20 @@ using namespace std; int main(int argc, char *argv[]) { - int x = 4; - int y = 4; + int x = 20; - std::vector alpha(x * y, 1 * pow(10, -9)); - std::vector input(x * y, 1 * std::pow(10,-6)); - input[x + 1] = 5 * std::pow(10, -6); + std::vector alpha(x, 1 * pow(10, -1)); + std::vector input(x, 1 * std::pow(10, -6)); + std::vector bc; + + bc.push_back(5. * std::pow(10, -6)); + bc.push_back(1. * std::pow(10, -6)); // input[x + 2] = 5.5556554 * std::pow(10, -6); // input[x + 3] = 5.234564213 * std::pow(10, -6); - BTCS2D(x, y, input, alpha, 1.); + for (int i = 0; i < 100; i++) { + BTCS1D(x, input, alpha, 1., bc); + } return 0; } From f2e80c2c4858141a35d61da67f3cab169cacab2e Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 23 Nov 2021 14:56:52 +0100 Subject: [PATCH 05/12] implement 1D diffusion as class --- src/BTCSDiffusion.cpp | 96 +++++++++++++++++++++++++++++++++++++++++++ src/BTCSDiffusion.hpp | 32 +++++++++++++++ src/CMakeLists.txt | 3 ++ 3 files changed, 131 insertions(+) create mode 100644 src/BTCSDiffusion.cpp create mode 100644 src/BTCSDiffusion.hpp diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp new file mode 100644 index 0000000..919a1e5 --- /dev/null +++ b/src/BTCSDiffusion.cpp @@ -0,0 +1,96 @@ +#include "BTCSDiffusion.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +const BCSide BTCSDiffusion::LEFT = 0; +const BCSide BTCSDiffusion::RIGHT = 1; + +BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { this->grid_dim = 1; } +BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) { + + this->grid_dim = 2; +} +BTCSDiffusion::BTCSDiffusion(int x, int y, int z) + : dim_x(x), dim_y(y), dim_z(z) { + + this->grid_dim = 3; +} + +void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, + double timestep) { + double dx = 1. / this->dim_x; + + int size = this->dim_x + 2; + + Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); + Eigen::VectorXd x_out(size); + std::vector tripletList; + tripletList.reserve(c.size() * 3 + bc.size()); + + int A_line = 0; + + 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++; + } + + tripletList.push_back(T(A_line, 0, 1)); + if (bc[0] == -1) + b[A_line] = c[0]; + 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[0] == -1) + b[A_line] = c[c.size() - 1]; + else + b[A_line] = this->bc[1]; + + // std::cout << b << std::endl; + + Eigen::SparseMatrix A(size, size); + A.setFromTriplets(tripletList.begin(), tripletList.end()); + + // std::cout << A << std::endl; + Eigen::SparseQR, Eigen::COLAMDOrdering> + solver; + + // Eigen::SparseLU, Eigen::COLAMDOrdering> + // 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]; + } +} diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp new file mode 100644 index 0000000..cee17c9 --- /dev/null +++ b/src/BTCSDiffusion.hpp @@ -0,0 +1,32 @@ +#ifndef BTCSDIFFUSION_H_ +#define BTCSDIFFUSION_H_ + +#include +#include + +typedef int BCSide; +typedef Eigen::Triplet T; + +class BTCSDiffusion { + +public: + static const BCSide LEFT; + static const BCSide RIGHT; + + BTCSDiffusion(int x); + BTCSDiffusion(int x, int y); + BTCSDiffusion(int x, int y, int z); + + void setBoundaryCondition(std::vector input, BCSide side); + void simulate(std::vector &c, std::vector &alpha, + double timestep); + +private: + std::vector bc; + int grid_dim; + int dim_x; + int dim_y; + int dim_z; +}; + +#endif // BTCSDIFFUSION_H_ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 33e2d28..97ce3f1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,8 @@ add_library(diffusion OBJECT diffusion.cpp diffusion.hpp) target_link_libraries(diffusion Eigen3::Eigen) +add_library(diffusion_class OBJECT BTCSDiffusion.cpp) +target_link_libraries(diffusion_class Eigen3::Eigen) + add_executable(test main.cpp) target_link_libraries(test PUBLIC diffusion) From 9461dad14724f2110e83e58b7688fc44905238b5 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 23 Nov 2021 15:10:51 +0100 Subject: [PATCH 06/12] added test executable --- src/BTCSDiffusion.cpp | 13 ++++++++++++- src/test_class.cpp | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 1 deletion(-) create mode 100644 src/test_class.cpp diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 919a1e5..9a92418 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -19,17 +19,28 @@ const BCSide BTCSDiffusion::LEFT = 0; const BCSide BTCSDiffusion::RIGHT = 1; -BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { this->grid_dim = 1; } +BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { + this->grid_dim = 1; + this->bc.reserve(2); +} BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) { this->grid_dim = 2; + this->bc.reserve(x * 2 + y * 2); } 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 input, + BCSide side) { + if (this->grid_dim == 1) { + bc[side] = input[0]; + } +} void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, double timestep) { double dx = 1. / this->dim_x; diff --git a/src/test_class.cpp b/src/test_class.cpp new file mode 100644 index 0000000..48d8b14 --- /dev/null +++ b/src/test_class.cpp @@ -0,0 +1,33 @@ +#include "BTCSDiffusion.hpp" +#include "diffusion.hpp" +#include +#include +#include + +using namespace std; + +int main(int argc, char *argv[]) { + + int x = 20; + + std::vector alpha(x, 1 * pow(10, -1)); + std::vector input(x, 1 * std::pow(10, -6)); + std::vector bc_left, bc_right; + + bc_left.push_back(5. * std::pow(10, -6)); + bc_right.push_back(1. * std::pow(10, -6)); + // input[x + 2] = 5.5556554 * std::pow(10, -6); + // input[x + 3] = 5.234564213 * std::pow(10, -6); + + BTCSDiffusion diffu(x); + + diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); + diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); + + for (int i = 0; i < 100; i++) { + diffu.simulate(input, alpha, 1.); + // BTCS1D(x, input, alpha, 1., bc); + } + + return 0; +} From 57a0e8a1a67f7d1b87c2293b17a992a91e28eecb Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 23 Nov 2021 15:22:46 +0100 Subject: [PATCH 07/12] fix some cmake bugs --- src/BTCSDiffusion.cpp | 2 +- src/BTCSDiffusion.hpp | 2 +- src/CMakeLists.txt | 5 ++++- src/test_class.cpp | 3 +-- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 9a92418..6c9da99 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -75,7 +75,7 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, A_line++; tripletList.push_back(T(A_line, size - 1, 1)); // b[A_line] = bc[1]; - if (bc[0] == -1) + if (bc[1] == -1) b[A_line] = c[c.size() - 1]; else b[A_line] = this->bc[1]; diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index cee17c9..693b8f7 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -2,7 +2,7 @@ #define BTCSDIFFUSION_H_ #include -#include +#include typedef int BCSide; typedef Eigen::Triplet T; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 97ce3f1..9678862 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,11 @@ add_library(diffusion OBJECT diffusion.cpp diffusion.hpp) target_link_libraries(diffusion Eigen3::Eigen) -add_library(diffusion_class OBJECT BTCSDiffusion.cpp) +add_library(diffusion_class OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp) target_link_libraries(diffusion_class Eigen3::Eigen) add_executable(test main.cpp) target_link_libraries(test PUBLIC diffusion) + +add_executable(test_class test_class.cpp) +target_link_libraries(test_class PUBLIC diffusion_class) diff --git a/src/test_class.cpp b/src/test_class.cpp index 48d8b14..fd28c85 100644 --- a/src/test_class.cpp +++ b/src/test_class.cpp @@ -1,5 +1,4 @@ #include "BTCSDiffusion.hpp" -#include "diffusion.hpp" #include #include #include @@ -15,7 +14,7 @@ int main(int argc, char *argv[]) { std::vector bc_left, bc_right; bc_left.push_back(5. * std::pow(10, -6)); - bc_right.push_back(1. * std::pow(10, -6)); + bc_right.push_back(-1); // input[x + 2] = 5.5556554 * std::pow(10, -6); // input[x + 3] = 5.234564213 * std::pow(10, -6); From 96a2d1cc5bb5c7e0c31770ab67c9a5b29f510193 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 1 Dec 2021 17:59:08 +0100 Subject: [PATCH 08/12] remove old library files --- src/CMakeLists.txt | 8 +-- src/diffusion.cpp | 132 --------------------------------------------- src/diffusion.hpp | 14 ----- src/main.cpp | 16 ++++-- src/test_class.cpp | 32 ----------- 5 files changed, 12 insertions(+), 190 deletions(-) delete mode 100644 src/diffusion.cpp delete mode 100644 src/diffusion.hpp delete mode 100644 src/test_class.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9678862..fab0508 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,11 +1,5 @@ -add_library(diffusion OBJECT diffusion.cpp diffusion.hpp) +add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp) target_link_libraries(diffusion Eigen3::Eigen) -add_library(diffusion_class OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp) -target_link_libraries(diffusion_class Eigen3::Eigen) - add_executable(test main.cpp) target_link_libraries(test PUBLIC diffusion) - -add_executable(test_class test_class.cpp) -target_link_libraries(test_class PUBLIC diffusion_class) diff --git a/src/diffusion.cpp b/src/diffusion.cpp deleted file mode 100644 index 5b398a2..0000000 --- a/src/diffusion.cpp +++ /dev/null @@ -1,132 +0,0 @@ -#include "diffusion.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -void BTCS1D(int x, std::vector &c, std::vector &alpha, - double timestep, std::vector &bc) { - double dx = 1. / x; - - int size = x + 2; - - Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); - Eigen::VectorXd x_out(size); - std::vector tripletList; - tripletList.reserve(c.size() * 3 + bc.size()); - - int A_line = 0; - - for (int i = 1; i < 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++; - } - - tripletList.push_back(T(A_line, 0, 1)); - b[A_line] = bc[0]; - - A_line++; - tripletList.push_back(T(A_line, size-1, 1)); - // b[A_line] = bc[1]; - b[A_line] = c[c.size()-1]; - - // std::cout << b << std::endl; - - Eigen::SparseMatrix A(size, size); - A.setFromTriplets(tripletList.begin(), tripletList.end()); - - // std::cout << A << std::endl; - Eigen::SparseQR, Eigen::COLAMDOrdering> - solver; - - // Eigen::SparseLU, Eigen::COLAMDOrdering> - // 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]; - } -} - -void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, - double timestep) { - - double dx = 1. / x; - double dy = 1. / y; - - int size = (x - 2) * (y - 2); - - Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); - Eigen::VectorXd x_out(x * y); - std::vector tripletList; - tripletList.reserve(size * 5); - - int A_line = 0; - - for (int i = 1; i < y - 1; i++) { - for (int j = 1; j < x - 1; j++) { - double sx = (alpha[i * x + j] * timestep) / (dx * dx); - double sy = (alpha[i * x + j] * timestep) / (dy * dy); - - tripletList.push_back(T(A_line, i * x + j, (1. + 2. * sx + 2. * sy))); - - std::cout << sx << std::endl; - - tripletList.push_back(T(A_line, (i - 1) * x + j, sy)); - tripletList.push_back(T(A_line, (i + 1) * x + j, sy)); - tripletList.push_back(T(A_line, i * x + (j + 1), sx)); - tripletList.push_back(T(A_line, i * x + (j - 1), sx)); - - b[A_line] = -c[i * x + j]; - A_line++; - } - } - - std::cout << b << std::endl; - - Eigen::SparseMatrix A(size, x * y); - A.setFromTriplets(tripletList.begin(), tripletList.end()); - - Eigen::SparseQR, Eigen::COLAMDOrdering> - solver; - - // Eigen::SparseLU, Eigen::COLAMDOrdering> - // 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; -} diff --git a/src/diffusion.hpp b/src/diffusion.hpp deleted file mode 100644 index 9d80e45..0000000 --- a/src/diffusion.hpp +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef DIFFUSION_H_ -#define DIFFUSION_H_ - -#include -#include - -typedef Eigen::Triplet T; - -extern void BTCS1D(int x, std::vector &c, std::vector &alpha, - double timestep, std::vector &bc); - -extern void BTCS2D(int x, int y, std::vector &c, - std::vector &alpha, double timestep); -#endif // DIFFUSION_H_ diff --git a/src/main.cpp b/src/main.cpp index 3613d17..fd28c85 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,4 +1,4 @@ -#include "diffusion.hpp" +#include "BTCSDiffusion.hpp" #include #include #include @@ -11,15 +11,21 @@ int main(int argc, char *argv[]) { std::vector alpha(x, 1 * pow(10, -1)); std::vector input(x, 1 * std::pow(10, -6)); - std::vector bc; + std::vector bc_left, bc_right; - bc.push_back(5. * std::pow(10, -6)); - bc.push_back(1. * std::pow(10, -6)); + bc_left.push_back(5. * std::pow(10, -6)); + bc_right.push_back(-1); // input[x + 2] = 5.5556554 * std::pow(10, -6); // input[x + 3] = 5.234564213 * std::pow(10, -6); + BTCSDiffusion diffu(x); + + diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); + diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); + for (int i = 0; i < 100; i++) { - BTCS1D(x, input, alpha, 1., bc); + diffu.simulate(input, alpha, 1.); + // BTCS1D(x, input, alpha, 1., bc); } return 0; diff --git a/src/test_class.cpp b/src/test_class.cpp deleted file mode 100644 index fd28c85..0000000 --- a/src/test_class.cpp +++ /dev/null @@ -1,32 +0,0 @@ -#include "BTCSDiffusion.hpp" -#include -#include -#include - -using namespace std; - -int main(int argc, char *argv[]) { - - int x = 20; - - std::vector alpha(x, 1 * pow(10, -1)); - std::vector input(x, 1 * std::pow(10, -6)); - std::vector bc_left, bc_right; - - bc_left.push_back(5. * std::pow(10, -6)); - bc_right.push_back(-1); - // input[x + 2] = 5.5556554 * std::pow(10, -6); - // input[x + 3] = 5.234564213 * std::pow(10, -6); - - BTCSDiffusion diffu(x); - - diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); - diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); - - for (int i = 0; i < 100; i++) { - diffu.simulate(input, alpha, 1.); - // BTCS1D(x, input, alpha, 1., bc); - } - - return 0; -} From 85278bcaff8fbb1fb6827085cc2d6b1f7d3d735d Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 1 Dec 2021 18:01:52 +0100 Subject: [PATCH 09/12] cleanup code --- src/BTCSDiffusion.cpp | 3 --- src/main.cpp | 3 --- 2 files changed, 6 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 6c9da99..4588a6f 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -80,12 +80,9 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, else b[A_line] = this->bc[1]; - // std::cout << b << std::endl; - Eigen::SparseMatrix A(size, size); A.setFromTriplets(tripletList.begin(), tripletList.end()); - // std::cout << A << std::endl; Eigen::SparseQR, Eigen::COLAMDOrdering> solver; diff --git a/src/main.cpp b/src/main.cpp index fd28c85..03b5152 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -15,8 +15,6 @@ int main(int argc, char *argv[]) { bc_left.push_back(5. * std::pow(10, -6)); bc_right.push_back(-1); - // input[x + 2] = 5.5556554 * std::pow(10, -6); - // input[x + 3] = 5.234564213 * std::pow(10, -6); BTCSDiffusion diffu(x); @@ -25,7 +23,6 @@ int main(int argc, char *argv[]) { for (int i = 0; i < 100; i++) { diffu.simulate(input, alpha, 1.); - // BTCS1D(x, input, alpha, 1., bc); } return 0; From eb595bc0a334183f165967a304f147cc6d50cbbd Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 1 Dec 2021 18:03:17 +0100 Subject: [PATCH 10/12] Use LU solver instead of QR --- src/BTCSDiffusion.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 4588a6f..d2baef1 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -83,11 +83,8 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, Eigen::SparseMatrix A(size, size); A.setFromTriplets(tripletList.begin(), tripletList.end()); - Eigen::SparseQR, Eigen::COLAMDOrdering> + Eigen::SparseLU, Eigen::COLAMDOrdering> solver; - - // Eigen::SparseLU, Eigen::COLAMDOrdering> - // solver; solver.analyzePattern(A); solver.factorize(A); From 971f8212afaba568a09fa0da92066e9306c4dbc2 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 2 Dec 2021 09:25:34 +0100 Subject: [PATCH 11/12] Added comments --- src/BTCSDiffusion.cpp | 32 +++++++++++++++++++++-- src/BTCSDiffusion.hpp | 59 ++++++++++++++++++++++++++++++++++++++++++- src/main.cpp | 4 ++- 3 files changed, 91 insertions(+), 4 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index d2baef1..287fcb7 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -13,6 +13,7 @@ #include #include +#include #include #include @@ -21,18 +22,23 @@ const BCSide BTCSDiffusion::RIGHT = 1; BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { this->grid_dim = 1; - this->bc.reserve(2); + + // 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 + // TODO: reserve memory for boundary conditions } void BTCSDiffusion::setBoundaryCondition(std::vector input, @@ -43,17 +49,27 @@ void BTCSDiffusion::setBoundaryCondition(std::vector input, } void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, double timestep) { + // calculate dx double dx = 1. / this->dim_x; + // 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 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); @@ -66,9 +82,14 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, 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]; @@ -80,6 +101,13 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, 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 A(size, size); A.setFromTriplets(tripletList.begin(), tripletList.end()); diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 693b8f7..cadc1e7 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -1,23 +1,80 @@ #ifndef BTCSDIFFUSION_H_ #define BTCSDIFFUSION_H_ -#include #include +#include +/*! + * 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 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 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 &c, std::vector &alpha, double timestep); diff --git a/src/main.cpp b/src/main.cpp index 03b5152..2836c7c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -19,7 +19,9 @@ int main(int argc, char *argv[]) { BTCSDiffusion diffu(x); diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); - diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); + // 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.); From 3562f30efaaa6e7c5abc5d46622e218e8d6d95f9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 2 Dec 2021 10:51:47 +0100 Subject: [PATCH 12/12] change dx to 1/n-1 --- src/BTCSDiffusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 287fcb7..fcebc92 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -50,7 +50,7 @@ void BTCSDiffusion::setBoundaryCondition(std::vector input, void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, double timestep) { // calculate dx - double dx = 1. / this->dim_x; + double dx = 1. / (this->dim_x - 1); // calculate size needed for A matrix and b,x vectors int size = this->dim_x + 2;