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);
232 const label procNoInSubset = myProcNo % nAgglomerationProcs_;
235 if (procNoInSubset != 0)
237 const label procNoSubsetMaster = myProcNo - procNoInSubset;
239 UOPstream toSubsetMaster(procNoSubsetMaster, pBufs);
240 toSubsetMaster <<
Rx;
245 pBufs.finishedSends();
248 if (procNoInSubset == 0)
255 label nbr = myProcNo + 1;
257 nbr < myProcNo + nAgglomerationProcs_
265 UIPstream fromNbr(nbr, pBufs);
269 if (recvMtrx.size() > 0)
271 Rx.resize(
Rx.m() + recvMtrx.m(),
Rx.n());
274 Rx.m() - recvMtrx.m(),
275 Rx.n() - recvMtrx.n()
284 QRMatrix<RMatrix> QRM
286 QRMatrix<RMatrix>::modes::ECONOMY,
287 QRMatrix<RMatrix>::outputs::ONLY_R,
288 QRMatrix<RMatrix>::pivoting::FALSE,
291 RMatrix&
R = QRM.R();
300 if (procNoInSubset == 0)
309 pBufs.finishedSends();
319 nbr += nAgglomerationProcs_
324 UIPstream fromSubsetMaster(nbr, pBufs);
325 fromSubsetMaster >> recvMtrx;
328 if (recvMtrx.size() > 0)
330 Rx.resize(
Rx.m() + recvMtrx.m(),
Rx.n());
333 Rx.m() - recvMtrx.m(),
334 Rx.n() - recvMtrx.n()
341 QRMatrix<RMatrix> QRM
343 QRMatrix<RMatrix>::modes::ECONOMY,
344 QRMatrix<RMatrix>::outputs::ONLY_R,
345 QRMatrix<RMatrix>::pivoting::FALSE,
348 RMatrix&
R = QRM.R();
364 SMatrix A2(Qupper_ & Qlower_);
365 reduce(A2, sumOp<SMatrix>());
369 SMatrix A3(Qupper_ & Qupper_);
370 reduce(A3, sumOp<SMatrix>());
375 return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
387 SMatrix A2(Qupper_ & Qlower_);
391 SMatrix A3(Qupper_ & Qupper_);
394 return SMatrix(RxInv_ & ((A2*(G_*A3))*A1));
399 bool Foam::DMDModels::STDMD::eigendecomposition(SMatrix& Atilde)
406 Info<<
tab <<
"Computing eigendecomposition" <<
endl;
409 const EigenMatrix<scalar> EM(Atilde);
410 const DiagonalMatrix<scalar>& evalsRe = EM.EValsRe();
411 const DiagonalMatrix<scalar>& evalsIm = EM.EValsIm();
413 evals_.resize(evalsRe.size());
414 evecs_ = RCMatrix(EM.complexEVecs());
418 for (
auto& eval : evals_)
420 eval =
complex(evalsRe[i], evalsIm[i]);
426 List<complex>
cp(evals_.size());
440 <<
"No eigenvalue with mag(eigenvalue) larger than " 441 <<
"minEval = " << minEval_ <<
" was found." <<
nl 442 <<
" Input field may contain only zero-value elements," <<
nl 443 <<
" e.g. no-slip velocity condition on a given patch." <<
nl 444 <<
" Otherwise, please decrease the value of minEval." <<
nl 445 <<
" Skipping further dynamics/mode computations." <<
nl 469 void Foam::DMDModels::STDMD::frequencies()
475 freqs_.resize(evals_.size());
478 auto frequencyEquation =
493 Info<<
tab <<
"Computing frequency indices" <<
endl;
496 auto margin = [&](
const scalar&
x){
return (fMin_ <=
x &&
x < fMax_); };
498 auto it =
std::find_if(freqs_.cbegin(), freqs_.cend(), margin);
500 while (it != freqs_.end())
503 it =
std::find_if(std::next(it), freqs_.cend(), margin);
511 void Foam::DMDModels::STDMD::amplitudes()
513 IOField<scalar> snapshot0
517 "snapshot0_" + name_ +
"_" + fieldName_,
528 List<scalar> tmp(Qupper_.n(),
Zero);
529 const label sz0 = snapshot0.size();
530 const label sz1 = tmp.size();
532 for (label i = 0; i < sz0; ++i)
534 for (label j = 0; j < sz1; ++j)
536 tmp[j] += Qupper_(i,j)*snapshot0[i];
542 List<scalar>
R(Qupper_.n(),
Zero);
543 for (label i = 0; i < sz1; ++i)
545 for (label j = 0; j <
R.size(); ++j)
547 R[j] += RxInv_(i,j)*tmp[i];
552 reduce(
R, sumOp<List<scalar>>());
558 amps_.resize(
R.size());
562 for (label i = 0; i < amps_.size(); ++i)
564 for (label j = 0; j <
R.size(); ++j)
566 amps_[i] += pEvecs(i,j)*
R[j];
574 void Foam::DMDModels::STDMD::magnitudes()
580 mags_.resize(amps_.size());
582 Info<<
tab <<
"Sorting modes with ";
586 case modeSorterType::FIRST_SNAPSHOT:
588 Info<<
"method of first snapshot" <<
endl;
600 case modeSorterType::KIEWAT:
602 Info<<
"modified weighted amplitude scaling method" <<
endl;
606 const scalar modeNorm = 1;
608 List<scalar> w(step_);
609 std::iota(w.begin(), w.end(), 1);
610 w =
sin(
twoPi/step_*(w - 1 - 0.25*step_))*pr + pr;
614 mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
620 case modeSorterType::KOU_ZHANG:
622 Info<<
"weighted amplitude scaling method" <<
endl;
624 const scalar modeNorm = 1;
625 const List<scalar> w(step_, 1);
629 mags_[i] = sorter(w, amps_[i], evals_[i], modeNorm);
639 Info<<
tab <<
"Computing magnitude indices" <<
endl;
644 [&](
const label i1,
const label i2)
646 return !(mags_[i1] < mags_[i2]);
649 std::sort(magsi_.begin(), magsi_.end(), descend);
656 Foam::scalar Foam::DMDModels::STDMD::sorter
658 const List<scalar>& weight,
661 const scalar modeNorm
665 if (!(
mag(eigenvalue) < GREAT &&
mag(eigenvalue) > VSMALL))
667 Info<<
" Returning zero magnitude for mag(eigenvalue) = " 674 if (
mag(eigenvalue)*step_ > sortLimiter_)
676 Info<<
" Returning zero magnitude for" 677 <<
" mag(eigenvalue) = " <<
mag(eigenvalue)
678 <<
" current index = " << step_
679 <<
" sortLimiter = " << sortLimiter_
685 scalar magnitude = 0;
687 for (label j = 0; j < step_; ++j)
689 magnitude += weight[j]*modeNorm*
mag(amplitude*
pow(eigenvalue, j + 1));
696 bool Foam::DMDModels::STDMD::dynamics()
698 SMatrix Atilde(reducedKoopmanOperator());
701 if(!eigendecomposition(Atilde))
718 bool Foam::DMDModels::STDMD::modes()
722 bool processed =
false;
723 processed = processed || modes<scalar>();
724 processed = processed || modes<vector>();
725 processed = processed || modes<sphericalTensor>();
726 processed = processed || modes<symmTensor>();
727 processed = processed || modes<tensor>();
740 Info<<
tab <<
"Writing objects of dynamics" <<
endl;
744 autoPtr<OFstream> osPtr =
747 fileName +
"_" + fieldName_,
748 mesh_.time().timeOutputValue()
750 OFstream&
os = osPtr.ref();
754 for (
const auto& i : labelRange(0, freqs_.size()))
756 os << freqs_[i] <<
tab 758 << amps_[i].real() <<
tab 759 << amps_[i].imag() <<
tab 760 << evals_[i].real() <<
tab 761 << evals_[i].imag() <<
endl;
765 Info<<
tab <<
"Ends output processing for field: " << fieldName_
770 void Foam::DMDModels::STDMD::writeFileHeader(Ostream&
os)
const 773 writeCommented(
os,
"Frequency");
774 writeTabbed(
os,
"Magnitude");
775 writeTabbed(
os,
"Amplitude (real)");
776 writeTabbed(
os,
"Amplitude (imag)");
777 writeTabbed(
os,
"Eigenvalue (real)");
778 writeTabbed(
os,
"Eigenvalue (imag)");
783 void Foam::DMDModels::STDMD::filter()
785 Info<<
tab <<
"Filtering objects of dynamics" <<
endl;
788 filterIndexed(evals_, magsi_);
789 filterIndexed(evecs_, magsi_);
790 filterIndexed(freqs_, magsi_);
791 filterIndexed(amps_, magsi_);
792 filterIndexed(mags_, magsi_);
795 if (freqs_.size() > nModes_)
797 evals_.resize(nModes_);
798 evecs_.resize(evecs_.m(), nModes_);
799 freqs_.resize(nModes_);
800 amps_.resize(nModes_);
801 mags_.resize(nModes_);
818 modeSorterTypeNames.getOrDefault
822 modeSorterType::FIRST_SNAPSHOT
857 nAgglomerationProcs_(20),
866 writeToFile_ =
dict.getOrDefault<
bool>(
"writeToFile",
true);
869 modeSorterTypeNames.getOrDefault
873 modeSorterType::FIRST_SNAPSHOT
883 dict.getCheckOrDefault<scalar>
888 *mesh_.time().deltaT().value()
900 dict.getCheckOrDefault<label>
908 dict.getCheckOrDefault<label>
918 dict.getCheckOrDefault<label>
925 nAgglomerationProcs_ =
926 dict.getCheckOrDefault<label>
928 "nAgglomerationProcs",
933 Info<<
tab <<
"Settings are read for:" <<
nl 934 <<
tab <<
" field: " << fieldName_ <<
nl 935 <<
tab <<
" modeSorter: " << modeSorterTypeNames[modeSorter_] <<
nl 936 <<
tab <<
" nModes: " << nModes_ <<
nl 937 <<
tab <<
" maxRank: " << maxRank_ <<
nl 938 <<
tab <<
" nGramSchmidt: " << nGramSchmidt_ <<
nl 939 <<
tab <<
" fMin: " << fMin_ <<
nl 940 <<
tab <<
" fMax: " << fMax_ <<
nl 941 <<
tab <<
" minBasis: " << minBasis_ <<
nl 942 <<
tab <<
" minEVal: " << minEval_ <<
nl 943 <<
tab <<
" sortLimiter: " << sortLimiter_ <<
nl 944 <<
tab <<
" nAgglomerationProcs: " << nAgglomerationProcs_ <<
nl 953 const scalar norm = L2norm(z);
960 const label nSnap = z.m()/2;
962 mesh_.time().timeName(mesh_.time().startTime().value());
973 "snapshot0_" + name_ +
"_" + fieldName_,
992 mesh_.time().writeFormat(),
993 mesh_.time().writeCompression()
996 fileHandler().writeObject(snapshot0, streamOpt,
true);
1001 G_(0,0) =
sqr(norm);
1017 const RMatrix ez(orthonormalise(z));
1020 const scalar ezNorm = L2norm(ez);
1021 const scalar zNorm = L2norm(z);
1023 if (ezNorm/zNorm > minBasis_)
1033 if (Q_.n() >= maxRank_)
1047 const label nSnap = Q_.m()/2;
1050 Qupper_ = RMatrix(Q_.subMatrix(0, 0,
max(nSnap, 1)));
1051 Qlower_ = RMatrix(Q_.subMatrix(nSnap, 0,
max(nSnap, 1)));
1063 writeToFile(
word(
"dynamics"));
1067 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.
static bool & parRun() noexcept
Test if this a parallel run.
autoPtr< fileOperation > fileHandler(std::nullptr_t)
Delete current file 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)
Number of this process (starting from masterNo() = 0)
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 processes in communicator.
#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) is 1 for serial run.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for 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
Process index of the master (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.
label m() const noexcept
The number of rows.
dimensionedScalar sin(const dimensionedScalar &ds)
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
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)
Am I the master rank.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
"nonBlocking" : (MPI_Isend, MPI_Irecv)
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...
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)