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

changed the solver, now it is running still no valid solution

Browse Source
This commit is contained in:
Max Lübke 2021-10-26 15:01:39 +02:00
parent 72f527a59b
commit 2f728caa66
2 changed files with 72 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
src/diffusion.cpp
View File

@ -4,10 +4,14 @@
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h> #include <Eigen/src/Core/util/Constants.h>
#include <Eigen/src/OrderingMethods/Ordering.h> #include <Eigen/src/OrderingMethods/Ordering.h>
#include <Eigen/src/SparseCholesky/SimplicialCholesky.h>
#include <Eigen/src/SparseCore/SparseMap.h> #include <Eigen/src/SparseCore/SparseMap.h>
#include <Eigen/src/SparseCore/SparseMatrix.h> #include <Eigen/src/SparseCore/SparseMatrix.h>
#include <Eigen/src/SparseCore/SparseMatrixBase.h> #include <Eigen/src/SparseCore/SparseMatrixBase.h>
#include <Eigen/src/SparseLU/SparseLU.h> #include <Eigen/src/SparseLU/SparseLU.h>
#include<Eigen/SparseCholesky>
#include<Eigen/SparseQR>
#include <Eigen/src/SparseQR/SparseQR.h>
#include <iostream> #include <iostream>
void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha, void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
@ -16,7 +20,7 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
double dx = 1. / x; double dx = 1. / x;
double dy = 1. / y; double dy = 1. / y;
int size = (x - 2) * (y - 2); int size = (x * y) - 4;
int local_x = x - 2; int local_x = x - 2;
@ -25,22 +29,62 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
std::vector<T> tripletList; std::vector<T> tripletList;
tripletList.reserve(size * 5); tripletList.reserve(size * 5);
for (int i = 0; i < (size-local_x); i++) { for (int i = x - 1 ; i < 2*x - 3 ; i++) {
int current = local_x + i; double sx = (alpha[i + 2] * timestep) / (dx * dx);
double sx = (alpha[current] * timestep) / (dx * dx); double sy = (alpha[i + 2] * timestep) / (dy * dy);
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; 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[x+i]; 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 i = 0; i < y; i++) {
// for (int j = 0; j < x; j++) { // for (int j = 0; j < x; j++) {
// double sx = (alpha[i * y + j] * timestep) / (dx * dx); // double sx = (alpha[i * y + j] * timestep) / (dx * dx);
@ -56,16 +100,23 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
// } // }
// } // }
Eigen::SparseMatrix<double> A(size, (x*y)-4); std::cout << b << std::endl;
Eigen::SparseMatrix<double> A(size, (x * y) - 4);
A.setFromTriplets(tripletList.begin(), tripletList.end()); A.setFromTriplets(tripletList.begin(), tripletList.end());
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
solver;
// Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
// solver;
solver.analyzePattern(A); solver.analyzePattern(A);
solver.factorize(A); solver.factorize(A);
std::cout << A << std::endl; std::cout << A << std::endl;
std::cout << solver.lastErrorMessage() << std::endl;
x_out = solver.solve(b); x_out = solver.solve(b);
std::cout << x_out << std::endl;
} }

10
src/main.cpp
View File

@ -7,12 +7,14 @@ using namespace std;
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
int x = 10; int x = 5;
int y = 10; int y = 5;
std::vector<double> alpha(x * y, 1.5 * pow(10, -9)); std::vector<double> alpha(x * y, -1.5 * pow(10, -2));
std::vector<double> input(x * y, 0); std::vector<double> input(x * y, 0);
input[x + 1] = 5 * std::pow(10, 6); 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);
BTCS2D(x, y, input, alpha, 10.); BTCS2D(x, y, input, alpha, 10.);