35 #include <boost/core/ignore_unused.hpp>
45 CompressibleFlowSystem::CompressibleFlowSystem(
59 for (
int i = 0; i <
m_fields.size(); i++)
70 "No UPWINDTYPE defined in session.");
102 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
104 std::string type =
m_fields[0]->GetBndConditions()[n]->GetUserDefined();
106 if (
m_fields[0]->GetBndConditions()[n]->GetBoundaryConditionType() ==
117 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
160 "Local time stepping requires CFL parameter.");
171 "Unsupported projection type.");
173 string advName, riemName;
174 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
191 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Average");
215 int nvariables = inarray.size();
232 for (
int i = 0; i < nvariables; ++i)
236 m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
248 for (
int i = 0; i < nvariables; ++i)
262 x->Apply(
m_fields, inarray, outarray, time);
267 int nElements =
m_fields[0]->GetExpSize();
276 for (
int n = 0; n < nElements; ++n)
278 nq =
m_fields[0]->GetExp(n)->GetTotPoints();
279 offset =
m_fields[0]->GetPhys_Offset(n);
281 for (
int i = 0; i < nvariables; ++i)
284 tmp = outarray[i] + offset, 1);
298 int nvariables = inarray.size();
307 for (
int i = 0; i < nvariables; ++i)
324 "compressible Navier-Stokes equations");
342 int nvariables = inarray.size();
365 int nvariables = physarray.size();
368 for (
int i = 0; i < nvariables; ++i)
371 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
379 x->Apply(Fwd, physarray, time);
409 auto nVariables = physfield.size();
410 auto nPts = physfield[0].size();
412 constexpr
unsigned short maxVel = 3;
413 constexpr
unsigned short maxFld = 5;
416 ASSERTL1(nVariables <= maxFld,
"GetFluxVector, hard coded max fields");
418 for (
size_t p = 0;
p < nPts; ++
p)
421 std::array<NekDouble, maxFld> fieldTmp;
422 std::array<NekDouble, maxVel> velocity;
425 for (
size_t f = 0; f < nVariables; ++f)
427 fieldTmp[f] = physfield[f][
p];
436 flux[0][d][
p] = fieldTmp[d + 1];
438 velocity[d] = fieldTmp[d + 1] * oneOrho;
448 flux[f + 1][d][
p] = velocity[d] * fieldTmp[f + 1];
472 int nq = physfield[0].size();
477 nq =
m_fields[0]->Get1DScaledTotPoints(OneDptscale);
485 for (i = 0; i < nVariables; ++i)
489 m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
490 physfield_interp[i]);
504 m_fields[0]->PhysGalerkinProjection1DScaled(
505 OneDptscale, physfield_interp[i + 1], flux[0][i]);
508 m_varConv->GetVelocityVector(physfield_interp, velocity);
516 Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i + 1], 1,
517 flux_interp[i + 1][j], 1);
522 flux_interp[i + 1][i], 1);
530 m_fields[0]->PhysGalerkinProjection1DScaled(
531 OneDptscale, flux_interp[i + 1][j], flux[i + 1][j]);
545 m_fields[0]->PhysGalerkinProjection1DScaled(
559 boost::ignore_unused(inarray);
561 int nElements =
m_fields[0]->GetExpSize();
574 for (
int n = 0; n < nElements; ++n)
589 int nElements =
m_fields[0]->GetExpSize();
623 boost::ignore_unused(pFwd, pBwd);
634 bool dumpInitialConditions,
637 boost::ignore_unused(domain);
643 int phystot =
m_fields[0]->GetTotPoints();
646 m_session->LoadParameter(
"Noise", Noise, 0.0);
647 int m_nConvectiveFields =
m_fields.size();
651 int seed = -
m_comm->GetRank() * m_nConvectiveFields;
652 for (
int i = 0; i < m_nConvectiveFields; i++)
678 int n_element =
m_fields[0]->GetExpSize();
679 int expdim =
m_fields[0]->GetGraph()->GetMeshDimension();
685 for (
int i = 0; i < nfields; ++i)
687 physfields[i] =
m_fields[i]->GetPhys();
706 m_varConv->GetVelocityVector(physfields, velocity);
707 m_varConv->GetSoundSpeed(physfields, soundspeed);
709 for (
int el = 0; el < n_element; ++el)
711 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
712 offset =
m_fields[0]->GetPhys_Offset(el);
713 int nq =
m_fields[0]->GetExp(el)->GetTotPoints();
716 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
722 ->GetDerivFactors(ptsKeys);
730 for (
int i = 0; i < expdim; ++i)
732 Vmath::Vmul(nq, gmat[i], 1, velocity[0] + offset, 1,
733 tmp = stdVelocity[i] + offset, 1);
734 Vmath::Vmul(nq, gmat[i], 1, soundspeed + offset, 1,
735 tmp = stdSoundSpeed[i] + offset, 1);
736 for (
int j = 1; j < expdim; ++j)
739 velocity[j] + offset, 1,
740 stdVelocity[i] + offset, 1,
741 tmp = stdVelocity[i] + offset, 1);
743 soundspeed + offset, 1,
744 stdSoundSpeed[i] + offset, 1,
745 tmp = stdSoundSpeed[i] + offset, 1);
751 for (
int i = 0; i < expdim; ++i)
753 Vmath::Smul(nq, gmat[i][0], velocity[0] + offset, 1,
754 tmp = stdVelocity[i] + offset, 1);
755 Vmath::Smul(nq, gmat[i][0], soundspeed + offset, 1,
756 tmp = stdSoundSpeed[i] + offset, 1);
757 for (
int j = 1; j < expdim; ++j)
760 velocity[j] + offset, 1,
761 stdVelocity[i] + offset, 1,
762 tmp = stdVelocity[i] + offset, 1);
764 soundspeed + offset, 1,
765 stdSoundSpeed[i] + offset, 1,
766 tmp = stdSoundSpeed[i] + offset, 1);
772 for (
int i = 0; i < nq; ++i)
775 for (
int j = 0; j < expdim; ++j)
778 vel =
std::abs(stdVelocity[j][offset + i]) +
779 SpeedSoundFactor *
std::abs(stdSoundSpeed[j][offset + i]);
780 pntVelocity += vel * vel;
782 pntVelocity =
sqrt(pntVelocity);
783 if (pntVelocity > stdV[el])
785 stdV[el] = pntVelocity;
802 ASSERTL0(n <= 20,
"Illegal modes dimension for CFL calculation "
803 "(P has to be less then 20)");
805 NekDouble CFLDG[21] = {2.0000, 6.0000, 11.8424, 19.1569, 27.8419,
806 37.8247, 49.0518, 61.4815, 75.0797, 89.8181,
807 105.6700, 122.6200, 140.6400, 159.7300, 179.8500,
808 201.0100, 223.1800, 246.3600, 270.5300, 295.6900,
819 "coefficients not introduced yet.");
837 for (i = 0; i <
m_fields[0]->GetExpSize(); i++)
846 std::vector<std::string> &variables)
849 m_session->MatchSolverInfo(
"OutputExtraFields",
"True", extraFields,
true);
852 const int nPhys =
m_fields[0]->GetNpoints();
853 const int nCoeffs =
m_fields[0]->GetNcoeffs();
856 for (
int i = 0; i <
m_fields.size(); ++i)
874 m_varConv->GetVelocityVector(tmp, velocity);
876 m_varConv->GetTemperature(tmp, temperature);
878 m_varConv->GetSoundSpeed(tmp, soundspeed);
879 m_varConv->GetMach(tmp, soundspeed, mach);
882 m_session->LoadParameter(
"SensorOffset", sensorOffset, 1);
891 string velNames[3] = {
"u",
"v",
"w"};
894 m_fields[0]->FwdTransLocalElmt(velocity[i], velFwd[i]);
895 variables.push_back(velNames[i]);
896 fieldcoeffs.push_back(velFwd[i]);
900 m_fields[0]->FwdTransLocalElmt(temperature, TFwd);
901 m_fields[0]->FwdTransLocalElmt(entropy, sFwd);
902 m_fields[0]->FwdTransLocalElmt(soundspeed, aFwd);
903 m_fields[0]->FwdTransLocalElmt(mach, mFwd);
904 m_fields[0]->FwdTransLocalElmt(sensor, sensFwd);
906 variables.push_back(
"p");
907 variables.push_back(
"T");
908 variables.push_back(
"s");
909 variables.push_back(
"a");
910 variables.push_back(
"Mach");
911 variables.push_back(
"Sensor");
912 fieldcoeffs.push_back(pFwd);
913 fieldcoeffs.push_back(TFwd);
914 fieldcoeffs.push_back(sFwd);
915 fieldcoeffs.push_back(aFwd);
916 fieldcoeffs.push_back(mFwd);
917 fieldcoeffs.push_back(sensFwd);
926 variables.push_back(
"ArtificialVisc");
927 fieldcoeffs.push_back(sensorFwd);
949 density = physfield[0];
959 m_varConv->GetVelocityVector(physfield, velocity);
965 boost::ignore_unused(step);
967 const int nFields =
m_fields.size();
970 for (
int i = 0; i < nFields; ++i)
973 inarray[i] =
m_fields[i]->UpdatePhys();
983 for (
int i = 0; i < nFields; ++i)
994 for (
int i = 0; i < nFields; ++i)
996 L2[i] =
sqrt(residual[i] * onPoints);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
void DoAdvection(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Compute the advection terms for the right-hand side.
virtual void GetVelocity(const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &velocity)
Extract array with velocity from physfield.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the right-hand side.
NekDouble m_bndEvaluateTime
void DoDiffusion(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Add the diffusions terms to the right-hand side.
std::string m_shockCaptureType
void GetFluxVector(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &flux)
Return the flux vector for the compressible Euler equations.
Array< OneD, NekDouble > GetStabilityLimitVector(const Array< OneD, int > &ExpOrder)
Function to calculate the stability limit for DG/CG (a vector of them).
void InitAdvection()
Create advection and diffusion objects for CFS.
virtual ~CompressibleFlowSystem()
Destructor for CompressibleFlowSystem class.
virtual void v_DoDiffusion(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Apply artificial diffusion (Laplacian operator)
virtual NekDouble v_GetTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray)
Calculate the maximum timestep subject to CFL restrictions.
virtual Array< OneD, NekDouble > v_GetMaxStdVelocity(const NekDouble SpeedSoundFactor)
Compute the advection velocity in the standard space for each element of the expansion.
void GetFluxVectorDeAlias(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &flux)
Return the flux vector for the compressible Euler equations by using the de-aliasing technique.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble >> &physarray, NekDouble time)
virtual void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Set up logic for residual calculation.
NekDouble m_filterExponent
NekDouble GetStabilityLimit(int n)
Function to calculate the stability limit for DG/CG.
std::vector< CFSBndCondSharedPtr > m_bndConds
void GetElmtTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &tstep)
Calculate the maximum timestep on each element subject to CFL restrictions.
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
virtual void v_SteadyStateResidual(int step, Array< OneD, NekDouble > &L2)
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...
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
VariableConverterSharedPtr m_varConv
void SetBoundaryConditionsBwdWeight()
Set up a weight on physical boundaries for boundary condition applications.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
ArtificialDiffusionSharedPtr m_artificialDiffusion
virtual void GetDensity(const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, NekDouble > &density)
Extract array with density from physfield.
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
void InitialiseParameters()
Load CFS parameters from the session file.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetElmtCFLVals(const bool FlagAcousticCFL=true)
int m_spacedim
Spatial dimension (>= expansion dim).
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
SOLVER_UTILS_EXPORT int GetExpSize()
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
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.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
int m_nchk
Number of checkpoints written so far.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT int GetTotPoints()
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.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator()
Get the maximum timestep estimator for cfl control.
NekDouble m_cflNonAcoustic
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
@ eMixed_CG_Discontinuous
std::shared_ptr< RiemannSolver > RiemannSolverSharedPtr
A shared pointer to an EquationSystem object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
The above copyright notice and this permission notice shall be included.
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
ArtificialDiffusionFactory & GetArtificialDiffusionFactory()
Declaration of the artificial diffusion factory singleton.
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.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
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
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)