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

use new indexing method for linear system

Browse Source
This commit is contained in:
Max Lübke 2021-11-05 10:07:54 +01:00
parent 2f728caa66
commit 00dad33b0f
1 changed files with 19 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
src/diffusion.cpp
View File

@ -1,6 +1,8 @@
#include "diffusion.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>
@ -9,8 +11,6 @@
#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>
@ -20,92 +20,38 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &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<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);
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<double> A(size, (x * y) - 4);
A.setFromTriplets(tripletList.begin(), tripletList.end());
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;
// Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
// solver;