preCICE v3.1.1
Loading...
Searching...
No Matches
Preconditioner.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Core>
4#include <numeric>
5#include <vector>
6
9#include "logging/Logger.hpp"
10#include "utils/assertion.hpp"
11
12namespace precice {
13namespace acceleration {
14namespace impl {
15
26public:
27 Preconditioner(int maxNonConstTimeWindows)
28 : _maxNonConstTimeWindows(maxNonConstTimeWindows)
29 {
30 }
31
33 virtual ~Preconditioner() {}
34
40 {
42
44 _subVectorSizes = svs;
45
46 size_t N = std::accumulate(_subVectorSizes.begin(), _subVectorSizes.end(), static_cast<std::size_t>(0));
47
48 // cannot do this already in the constructor as the size is unknown at that point
49 _weights.resize(N, 1.0);
50 _invWeights.resize(N, 1.0);
51 }
52
57 void apply(Eigen::MatrixXd &M, bool transpose)
58 {
60 if (transpose) {
61 PRECICE_ASSERT(M.cols() == (int) _weights.size(), M.cols(), _weights.size());
62 for (int i = 0; i < M.cols(); i++) {
63 for (int j = 0; j < M.rows(); j++) {
64 M(j, i) *= _weights[i];
65 }
66 }
67 } else {
68 PRECICE_ASSERT(M.rows() == (int) _weights.size(), M.rows(), (int) _weights.size());
69 for (int i = 0; i < M.cols(); i++) {
70 for (int j = 0; j < M.rows(); j++) {
71 M(j, i) *= _weights[j];
72 }
73 }
74 }
75 }
76
81 void revert(Eigen::MatrixXd &M, bool transpose)
82 {
84 //PRECICE_ASSERT(_needsGlobalWeights);
85 if (transpose) {
86 PRECICE_ASSERT(M.cols() == (int) _invWeights.size());
87 for (int i = 0; i < M.cols(); i++) {
88 for (int j = 0; j < M.rows(); j++) {
89 M(j, i) *= _invWeights[i];
90 }
91 }
92 } else {
93 PRECICE_ASSERT(M.rows() == (int) _invWeights.size(), M.rows(), (int) _invWeights.size());
94 for (int i = 0; i < M.cols(); i++) {
95 for (int j = 0; j < M.rows(); j++) {
96 M(j, i) *= _invWeights[j];
97 }
98 }
99 }
100 }
101
103 void apply(Eigen::MatrixXd &M)
104 {
106 PRECICE_ASSERT(M.rows() == (int) _weights.size(), M.rows(), (int) _weights.size());
107
108 // scale matrix M
109 for (int i = 0; i < M.cols(); i++) {
110 for (int j = 0; j < M.rows(); j++) {
111 M(j, i) *= _weights[j];
112 }
113 }
114 }
115
117 void apply(Eigen::VectorXd &v)
118 {
120
121 PRECICE_ASSERT(v.size() == (int) _weights.size());
122
123 // scale residual
124 for (int j = 0; j < v.size(); j++) {
125 v[j] *= _weights[j];
126 }
127 }
128
130 void revert(Eigen::MatrixXd &M)
131 {
133
134 PRECICE_ASSERT(M.rows() == (int) _weights.size());
135
136 // scale matrix M
137 for (int i = 0; i < M.cols(); i++) {
138 for (int j = 0; j < M.rows(); j++) {
139 M(j, i) *= _invWeights[j];
140 }
141 }
142 }
143
145 void revert(Eigen::VectorXd &v)
146 {
148
149 PRECICE_ASSERT(v.size() == (int) _weights.size());
150
151 // scale residual
152 for (int j = 0; j < v.size(); j++) {
153 v[j] *= _invWeights[j];
154 }
155 }
156
162 void update(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res)
163 {
165
166 // if number of allowed non-const time windows is exceeded, do not update weights
167 if (_frozen)
168 return;
169
170 // increment number of time windows that has been scaled with changing preconditioning weights
171 if (timeWindowComplete) {
174 _frozen = true;
175 }
176
177 // type specific update functionality
178 _update_(timeWindowComplete, oldValues, res);
179 }
180
183 {
185 return _requireNewQR;
186 }
187
190 {
191 _requireNewQR = false;
192 }
193
195 {
196 return _weights;
197 }
198
199 bool isConst()
200 {
201 return _frozen;
202 }
203
204protected:
207
210
213
218
221
223 bool _requireNewQR = false;
224
226 bool _frozen = false;
227
233 virtual void _update_(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res) = 0;
234
235private:
236 logging::Logger _log{"acceleration::Preconditioner"};
237};
238
239} // namespace impl
240} // namespace acceleration
241} // namespace precice
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
T accumulate(T... args)
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
Interface for preconditioner variants that can be applied to quasi-Newton acceleration schemes.
void apply(Eigen::MatrixXd &M, bool transpose)
Apply preconditioner to matrix.
void revert(Eigen::MatrixXd &M, bool transpose)
Apply inverse preconditioner to matrix.
int _nbNonConstTimeWindows
Counts the number of completed time windows with a non-const weighting.
std::vector< size_t > _subVectorSizes
Sizes of each sub-vector, i.e. each coupling data.
void newQRfulfilled()
to tell the preconditioner that QR-decomposition has been recomputed
void update(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res)
Update the scaling after every FSI iteration and require a new QR decomposition (if necessary)
virtual void _update_(bool timeWindowComplete, const Eigen::VectorXd &oldValues, const Eigen::VectorXd &res)=0
Update the scaling after every FSI iteration and require a new QR decomposition (if necessary)
virtual void initialize(std::vector< size_t > &svs)
initialize the preconditioner
virtual ~Preconditioner()
Destructor, empty.
void revert(Eigen::MatrixXd &M)
To transform balanced values back to physical values. Matrix version.
std::vector< double > _invWeights
Inverse weights (for efficiency reasons)
bool requireNewQR()
returns true if a QR decomposition from scratch is necessary
int _maxNonConstTimeWindows
maximum number of non-const time windows, i.e., after this number of time windows,...
std::vector< double > _weights
Weights used to scale the matrix V and the residual.
void revert(Eigen::VectorXd &v)
To transform balanced values back to physical values. Vector version.
bool _requireNewQR
True if a QR decomposition from scratch is necessary.
void apply(Eigen::VectorXd &v)
To transform physical values to balanced values. Vector version.
void apply(Eigen::MatrixXd &M)
To transform physical values to balanced values. Matrix version.
bool _frozen
True if _nbNonConstTimeWindows >= _maxNonConstTimeWindows, i.e., preconditioner is not updated any mo...
This class provides a lightweight logger.
Definition Logger.hpp:16
T empty(T... args)
T end(T... args)
Main namespace of the precice library.
static std::unique_ptr< precice::Participant > impl
Definition preciceC.cpp:21
T resize(T... args)
T size(T... args)