55 "APE1/APE4 (Acoustic Perturbation Equations)");
77 "Only Projection=DisContinuous supported by the APE class.");
120 ASSERTL0(
false,
"Expansion dimension not recognised");
136 for (
int i = 0; i <
m_bf.num_elements(); ++i)
163 m_vecLocs[0][i] = 1 + i;
167 m_session->LoadSolverInfo(
"UpwindType", riemName,
"APEUpwind");
177 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
191 ASSERTL0(
false,
"Implicit APE not set up.");
207 int nElm =
m_fields[0]->GetExpSize();
221 for (
int el = 0; el < nElm; ++el)
223 NekDouble lambdaMax = stdVelocity[el] * cLambda
224 * (expOrder[el] - 1) * (expOrder[el] - 1);
246 int nq = physfield[0].num_elements();
251 "Dimension of flux array and velocity array do not match");
260 Vmath::Vmul(nq, tmp1, 1, physfield[j+1], 1, tmp1, 1);
269 for (
int i = 1; i < flux.num_elements(); ++i)
271 ASSERTL1(flux[i].num_elements() == m_spacedim,
272 "Dimension of flux array and velocity array do not match");
292 Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
307 std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
310 for (
int i = 0; i < (*x)->GetForces().num_elements(); ++i)
312 m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
313 m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
316 m_bfField->BwdTrans(tmpC, (*x)->UpdateForces()[i]);
324 m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
342 if (
m_comm->GetRank() == 0)
344 cout <<
"CFL: " << cfl << endl;
359 int nVariables = inarray.num_elements();
367 for (
int i = 0; i < nVariables; ++i)
372 std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
388 int nvariables = inarray.num_elements();
392 for (
int i = 0; i < nvariables; ++i)
408 int nvariables =
m_fields.num_elements();
415 for (
int i = 0; i < nvariables; ++i)
418 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
422 for(
int n = 0; n <
m_fields[0]->GetBndConditions().num_elements(); ++n)
425 if (boost::iequals(
m_fields[0]->GetBndConditions()[n]->GetUserDefined(),
"Wall"))
427 WallBC(n, cnt, Fwd, inarray);
431 if (
m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
433 for (
int i = 0; i < nvariables; ++i)
436 m_fields[i]->EvaluateBoundaryConditions(time, varName);
439 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
451 int nVariables = physarray.num_elements();
457 int id1, id2, nBCEdgePts;
458 int eMax =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
460 for (
int e = 0; e < eMax; ++e)
462 nBCEdgePts =
m_fields[0]->GetBndCondExpansions()[bcRegion]->
463 GetExp(e)->GetTotPoints();
464 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
465 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
481 Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
494 for (
int i = 0; i < nVariables; ++i)
498 &(
m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1], 1);
512 int nElm =
m_fields[0]->GetExpSize();
514 ASSERTL1(stdV.num_elements() == nElm,
"stdV malformed");
525 for (
int el = 0; el < nElm; ++el)
527 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
531 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
533 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()
534 ->GetDerivFactors(ptsKeys);
536 int nq =
m_fields[0]->GetExp(el)->GetTotPoints();
543 for (
int j = 0; j < nq; ++j)
548 velocity[i][j] =
m_bf[i+2][cnt+j] + c;
558 Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1, stdVelocity[i], 1);
561 Vmath::Vvtvp(nq, gmat[m_spacedim * j + i], 1, velocity[j],
562 1, stdVelocity[i], 1, stdVelocity[i], 1);
570 Vmath::Smul(nq, gmat[i][0], velocity[0], 1, stdVelocity[i], 1);
573 Vmath::Svtvp(nq, gmat[m_spacedim * j + i][0], velocity[j],
574 1, stdVelocity[i], 1, stdVelocity[i], 1);
580 for (
int i = 0; i < nq; ++i)
585 pntVelocity += stdVelocity[j][i] * stdVelocity[j][i];
587 pntVelocity = sqrt(pntVelocity);
589 if (pntVelocity > stdV[el])
591 stdV[el] = pntVelocity;
602 std::vector<std::string> &variables)
610 m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
615 fieldcoeffs.push_back(tmpC);
619 std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
622 for (
int i = 0; i < (*x)->GetForces().num_elements(); ++i)
626 m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
627 m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
631 variables.push_back(
"F_" + boost::lexical_cast<string>(f) +
633 fieldcoeffs.push_back(tmpC);
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the APE equations.
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Get the normal vectors.
virtual bool v_PostIntegrate(int step)
v_PostIntegrate
#define ASSERTL0(condition, msg)
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
SolverUtils::AdvectionSharedPtr m_advection
std::vector< PointsKey > PointsKeyVector
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
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.
MultiRegions::ExpListSharedPtr m_bfField
NekDouble m_time
Current time of simulation.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp()
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
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
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
virtual bool v_PreIntegrate(int step)
v_PostIntegrate
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the APE equations.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
LibUtilities::CommSharedPtr m_comm
Communicator.
SOLVER_UTILS_EXPORT int GetTotPoints()
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
NekDouble GetGamma()
Get the heat capacity ratio.
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
void GetStdVelocity(Array< OneD, NekDouble > &stdV)
Compute the advection velocity in the standard space for each element of the expansion.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Get the locations of the components of the directed fields within the fields array.
Array< OneD, Array< OneD, NekDouble > > m_traceBasefield
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate(int step)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual ~APE()
Destructor.
Base class for unsteady solvers.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the right-hand side.
RiemannSolverFactory & GetRiemannSolverFactory()
int m_cflsteps
dump cfl estimate
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
void Neg(int n, T *x, const int incx)
Negate x = -x.
virtual void v_InitObject()
Initialization object for the APE class.
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Apply the Boundary Conditions to the APE equations.
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
NekDouble m_gamma
Isentropic coefficient, Ratio of specific heats (APE)
const Array< OneD, const Array< OneD, NekDouble > > & GetBasefield()
Get the baseflow field.
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
NekDouble GetCFLEstimate()
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
void Zero(int n, T *x, const int incx)
Zero vector.
#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)
Geometry is curved or has non-constant factors.
Array< OneD, Array< OneD, NekDouble > > m_bf
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
std::vector< std::string > m_bfNames