53 string LinearisedAdvection::className
56 LinearisedAdvection::create);
65 LinearisedAdvection::LinearisedAdvection():
80 m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
81 m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
91 if(
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
93 std::string HomoStr =
m_session->GetSolverInfo(
"HOMOGENEOUS");
96 if((HomoStr ==
"HOMOGENEOUS1D")||(HomoStr ==
"Homogeneous1D")||
97 (HomoStr ==
"1D")||(HomoStr ==
"Homo1D"))
104 "Need to specify ModeType as HalfMode,SingleMode or "
107 m_session->MatchSolverInfo(
"ModeType",
"SingleMode",
109 m_session->MatchSolverInfo(
"ModeType",
"HalfMode",
111 m_session->MatchSolverInfo(
"ModeType",
"MultipleModes",
128 if((HomoStr ==
"HOMOGENEOUS2D")||(HomoStr ==
"Homogeneous2D")||
129 (HomoStr ==
"2D")||(HomoStr ==
"Homo2D"))
139 if((HomoStr ==
"HOMOGENEOUS3D")||(HomoStr ==
"Homogeneous3D")||
140 (HomoStr ==
"3D")||(HomoStr ==
"Homo3D"))
152 if(
m_session->DefinesSolverInfo(
"USEFFT"))
162 int nvar =
m_session->GetVariables().size();
164 for (
int i = 0; i < nvar; ++i)
170 "Base flow must be defined for linearised forms.");
171 string file =
m_session->GetFunctionFilename(
"BaseFlow", 0);
175 if(
m_session->DefinesParameter(
"N_slices"))
184 ASSERTL0(
false,
"Number of slices must be a positive number "
203 int nq = pFields[0]->GetNpoints();
210 pFields[0]->GetCoords(x0,x1,x2);
212 for(
unsigned int i = 0 ; i < pFields.num_elements(); i++)
222 if(
m_session->DefinesParameter(
"period"))
241 const int nConvectiveFields,
248 ASSERTL1(nConvectiveFields == inarray.num_elements(),
249 "Number of convective fields and Inarray are not compatible");
251 int nPointsTot = fields[0]->GetNpoints();
252 int ndim = advVel.num_elements();
279 "Base flow should be a sequence of files.");
281 for (
int i = 0; i < ndim; ++i)
293 ASSERTL0(
false,
"Not set up for 1D");
302 fields[0]-> PhysDeriv(
m_baseflow[0], grad_base_u0, grad_base_u1);
303 fields[0]-> PhysDeriv(
m_baseflow[1], grad_base_v0, grad_base_v1);
306 fields[0]->PhysDeriv(inarray[0],grad0,grad1);
313 Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[0],1,
316 Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[0],1,
321 fields[0]->PhysDeriv(inarray[1],grad0,grad1);
328 Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[1],1,
331 Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[1],1,
354 grad_base_u0, grad_base_u1);
356 grad_base_v0, grad_base_v1);
360 fields[0]->PhysDeriv(
m_baseflow[0], grad_base_u0,
362 fields[0]->PhysDeriv(
m_baseflow[1], grad_base_v0,
364 fields[0]->PhysDeriv(
m_baseflow[2], grad_base_w0,
370 bool oldwavespace = fields[0]->GetWaveSpace();
371 fields[0]->SetWaveSpace(
false);
372 fields[0]->PhysDeriv(
m_baseflow[0], grad_base_u0,
373 grad_base_u1, grad_base_u2);
374 fields[0]->PhysDeriv(
m_baseflow[1], grad_base_v0,
375 grad_base_v1, grad_base_v2);
376 fields[0]->PhysDeriv(
m_baseflow[2], grad_base_w0,
377 grad_base_w1, grad_base_w2);
378 fields[0]->SetWaveSpace(oldwavespace);
382 ASSERTL0(
false,
"ERROR: Must be one of half, single or multiple modes");
391 fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
393 fields[0]->HomogeneousBwdTrans(grad0, grad0);
394 fields[0]->HomogeneousBwdTrans(grad1, grad1);
395 fields[0]->HomogeneousBwdTrans(grad2, grad2);
401 fields[0]->PhysDeriv(inarray[0], grad0, grad1);
405 fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
414 outarray[0], 1, outarray[0], 1);
417 outarray[0],1,outarray[0],1);
420 outarray[0],1,outarray[0],1);
426 1,outarray[0], 1, outarray[0], 1);
433 advVel[2],1,outarray[0],1,outarray[0],1);
434 fields[0]->HomogeneousFwdTrans(outarray[0],outarray[0]);
441 fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
443 fields[0]->HomogeneousBwdTrans(grad0, grad0);
444 fields[0]->HomogeneousBwdTrans(grad1, grad1);
445 fields[0]->HomogeneousBwdTrans(grad2, grad2);
451 fields[0]->PhysDeriv(inarray[1], grad0, grad1);
455 fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
464 outarray[1], 1, outarray[1], 1);
467 ,outarray[1],1,outarray[1],1);
470 outarray[1],1,outarray[1],1);
476 outarray[1], 1, outarray[1], 1);
483 outarray[1],1,outarray[1],1);
484 fields[0]->HomogeneousFwdTrans(outarray[1],outarray[1]);
492 fields[0]->PhysDeriv(inarray[2], grad0, grad1, grad2);
494 fields[0]->HomogeneousBwdTrans(grad0, grad0);
495 fields[0]->HomogeneousBwdTrans(grad1, grad1);
496 fields[0]->HomogeneousBwdTrans(grad2, grad2);
502 fields[0]->PhysDeriv(inarray[2], grad0, grad1);
506 fields[0]->PhysDeriv(inarray[2], grad0, grad1, grad2);
515 outarray[2], 1, outarray[2], 1);
521 outarray[2],1,outarray[2],1);
524 outarray[2],1,outarray[2],1);
528 outarray[2], 1, outarray[2], 1);
535 outarray[2],1,outarray[2],1);
536 fields[0]->HomogeneousFwdTrans(outarray[2],outarray[2]);
546 if (
m_session->GetSolverInfo(
"EqType") ==
"UnsteadyNavierStokes")
551 "Number of base flow variables does not match what is "
557 "Number of base flow variables does not match what is expected.");
560 int npts = inarray[0].num_elements();
562 for (
int i = 0; i < inarray.num_elements(); ++i)
565 "Size of base flow array does not match expected.");
583 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
584 std::vector<std::vector<NekDouble> > FieldData;
589 int numexp = pFields[0]->GetExpSize();
593 for(
int i = 0; i < numexp; ++i)
595 ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
601 fld->Import(pInfile, FieldDef, FieldData,
605 int nSessionVar =
m_session->GetVariables().size();
606 int nSessionConvVar = nSessionVar - 1;
607 int nFileVar = FieldDef[0]->m_fields.size();
608 int nFileConvVar = nFileVar - 1;
611 ASSERTL0(nFileVar == 3,
"For half mode, expect 2D2C base flow.");
615 for(
int j = 0; j < nFileConvVar; ++j)
617 for(
int i = 0; i < FieldDef.size(); ++i)
619 bool flag = FieldDef[i]->m_fields[j] ==
622 ASSERTL0(flag, (std::string(
"Order of ") + pInfile
623 + std::string(
" data and that defined in "
624 "the session file differs")).c_str());
626 pFields[j]->ExtractDataToCoeffs(
629 FieldDef[i]->m_fields[j],
635 pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff,
m_baseflow[j]);
641 int ncplane=(pFields[0]->GetNpoints())/
m_npointsZ;
647 bool oldwavespace = pFields[j]->GetWaveSpace();
648 pFields[j]->SetWaveSpace(
false);
649 pFields[j]->BwdTrans(tmp_coeff,
m_baseflow[j]);
650 pFields[j]->SetWaveSpace(oldwavespace);
655 for (
int j = nFileConvVar; j < nSessionConvVar; ++j) {
660 if(
m_session->DefinesParameter(
"N_slices"))
662 for(
int i = 0; i < nSessionConvVar; ++i)
683 Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
685 for (
int i = 2; i <
m_slices; i += 2)
687 phase = (i>>1) * BetaT;
689 Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
690 Vmath::Svtvp(npoints, -sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
725 for(
int i = 0; i < n_exp; ++i)
727 BlkMatrix->SetBlock(i,i,loc_mat);
738 int ConvectedFields =
m_baseflow.num_elements()-1;
742 for (
int i = 0; i < ConvectedFields; ++i)
750 size_t found = file.find(
"%d");
751 ASSERTL0(found != string::npos && file.find(
"%d", found+1) == string::npos,
752 "Since N_slices is specified, the filename provided for function "
753 "'BaseFlow' must include exactly one instance of the format "
754 "specifier '%d', to index the time-slices.");
755 char*
buffer =
new char[file.length() + 8];
758 sprintf(buffer, file.c_str(), i);
765 for(
int k=0; k< ConvectedFields;++k)
767 #ifdef NEKTAR_USING_FFTW
785 for(
int i=0; i<npoints; i++)
787 m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
805 int nrows = blkmat->GetRows();
806 int ncols = blkmat->GetColumns();
830 for(
int r=0; r<sortedinarray.num_elements();++r)
840 for (
int i = 2; i <
m_slices; i += 2)
enum HomogeneousType m_HomogeneousType
virtual ~LinearisedAdvection()
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
#define ASSERTL0(condition, msg)
Array< OneD, NekDouble > m_tmpIN
int m_HomoDirec
number of homogenous directions
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
LibUtilities::NektarFFTSharedPtr m_FFT
auxiliary variables
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
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 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
int m_slices
number of slices
Array< OneD, NekDouble > m_tmpOUT
boost::shared_ptr< DNekMat > DNekMatSharedPtr
NektarFFTFactory & GetNektarFFTFactory()
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Overrides the base flow used during linearised advection.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
Class representing a segment element in reference space.
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
1D Evenly-spaced points using Fourier Fit
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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)
Advects a vector field.
int m_npointsZ
number of points in Z direction (if homogeneous)
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Defines a specification for a set of points.
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< FieldIO > FieldIOSharedPtr
NekDouble m_period
period length
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
int m_npointsY
number of points in Y direction (if homogeneous)
boost::shared_ptr< Equation > EquationSharedPtr
bool m_multipleModes
flag to determine if use multiple mode or not
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
int m_npointsX
number of points in X direction (if homogeneous)
NekDouble m_LhomX
physical length in X direction (if homogeneous)
bool m_halfMode
flag to determine if use half mode or not
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
void Zero(int n, T *x, const int incx)
Zero vector.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
static FieldMetaDataMap NullFieldMetaDataMap
bool m_useFFT
flag to determine if use or not the FFT for transformations
Describes the specification for a Basis.
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
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.
Defines a callback function which evaluates the flux vector.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.