56 AdjointAdvection::create);
61 AdjointAdvection::AdjointAdvection()
77 m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
78 m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
90 if(
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
92 std::string HomoStr =
m_session->GetSolverInfo(
"HOMOGENEOUS");
95 if((HomoStr ==
"HOMOGENEOUS1D")||(HomoStr ==
"Homogeneous1D")||
96 (HomoStr ==
"1D")||(HomoStr ==
"Homo1D"))
104 if(
m_session->DefinesSolverInfo(
"ModeType"))
110 if(
m_session->DefinesSolverInfo(
"ModeType"))
126 ASSERTL0(
false,
"SolverInfo ModeType not valid");
137 if((HomoStr ==
"HOMOGENEOUS2D")||(HomoStr ==
"Homogeneous2D")||
138 (HomoStr ==
"2D")||(HomoStr ==
"Homo2D"))
148 if((HomoStr ==
"HOMOGENEOUS3D")||(HomoStr ==
"Homogeneous3D")||
149 (HomoStr ==
"3D")||(HomoStr ==
"Homo3D"))
161 if(
m_session->DefinesSolverInfo(
"USEFFT"))
171 if(
m_session->DefinesSolverInfo(
"PROJECTION"))
173 std::string ProjectStr
174 =
m_session->GetSolverInfo(
"PROJECTION");
176 if((ProjectStr ==
"Continuous")||(ProjectStr ==
"Galerkin")||
177 (ProjectStr ==
"CONTINUOUS")||(ProjectStr ==
"GALERKIN"))
181 else if(ProjectStr ==
"DisContinuous")
187 ASSERTL0(
false,
"PROJECTION value not recognised");
192 cerr <<
"Projection type not specified in SOLVERINFO,"
193 "defaulting to continuous Galerkin" << endl;
197 int nvar =
m_session->GetVariables().size();
199 for (
int i = 0; i < nvar; ++i)
205 "Base flow must be defined for linearised forms.");
206 string file =
m_session->GetFunctionFilename(
"BaseFlow", 0);
210 if(
m_session->DefinesParameter(
"N_slices"))
219 ASSERTL0(
false,
"Number of slices must be a positive number greater than 1");
237 int nq = pFields[0]->GetNpoints();
244 pFields[0]->GetCoords(x0,x1,x2);
245 for(
unsigned int i = 0 ; i <
m_baseflow.num_elements(); i++)
265 const int nConvectiveFields,
272 int nqtot = fields[0]->GetTotPoints();
273 ASSERTL1(nConvectiveFields == inarray.num_elements(),
"Number of convective fields and Inarray are not compatible");
277 for(
int n = 0; n < nConvectiveFields; ++n)
280 int ndim = advVel.num_elements();
281 int nPointsTot = fields[0]->GetNpoints();
303 if (
m_session->GetFunctionType(
"BaseFlow", 0)
306 for(
int i=0; i<ndim;++i)
313 ASSERTL0(
false,
"Periodic Base flow requires .fld files");
322 fields[0]->PhysDeriv(inarray[n],grad0);
323 fields[0]->PhysDeriv(
m_baseflow[0],grad_base_u0);
327 Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
337 fields[0]->PhysDeriv(inarray[n],grad0,grad1);
340 fields[0]-> PhysDeriv(
m_baseflow[0], grad_base_u0, grad_base_u1);
341 fields[0]-> PhysDeriv(
m_baseflow[1], grad_base_v0, grad_base_v1);
356 Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
358 Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
370 Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
372 Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
391 fields[0]->PhysDeriv(
m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
392 fields[0]->PhysDeriv(
m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
393 fields[0]->PhysDeriv(
m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
398 for(
int i=0; i<grad_base_u2.num_elements();++i)
407 fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
422 Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
424 Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
426 Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[2],1,outarray[n],1,outarray[n],1);
439 Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
441 Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
443 Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[2],1,outarray[n],1,outarray[n],1);
457 Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[0],1,outarray[n],1,outarray[n],1);
459 Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[1],1,outarray[n],1,outarray[n],1);
461 Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
466 ASSERTL0(
false,
"dimension unknown");
479 "Number of base flow variables does not match what is"
482 int npts = inarray[0].num_elements();
483 for (
int i = 0; i < inarray.num_elements(); ++i)
486 "Size of base flow array does not match expected.");
503 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
504 std::vector<std::vector<NekDouble> > FieldData;
509 int numexp = pFields[0]->GetExpSize();
513 for(
int i = 0; i < numexp; ++i)
515 ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
521 fld->Import(pInfile, FieldDef, FieldData,
525 int nSessionVar =
m_session->GetVariables().size();
526 int nSessionConvVar = nSessionVar - 1;
527 int nFileVar = FieldDef[0]->m_fields.size();
528 int nFileConvVar = nFileVar - 1;
531 ASSERTL0(nFileVar == 3,
"For half mode, expect 2D2C base flow.");
535 for(
int j = 0; j < nFileConvVar; ++j)
537 for(
int i = 0; i < FieldDef.size(); ++i)
539 bool flag = FieldDef[i]->m_fields[j] ==
542 ASSERTL0(flag, (std::string(
"Order of ") + pInfile
543 + std::string(
" data and that defined in "
544 "the session file differs")).c_str());
546 pFields[j]->ExtractDataToCoeffs(
549 FieldDef[i]->m_fields[j],
555 pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff,
m_baseflow[j]);
561 int ncplane=(pFields[0]->GetNpoints())/
m_npointsZ;
567 bool oldwavespace = pFields[j]->GetWaveSpace();
568 pFields[j]->SetWaveSpace(
false);
569 pFields[j]->BwdTrans(tmp_coeff,
m_baseflow[j]);
570 pFields[j]->SetWaveSpace(oldwavespace);
575 for (
int j = nFileConvVar; j < nSessionConvVar; ++j) {
580 if(
m_session->DefinesParameter(
"N_slices"))
582 for(
int i = 0; i < nSessionConvVar; ++i)
605 Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
607 for (
int i = 2; i <
m_slices; i += 2)
609 phase = (i>>1) * BetaT;
611 Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
612 Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
647 for(
int i = 0; i < n_exp; ++i)
649 BlkMatrix->SetBlock(i,i,loc_mat);
663 int ConvectedFields=
m_baseflow.num_elements()-1;
666 for(
int i=0; i<ConvectedFields;++i)
674 size_t found = file.find(
"%d");
675 ASSERTL0(found != string::npos && file.find(
"%d", found+1) == string::npos,
676 "Since N_slices is specified, the filename provided for function "
677 "'BaseFlow' must include exactly one instance of the format "
678 "specifier '%d', to index the time-slices.");
679 char*
buffer =
new char[file.length() + 8];
682 sprintf(buffer, file.c_str(), i);
689 for(
int k=0; k< ConvectedFields;++k)
691 #ifdef NEKTAR_USING_FFTW
711 for(
int i=0; i<npoints; i++)
713 m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
733 int nrows = blkmat->GetRows();
734 int ncols = blkmat->GetColumns();
758 Vmath::Zero(sortedinarray.num_elements(),&sortedinarray[0],1);
759 Vmath::Zero(sortedoutarray.num_elements(),&sortedoutarray[0],1);
765 for (
int i = 2; i <
m_slices; i += 2)
774 if(
m_session->DefinesParameter(
"period"))
Array< OneD, NekDouble > m_tmpOUT
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
#define ASSERTL0(condition, msg)
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.
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
LibUtilities::SessionReaderSharedPtr m_session
int m_npointsX
number of points in X direction (if homogeneous)
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
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
bool m_homogen_dealiasing
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
LibUtilities::NektarFFTSharedPtr m_FFT
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
NektarFFTFactory & GetNektarFFTFactory()
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
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.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Class representing a segment element in reference space.
Array< OneD, Array< OneD, NekDouble > > m_interp
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
1D Evenly-spaced points using Fourier Fit
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
virtual ~AdjointAdvection()
int m_npointsY
number of points in Y direction (if homogeneous)
bool m_HalfMode
flag to determine if use half mode or not
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
int m_HomoDirec
number of homogenous directions
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
boost::shared_ptr< Equation > EquationSharedPtr
NekDouble m_LhomX
physical length in X direction (if homogeneous)
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
MultiRegions::CoeffState m_CoeffState
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
enum HomogeneousType m_HomogeneousType
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Overrides the base flow used during linearised advection.
bool m_SingleMode
flag to determine if use single mode or not
int m_npointsZ
number of points in Z direction (if homogeneous)
Array< OneD, NekDouble > m_tmpIN
MultiRegions::ProjectionType m_projectionType
void Zero(int n, T *x, const int incx)
Zero vector.
bool m_useFFT
flag to determine if use or not the FFT for transformations
#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
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.
bool m_MultipleModes
flag to determine if use multiple mode or not
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.