preCICE v3.1.2
Loading...
Searching...
No Matches
IQNILSAcceleration.cpp
Go to the documentation of this file.
2#include <Eigen/Core>
3#include <algorithm>
4#include <deque>
5#include <memory>
6#include <utility>
7
11#include "com/Communication.hpp"
12#include "com/SharedPointer.hpp"
15#include "logging/LogMacros.hpp"
17#include "utils/Helpers.hpp"
18#include "utils/IntraComm.hpp"
19#include "utils/assertion.hpp"
20
21//#include "utils/NumericalCompare.hpp"
22
24
25namespace precice::acceleration {
26
28 double initialRelaxation,
29 bool forceInitialRelaxation,
30 int maxIterationsUsed,
31 int pastTimeWindowsReused,
32 int filter,
33 double singularityLimit,
34 std::vector<int> dataIDs,
35 impl::PtrPreconditioner preconditioner)
36 : BaseQNAcceleration(initialRelaxation, forceInitialRelaxation, maxIterationsUsed, pastTimeWindowsReused,
37 filter, singularityLimit, std::move(dataIDs), std::move(preconditioner))
38{
39}
40
42 const DataMap &cplData)
43{
44 // do common QN acceleration initialization
46
47 // Fetch secondary data IDs, to be relaxed with same coefficients from IQN-ILS
48 for (const DataMap::value_type &pair : cplData) {
49 if (not utils::contained(pair.first, _dataIDs)) {
50 int secondaryEntries = pair.second->getSize();
51 utils::append(_secondaryOldXTildes[pair.first], Eigen::VectorXd(Eigen::VectorXd::Zero(secondaryEntries)));
52 }
53 }
54}
55
57 const DataMap &cplData)
58{
59 // Compute residuals of secondary data
60 for (int id : _secondaryDataIDs) {
61 Eigen::VectorXd &secResiduals = _secondaryResiduals[id];
62 PtrCouplingData data = cplData.at(id);
63 PRECICE_ASSERT(secResiduals.size() == data->getSize(),
64 secResiduals.size(), data->getSize());
65 secResiduals = data->values();
66 secResiduals -= data->previousIteration();
67 }
68
70 // constant relaxation: for secondary data called from base class
71 } else {
72 if (not _firstIteration) {
73 bool columnLimitReached = getLSSystemCols() == _maxIterationsUsed;
74 bool overdetermined = getLSSystemCols() <= getLSSystemRows();
75 if (not columnLimitReached && overdetermined) {
76
77 // Append column for secondary W matrices
78 for (int id : _secondaryDataIDs) {
80 }
81 } else {
82 // Shift column for secondary W matrices
83 for (int id : _secondaryDataIDs) {
85 }
86 }
87
88 // Compute delta_x_tilde for secondary data
89 for (int id : _secondaryDataIDs) {
90 Eigen::MatrixXd &secW = _secondaryMatricesW[id];
91 PRECICE_ASSERT(secW.rows() == cplData.at(id)->getSize(), secW.rows(), cplData.at(id)->getSize());
92 secW.col(0) = cplData.at(id)->values();
93 secW.col(0) -= _secondaryOldXTildes[id];
94 }
95 }
96
97 // Store x_tildes for secondary data
98 for (int id : _secondaryDataIDs) {
99 PRECICE_ASSERT(_secondaryOldXTildes[id].size() == cplData.at(id)->getSize(),
100 _secondaryOldXTildes[id].size(), cplData.at(id)->getSize());
101 _secondaryOldXTildes[id] = cplData.at(id)->values();
102 }
103 }
104
105 // call the base method for common update of V, W matrices
107}
108
110 const DataMap &cplData)
111{
112 //Store x_tildes for secondary data
113 for (int id : _secondaryDataIDs) {
114 PRECICE_ASSERT(_secondaryOldXTildes.at(id).size() == cplData.at(id)->getSize(),
115 _secondaryOldXTildes.at(id).size(), cplData.at(id)->getSize());
116 _secondaryOldXTildes[id] = cplData.at(id)->values();
117 }
118
119 // Perform underrelaxation with initial relaxation factor for secondary data
120 for (int id : _secondaryDataIDs) {
121 PtrCouplingData data = cplData.at(id);
122 Eigen::VectorXd &values = data->values();
123 values *= _initialRelaxation; // new * omg
124 Eigen::VectorXd &secResiduals = _secondaryResiduals[id];
125 secResiduals = data->previousIteration(); // old
126 secResiduals *= 1.0 - _initialRelaxation; // (1-omg) * old
127 values += secResiduals; // (1-omg) * old + new * omg
128 }
129}
130
131void IQNILSAcceleration::computeQNUpdate(const DataMap &cplData, Eigen::VectorXd &xUpdate)
132{
134 PRECICE_DEBUG(" Compute Newton factors");
135
136 // Calculate QR decomposition of matrix V and solve Rc = -Qr
137 Eigen::VectorXd c;
138
139 // for procs with no vertices,
140 // qrV.cols() = getLSSystemCols() and _qrV.rows() = 0
141 auto Q = _qrV.matrixQ();
142 auto R = _qrV.matrixR();
143
146 PRECICE_ASSERT(_qrV.rows() == 0, _qrV.rows());
147 PRECICE_ASSERT(Q.size() == 0, Q.size());
148 }
149
150 Eigen::VectorXd _local_b = Eigen::VectorXd::Zero(_qrV.cols());
151 Eigen::VectorXd _global_b;
152
153 // need to scale the residual to compensate for the scaling in c = R^-1 * Q^T * P^-1 * residual'
154 // it is also possible to apply the inverse scaling weights from the right to the vector c
156 _local_b = Q.transpose() * _residuals;
158 _local_b *= -1.0; // = -Qr
159
160 PRECICE_ASSERT(c.size() == 0, c.size());
161 // reserve memory for c
162 utils::append(c, Eigen::VectorXd(Eigen::VectorXd::Zero(_local_b.size())));
163
164 // compute rhs Q^T*res in parallel
166 PRECICE_ASSERT(Q.cols() == getLSSystemCols(), Q.cols(), getLSSystemCols());
167 // back substitution
168 c = R.triangularView<Eigen::Upper>().solve<Eigen::OnTheLeft>(_local_b);
169 } else {
173 PRECICE_ASSERT(Q.cols() == getLSSystemCols(), Q.cols(), getLSSystemCols());
174 }
175 PRECICE_ASSERT(_local_b.size() == getLSSystemCols(), _local_b.size(), getLSSystemCols());
176
178 PRECICE_ASSERT(_global_b.size() == 0, _global_b.size());
179 }
180 utils::append(_global_b, Eigen::VectorXd(Eigen::VectorXd::Zero(_local_b.size())));
181
182 // do a reduce operation to sum up all the _local_b vectors
183 utils::IntraComm::reduceSum(_local_b, _global_b);
184
185 // back substitution R*c = b only on the primary rank
187 c = R.triangularView<Eigen::Upper>().solve<Eigen::OnTheLeft>(_global_b);
188 }
189
190 // broadcast coefficients c to all secondary ranks
192 }
193
194 PRECICE_DEBUG(" Apply Newton factors");
195 // compute x updates from W and coefficients c, i.e, xUpdate = c*W
196 xUpdate = _matrixW * c;
197
202 // If the previous time window converged within one single iteration, nothing was added
203 // to the LS system matrices and they need to be restored from the backup at time T-2
205 PRECICE_DEBUG(" Last time window converged after one iteration. Need to restore the secondaryMatricesW from backup.");
207 }
208
209 // Perform QN relaxation for secondary data
210 for (int id : _secondaryDataIDs) {
211 PtrCouplingData data = cplData.at(id);
212 auto & values = data->values();
213 PRECICE_ASSERT(_secondaryMatricesW[id].cols() == c.size(), _secondaryMatricesW[id].cols(), c.size());
214 values = _secondaryMatricesW[id] * c;
215 PRECICE_ASSERT(data->getSize() == data->getPreviousIterationSize(), data->getSize(), data->getPreviousIterationSize());
216 values += data->previousIteration();
217 PRECICE_ASSERT(data->getSize() == _secondaryResiduals[id].size(), data->getSize(), _secondaryResiduals[id].size());
218 values += _secondaryResiduals[id];
219 }
220
221 // pending deletion: delete old secondaryMatricesW
223 // save current secondaryMatrix data in case the coupling for the next time window will terminate
224 // after the first iteration (no new data, i.e., V = W = 0)
225 if (getLSSystemCols() > 0) {
227 }
228 for (int id : _secondaryDataIDs) {
229 _secondaryMatricesW[id].resize(0, 0);
230 }
231 }
232}
233
235 const DataMap &cplData)
236{
238 if (_matrixCols.empty()) {
239 PRECICE_WARN("The IQN matrix has no columns.");
240 } else {
241 if (_matrixCols.front() == 0) { // Did only one iteration
243 }
244
245 if (_timeWindowsReused == 0) {
247 for (int id : _secondaryDataIDs) {
248 _secondaryMatricesW[id].resize(0, 0);
249 }
250 } else {
256 }
257 } else if (static_cast<int>(_matrixCols.size()) > _timeWindowsReused) {
258 int toRemove = _matrixCols.back();
259 for (int id : _secondaryDataIDs) {
260 Eigen::MatrixXd &secW = _secondaryMatricesW[id];
261 PRECICE_ASSERT(secW.cols() > toRemove, secW, toRemove, id);
262 for (int i = 0; i < toRemove; i++) {
263 utils::removeColumnFromMatrix(secW, secW.cols() - 1);
264 }
265 }
266 }
267 }
268}
269
271 int columnIndex)
272{
273 PRECICE_ASSERT(_matrixV.cols() > 1);
274 // remove column from secondary Data Matrix W
275 for (int id : _secondaryDataIDs) {
277 }
278
280}
281} // namespace precice::acceleration
#define PRECICE_WARN(...)
Definition LogMacros.hpp:11
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T at(T... args)
T back(T... args)
Base Class for quasi-Newton acceleration schemes.
virtual void updateDifferenceMatrices(const DataMap &cplData)
Updates the V, W matrices (as well as the matrices for the secondary data)
double _initialRelaxation
Constant relaxation factor used for first iteration.
int _timeWindowsReused
Maximum number of old time windows (with data values) kept.
virtual int getLSSystemCols() const
: computes number of cols in least squares system, i.e, number of cols in _matrixV,...
Eigen::MatrixXd _matrixV
Stores residual deltas.
Eigen::VectorXd _residuals
Current iteration residuals of IQN data. Temporary.
virtual void removeMatrixColumn(int columnIndex)
Removes one iteration from V,W matrices and adapts _matrixCols.
Eigen::MatrixXd _matrixW
Stores x tilde deltas, where x tilde are values computed by solvers.
std::deque< int > _matrixCols
Indices (of columns in W, V matrices) of 1st iterations of time windows.
bool _firstIteration
Indicates the first iteration, where constant relaxation is used.
std::map< int, Eigen::VectorXd > _secondaryResiduals
Current iteration residuals of secondary data.
impl::PtrPreconditioner _preconditioner
Preconditioner for least-squares system if vectorial system is used.
impl::QRFactorization _qrV
Stores the current QR decomposition ov _matrixV, can be updated via deletion/insertion of columns.
int _maxIterationsUsed
Maximum number of old data iterations kept.
std::vector< int > _secondaryDataIDs
Data IDs of data not involved in IQN coefficient computation.
std::vector< int > _dataIDs
Data IDs of data to be involved in the IQN algorithm.
virtual void initialize(const DataMap &cplData)
Initializes the acceleration.
std::map< int, Eigen::MatrixXd > _secondaryMatricesWBackup
virtual void initialize(const DataMap &cplData)
Initializes the acceleration.
virtual void updateDifferenceMatrices(const DataMap &cplData)
updates the V, W matrices (as well as the matrices for the secondary data)
virtual void removeMatrixColumn(int columnIndex)
Removes one iteration from V,W matrices and adapts _matrixCols.
virtual void computeUnderrelaxationSecondaryData(const DataMap &cplData)
computes underrelaxation for the secondary data
IQNILSAcceleration(double initialRelaxation, bool forceInitialRelaxation, int maxIterationsUsed, int pastTimeWindowsReused, int filter, double singularityLimit, std::vector< int > dataIDs, impl::PtrPreconditioner preconditioner)
std::map< int, Eigen::MatrixXd > _secondaryMatricesW
virtual void computeQNUpdate(const DataMap &cplData, Eigen::VectorXd &xUpdate)
computes the IQN-ILS update using QR decomposition
std::map< int, Eigen::VectorXd > _secondaryOldXTildes
Secondary data solver output from last iteration.
virtual void specializedIterationsConverged(const DataMap &cplData)
Marks a iteration sequence as converged.
Eigen::MatrixXd & matrixR()
returns a matrix representation of the upper triangular matrix R
Eigen::MatrixXd & matrixQ()
returns a matrix representation of the orthogonal matrix Q
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static void broadcast(bool &value)
static bool isParallel()
True if this process is running in parallel.
Definition IntraComm.cpp:62
static com::PtrCommunication & getCommunication()
Intra-participant communication.
Definition IntraComm.hpp:31
static void reduceSum(precice::span< const double > sendData, precice::span< double > rcvData)
T empty(T... args)
T front(T... args)
contains implementations of acceleration schemes.
void removeColumnFromMatrix(Eigen::MatrixXd &A, int col)
bool contained(const ELEMENT_T &element, const std::vector< ELEMENT_T > &vec)
Returns true, if given element is in vector, otherwise false.
Definition Helpers.hpp:39
void appendFront(Eigen::MatrixXd &A, Eigen::VectorXd &v)
void append(Eigen::VectorXd &v, double value)
void shiftSetFirst(Eigen::MatrixXd &A, const Eigen::VectorXd &v)
STL namespace.
T pop_front(T... args)
T size(T... args)