48 Foam::DMDModels::STDMD::modeSorterType
50Foam::DMDModels::STDMD::modeSorterTypeNames
52 { modeSorterType::FIRST_SNAPSHOT,
"firstSnapshot" },
53 { modeSorterType::KIEWAT ,
"kiewat" },
54 { modeSorterType::KOU_ZHANG ,
"kouZhang" }
60Foam::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);
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);
104 for (label i = 0; i < sz0; ++i)
107 for (label j = 0; j < sz1; ++j)
119void 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);
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);
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];
161void Foam::DMDModels::STDMD::compress()
163 Info<<
tab <<
"Compressing orthonormal basis for field: " << fieldName_
166 RMatrix q(1, 1,
Zero);
170 const bool symmetric =
true;
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]);
200reducedKoopmanOperator()
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_
269 if (recvMtrx.size() > 0)
271 Rx.resize(
Rx.
m() + recvMtrx.m(),
Rx.
n());
274 Rx.
m() - recvMtrx.m(),
275 Rx.
n() - recvMtrx.n()
291 RMatrix&
R = QRM.R();
300 if (procNoInSubset == 0)
309 pBufs.finishedSends();
319 nbr += nAgglomerationProcs_
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()
348 RMatrix&
R = QRM.R();
364 SMatrix A2(Qupper_ & Qlower_);
369 SMatrix A3(Qupper_ & Qupper_);
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));
399bool Foam::DMDModels::STDMD::eigendecomposition(SMatrix& Atilde)
406 Info<<
tab <<
"Computing eigendecomposition" <<
endl;
414 evecs_ = RCMatrix(EM.complexEVecs());
418 for (
auto& eval : evals_)
420 eval =
complex(evalsRe[i], evalsIm[i]);
433 [&](
const complex&
x){ return mag(x) > minEval_; }
435 cp.resize(std::distance(
cp.begin(), it));
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
469void 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())
502 freqsi_.append(std::distance(freqs_.cbegin(), it));
503 it = std::find_if(std::next(it), freqs_.cend(), margin);
511void Foam::DMDModels::STDMD::amplitudes()
517 "snapshot0_" + name_ +
"_" + fieldName_,
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];
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];
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];
574void 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;
595 [&](
const complex& val){ return mag(val); }
600 case modeSorterType::KIEWAT:
602 Info<<
"modified weighted amplitude scaling method" <<
endl;
606 const scalar modeNorm = 1;
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;
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);
656Foam::scalar Foam::DMDModels::STDMD::sorter
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));
696bool Foam::DMDModels::STDMD::dynamics()
698 SMatrix Atilde(reducedKoopmanOperator());
701 if(!eigendecomposition(Atilde))
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;
748 mesh_.time().timeOutputValue()
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_
770void 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)");
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
845 fieldName_(
dict.get<
word>(
"field")),
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()
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"));
#define R(A, B, C, D, E, F, K, M)
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Abstract base class for DMD models to handle DMD characteristics for the DMD function object.
Streaming Total Dynamic Mode Decomposition (i.e. STDMD) is a variant of dynamic mode decomposition.
A templated (N x N) diagonal matrix of objects of <Type>, effectively containing N elements,...
EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes a diagonalisable nonsymmet...
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
A primitive field of type <T> with automated input and output.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
The IOstreamOption is a simple container for options an IOstream can normally have.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void resize(const label len)
Adjust allocated size of list.
static direction m() noexcept
The number of rows.
static direction n() noexcept
The number of columns.
const_iterator cbegin() const noexcept
Return const_iterator to begin traversing a constant Matrix.
label m() const noexcept
The number of rows.
ConstMatrixBlock< mType > subColumn(const label colIndex, const label rowIndex=0, label len=-1) const
Return const column or column's subset of Matrix.
static MinMax< scalar > ge(const scalar &minVal)
A semi-infinite range from minVal to the type max.
Output to file stream, using an OSstream.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Buffers for inter-processor communications streams (UOPstream, UIPstream).
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void scatter(const List< commsStruct > &comms, T &value, const int tag, const label comm)
Broadcast data: Distribute without modification.
QRMatrix computes QR decomposition of a given scalar/complex matrix A into the following:
modes
Options for the decomposition mode.
virtual bool read()
Re-read model coefficients if they have changed.
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
iterator begin() noexcept
Return an iterator to begin traversing the UList.
void size(const label n)
Older name for setAddressableSize.
@ nonBlocking
"nonBlocking"
static constexpr int masterNo() noexcept
Process index of the master (always 0)
static bool & parRun() noexcept
Test if this a parallel run.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
T & ref()
Return reference to the managed object without nullptr checking.
A complex number, similar to the C++ complex type.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
A class for handling file names.
virtual bool writeObject(const regIOobject &io, IOstreamOption streamOpt=IOstreamOption(), const bool valid=true) const
Writes a regIOobject (so header, contents and divider).
virtual bool writeToFile()
Write tracks to file.
Mesh data needed to do the Finite Volume discretisation.
virtual bool update()
Update the mesh for both mesh motion and topology change.
A range or interval of labels defined by a start and a size.
A traits class, which is primarily used for primitives.
int myProcNo() const noexcept
Return processor number.
void updateG()
Update G and calculate total heat flux on boundary.
splitCell * master() const
A class for managing temporary objects.
void clear() const noexcept
void initialise()
Initialise integral-scale box properties.
A List of wordRe with additional matching capabilities.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar twoPi(2 *M_PI)
const fileOperation & fileHandler()
Get current file handler.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
tensor Rx(const scalar omega)
Rotational transformation tensor about the x-axis by omega radians.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sin(const dimensionedScalar &ds)
static void writeHeader(Ostream &os, const word &fieldName)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
constexpr char tab
The tab '\t' character(0x09)
const volScalarField & cp
#define forAll(list, i)
Loop across all elements in list.
Functor wrapper of allow/deny lists of wordRe for filtering.