48 Foam::DMDModels::STDMD::modeSorterType
50 Foam::DMDModels::STDMD::modeSorterTypeNames
52 { modeSorterType::FIRST_SNAPSHOT,
"firstSnapshot" },
53 { modeSorterType::KIEWAT ,
"kiewat" },
54 { modeSorterType::KOU_ZHANG ,
"kouZhang" }
60 Foam::scalar Foam::DMDModels::STDMD::L2norm(
const RMatrix& z)
const 67 <<
"Input matrix is not a column vector." 71 const bool noSqrt =
true;
72 scalar result = z.columnNorm(0, noSqrt);
73 reduce(result, sumOp<scalar>());
85 List<scalar> dz(Q_.n());
86 const label sz0 = ez.
m();
87 const label sz1 = dz.size();
89 for (label i = 0; i < nGramSchmidt_; ++i)
93 for (label i = 0; i < sz0; ++i)
95 for (label j = 0; j < sz1; ++j)
97 dz[j] += Q_(i,j)*ez(i,0);
101 reduce(dz, sumOp<List<scalar>>());
104 for (label i = 0; i < sz0; ++i)
107 for (label j = 0; j < sz1; ++j)
119 void Foam::DMDModels::STDMD::expand(
const RMatrix& ez,
const scalar ezNorm)
121 Info<<
tab <<
"Expanding orthonormal basis for field: " << fieldName_
125 Q_.resize(Q_.m(), Q_.n() + 1);
126 Q_.subColumn(Q_.n() - 1) = ez/ezNorm;
129 G_.resize(G_.m() + 1);
133 void Foam::DMDModels::STDMD::updateG(
const RMatrix& z)
135 List<scalar> zTilde(Q_.n(),
Zero);
136 const label sz0 = z.m();
137 const label sz1 = zTilde.size();
140 for (label i = 0; i < sz0; ++i)
142 for (label j = 0; j < sz1; ++j)
144 zTilde[j] += Q_(i,j)*z(i,0);
148 reduce(zTilde, sumOp<List<scalar>>());
151 for (label i = 0; i < G_.m(); ++i)
153 for (label j = 0; j < G_.n(); ++j)
155 G_(i, j) += zTilde[i]*zTilde[j];
161 void Foam::DMDModels::STDMD::compress()
163 Info<<
tab <<
"Compressing orthonormal basis for field: " << fieldName_
166 RMatrix q(1, 1,
Zero);
170 const bool symmetric =
true;
171 const EigenMatrix<scalar> EM(G_, symmetric);
172 const SquareMatrix<scalar>& EVecs = EM.EVecs();
173 DiagonalMatrix<scalar> EVals(EM.EValsRe());
176 const auto descend = [&](scalar a, scalar
b){
return a >
b; };
177 const List<label> permutation(EVals.sortPermutation(descend));
178 EVals.applyPermutation(permutation);
179 EVals.resize(EVals.size() - 1);
182 G_ = SMatrix(maxRank_,
Zero);
185 q.resize(Q_.n(), maxRank_);
186 for (label i = 0; i < maxRank_; ++i)
188 q.subColumn(i) = EVecs.subColumn(permutation[i]);
200 reducedKoopmanOperator()
208 QRMatrix<RMatrix> QRM
211 QRMatrix<RMatrix>::modes::ECONOMY,
212 QRMatrix<RMatrix>::outputs::ONLY_R
214 RMatrix&
R = QRM.R();
220 RMatrix A1(1, 1,
Zero);
230 PstreamBuffers pBufs;
231 pBufs.allowClearRecv(
false);
234 const label procNoInSubset = myProcNo % nAgglomerationProcs_;
237 if (procNoInSubset != 0)
239 const label procNoSubsetMaster = myProcNo - procNoInSubset;
241 UOPstream toSubsetMaster(procNoSubsetMaster, pBufs);
242 toSubsetMaster <<
Rx;
247 pBufs.finishedSends();
250 if (procNoInSubset == 0)
257 label nbr = myProcNo + 1;
259 nbr < myProcNo + nAgglomerationProcs_
267 UIPstream fromNbr(nbr, pBufs);
271 if (recvMtrx.size() > 0)
273 Rx.resize(
Rx.m() + recvMtrx.m(),
Rx.n());
276 Rx.m() - recvMtrx.m(),
277 Rx.n() - recvMtrx.n()
286 QRMatrix<RMatrix> QRM
288 QRMatrix<RMatrix>::modes::ECONOMY,
289 QRMatrix<RMatrix>::outputs::ONLY_R,
290 QRMatrix<RMatrix>::pivoting::FALSE,
293 RMatrix&
R = QRM.R();
302 if (procNoInSubset == 0)
311 pBufs.finishedSends();
321 nbr += nAgglomerationProcs_
326 UIPstream fromSubsetMaster(nbr, pBufs);
327 fromSubsetMaster >> recvMtrx;
330 if (recvMtrx.size() > 0)
332 Rx.resize(
Rx.m() + recvMtrx.m(),
Rx.n());
335 Rx.m() - recvMtrx.m(),
336 Rx.n() - recvMtrx.n()
343 QRMatrix<RMatrix> QRM
345 QRMatrix<RMatrix>::modes::ECONOMY,
346 QRMatrix<RMatrix>::outputs::ONLY_R,
347 QRMatrix<RMatrix>::pivoting::FALSE,
350 RMatrix&
R = QRM.R();
366 SMatrix A2(Qupper_ & Qlower_);
367 reduce(A2, sumOp<SMatrix>());
371 SMatrix A3(Qupper_ & Qupper_);
372 reduce(A3, sumOp<SMatrix>());
377 return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
389 SMatrix A2(Qupper_ & Qlower_);
393 SMatrix A3(Qupper_ & Qupper_);
396 return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
401 bool Foam::DMDModels::STDMD::eigendecomposition(SMatrix& Atilde)
408 Info<<
tab <<
"Computing eigendecomposition" <<
endl;
411 const EigenMatrix<scalar> EM(Atilde);
412 const DiagonalMatrix<scalar>& evalsRe = EM.EValsRe();
413 const DiagonalMatrix<scalar>& evalsIm = EM.EValsIm();
415 evals_.resize(evalsRe.size());
416 evecs_ = RCMatrix(EM.complexEVecs());
420 for (
auto& eval : evals_)
422 eval =
complex(evalsRe[i], evalsIm[i]);
428 List<complex>
cp(evals_.size());
442 <<
"No eigenvalue with mag(eigenvalue) larger than " 443 <<
"minEval = " << minEval_ <<
" was found." <<
nl 444 <<
" Input field may contain only zero-value elements," <<
nl 445 <<
" e.g. no-slip velocity condition on a given patch." <<
nl 446 <<
" Otherwise, please decrease the value of minEval." <<
nl 447 <<
" Skipping further dynamics/mode computations." <<
nl 471 void Foam::DMDModels::STDMD::frequencies()
477 freqs_.resize(evals_.size());
480 auto frequencyEquation =
495 Info<<
tab <<
"Computing frequency indices" <<
endl;
498 auto margin = [&](
const scalar&
x){
return (fMin_ <=
x &&
x < fMax_); };
500 auto it =
std::find_if(freqs_.cbegin(), freqs_.cend(), margin);
502 while (it != freqs_.end())
505 it =
std::find_if(std::next(it), freqs_.cend(), margin);
513 void Foam::DMDModels::STDMD::amplitudes()
515 IOField<scalar> snapshot0
519 "snapshot0_" + name_ +
"_" + fieldName_,
530 List<scalar> tmp(Qupper_.n(),
Zero);
531 const label sz0 = snapshot0.size();
532 const label sz1 = tmp.size();
534 for (label i = 0; i < sz0; ++i)
536 for (label j = 0; j < sz1; ++j)
538 tmp[j] += Qupper_(i,j)*snapshot0[i];
544 List<scalar>
R(Qupper_.n(),
Zero);
545 for (label i = 0; i < sz1; ++i)
547 for (label j = 0; j <
R.size(); ++j)
549 R[j] += RxInv_(i,j)*tmp[i];
554 reduce(
R, sumOp<List<scalar>>());
560 amps_.resize(
R.size());
564 for (label i = 0; i < amps_.size(); ++i)
566 for (label j = 0; j <
R.size(); ++j)
568 amps_[i] += pEvecs(i,j)*
R[j];
576 void Foam::DMDModels::STDMD::magnitudes()
582 mags_.resize(amps_.size());
584 Info<<
tab <<
"Sorting modes with ";
588 case modeSorterType::FIRST_SNAPSHOT:
590 Info<<
"method of first snapshot" <<
endl;
602 case modeSorterType::KIEWAT:
604 Info<<
"modified weighted amplitude scaling method" <<
endl;
608 const scalar modeNorm = 1;
610 List<scalar> w(step_);
611 std::iota(w.begin(), w.end(), 1);
612 w =
sin(
twoPi/step_*(w - 1 - 0.25*step_))*pr + pr;
616 mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
622 case modeSorterType::KOU_ZHANG:
624 Info<<
"weighted amplitude scaling method" <<
endl;
626 const scalar modeNorm = 1;
627 const List<scalar> w(step_, 1);
631 mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
641 Info<<
tab <<
"Computing magnitude indices" <<
endl;
646 [&](
const label i1,
const label i2)
648 return !(mags_[i1] < mags_[i2]);
651 std::sort(magsi_.begin(), magsi_.end(), descend);
658 Foam::scalar Foam::DMDModels::STDMD::sorter
660 const List<scalar>& weight,
663 const scalar modeNorm
667 if (!(
mag(eigenvalue) < GREAT &&
mag(eigenvalue) > VSMALL))
669 Info<<
" Returning zero magnitude for mag(eigenvalue) = " 676 if (
mag(eigenvalue)*step_ > sortLimiter_)
678 Info<<
" Returning zero magnitude for" 679 <<
" mag(eigenvalue) = " <<
mag(eigenvalue)
680 <<
" current index = " << step_
681 <<
" sortLimiter = " << sortLimiter_
687 scalar magnitude = 0;
689 for (label j = 0; j < step_; ++j)
691 magnitude += weight[j]*modeNorm*
mag(amplitude*
pow(eigenvalue, j + 1));
698 bool Foam::DMDModels::STDMD::dynamics()
700 SMatrix Atilde(reducedKoopmanOperator());
703 if(!eigendecomposition(Atilde))
720 bool Foam::DMDModels::STDMD::modes()
724 bool processed =
false;
725 processed = processed || modes<scalar>();
726 processed = processed || modes<vector>();
727 processed = processed || modes<sphericalTensor>();
728 processed = processed || modes<symmTensor>();
729 processed = processed || modes<tensor>();
742 Info<<
tab <<
"Writing objects of dynamics" <<
endl;
746 autoPtr<OFstream> osPtr =
749 fileName +
"_" + fieldName_,
750 mesh_.time().timeOutputValue()
752 OFstream&
os = osPtr.ref();
756 for (
const auto& i : labelRange(0, freqs_.size()))
758 os << freqs_[i] <<
tab 760 << amps_[i].real() <<
tab 761 << amps_[i].imag() <<
tab 762 << evals_[i].real() <<
tab 763 << evals_[i].imag() <<
endl;
767 Info<<
tab <<
"Ends output processing for field: " << fieldName_
772 void Foam::DMDModels::STDMD::writeFileHeader(Ostream&
os)
const 775 writeCommented(
os,
"Frequency");
776 writeTabbed(
os,
"Magnitude");
777 writeTabbed(
os,
"Amplitude (real)");
778 writeTabbed(
os,
"Amplitude (imag)");
779 writeTabbed(
os,
"Eigenvalue (real)");
780 writeTabbed(
os,
"Eigenvalue (imag)");
785 void Foam::DMDModels::STDMD::filter()
787 Info<<
tab <<
"Filtering objects of dynamics" <<
endl;
790 filterIndexed(evals_, magsi_);
791 filterIndexed(evecs_, magsi_);
792 filterIndexed(freqs_, magsi_);
793 filterIndexed(amps_, magsi_);
794 filterIndexed(mags_, magsi_);
797 if (freqs_.size() > nModes_)
799 evals_.resize(nModes_);
800 evecs_.resize(evecs_.m(), nModes_);
801 freqs_.resize(nModes_);
802 amps_.resize(nModes_);
803 mags_.resize(nModes_);
820 modeSorterTypeNames.getOrDefault
824 modeSorterType::FIRST_SNAPSHOT
859 nAgglomerationProcs_(20),
868 writeToFile_ =
dict.getOrDefault<
bool>(
"writeToFile",
true);
871 modeSorterTypeNames.getOrDefault
875 modeSorterType::FIRST_SNAPSHOT
885 dict.getCheckOrDefault<scalar>
890 *mesh_.time().deltaTValue()
902 dict.getCheckOrDefault<label>
910 dict.getCheckOrDefault<label>
920 dict.getCheckOrDefault<label>
927 nAgglomerationProcs_ =
928 dict.getCheckOrDefault<label>
930 "nAgglomerationProcs",
935 Info<<
tab <<
"Settings are read for:" <<
nl 936 <<
tab <<
" field: " << fieldName_ <<
nl 937 <<
tab <<
" modeSorter: " << modeSorterTypeNames[modeSorter_] <<
nl 938 <<
tab <<
" nModes: " << nModes_ <<
nl 939 <<
tab <<
" maxRank: " << maxRank_ <<
nl 940 <<
tab <<
" nGramSchmidt: " << nGramSchmidt_ <<
nl 941 <<
tab <<
" fMin: " << fMin_ <<
nl 942 <<
tab <<
" fMax: " << fMax_ <<
nl 943 <<
tab <<
" minBasis: " << minBasis_ <<
nl 944 <<
tab <<
" minEVal: " << minEval_ <<
nl 945 <<
tab <<
" sortLimiter: " << sortLimiter_ <<
nl 946 <<
tab <<
" nAgglomerationProcs: " << nAgglomerationProcs_ <<
nl 955 const scalar norm = L2norm(z);
962 const label nSnap = z.m()/2;
964 mesh_.time().timeName(mesh_.time().startTime().value());
975 "snapshot0_" + name_ +
"_" + fieldName_,
994 mesh_.time().writeFormat(),
995 mesh_.time().writeCompression()
998 fileHandler().writeObject(snapshot0, streamOpt,
true);
1003 G_(0,0) =
sqr(norm);
1019 const RMatrix ez(orthonormalise(z));
1022 const scalar ezNorm = L2norm(ez);
1023 const scalar zNorm = L2norm(z);
1025 if (ezNorm/zNorm > minBasis_)
1035 if (Q_.n() >= maxRank_)
1049 const label nSnap = Q_.m()/2;
1052 Qupper_ = RMatrix(Q_.subMatrix(0, 0,
max(nSnap, 1)));
1053 Qlower_ = RMatrix(Q_.subMatrix(nSnap, 0,
max(nSnap, 1)));
1065 writeToFile(
word(
"dynamics"));
1069 writeToFile(
word(
"filtered_dynamics"));
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
static void writeHeader(Ostream &os, const word &fieldName)
dimensionedScalar log(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
virtual bool update(const RMatrix &z)
Incremental orthonormal basis update (K:Fig. 15)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Streaming Total Dynamic Mode Decomposition (i.e. STDMD) is a variant of dynamic mode decomposition...
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual bool fit()
Compute and write modes and mode dynamics of model data members.
constexpr char nl
The newline '\n' character (0x0a)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool cp(const fileName &src, const fileName &dst, const bool followLink=true)
Copy the source to the destination (recursively if necessary).
A traits class, which is primarily used for primitives and vector-space.
static bool & parRun() noexcept
Test if this a parallel run.
refPtr< fileOperation > fileHandler(std::nullptr_t)
Delete current file handler - forwards to fileOperation::handler()
scalar distance(const vector &p1, const vector &p2)
constexpr char tab
The tab '\t' character(0x09)
A simple container for options an IOstream can normally have.
virtual bool initialise(const RMatrix &z)
Initialise 'Q' and 'G' (both require the first two snapshots)
Ignore writing from objectRegistry::writeObject()
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
Macros for easy insertion into run-time selection tables.
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all communicator ranks. Does nothing in non-paral...
#define forAll(list, i)
Loop across all elements in list.
virtual bool read(const dictionary &dict)
Read STDMD settings.
STDMD(const fvMesh &mesh, const word &name, const dictionary &dict)
Construct from components.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
constexpr scalar twoPi(2 *M_PI)
A class for handling words, derived from Foam::string.
void sort(UList< T > &list)
Sort the list.
static constexpr int masterNo() noexcept
Relative rank for the master process - is always 0.
static MinMax< scalar > ge(const scalar &minVal)
A semi-infinite range from minVal to the type max.
A List of wordRe with additional matching capabilities.
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements, derived from Matrix.
virtual bool writeToFile() const
Flag to allow writing to file.
iterator begin() noexcept
Return an iterator to begin traversing the UList.
dimensionedScalar sin(const dimensionedScalar &ds)
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
label m() const noexcept
The number of rows.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Mesh data needed to do the Finite Volume discretisation.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
messageStream Info
Information stream (stdout output on master, null elsewhere)
label find_if(const ListType &input, const UnaryPredicate &pred, const label start=0)
Find index of the first occurrence that satisfies the predicate.
Abstract base class for DMD models to handle DMD characteristics for the DMD function object...
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
tensor Rx(const scalar omega)
Rotational transformation tensor about the x-axis by omega radians.
Defines the attributes of an object for which implicit objectRegistry management is supported...
A primitive field of type <T> with automated input and output.
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix...
Do not request registration (bool: false)
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
string expand(const std::string &s, const HashTable< string > &mapping, const char sigil='$')
Expand occurrences of variables according to the mapping and return the expanded string.
static constexpr const zero Zero
Global zero (0)