46 namespace MultiRegions
51 m_lhom(1), m_homogeneous1DBlockMat(
60 const bool useFFT,
const bool dealiasing)
61 :
ExpList(type), m_useFFT(useFFT), m_lhom(lhom),
62 m_homogeneous1DBlockMat(
64 m_dealiasing(dealiasing)
67 m_comm = pSession->GetComm();
70 "Homogeneous Basis is a null basis");
75 ?
m_comm->GetColumnComm()->GetColumnComm()
79 "HomModesZ should be a multiple of npz.");
90 "HomModesZ/npz should be an even integer.");
118 ASSERTL0(
false,
"Dealiasing available just in combination "
128 :
ExpList(In, false), m_transposition(In.m_transposition),
129 m_StripZcomm(In.m_StripZcomm), m_useFFT(In.m_useFFT), m_FFT(In.m_FFT),
130 m_tmpIN(In.m_tmpIN), m_tmpOUT(In.m_tmpOUT),
131 m_homogeneousBasis(In.m_homogeneousBasis), m_lhom(In.m_lhom),
132 m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
133 m_dealiasing(In.m_dealiasing), m_padsize(In.m_padsize)
141 :
ExpList(In, eIDs, false, ImpType), m_transposition(In.m_transposition),
142 m_useFFT(In.m_useFFT), m_FFT(In.m_FFT), m_tmpIN(In.m_tmpIN),
143 m_tmpOUT(In.m_tmpOUT), m_homogeneousBasis(In.m_homogeneousBasis),
145 m_homogeneous1DBlockMat(
147 m_dealiasing(In.m_dealiasing), m_padsize(In.m_padsize)
186 int num_dofs = inarray1.size();
204 int num_points_per_plane = num_dofs /
m_planes.size();
206 if (!
m_session->DefinesSolverInfo(
"HomoStrip"))
208 num_proc =
m_comm->GetColumnComm()->GetSize();
214 int num_dfts_per_proc =
215 num_points_per_plane / num_proc + (num_points_per_plane % num_proc > 0);
233 for (
int i = 0; i < num_dfts_per_proc; i++)
237 Vmath::Vcopy(N, &(ShufV1[i * N]), 1, &(ShufV1_PAD_coef[0]), 1);
238 Vmath::Vcopy(N, &(ShufV2[i * N]), 1, &(ShufV2_PAD_coef[0]), 1);
241 m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
242 m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);
247 ShufV1V2_PAD_phys, 1);
251 m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
255 Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1, &(ShufV1V2[i * N]), 1);
285 int ndim = inarray1.size();
286 ASSERTL1(inarray2.size() % ndim == 0,
287 "Wrong dimensions for DealiasedDotProd.");
288 int nvec = inarray2.size() / ndim;
290 int num_dofs = inarray1[0].size();
293 int num_points_per_plane = num_dofs /
m_planes.size();
295 if (!
m_session->DefinesSolverInfo(
"HomoStrip"))
297 num_proc =
m_comm->GetColumnComm()->GetSize();
303 int num_dfts_per_proc =
304 num_points_per_plane / num_proc + (num_points_per_plane % num_proc > 0);
311 for (
int i = 0; i < ndim; i++)
315 for (
int i = 0; i < ndim * nvec; i++)
322 for (
int i = 0; i < ndim; i++)
327 for (
int i = 0; i < ndim * nvec; i++)
338 for (
int i = 0; i < ndim; i++)
347 for (
int i = 0; i < ndim * nvec; i++)
356 for (
int i = 0; i < nvec; i++)
362 for (
int i = 0; i < ndim; i++)
367 for (
int i = 0; i < ndim * nvec; i++)
374 for (
int i = 0; i < num_dfts_per_proc; i++)
376 for (
int j = 0; j < ndim; j++)
380 Vmath::Vcopy(N, &(ShufV1[j][i * N]), 1, &(ShufV1_PAD_coef[0]), 1);
382 m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys[j]);
384 for (
int j = 0; j < ndim * nvec; j++)
386 Vmath::Vcopy(N, &(ShufV2[j][i * N]), 1, &(ShufV2_PAD_coef[0]), 1);
387 m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys[j]);
392 for (
int j = 0; j < nvec; j++)
395 for (
int k = 0; k < ndim; k++)
398 ShufV2_PAD_phys[j * ndim + k], 1,
399 ShufV1V2_PAD_phys, 1, ShufV1V2_PAD_phys, 1);
403 m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
406 Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1, &(ShufV1V2[j][i * N]),
414 for (
int j = 0; j < nvec; j++)
423 for (
int j = 0; j < nvec; j++)
439 int cnt = 0, cnt1 = 0;
442 for (
int n = 0; n <
m_planes.size(); ++n)
444 m_planes[n]->FwdTrans(inarray + cnt, tmparray = outarray + cnt1);
462 int cnt = 0, cnt1 = 0;
466 for (
int n = 0; n <
m_planes.size(); ++n)
468 m_planes[n]->FwdTransLocalElmt(inarray + cnt,
469 tmparray = outarray + cnt1);
486 int cnt = 0, cnt1 = 0;
490 for (
int n = 0; n <
m_planes.size(); ++n)
492 m_planes[n]->FwdTransBndConstrained(inarray + cnt,
493 tmparray = outarray + cnt1);
510 int cnt = 0, cnt1 = 0;
513 for (
int n = 0; n <
m_planes.size(); ++n)
515 m_planes[n]->BwdTrans(inarray + cnt, tmparray = outarray + cnt1);
517 cnt1 +=
m_planes[n]->GetTotPoints();
532 int cnt = 0, cnt1 = 0;
545 for (
int n = 0; n <
m_planes.size(); ++n)
547 m_planes[n]->IProductWRTBase(tmpIn + cnt, tmparray = outarray + cnt1);
565 num_dofs = inarray.size();
569 num_dofs = outarray.size();
574 int num_points_per_plane = num_dofs /
m_planes.size();
575 int num_dfts_per_proc;
576 if (!
m_session->DefinesSolverInfo(
"HomoStrip"))
578 int nP =
m_comm->GetColumnComm()->GetSize();
580 num_points_per_plane / nP + (num_points_per_plane % nP > 0);
586 num_points_per_plane / nP + (num_points_per_plane % nP > 0);
602 inarray, 1, fft_in, 1);
607 for (
int i = 0; i < num_dfts_per_proc; i++)
617 for (
int i = 0; i < num_dfts_per_proc; i++)
634 fft_out, 1, outarray, 1);
664 int nrows = blkmat->GetRows();
665 int ncols = blkmat->GetColumns();
685 out = (*blkmat) * in;
711 return matrixIter->second;
721 int num_trans_per_proc = 0;
730 n_exp =
m_planes[0]->GetTotPoints();
733 num_trans_per_proc = n_exp /
m_comm->GetColumnComm()->GetSize() +
734 (n_exp %
m_comm->GetColumnComm()->GetSize() > 0);
804 for (
int i = 0; i < num_trans_per_proc; ++i)
806 BlkMatrix->SetBlock(i, i, loc_mat);
815 std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
820 std::vector<NekDouble> HomoLen;
821 HomoLen.push_back(
m_lhom);
823 std::vector<unsigned int> StripsIDs;
826 m_session->MatchSolverInfo(
"HomoStrip",
"True", strips,
false);
832 std::vector<unsigned int> PlanesIDs;
842 for (
int i = 0; i <
m_planes.size(); i++)
847 m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, HomoBasis, HomoLen,
848 strips, StripsIDs, PlanesIDs);
853 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
858 std::vector<NekDouble> HomoLen;
859 HomoLen.push_back(
m_lhom);
861 std::vector<unsigned int> StripsIDs;
864 m_session->MatchSolverInfo(
"HomoStrip",
"True", strips,
false);
870 std::vector<unsigned int> PlanesIDs;
878 for (
int i = 0; i <
m_planes.size(); i++)
884 m_planes[0]->GeneralGetFieldDefinitions(fielddef, 1, HomoBasis, HomoLen,
885 strips, StripsIDs, PlanesIDs);
899 int ncoeffs_per_plane =
m_planes[0]->GetNcoeffs();
905 for (i = 0; i <
m_planes[0]->GetExpSize(); ++i)
911 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
914 int datalen = (*m_exp)[eid]->GetNcoeffs();
916 for (n = 0; n <
m_planes.size(); ++n)
928 std::vector<NekDouble> &fielddata)
936 std::vector<NekDouble> &fielddata, std::string &field,
942 int datalen = fielddata.size() / fielddef->m_fields.size();
943 std::vector<unsigned int> fieldDefHomoZids;
946 for (i = 0; i < fielddef->m_fields.size(); ++i)
948 if (fielddef->m_fields[i] == field)
955 if (i == fielddef->m_fields.size())
957 cout <<
"Field " << field <<
"not found in data file. " << endl;
962 int modes_offset = 0;
963 int planes_offset = 0;
970 for (i = 0; i <
m_planes.size(); ++i)
975 for (i = 0; i <
m_planes[0]->GetExpSize(); ++i)
981 if (fielddef->m_numHomogeneousDir)
983 nzmodes = fielddef->m_homogeneousZIDs.size();
984 fieldDefHomoZids = fielddef->m_homogeneousZIDs;
989 fieldDefHomoZids.push_back(0);
993 int ncoeffs_per_plane =
m_planes[0]->GetNcoeffs();
995 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
997 if (fielddef->m_uniOrder ==
true)
1003 fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
1011 offset += datalen * nzmodes;
1013 (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
1017 int eid = it->second;
1018 bool sameBasis =
true;
1019 for (
int j = 0; j < fielddef->m_basis.size() - 1; ++j)
1021 if (fielddef->m_basis[j] != (*
m_exp)[eid]->GetBasisType(j))
1028 for (n = 0; n < nzmodes; ++n, offset += datalen)
1039 planes_offset = it->second;
1044 planes_offset * ncoeffs_per_plane],
1049 (*m_exp)[eid]->ExtractDataToCoeffs(
1050 &fielddata[offset], fielddef->m_numModes, modes_offset,
1052 planes_offset * ncoeffs_per_plane],
1057 (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
1064 const std::shared_ptr<ExpList> &fromExpList,
1069 int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
1070 int toNcoeffs_per_plane =
m_planes[0]->GetNcoeffs();
1073 for (i = 0; i <
m_planes.size(); ++i)
1075 m_planes[i]->ExtractCoeffsToCoeffs(
1076 fromExpList->GetPlane(i),
1077 fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane * i,
1078 tocoeffs_tmp = toCoeffs + toNcoeffs_per_plane * i);
1083 int expansion, std::string var)
1088 m_planes[0]->WriteVtkPieceData(outfile, expansion, var);
1093 int nq = (*m_exp)[expansion]->GetTotPoints();
1094 int npoints_per_plane =
m_planes[0]->GetTotPoints();
1097 int outputExtraPlane = 0;
1103 outputExtraPlane = 1;
1114 int fromRank = (rank + 1) % size;
1115 int toRank = (rank == 0) ? size - 1 : rank - 1;
1119 m_StripZcomm->SendRecv(toRank, send, fromRank, extraPlane);
1124 outfile <<
" <DataArray type=\"Float64\" Name=\"" << var <<
"\">"
1127 for (
int n = 0; n <
m_planes.size(); ++n)
1131 for (i = 0; i < nq; ++i)
1137 if (outputExtraPlane)
1139 for (i = 0; i < nq; ++i)
1148 outfile <<
" </DataArray>" << endl;
1158 cnt1 =
m_planes[0]->Get1DScaledTotPoints(scale);
1161 "size of outarray does not match internal estimage");
1163 for (
int i = 0; i <
m_planes.size(); i++)
1166 m_planes[i]->PhysInterp1DScaled(scale, inarray + i * cnt,
1167 tmparray = outarray + i * cnt1);
1177 cnt =
m_planes[0]->Get1DScaledTotPoints(scale);
1178 cnt1 =
m_planes[0]->GetTotPoints();
1180 "size of outarray does not match internal estimage");
1182 for (
int i = 0; i <
m_planes.size(); i++)
1184 m_planes[i]->PhysGalerkinProjection1DScaled(
1185 scale, inarray + i * cnt, tmparray = outarray + i * cnt1);
1192 int nT_pts = inarray.size();
1204 for (
int i = 0; i <
m_planes.size(); i++)
1206 m_planes[i]->PhysDeriv(inarray + i * nP_pts, tmp2 = out_d0 + i * nP_pts,
1207 tmp3 = out_d1 + i * nP_pts);
1222 temparray = inarray;
1251 for (
int i = 0; i <
m_planes.size(); i++)
1257 tmp2 = outarray + (i -
int(
sign)) * nP_pts, 1);
1274 if (!
m_session->DefinesSolverInfo(
"HomoStrip"))
1277 "Parallelisation in the homogeneous direction "
1278 "implemented just for Fourier basis");
1283 "Parallelisation in the homogeneous direction "
1284 "implemented just for Fourier basis");
1289 ASSERTL0(
false,
"Semi-phyisical time-stepping not "
1290 "implemented yet for non-Fourier "
1300 for (
int i = 0; i < nP_pts; i++)
1303 tmp2 = outarray + i *
m_planes.size());
1320 int nT_pts = inarray.size();
1326 int dir = (int)edir;
1335 for (
int i = 0; i <
m_planes.size(); i++)
1337 m_planes[i]->PhysDeriv(edir, inarray + i * nP_pts,
1338 tmp2 = out_d + i * nP_pts);
1353 temparray = inarray;
1381 for (
int i = 0; i <
m_planes.size(); i++)
1387 tmp2 = outarray + (i -
int(
sign)) * nP_pts, 1);
1403 if (!
m_session->DefinesSolverInfo(
"HomoStrip"))
1406 "Parallelisation in the homogeneous direction "
1407 "implemented just for Fourier basis");
1412 "Parallelisation in the homogeneous direction "
1413 "implemented just for Fourier basis");
1418 ASSERTL0(
false,
"Semi-phyisical time-stepping not implemented "
1419 "yet for non-Fourier basis");
1428 for (
int i = 0; i < nP_pts; i++)
1431 tmp2 = outarray + i *
m_planes.size());
1485 for (i = 0; i < (*m_exp).size(); ++i)
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
#define sign(a, b)
return the sign(b)*a
Describes the specification for a Basis.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Array< OneD, NekDouble > m_tmpIN
NekDouble m_lhom
Width of homogeneous direction.
LibUtilities::TranspositionSharedPtr m_transposition
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
virtual NekDouble v_GetHomoLen(void)
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Array< OneD, ExpListSharedPtr > m_planes
std::unordered_map< int, int > m_zIdToPlane
LibUtilities::CommSharedPtr m_StripZcomm
Homo1DBlockMatrixMapShPtr m_homogeneous1DBlockMat
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble >> &inarray1, const Array< OneD, Array< OneD, NekDouble >> &inarray2, Array< OneD, Array< OneD, NekDouble >> &outarray)
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
LibUtilities::NektarFFTSharedPtr m_FFT_deal
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
void Homogeneous1DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, bool Shuff=true, bool UnShuff=true)
Array< OneD, NekDouble > m_tmpOUT
virtual void v_SetHomoLen(const NekDouble lhom)
LibUtilities::NektarFFTSharedPtr m_FFT
virtual ~ExpListHomogeneous1D()
Destructor.
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
DNekBlkMatSharedPtr GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype) const
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekBlkMatSharedPtr GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype) const
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
bool m_useFFT
FFT variables.
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
ExpListHomogeneous1D(const ExpansionType type)
Default constructor.
Base class for all multi-elemental spectral/hp expansions.
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
LibUtilities::CommSharedPtr m_comm
Communicator.
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
LibUtilities::SessionReaderSharedPtr m_session
Session.
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Class representing a segment element in reference space.
static BasisSharedPtr NullBasisSharedPtr
BasisManagerT & BasisManager(void)
NektarFFTFactory & GetNektarFFTFactory()
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset=0)
std::shared_ptr< Transposition > TranspositionSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
@ eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
@ eFourier
Fourier Expansion .
std::map< Homogeneous1DMatType, DNekBlkMatSharedPtr > Homo1DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)