42 string LinearisedAdvection::className =
44 "Linearised", LinearisedAdvection::create,
45 "Linearised Non-Conservative");
54 LinearisedAdvection::LinearisedAdvection() :
Advection()
68 m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
69 m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
79 if (
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
81 std::string HomoStr =
m_session->GetSolverInfo(
"HOMOGENEOUS");
84 if ((HomoStr ==
"HOMOGENEOUS1D") || (HomoStr ==
"Homogeneous1D") ||
85 (HomoStr ==
"1D") || (HomoStr ==
"Homo1D"))
92 "Need to specify ModeType as HalfMode,SingleMode or "
99 m_session->MatchSolverInfo(
"ModeType",
"MultipleModes",
116 if ((HomoStr ==
"HOMOGENEOUS2D") || (HomoStr ==
"Homogeneous2D") ||
117 (HomoStr ==
"2D") || (HomoStr ==
"Homo2D"))
127 if ((HomoStr ==
"HOMOGENEOUS3D") || (HomoStr ==
"Homogeneous3D") ||
128 (HomoStr ==
"3D") || (HomoStr ==
"Homo3D"))
140 if (
m_session->DefinesSolverInfo(
"USEFFT"))
151 int nvar =
m_session->GetVariables().size();
153 for (
int i = 0; i < nvar; ++i)
160 for (
int i = 0; i < nvar; ++i)
162 for (
int j = 0; j < nBaseDerivs; ++j)
170 "Base flow must be defined for linearised forms.");
171 string file =
m_session->GetFunctionFilename(
"BaseFlow", 0);
174 if (
m_session->DefinesParameter(
"N_slices"))
181 "Base flow should be a sequence of files.");
184 "BaseFlow_interporder should be smaller than or equal to "
193 ASSERTL0(
false,
"Number of slices must be a positive number "
211 int nq = pFields[0]->GetNpoints();
218 pFields[0]->GetCoords(x0, x1, x2);
220 for (
unsigned int i = 0; i < pFields.size(); i++)
230 for (
int i = 0; i < nvar; ++i)
235 if (
m_session->DefinesParameter(
"period"))
244 if (
m_session->GetComm()->GetRank() == 0)
246 cout <<
"baseflow info : interpolation order " <<
m_interporder
247 <<
", period " <<
m_period <<
", periodicity ";
256 cout <<
"baseflow info : files from " <<
m_start <<
" to "
259 <<
" time intervals" << endl;
270 const int nConvectiveFields,
278 ASSERTL1(nConvectiveFields == inarray.size(),
279 "Number of convective fields and Inarray are not compatible");
281 int nPointsTot = fields[0]->GetNpoints();
282 int ndim = advVel.size();
287 for (
int i = 0; i < ndim; ++i)
292 fields[i]->HomogeneousBwdTrans(advVel[i], velocity[i]);
296 velocity[i] = advVel[i];
301 for (
int i = 0; i < nDerivs; ++i)
309 for (
int i = 0; i < ndim; ++i)
317 for (
int i = 0; i < nConvectiveFields; ++i)
324 fields[i]->PhysDeriv(inarray[i], grad[0]);
329 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
334 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1], grad[2]);
338 fields[i]->HomogeneousBwdTrans(grad[0], grad[0]);
339 fields[i]->HomogeneousBwdTrans(grad[1], grad[1]);
340 fields[i]->HomogeneousBwdTrans(grad[2], grad[2]);
348 for (
int j = 1; j < nDerivs; ++j)
360 for (
int j = 0; j < lim; ++j)
363 velocity[j], 1, outarray[i], 1, outarray[i], 1);
368 fields[i]->HomogeneousFwdTrans(outarray[i], outarray[i]);
378 if (
m_session->GetSolverInfo(
"EqType") ==
"UnsteadyNavierStokes")
383 "Number of base flow variables does not match what is "
390 "Number of base flow variables does not match what is expected.");
393 int npts = inarray[0].size();
395 for (
int i = 0; i < inarray.size(); ++i)
398 "Size of base flow array does not match expected.");
415 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
416 std::vector<std::vector<NekDouble>> FieldData;
421 int numexp = pFields[0]->GetExpSize();
425 for (
int i = 0; i < numexp; ++i)
427 ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
433 fld->Import(pInfile, FieldDef, FieldData,
436 int nSessionVar =
m_session->GetVariables().size();
437 int nSessionConvVar = nSessionVar - 1;
438 int nFileVar = FieldDef[0]->m_fields.size();
439 int nFileConvVar = nFileVar - 1;
442 ASSERTL0(nFileVar == 3,
"For half mode, expect 2D2C base flow.");
446 for (
int j = 0; j < nFileConvVar; ++j)
448 for (
int i = 0; i < FieldDef.size(); ++i)
450 bool flag = FieldDef[i]->m_fields[j] ==
m_session->GetVariable(j);
452 ASSERTL0(flag, (std::string(
"Order of ") + pInfile +
453 std::string(
" data and that defined in "
454 "the session file differs"))
457 pFields[j]->ExtractDataToCoeffs(
458 FieldDef[i], FieldData[i], FieldDef[i]->m_fields[j], tmp_coeff);
463 pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff,
m_baseflow[j]);
469 int ncplane = (pFields[0]->GetNpoints()) /
m_npointsZ;
476 bool oldwavespace = pFields[j]->GetWaveSpace();
477 pFields[j]->SetWaveSpace(
false);
478 pFields[j]->BwdTrans(tmp_coeff,
m_baseflow[j]);
479 pFields[j]->SetWaveSpace(oldwavespace);
484 for (
int j = nFileConvVar; j < nSessionConvVar; ++j)
492 for (
int i = 0; i < nSessionConvVar; ++i)
514 &outarray[0], 1, &outarray[0], 1);
516 for (
int i = 2; i <
m_slices; i += 2)
518 phase = (i >> 1) * BetaT;
520 Vmath::Svtvp(npoints, cos(phase), &inarray[i * npoints], 1,
521 &outarray[0], 1, &outarray[0], 1);
522 Vmath::Svtvp(npoints, -sin(phase), &inarray[(i + 1) * npoints], 1,
523 &outarray[0], 1, &outarray[0], 1);
557 coeff[i] *= (x - ix + padleft - (
NekDouble)j) /
563 for (
int i = ix - padleft; i < ix + padright + 1; ++i)
566 &inarray[i * npoints], 1, &outarray[0], 1,
610 bool oldwavespace = field->GetWaveSpace();
611 field->SetWaveSpace(
false);
616 field->SetWaveSpace(oldwavespace);
653 for (
int i = 0; i < n_exp; ++i)
655 BlkMatrix->SetBlock(i, i, loc_mat);
670 for (
int i = 0; i < ConvectedFields; ++i)
678 size_t found = file.find(
"%d");
680 file.find(
"%d", found + 1) == string::npos,
681 "Since N_slices is specified, the filename provided for function "
682 "'BaseFlow' must include exactly one instance of the format "
683 "specifier '%d', to index the time-slices.");
684 char *
buffer =
new char[file.length() + 8];
688 sprintf(
buffer, file.c_str(), i);
690 if (
m_session->GetComm()->GetRank() == 0)
692 cout <<
"read base flow file " <<
buffer << endl;
702 for (
int k = 0; k < ConvectedFields; ++k)
704 #ifdef NEKTAR_USING_FFTW
724 for (
int i = 0; i < npoints; i++)
744 int nrows = blkmat->GetRows();
745 int ncols = blkmat->GetColumns();
762 out = (*blkmat) * in;
771 for (
int r = 0; r < sortedinarray.size(); ++r)
773 sortedinarray[0] = 0;
774 sortedoutarray[0] = 0;
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Describes the specification for a Basis.
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Defines a specification for a set of points.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
int m_npointsZ
number of points in Z direction (if homogeneous)
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
enum HomogeneousType m_HomogeneousType
NekDouble m_period
period length
int m_npointsX
number of points in X direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
int m_HomoDirec
number of homogenous directions
virtual void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &advVel, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray)
Advects a vector field.
bool m_useFFT
flag to determine if use or not the FFT for transformations
LibUtilities::SessionReaderSharedPtr m_session
Array< OneD, NekDouble > m_tmpIN
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
int m_start
number of slices
void UpdateGradBase(const int var, const MultiRegions::ExpListSharedPtr &field)
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
virtual ~LinearisedAdvection()
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
Array< OneD, NekDouble > m_tmpOUT
Array< OneD, Array< OneD, NekDouble > > m_gradBase
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
bool m_multipleModes
flag to determine if use multiple mode or not
LibUtilities::NektarFFTSharedPtr m_FFT
auxiliary variables
bool m_halfMode
flag to determine if use half mode or not
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
NekDouble m_LhomX
physical length in X direction (if homogeneous)
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble >> &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Overrides the base flow used during linearised advection.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
An abstract base class encapsulating the concept of advection of a vector field.
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Class representing a segment element in reference space.
NektarFFTFactory & GetNektarFFTFactory()
std::shared_ptr< FieldIO > FieldIOSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
static FieldMetaDataMap NullFieldMetaDataMap
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
@ eFourier
Fourier Expansion .
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
void Neg(int n, T *x, const int incx)
Negate x = -x.
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 Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)