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

Use local matrices and vectors.

Browse Source 
- Subsitute private members of Eigen matrices and vectors by local
variables in `simulate_base`
This commit is contained in:
Max Luebke 2022-03-09 11:22:53 +01:00
parent 1cc6b247b9
commit c2211c8a6f
2 changed files with 40 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
src/BTCSDiffusion.cpp
View File

@ -5,6 +5,8 @@
#include <Eigen/src/Core/Map.h> #include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/SparseCore/SparseMatrix.h>
#include <Eigen/src/SparseCore/SparseMatrixBase.h>
#include <algorithm> #include <algorithm>
#include <array> #include <array>
#include <cassert> #include <cassert>
@ -86,28 +88,32 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
int size, int size,
const DVectorRowMajor &t0_c) { const DVectorRowMajor &t0_c) {
reserveMemory(size, BTCS_MAX_DEP_PER_CELL); Eigen::SparseMatrix<double> A_matrix;
Eigen::VectorXd b_vector;
Eigen::VectorXd x_vector;
fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step); A_matrix.resize(size + 2, size + 2);
fillVectorFromRow(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, dx, A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL));
time_step);
solveLES(); b_vector.resize(size + 2);
x_vector.resize(size + 2);
fillMatrixFromRow(A_matrix, alpha.row(0), bc.row(0), size, dx, time_step);
fillVectorFromRow(b_vector, c, alpha, bc, Eigen::VectorXd::Constant(size, 0),
size, dx, time_step);
// start to solve
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;
solver.analyzePattern(A_matrix);
solver.factorize(A_matrix);
x_vector = solver.solve(b_vector);
c = x_vector.segment(1, size); c = x_vector.segment(1, size);
} }
inline void Diffusion::BTCSDiffusion::reserveMemory(int size,
int max_count_per_line) {
size += 2;
A_matrix.resize(size, size);
A_matrix.reserve(Eigen::VectorXi::Constant(size, max_count_per_line));
b_vector.resize(size);
x_vector.resize(size);
}
void Diffusion::BTCSDiffusion::simulate1D( void Diffusion::BTCSDiffusion::simulate1D(
Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha, Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
Eigen::Map<const BCVectorRowMajor> &bc) { Eigen::Map<const BCVectorRowMajor> &bc) {
@ -208,9 +214,9 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
return t0_c; return t0_c;
} }
inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( void Diffusion::BTCSDiffusion::fillMatrixFromRow(
const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, Eigen::SparseMatrix<double> &A_matrix, const DVectorRowMajor &alpha,
double dx, double time_step) { const BCVectorRowMajor &bc, int size, double dx, double time_step) {
Diffusion::boundary_condition left = bc[0]; Diffusion::boundary_condition left = bc[0];
Diffusion::boundary_condition right = bc[size - 1]; Diffusion::boundary_condition right = bc[size - 1];
@ -247,10 +253,10 @@ inline void Diffusion::BTCSDiffusion::fillMatrixFromRow(
} }
} }
inline void Diffusion::BTCSDiffusion::fillVectorFromRow( void Diffusion::BTCSDiffusion::fillVectorFromRow(
const DVectorRowMajor &c, const DVectorRowMajor &alpha, Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, const DVectorRowMajor &alpha, const BCVectorRowMajor &bc,
double dx, double time_step) { const DVectorRowMajor &t0_c, int size, double dx, double time_step) {
Diffusion::boundary_condition left = bc[0]; Diffusion::boundary_condition left = bc[0];
Diffusion::boundary_condition right = bc[size - 1]; Diffusion::boundary_condition right = bc[size - 1];
@ -341,14 +347,3 @@ inline auto Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc,
return val; return val;
} }
inline void Diffusion::BTCSDiffusion::solveLES() {
// start to solve
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;
solver.analyzePattern(A_matrix);
solver.factorize(A_matrix);
x_vector = solver.solve(b_vector);
}

25
src/BTCSDiffusion.hpp
View File

@ -7,6 +7,7 @@
#include <Eigen/src/Core/Map.h> #include <Eigen/src/Core/Map.h>
#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/SparseCore/SparseMatrix.h>
#include <cstddef> #include <cstddef>
#include <tuple> #include <tuple>
#include <type_traits> #include <type_traits>
@ -134,28 +135,24 @@ private:
const BCMatrixRowMajor &bc, double time_step, double dx) const BCMatrixRowMajor &bc, double time_step, double dx)
-> DMatrixRowMajor; -> DMatrixRowMajor;
inline void fillMatrixFromRow(const DVectorRowMajor &alpha, void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
const BCVectorRowMajor &bc, int size, double dx, const DVectorRowMajor &alpha,
double time_step); const BCVectorRowMajor &bc, int size, double dx,
inline void fillVectorFromRow(const DVectorRowMajor &c, double time_step);
const DVectorRowMajor &alpha,
const BCVectorRowMajor &bc, void fillVectorFromRow(Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
const DVectorRowMajor &t0_c, int size, const DVectorRowMajor &alpha,
double dx, double time_step); const BCVectorRowMajor &bc,
const DVectorRowMajor &t0_c, int size, double dx,
double time_step);
void simulate3D(std::vector<double> &c); void simulate3D(std::vector<double> &c);
inline void reserveMemory(int size, int max_count_per_line);
inline static auto getBCFromFlux(Diffusion::boundary_condition bc, inline static auto getBCFromFlux(Diffusion::boundary_condition bc,
double neighbor_c, double neighbor_alpha) double neighbor_c, double neighbor_alpha)
-> double; -> double;
void solveLES();
void updateInternals(); void updateInternals();
Eigen::SparseMatrix<double> A_matrix;
Eigen::VectorXd b_vector;
Eigen::VectorXd x_vector;
double time_step; double time_step;
unsigned int grid_dim; unsigned int grid_dim;