37 namespace functionObjects
47 Foam::functionObjects::STDMD::modeSorterType
49 Foam::functionObjects::STDMD::modeSorterTypeNames
51 { modeSorterType::KIEWAT ,
"kiewat" },
52 { modeSorterType::KOU_ZHANG ,
"kouZhang" },
53 { modeSorterType::FIRST_SNAPSHOT,
"firstSnap" }
59 Foam::scalar Foam::functionObjects::STDMD::parnorm
61 const RMatrix& colVector
66 if (colVector.n() != 1)
69 <<
" # Input matrix is not a column vector. #"
73 const bool noSqrt =
true;
74 scalar result(colVector.columnNorm(0, noSqrt));
75 reduce(result, sumOp<scalar>());
82 void Foam::functionObjects::STDMD::snapshot()
84 bool processed =
false;
85 processed = processed || getSnapshotType<scalar>();
86 processed = processed || getSnapshotType<vector>();
87 processed = processed || getSnapshotType<sphericalTensor>();
88 processed = processed || getSnapshotType<symmTensor>();
89 processed = processed || getSnapshotType<tensor>();
94 <<
" # Unknown type of input field during snapshot loading = "
95 << fieldName_ <<
" #" <<
nl
96 <<
" # Do you execute required functionObjects "
97 <<
"before executing STDMD, e.g. mapFields?"
103 void Foam::functionObjects::STDMD::init()
105 bool processed =
false;
106 processed = processed || getComps<scalar>();
107 processed = processed || getComps<vector>();
108 processed = processed || getComps<sphericalTensor>();
109 processed = processed || getComps<symmTensor>();
110 processed = processed || getComps<tensor>();
115 <<
" # Unknown type of input field during initialisation = "
116 << fieldName_ <<
" #" <<
nl
117 <<
" # Do you execute required functionObjects "
118 <<
"before executing STDMD, e.g. mapFields?"
122 nSnap_ = nComps_*mesh_.nCells();
127 <<
" # Zero-size input field = " << fieldName_ <<
" #"
135 z_ = RMatrix(2*nSnap_, 1,
Zero);
137 X1_ = RMatrix(nSnap_, 1,
Zero);
155 void Foam::functionObjects::STDMD::initBasis()
157 std::copy(z_.cbegin(), z_.cbegin() + nSnap_, X1_.begin());
159 zNorm_ = parnorm(z_);
161 Gz_(0,0) =
sqr(zNorm_);
165 void Foam::functionObjects::STDMD::GramSchmidt()
168 RMatrix dz(Qz_.n(), 1,
Zero);
170 for (label i = 0; i < nGramSchmidt_; ++i)
173 reduce(dz, sumOp<RMatrix>());
179 void Foam::functionObjects::STDMD::expandBasis()
182 <<
" Expanding orthonormal basis for field = " << fieldName_ <<
" #"
186 Qz_.resize(Qz_.m(), Qz_.n() + 1);
187 Qz_.subColumn(Qz_.n() - 1) = ez_/ezNorm_;
190 Gz_.resize(Gz_.m() + 1);
194 void Foam::functionObjects::STDMD::updateGz()
196 RMatrix zTilde(Qz_ & z_);
197 reduce(zTilde, sumOp<RMatrix>());
199 const SMatrix zTildes(zTilde^zTilde);
205 void Foam::functionObjects::STDMD::compressBasis()
208 <<
" Compressing orthonormal basis for field = " << fieldName_ <<
" #"
215 const bool symmetric =
true;
216 const EigenMatrix<scalar> EM(Gz_, symmetric);
217 const SquareMatrix<scalar>& EVecs = EM.EVecs();
218 DiagonalMatrix<scalar> EVals(EM.EValsRe());
222 testEigenvalues(Gz_, EVals);
223 testEigenvectors(Gz_, EVals, EVecs);
227 const auto descend = [&](scalar a, scalar
b){
return a >
b; };
228 const List<label> permut(EVals.sortPermutation(descend));
229 EVals.applyPermutation(permut);
230 EVals.resize(EVals.size() - 1);
233 Gz_ = SMatrix(maxRank_,
Zero);
236 qz.resize(Qz_.n(), maxRank_);
237 for (label i = 0; i < maxRank_; ++i)
239 qz.subColumn(i) = EVecs.subColumn(permut[i]);
250 void Foam::functionObjects::STDMD::calcAp()
257 QRMatrix<RMatrix> QRM
260 QRMatrix<RMatrix>::outputTypes::REDUCED_R,
261 QRMatrix<RMatrix>::colPivoting::FALSE
275 Log<<
tab <<
"# " <<
name() <<
": Gathering all local Rx #" <<
endl;
288 mRowSum += RxList[i].m();
291 Rx.resize(mRowSum,
Rx.n());
294 <<
": Populating the global Rx #" <<
endl;
298 label mRows = RxList[i].m();
300 Rx.subMatrix(m, 0, mRows) = RxList[i];
306 <<
": Computing the parallel-direct tall-skinny QR decomp. #"
308 QRMatrix<RMatrix> QRM
311 QRMatrix<RMatrix>::outputTypes::REDUCED_R,
312 QRMatrix<RMatrix>::storeMethods::IN_PLACE,
313 QRMatrix<RMatrix>::colPivoting::FALSE
318 <<
": Computing Moore-Penrose pseudo-inverse of Rx #"
323 const RMatrix Gx(
Rx*(Gz_^
Rx));
326 <<
": Computing Moore-Penrose pseudo-inverse of Gx #"
328 const RMatrix GxInv(
pinv(Gx));
337 SMatrix A2(QzUH_ & QzUH_);
338 reduce(A2, sumOp<SMatrix>());
341 SMatrix A3(QzUH_ & QzLH_);
342 reduce(A3, sumOp<SMatrix>());
347 Ap_ = SMatrix(RxInv_ & ((A3*(Gz_*A2))*A1));
353 <<
": Computing Moore-Penrose pseudo-inverse of Rx #"
358 const RMatrix Gx(
Rx*(Gz_^
Rx));
361 <<
": Computing Moore-Penrose pseudo-inverse of Gx #"
363 const RMatrix GxInv(
pinv(Gx));
370 SMatrix A2(QzUH_ & QzUH_);
373 SMatrix A3(QzUH_ & QzLH_);
378 Ap_ = SMatrix(RxInv_ & ((A3*(Gz_*A2))*A1));
383 void Foam::functionObjects::STDMD::calcEigen()
385 Log<<
tab <<
"# " <<
name() <<
": Computing eigendecomposition #" <<
endl;
391 const EigenMatrix<scalar> EM(Ap_);
392 const DiagonalMatrix<scalar>& evalsRe = EM.EValsRe();
393 const DiagonalMatrix<scalar>& evalsIm = EM.EValsIm();
395 oEVals_.resize(evalsRe.size());
396 oEVecs_ = RectangularMatrix<complex>(EM.complexEVecs());
400 for (
auto& eval : oEVals_)
402 eval =
complex(evalsRe[i], evalsIm[i]);
408 testEigenvalues(Ap_, evalsRe);
409 testEigenvectors(Ap_, oEVals_, oEVecs_);
415 bool Foam::functionObjects::STDMD::close
423 if (s1 == s2)
return true;
429 void Foam::functionObjects::STDMD::testEigenvalues
431 const SquareMatrix<scalar>&
A,
432 const DiagonalMatrix<scalar>& EValsRe
435 const scalar trace =
A.trace();
437 const scalar EValsSum =
sum(EValsRe);
439 const bool clse = close(EValsSum - trace, 0, absTol_, relTol_);
444 <<
" ## PASS: Eigenvalues ##" <<
endl;
449 <<
" ## INCONCLUSIVE: Eigenvalues ##" <<
nl
450 <<
" # sum(eigenvalues) = " << EValsSum <<
nl
451 <<
" # trace(A) = " << trace <<
nl
452 <<
" # sum(eigenvalues) - trace(A) ~ 0 ?= "
453 << EValsSum - trace <<
nl
454 <<
" #######################"
460 <<
" ## Operands ##" <<
nl
461 <<
" # eigenvalues:" <<
nl << EValsRe <<
nl
462 <<
" # input matrix A:" <<
nl <<
A <<
nl
470 void Foam::functionObjects::STDMD::testEigenvectors
472 const SquareMatrix<scalar>&
A,
473 const DiagonalMatrix<scalar>& EValsRe,
474 const SquareMatrix<scalar>& EVecs
477 unsigned nInconclusive = 0;
479 for (label i = 0; i <
A.n(); ++i)
481 const RectangularMatrix<scalar> EVec(EVecs.subColumn(i));
482 const scalar EVal = EValsRe[i];
483 const RectangularMatrix<scalar> leftSide(
A*EVec);
484 const RectangularMatrix<scalar> rightSide(EVal*EVec);
485 const scalar rslt = (leftSide - rightSide).norm();
487 const bool clse = close(rslt, 0, absTol_, relTol_);
492 <<
" ## INCONCLUSIVE: Eigenvector ##" <<
nl
493 <<
" # (A & EVec - EVal*EVec).norm() ~ 0 ?= " << rslt <<
nl
494 <<
" ##################################"
500 <<
" ## Operands ##" <<
nl
501 <<
" # eigenvalue:" <<
nl << EVal <<
nl
502 <<
" # input matrix A:" <<
nl <<
A <<
nl
503 <<
" # eigenvector:" <<
nl << EVec <<
nl
504 <<
" # (A & EVec):" <<
nl << leftSide <<
nl
505 <<
" # (EVal*EVec):" <<
nl << rightSide <<
nl
517 <<
" ## PASS: Eigenvectors ##" <<
endl;
522 void Foam::functionObjects::STDMD::testEigenvectors
524 const SquareMatrix<scalar>&
A,
525 const List<complex>& EVals,
526 const RectangularMatrix<complex>& EVecs
529 SquareMatrix<complex>
B(
A.m());
530 auto convertToComplex = [&](
const scalar& val) {
return complex(val); };
539 unsigned nInconclusive = 0;
541 for (label i = 0; i <
B.n(); ++i)
543 const RectangularMatrix<complex> EVec(EVecs.subColumn(i));
545 const RectangularMatrix<complex> leftSide(
B*EVec);
546 const RectangularMatrix<complex> rightSide(EVal*EVec);
547 const scalar rslt = (leftSide - rightSide).norm();
549 const bool clse = close(rslt, 0, absTol_, relTol_);
554 <<
" ## INCONCLUSIVE: Eigenvector ##" <<
nl
555 <<
" # (A & EVec - EVal*EVec).norm() ~ 0 ?= " << rslt <<
nl
556 <<
" ##################################"
562 <<
" ## Operands ##" <<
nl
563 <<
" # eigenvalue:" <<
nl << EVal <<
nl
564 <<
" # input matrix A:" <<
nl <<
A <<
nl
565 <<
" # eigenvector:" <<
nl << EVec <<
nl
566 <<
" # (A & EVec):" <<
nl << leftSide <<
nl
567 <<
" # (EVal*EVec):" <<
nl << rightSide <<
nl
579 <<
" ## PASS: Eigenvectors ##" <<
endl;
584 void Foam::functionObjects::STDMD::filterEVals()
586 Log<<
tab <<
"# " <<
name() <<
": Filtering eigenvalues #" <<
endl;
590 List<complex> cpEVals(oEVals_.size());
601 if (cpEVals.size() == 0)
604 <<
"No eigenvalue with mag(eigenvalue) larger than "
605 <<
"minMagEVal_ = " << minMagEVal_ <<
" was found."
618 void Foam::functionObjects::STDMD::calcFreqs()
620 Log<<
tab <<
"# " <<
name() <<
": Computing frequencies #" <<
endl;
624 oFreqs_.resize(oEVals_.size());
634 std::transform(oEVals_.cbegin(), oEVals_.cend(), oFreqs_.begin(), fEq);
639 void Foam::functionObjects::STDMD::calcFreqI()
641 Log<<
tab <<
"# " <<
name() <<
": Computing frequency indices #" <<
endl;
646 auto margin = [&](
const scalar&
x){
return (fMin_ <=
x &&
x < fMax_); };
648 auto it = std::find_if(oFreqs_.cbegin(), oFreqs_.cend(), margin);
650 while (it != oFreqs_.end())
653 it = std::find_if(std::next(it), oFreqs_.cend(), margin);
661 void Foam::functionObjects::STDMD::calcAmps()
663 Log<<
tab <<
"# " <<
name() <<
": Computing amplitudes #" <<
endl;
665 RMatrix temp((RxInv_.T()^QzUH_)*X1_);
666 reduce(temp, sumOp<RMatrix>());
670 oAmps_.resize(temp.m());
671 const RCMatrix pinvEVecs(
pinv(oEVecs_));
674 for (label i = 0; i < oAmps_.size(); ++i)
676 for (label
k = 0;
k < temp.m(); ++
k)
678 oAmps_[i] += pinvEVecs(i,
k)*temp(
k, 0);
686 void Foam::functionObjects::STDMD::calcMags()
688 Log<<
tab <<
"# " <<
name() <<
": Computing magnitudes #" <<
endl;
692 oMags_.resize(oAmps_.size());
694 Log<<
tab <<
"# " <<
name() <<
": Sorting modes with ";
698 case modeSorterType::FIRST_SNAPSHOT:
700 Log<<
"method of first snapshot #" <<
endl;
712 case modeSorterType::KIEWAT:
714 Log<<
"modified weighted amplitude scaling method #" <<
endl;
718 const scalar modeNorm = 1;
720 List<scalar> w(currIndex_);
721 std::iota(w.begin(), w.end(), 1);
722 w =
sin(
twoPi/currIndex_*(w - 1 - 0.25*currIndex_))*pr + pr;
725 for (
const auto& amp : oAmps_)
727 oMags_[i] = sorter(w, amp, oEVals_[i], modeNorm);
733 case modeSorterType::KOU_ZHANG:
735 Log<<
"weighted amplitude scaling method #" <<
endl;
737 const scalar modeNorm = 1;
738 const List<scalar> w(currIndex_, 1.0);
740 for (
const auto& amp : oAmps_)
742 oMags_[i] = sorter(w, amp, oEVals_[i], modeNorm);
755 void Foam::functionObjects::STDMD::calcMagI()
757 Log<<
tab <<
"# " <<
name() <<
": Computing magnitude indices #" <<
endl;
764 [&](
const label i1,
const label i2)
766 return !(oMags_[i1] < oMags_[i2]);
769 std::sort(iMags_.begin(), iMags_.end(), descend);
776 void Foam::functionObjects::STDMD::calcModes()
781 if (nModes_ < iMags_.size())
783 iMags_.resize(nModes_);
788 const RMatrix temp(QzUH_*RxInv_);
790 for (
const auto& i : iMags_)
795 for (label
p = 0;
p <
mode.size(); ++
p)
797 for (label q = 0; q < oEVecs_.m(); ++q)
799 mode[
p] += temp(
p, q)*oEVecs_(q, i);
807 "modeReal" + std::to_string(oModeI) + fieldName_ +
"_" +
name(),
821 "modeImag" + std::to_string(oModeI) + fieldName_ +
"_" +
name(),
832 for (label
k = 0;
k < nComps_; ++
k)
834 for (label j = 0; j < modeReal.size(); ++j)
837 modeReal[j].replace(
k, m.real());
838 modeImag[j].replace(
k, m.imag());
851 Foam::scalar Foam::functionObjects::STDMD::sorter
853 const List<scalar>& weight,
856 const scalar modeNorm
860 if (!(
mag(eval) < GREAT &&
mag(eval) > VSMALL))
862 Info<<
" Returning zero magnitude for mag(eval) = " <<
mag(eval)
869 if (
mag(eval)*currIndex_ > sortLimiter_)
871 Info<<
" Returning zero magnitude for"
872 <<
" mag(eval) = " <<
mag(eval)
873 <<
" currIndex = " << currIndex_
874 <<
" sortLimiter = " << sortLimiter_
880 scalar magnitude = 0;
882 for (label j = 0; j < currIndex_; ++j)
884 magnitude += weight[j]*modeNorm*
mag(amplitude*
pow(eval, j + 1));
891 void Foam::functionObjects::STDMD::writeFileHeader(Ostream& os)
const
894 writeCommented(os,
"Frequency");
895 writeTabbed(os,
"Magnitude");
896 writeTabbed(os,
"Amplitude (real)");
897 writeTabbed(os,
"Amplitude (imag)");
898 writeTabbed(os,
"Eigenvalue (real)");
899 writeTabbed(os,
"Eigenvalue (imag)");
904 void Foam::functionObjects::STDMD::filterOutput()
906 Log<<
tab <<
"# " <<
name() <<
": Filtering text output #" <<
endl;
911 filterIndexed(oEVals_, iMags_);
912 filterIndexed(oEVecs_, iMags_);
913 filterIndexed(oFreqs_, iMags_);
914 filterIndexed(oAmps_, iMags_);
915 filterIndexed(oMags_, iMags_);
918 if (nModes_ < oFreqs_.size())
920 oEVals_.resize(nModes_);
921 oEVecs_.resize(oEVecs_.m(), nModes_);
922 oFreqs_.resize(nModes_);
923 oAmps_.resize(nModes_);
924 oMags_.resize(nModes_);
930 void Foam::functionObjects::STDMD::writeOutput(OFstream& os)
const
936 for (
const auto& i : labelRange(0, oFreqs_.size()))
938 os << oFreqs_[i] <<
tab
940 << oAmps_[i].real() <<
tab
941 << oAmps_[i].imag() <<
tab
942 << oEVals_[i].real() <<
tab
943 << oEVals_[i].imag() <<
endl;
948 void Foam::functionObjects::STDMD::calcOutput()
951 <<
" Starts output processing for field = " << fieldName_ <<
" #"
955 QzUH_ = Qz_.subMatrix(0, 0, nSnap_);
956 QzLH_ = Qz_.subMatrix(nSnap_, 0, nSnap_);
985 autoPtr<OFstream> osPtr =
986 createFile(
"uSTDMD" + fieldName_, mesh_.time().timeOutputValue());
987 OFstream& os = osPtr.ref();
997 autoPtr<OFstream> osPtr =
998 createFile(
"STDMD" + fieldName_, mesh_.time().timeOutputValue());
999 OFstream& os = osPtr.ref();
1005 <<
" Ends output processing for field = " << fieldName_ <<
" #"
1023 modeSorterTypeNames.getOrDefault
1027 modeSorterType::FIRST_SNAPSHOT
1030 fieldName_(
dict.get<
word>(
"field")),
1031 initialised_(
false),
1036 dict.getCheckOrDefault<label>
1045 dict.getCheckOrDefault<label>
1054 dict.getCheckOrDefault<label>
1063 dict.getCheckOrDefault<label>
1072 dict.getCheckOrDefault<label>
1108 <<
" ### This STDMD release is the beta release ###" <<
nl
1109 <<
" Therefore small-to-medium changes in input/output interfaces "
1110 <<
"and internal structures should be expected in the next versions."
1114 if (
runTime.isAdjustTimeStep())
1117 <<
" # STDMD: Viable only for fixed time-step computations. #"
1122 if (mesh_.topoChanging())
1125 <<
" # STDMD: Available only for non-changing mesh topology. #"
1139 testEigen_ =
dict.getOrDefault<
bool>(
"testEigen",
false);
1143 dumpEigen_ =
dict.getOrDefault<
bool>(
"dumpEigen",
false);
1147 dict.getCheckOrDefault<scalar>
1155 dict.getCheckOrDefault<scalar>
1163 dict.getCheckOrDefault<scalar>
1171 dict.getCheckOrDefault<scalar>
1179 dict.getCheckOrDefault<scalar>
1184 *mesh_.time().deltaT().value()
1190 dict.getCheckOrDefault<scalar>
1197 writeToFile_ =
dict.getOrDefault<
bool>(
"writeToFile",
true);
1200 <<
tab <<
"# Input is read for field = "
1201 << fieldName_ <<
" #" <<
nl <<
endl;
1214 if (currIndex_ == 1)
1224 zNorm_ = parnorm(z_);
1225 ezNorm_ = parnorm(ez_);
1227 if (minBasis_ < ezNorm_/zNorm_)
1235 if (maxRank_ < Qz_.n())
1243 <<
" Execution index = " << currIndex_
1244 <<
" for field = " << fieldName_ <<
" #"
1258 <<
" # STDMD needs at least three snapshots to produce output #"
1260 <<
" # Only " << currIndex_ + 1 <<
" snapshots are available #"
1262 <<
" # Skipping STDMD output calculation and write #"
1272 initialised_ =
false;
1274 mesh_.time().printExecutionTime(
Info);