Skip to content

Commit fb9f783

Browse files
giovastabilegiovastabile
authored
Merge pull request ITHACA-FV#580 from Nkana-valentin/ROMs4MovingDomainsUpdated
Updating all dependencies files for the tutorial of the flow around moving cylinder
2 parents d3083b0 + eceaafa commit fb9f783

File tree

171 files changed

+1136495
-80
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

171 files changed

+1136495
-80
lines changed

‎src/ITHACA_CORE/Containers/Modes.C‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -491,3 +491,4 @@ void Modes<Type, PatchField, GeoMesh>::operator=(const
491491
template class Modes<scalar, fvPatchField, volMesh>;
492492
template class Modes<vector, fvPatchField, volMesh>;
493493
template class Modes<scalar, fvsPatchField, surfaceMesh>;
494+
template class Modes<vector, pointPatchField, pointMesh>;

‎src/ITHACA_CORE/Containers/Modes.H‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,7 @@ class Modes : public PtrList<GeometricField<Type, PatchField, GeoMesh >>
255255
typedef Modes<scalar, fvPatchField, volMesh> volScalarModes;
256256
typedef Modes<vector, fvPatchField, volMesh> volVectorModes;
257257
typedef Modes<scalar, fvsPatchField, surfaceMesh> surfaceScalarModes;
258+
typedef Modes<vector, pointPatchField, pointMesh> pointVectorModes;
258259

259260
#endif
260261

‎src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.C‎

Lines changed: 168 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -220,16 +220,34 @@ List<Eigen::VectorXd> Foam2Eigen::field2EigenBC(
220220
label size = field.boundaryField().size();
221221
Out.resize(size);
222222

223-
for (label i = 0; i < size; i++)
224-
{
225-
label sizei = field.boundaryField()[i].size();
226-
Out[i].resize(sizei * 3);
223+
constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value || std::is_same<surfaceMesh, GeoMesh>::value;
224+
if constexpr(check_vol){
225+
for (label i = 0; i < size; i++ )
226+
{
227+
label sizei = field.boundaryField()[i].size();
228+
Out[i].resize(sizei * 3);
227229

228-
for (label k = 0; k < sizei; k++)
230+
for (label k = 0; k < sizei ; k++)
231+
{
232+
Out[i](k) = field.boundaryField()[i][k][0];
233+
Out[i](k + sizei) = field.boundaryField()[i][k][1];
234+
Out[i](k + 2 * sizei) = field.boundaryField()[i][k][2];
235+
}
236+
}
237+
}
238+
239+
else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
240+
{
241+
for (label i = 0; i < size; i++ )
229242
{
230-
for (label j = 0; j < 3; j++)
243+
label sizei = field.boundaryField()[i].size();
244+
Out[i].resize(sizei * 3);
245+
246+
for (label k = 0; k < sizei ; k++)
231247
{
232-
Out[i](k + j * sizei) = field.boundaryField()[i][k][j];
248+
Out[i](k) = field.boundaryField()[i].patchInternalField()()[k][0];
249+
Out[i](k + sizei) = field.boundaryField()[i].patchInternalField()()[k][1];
250+
Out[i](k + 2 * sizei) = field.boundaryField()[i].patchInternalField()()[k][2];
233251
}
234252
}
235253
}
@@ -355,6 +373,8 @@ List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
355373
template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
356374
PtrList<volVectorField>& fields, label Nfields);
357375

376+
template List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
377+
PtrList<pointVectorField>& fields, label Nfields);
358378

359379
template <template <class> class PatchField, class GeoMesh>
360380
List<Eigen::MatrixXd> Foam2Eigen::PtrList2EigenBC(
@@ -509,6 +529,8 @@ GeometricField<vector, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
509529

510530
template volVectorField Foam2Eigen::Eigen2field(
511531
volVectorField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC);
532+
template pointVectorField Foam2Eigen::Eigen2field(
533+
pointVectorField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC);
512534

513535
template <template <class> class PatchField, class GeoMesh>
514536
GeometricField<scalar, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
@@ -679,16 +701,20 @@ template Eigen::MatrixXd Foam2Eigen::PtrList2Eigen(PtrList<surfaceScalarField>&
679701
fields,
680702
label Nfields);
681703
template Eigen::MatrixXd
682-
Foam2Eigen::PtrList2Eigen<vector, fvPatchField, volMesh>
683-
(PtrList<volVectorField>&
684-
fields,
685-
label Nfields);
704+
Foam2Eigen::PtrList2Eigen<vector, fvPatchField, volMesh>(PtrList<volVectorField>&
705+
fields,
706+
label Nfields);
707+
template Eigen::MatrixXd
708+
Foam2Eigen::PtrList2Eigen<vector, pointPatchField, pointMesh>(PtrList<pointVectorField>&
709+
fields, label Nfields);
710+
686711
template Eigen::MatrixXd
687712
Foam2Eigen::PtrList2Eigen<tensor, fvPatchField, volMesh>
688713
(PtrList<volTensorField>&
689714
fields,
690715
label Nfields);
691716

717+
692718
template <>
693719
void Foam2Eigen::fvMatrix2Eigen(fvMatrix<scalar> foam_matrix,
694720
Eigen::MatrixXd& A,
@@ -774,6 +800,137 @@ void Foam2Eigen::fvMatrix2Eigen(fvMatrix<vector> foam_matrix,
774800
}
775801
}
776802

803+
804+
template <typename SparseMatType, typename VecType>
805+
void Foam2Eigen::fvMat2Eigen(fvMatrix<scalar> foam_matrix,
806+
SparseMatType& A,
807+
VecType& b)
808+
{
809+
//using Trip = Eigen::Triplet<typename SparseMatType::Scalar>;
810+
label sizeA = foam_matrix.diag().size();
811+
label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
812+
foam_matrix.lower().size();
813+
A.resize(sizeA, sizeA);
814+
b.resize(sizeA);
815+
A.reserve(nel);
816+
typedef Eigen::Triplet<double> Trip;
817+
std::vector<Trip> tripletList;
818+
tripletList.reserve(nel);
819+
820+
for (label i = 0; i < sizeA; ++i)
821+
{
822+
tripletList.emplace_back(i, i, foam_matrix.diag()[i]);
823+
b(i) = foam_matrix.source()[i];
824+
}
825+
826+
const lduAddressing& addr = foam_matrix.lduAddr();
827+
const labelList& lowerAddr = addr.lowerAddr();
828+
const labelList& upperAddr = addr.upperAddr();
829+
forAll(lowerAddr, i)
830+
{
831+
tripletList.emplace_back(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]);
832+
tripletList.emplace_back(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]);
833+
}
834+
835+
forAll(foam_matrix.psi().boundaryField(), I)
836+
{
837+
const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
838+
forAll(ptch, J)
839+
{
840+
label w = ptch.faceCells()[J];
841+
tripletList.emplace_back(w, w, foam_matrix.internalCoeffs()[I][J]);
842+
b(w) += foam_matrix.boundaryCoeffs()[I][J];
843+
}
844+
}
845+
846+
A.setFromTriplets(tripletList.begin(), tripletList.end());
847+
}
848+
849+
// Explicit instantiations
850+
template void Foam2Eigen::fvMat2Eigen<Eigen::SparseMatrix<double, Eigen::RowMajor>, Eigen::VectorXd>(
851+
fvMatrix<scalar> foam_matrix,
852+
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
853+
Eigen::VectorXd& b);
854+
template void Foam2Eigen::fvMat2Eigen<Eigen::SparseMatrix<double, Eigen::ColMajor>, Eigen::VectorXd>(
855+
fvMatrix<scalar> foam_matrix,
856+
Eigen::SparseMatrix<double, Eigen::ColMajor>& A,
857+
Eigen::VectorXd& b);
858+
859+
860+
template<typename SparseMatType, typename VecType>
861+
void Foam2Eigen::fvMat2Eigen(fvMatrix<vector> foam_matrix,
862+
SparseMatType& A,
863+
VecType& b)
864+
{
865+
label sizeA = foam_matrix.diag().size();
866+
label nel = foam_matrix.diag().size() + foam_matrix.upper().size() +
867+
foam_matrix.lower().size();
868+
A.resize(sizeA * 3, sizeA * 3);
869+
A.reserve(nel * 3);
870+
b.resize(sizeA * 3);
871+
typedef Eigen::Triplet<double> Trip;
872+
std::vector<Trip> tripletList;
873+
tripletList.reserve(nel * 3);
874+
875+
for (auto i = 0; i < sizeA; i++)
876+
{
877+
tripletList.push_back(Trip(i, i, foam_matrix.diag()[i]));
878+
tripletList.push_back(Trip(sizeA + i, sizeA + i, foam_matrix.diag()[i]));
879+
tripletList.push_back(Trip(2 * sizeA + i, 2 * sizeA + i,
880+
foam_matrix.diag()[i]));
881+
b(i) = foam_matrix.source()[i][0];
882+
b(sizeA + i) = foam_matrix.source()[i][1];
883+
b(2 * sizeA + i) = foam_matrix.source()[i][2];
884+
}
885+
886+
const lduAddressing& addr = foam_matrix.lduAddr();
887+
const labelList& lowerAddr = addr.lowerAddr();
888+
const labelList& upperAddr = addr.upperAddr();
889+
forAll(lowerAddr, i)
890+
{
891+
tripletList.push_back(Trip(lowerAddr[i], upperAddr[i], foam_matrix.upper()[i]));
892+
tripletList.push_back(Trip(lowerAddr[i] + sizeA, upperAddr[i] + sizeA,
893+
foam_matrix.upper()[i]));
894+
tripletList.push_back(Trip(lowerAddr[i] + sizeA * 2, upperAddr[i] + sizeA * 2,
895+
foam_matrix.upper()[i]));
896+
tripletList.push_back(Trip(upperAddr[i], lowerAddr[i], foam_matrix.lower()[i]));
897+
tripletList.push_back(Trip(upperAddr[i] + sizeA, lowerAddr[i] + sizeA,
898+
foam_matrix.lower()[i]));
899+
tripletList.push_back(Trip(upperAddr[i] + sizeA * 2, lowerAddr[i] + sizeA * 2,
900+
foam_matrix.lower()[i]));
901+
}
902+
903+
forAll(foam_matrix.psi().boundaryField(), I)
904+
{
905+
const fvPatch& ptch = foam_matrix.psi().boundaryField()[I].patch();
906+
forAll(ptch, J)
907+
{
908+
label w = ptch.faceCells()[J];
909+
tripletList.push_back(Trip(w, w, foam_matrix.internalCoeffs()[I][J][0]));
910+
tripletList.push_back(Trip(w + sizeA, w + sizeA,
911+
foam_matrix.internalCoeffs()[I][J][1]));
912+
tripletList.push_back(Trip(w + sizeA * 2, w + sizeA * 2,
913+
foam_matrix.internalCoeffs()[I][J][2]));
914+
b(w) += foam_matrix.boundaryCoeffs()[I][J][0];
915+
b(w + sizeA) += foam_matrix.boundaryCoeffs()[I][J][1];
916+
b(w + sizeA * 2) += foam_matrix.boundaryCoeffs()[I][J][2];
917+
}
918+
}
919+
920+
A.setFromTriplets(tripletList.begin(), tripletList.end());
921+
}
922+
923+
template void Foam2Eigen::fvMat2Eigen<Eigen::SparseMatrix<double, Eigen::RowMajor>, Eigen::VectorXd>(
924+
fvMatrix<vector> foam_matrix,
925+
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
926+
Eigen::VectorXd& b);
927+
928+
template void Foam2Eigen::fvMat2Eigen<Eigen::SparseMatrix<double, Eigen::ColMajor>, Eigen::VectorXd>(
929+
fvMatrix<vector> foam_matrix,
930+
Eigen::SparseMatrix<double, Eigen::ColMajor>& A,
931+
Eigen::VectorXd& b);
932+
933+
/////////////////////////////////////////////////////////////////////////////////////////////
777934
template <>
778935
void Foam2Eigen::fvMatrix2Eigen(fvMatrix<scalar> foam_matrix,
779936
Eigen::SparseMatrix<double>& A, Eigen::VectorXd& b)

‎src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.H‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,16 @@
3939

4040
#include "fvCFD.H"
4141
#include "IOmanip.H"
42+
#include "pointMesh.H"
43+
#include "pointPatchField.H"
4244
#include "ITHACAassert.H"
4345
#include "ITHACAutilities.H"
4446
#include <tuple>
4547
#include <sys/stat.h>
4648
#pragma GCC diagnostic push
4749
#pragma GCC diagnostic ignored "-Wold-style-cast"
4850
#include <Eigen/Eigen>
51+
#include <type_traits>
4952
#pragma GCC diagnostic pop
5053
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
5154

@@ -86,6 +89,20 @@ class Foam2Eigen
8689
template <class type_foam_matrix, class type_A, class type_B>
8790
static void fvMatrix2Eigen(fvMatrix<type_foam_matrix> foam_matrix, type_A& A,
8891
type_B& b);
92+
/*
93+
template <typename ScalarType, typename SparseMatType, typename VecType>
94+
static void fvMat2Eigen(fvMatrix<ScalarType> foam_matrix,
95+
SparseMatType& A,
96+
VecType& b);*/
97+
template<typename SparseMatType, typename VecType>
98+
static void fvMat2Eigen(Foam::fvMatrix<scalar> foam_matrix,
99+
SparseMatType& A,
100+
VecType& b);
101+
102+
template<typename SparseMatType, typename VecType>
103+
static void fvMat2Eigen(Foam::fvMatrix<vector> foam_matrix,
104+
SparseMatType& A,
105+
VecType& b);
89106

90107
//----------------------------------------------------------------------
91108
/// @brief Convert a ldu OpenFOAM matrix into a Eigen Matrix A

‎src/ITHACA_CORE/ITHACAPOD/ITHACAPOD.C‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,7 @@ void getModes(
9494
bool correctBC)
9595
{
9696
ITHACAparameters* para(ITHACAparameters::getInstance());
97+
constexpr bool check_vol = std::is_same<volMesh, GeoMesh>::value || std::is_same<surfaceMesh, GeoMesh>::value;
9798
word PODkey = "POD_" + fieldName;
9899
word PODnorm = para->ITHACAdict->lookupOrDefault<word>(PODkey, "L2");
99100
M_Assert(PODnorm == "L2" ||
@@ -197,6 +198,15 @@ void getModes(
197198
{
198199
normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
199200
modesEig.col(i))(0, 0));
201+
if constexpr(check_vol){
202+
normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * V.asDiagonal() *
203+
modesEig.col(i))(0, 0));
204+
}
205+
206+
else if constexpr(std::is_same<pointMesh, GeoMesh>::value){
207+
normFact(i, 0) = std::sqrt((modesEig.col(i).transpose() * modesEig.col(i))(0, 0));
208+
}
209+
200210
if (Pstream::parRun())
201211
{
202212
normFact(i, 0) = (modesEig.col(i).transpose() * V.asDiagonal() *
@@ -312,6 +322,10 @@ template void getModes(
312322
PtrList<surfaceScalarField>& snapshots, PtrList<surfaceScalarField>& modes,
313323
word fieldName, bool podex, bool supex, bool sup, label nmodes,
314324
bool correctBC);
325+
template void getModes(
326+
PtrList<pointVectorField>& snapshots, PtrList<pointVectorField>& modes,
327+
word fieldName, bool podex, bool supex, bool sup, label nmodes,
328+
bool correctBC);
315329

316330
template<class Type, template<class> class PatchField, class GeoMesh>
317331
void getModesMemoryEfficient(

0 commit comments

Comments
 (0)