39#include <boost/version.hpp>
40#if BOOST_VERSION >= 106900 && BOOST_VERSION < 107000
41#define BOOST_ALLOW_DEPRECATED_HEADERS
44#include <boost/random/normal_distribution.hpp>
45#include <boost/random/variate_generator.hpp>
58 m_ip(-1), m_irho(-1), m_iu(1), m_conservative(false)
71 "Only Projection=DisContinuous supported by the AcousticSystem class.");
97 if (
m_session->DefinesElement(
"Nektar/Coupling"))
99 TiXmlElement *vCoupling =
m_session->GetElement(
"Nektar/Coupling");
101 ASSERTL0(vCoupling->Attribute(
"TYPE"),
102 "Missing TYPE attribute in Coupling");
103 std::string vType = vCoupling->Attribute(
"TYPE");
105 "TYPE attribute must be non-empty in Coupling");
124 int numForceFields = 0;
127 numForceFields += x->GetForces().size();
129 std::vector<std::string> varNames;
132 for (
int i = 0; i <
m_fields.size(); ++i)
134 varNames.push_back(
m_session->GetVariable(i));
135 phys[i] =
m_fields[i]->UpdatePhys();
137 for (
int i = 0; i <
m_bfNames.size(); ++i)
146 for (
int i = 0; i < x->GetForces().size(); ++i)
150 varNames.push_back(
"F_" + std::to_string(f) +
"_" +
180 int nVariables = inarray.size();
188 for (
int i = 0; i < nVariables; ++i)
209 int nvariables = inarray.size();
213 if (inarray != outarray)
215 for (
int i = 0; i < nvariables; ++i)
240 for (
int i = 0; i < nvariables; ++i)
243 m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
248 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
250 std::string userDefStr =
251 m_fields[0]->GetBndConditions()[n]->GetUserDefined();
253 if (!userDefStr.empty())
256 if (boost::iequals(userDefStr,
"Wall"))
258 WallBC(n, cnt, Fwd, inarray);
260 else if (boost::iequals(userDefStr,
"WhiteNoise"))
264 else if (boost::iequals(userDefStr,
"RiemannInvariantBC"))
268 else if (boost::iequals(userDefStr,
"TimeDependent"))
270 for (
int i = 0; i < nvariables; ++i)
273 m_fields[i]->EvaluateBoundaryConditions(time, varName);
278 std::string errmsg =
"Unrecognised boundary condition: ";
279 errmsg += userDefStr;
285 for (
int i = 0; i < nvariables; ++i)
288 m_fields[i]->EvaluateBoundaryConditions(time, varName);
292 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
303 int nVariables = physarray.size();
309 int id1, id2, nBCEdgePts;
310 int eMax =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
312 for (
int e = 0; e < eMax; ++e)
315 ->GetBndCondExpansions()[bcRegion]
318 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
319 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
332 Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
338 &Fwd[
m_iu + i][id2], 1, &Fwd[
m_iu + i][id2], 1);
342 for (
int i = 0; i < nVariables; ++i)
346 ->GetBndCondExpansions()[bcRegion]
347 ->UpdatePhys())[id1],
357 int bcRegion,
int cnt,
362 int id1, id2, nBCEdgePts;
363 int nVariables = physarray.size();
367 if (
m_rng.count(bcRegion) == 0)
369 m_rng[bcRegion] = boost::mt19937(bcRegion);
373 m_fields[0]->GetBndConditions()[bcRegion]->GetBoundaryConditionType() ==
375 "WhiteNoise BCs must be Dirichlet type BCs");
378 std::static_pointer_cast<SpatialDomains::DirichletBoundaryCondition>(
379 m_fields[0]->GetBndConditions()[bcRegion])
380 ->m_dirichletCondition;
384 "sigma must be greater than zero");
391 boost::normal_distribution<> dist(0, sigma);
395 int eMax =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
396 for (
int e = 0; e < eMax; ++e)
399 ->GetBndCondExpansions()[bcRegion]
402 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
403 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
406 for (
int i = 0; i < nVariables; ++i)
416 for (
int i = 0; i < nBCEdgePts; ++i)
433 for (
int i = 0; i < nBCEdgePts; ++i)
437 (
sqrt(BfFwd[0][id2 + i]) * BfFwd[1][id2 + i]);
447 for (
int i = 0; i < nVariables; ++i)
451 ->GetBndCondExpansions()[bcRegion]
452 ->UpdatePhys())[id1],
465 [[maybe_unused]]
const NekDouble SpeedSoundFactor)
467 int nElm =
m_fields[0]->GetExpSize();
477 for (
int el = 0; el < nElm; ++el)
479 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
483 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
489 ->GetDerivFactors(ptsKeys);
491 int nq =
m_fields[0]->GetExp(el)->GetTotPoints();
498 for (
int j = 0; j < nq; ++j)
503 velocity[i][j] =
m_bf[i + 2][cnt + j] + c;
513 Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1, stdVelocity[i], 1);
517 1, stdVelocity[i], 1, stdVelocity[i], 1);
525 Vmath::Smul(nq, gmat[i][0], velocity[0], 1, stdVelocity[i], 1);
529 1, stdVelocity[i], 1, stdVelocity[i], 1);
535 for (
int i = 0; i < nq; ++i)
540 pntVelocity += stdVelocity[j][i] * stdVelocity[j][i];
542 pntVelocity =
sqrt(pntVelocity);
544 if (pntVelocity > stdV[el])
546 stdV[el] = pntVelocity;
558 std::vector<std::string> &variables)
560 for (
int i = 0; i <
m_bfNames.size(); i++)
565 fieldcoeffs.push_back(tmpC);
571 for (
int i = 0; i < x->GetForces().size(); ++i)
573 variables.push_back(
"F_" + std::to_string(f) +
"_" +
576 m_fields[0]->FwdTrans(x->GetForces()[i], tmpC);
577 fieldcoeffs.push_back(tmpC);
585 for (
int i = 0; i <
m_bfNames.size(); i++)
598 for (
int bcRegion = 0; bcRegion <
m_fields[0]->GetBndConditions().size();
610 ->GetBndCondExpansions()[bcRegion]
613 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
614 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt +
620 cnt +=
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
#define ASSERTL0(condition, msg)
SolverUtils::CouplingSharedPtr m_coupling
NekDouble m_whiteNoiseBC_p
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the right-hand side.
Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor) override
Compute the advection velocity in the standard space for each element of the expansion.
std::vector< std::string > m_bfNames
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
void CopyBoundaryTrace(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
bool m_conservative
we are dealing with a conservative formualtion
void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables) override
virtual void v_RiemannInvariantBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &BfFwd, Array< OneD, Array< OneD, NekDouble > > &physarray)=0
SolverUtils::AdvectionSharedPtr m_advection
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the AcousticSystem equations.
AcousticSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
void WhiteNoiseBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &BfFwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the AcousticSystem equations.
int m_ip
indices of the fields
std::map< int, boost::mt19937 > m_rng
Array< OneD, Array< OneD, NekDouble > > m_bfFwdBwd
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...
Array< OneD, Array< OneD, NekDouble > > m_bf
void v_InitObject(bool DeclareFields=true) override
Initialization object for the AcousticSystem class.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
void UpdateBasefieldFwdBwd()
const Array< OneD, const Array< OneD, NekDouble > > & GetBasefieldFwdBwd()
virtual void v_AddLinTerm(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)=0
NekDouble m_whiteNoiseBC_lastUpdate
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Apply the Boundary Conditions to the AcousticSystem equations.
bool v_PreIntegrate(int step) override
v_PreIntegrate
NekDouble Evaluate() const
A base class for PDEs which include an advection component.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT int GetExpSize()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
SOLVER_UTILS_EXPORT int GetTotPoints()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
virtual SOLVER_UTILS_EXPORT void v_Output(void)
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Base class for unsteady solvers.
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate(int step)
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
static const NekDouble kNekZeroTol
CouplingFactory & GetCouplingFactory()
Declaration of the Coupling factory singleton.
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
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)
scalarT< T > sqrt(scalarT< T > in)