41 string NavierStokesCFE::className =
43 "NavierStokesCFE", NavierStokesCFE::create,
44 "NavierStokes equations in conservative variables.");
46 NavierStokesCFE::NavierStokesCFE(
73 m_session->LoadParameter(
"GasConstant", gasConstant, 287.058);
88 if (
m_session->DefinesParameter(
"thermalConductivity"))
91 "Cannot define both Pr and thermalConductivity.");
93 m_session->LoadParameter(
"thermalConductivity",
109 m_session->LoadSolverInfo(
"DiffusionType", diffName,
"LDGNS");
113 if (
"InteriorPenalty" == diffName)
119 if (
"LDGNS" == diffName ||
"LDGNS3DHomogeneous1D" == diffName)
136 &NavierStokesCFE::GetViscousFluxVectorConservVar<false>,
this);
139 &NavierStokesCFE::GetViscousFluxVectorConservVar<true>,
this);
156 size_t nvariables = inarray.size();
162 for (
int i = 0; i < nvariables; ++i)
185 for (
int i = 0; i < nvariables; ++i)
187 Vmath::Vadd(npoints, outarrayDiff[i], 1, outarray[i], 1,
198 for (
int i = 0; i < nvariables - 1; ++i)
207 m_varConv->GetPressure(inarray, inarrayDiff[0]);
210 m_varConv->GetTemperature(inarray, inarrayDiff[nvariables - 2]);
213 m_varConv->GetVelocityVector(inarray, inarrayDiff);
227 m_varConv->GetTemperature(pFwd, inFwd[nvariables - 2]);
228 m_varConv->GetTemperature(pBwd, inBwd[nvariables - 2]);
230 m_varConv->GetVelocityVector(pFwd, inFwd);
231 m_varConv->GetVelocityVector(pBwd, inBwd);
238 for (
int i = 0; i < nvariables; ++i)
240 Vmath::Vadd(npoints, outarrayDiff[i], 1, outarray[i], 1,
256 size_t nScalar = physfield.size();
257 size_t nPts = physfield[0].size();
267 thermalConductivity);
272 Vmath::Vadd(nPts, divVel, 1, derivativesO1[j][j], 1, divVel, 1);
290 Vmath::Vadd(nPts, derivativesO1[i][j], 1, derivativesO1[j][i], 1,
291 viscousTensor[i][j + 1], 1);
293 Vmath::Vmul(nPts, mu, 1, viscousTensor[i][j + 1], 1,
294 viscousTensor[i][j + 1], 1);
299 Vmath::Vadd(nPts, viscousTensor[i][j + 1], 1, divVel, 1,
300 viscousTensor[i][j + 1], 1);
306 viscousTensor[j][i + 1], 1);
318 Vmath::Vvtvp(nPts, physfield[j], 1, viscousTensor[i][j + 1], 1,
341 size_t nScalar = physfield.size();
342 int nPts =
m_fields[0]->Get1DScaledTotPoints(OneDptscale);
343 size_t nPts_orig = physfield[0].size();
355 thermalConductivity);
365 m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
373 m_fields[0]->PhysInterp1DScaled(OneDptscale, derivativesO1[i][j],
388 Vmath::Vadd(nPts, divVel, 1, deriv_interp[j][j], 1, divVel, 1);
406 Vmath::Vadd(nPts, deriv_interp[i][j], 1, deriv_interp[j][i], 1,
407 out_interp[i][j + 1], 1);
410 out_interp[i][j + 1], 1);
415 Vmath::Vadd(nPts, out_interp[i][j + 1], 1, divVel, 1,
416 out_interp[i][j + 1], 1);
421 out_interp[j][i + 1] = out_interp[i][j + 1];
433 Vmath::Vvtvp(nPts, vel_interp[j], 1, out_interp[i][j + 1], 1,
448 m_fields[0]->PhysGalerkinProjection1DScaled(
449 OneDptscale, out_interp[i][j], viscousTensor[i][j]);
466 size_t nConvectiveFields = consvar.size();
468 int nengy = nConvectiveFields - 1;
475 int nLengthArray = 0;
478 int nBndRegions =
m_fields[nengy]->GetBndCondExpansions().size();
479 for (
int j = 0; j < nBndRegions; ++j)
482 ->GetBndConditions()[j]
489 m_fields[nengy]->GetBndCondExpansions()[j]->GetExpSize();
490 for (
int e = 0; e < nBndEdges; ++e)
492 size_t nBndEdgePts =
m_fields[nengy]
493 ->GetBndCondExpansions()[j]
497 int id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
498 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
502 m_fields[nengy]->GetBndConditions()[j]->GetUserDefined(),
505 if (nBndEdgePts != nLengthArray)
507 for (
int i = 0; i < nConvectiveFields; ++i)
515 nLengthArray = nBndEdgePts;
526 for (
int k = 0; k < nConvectiveFields; ++k)
532 m_varConv->GetPressure(bndCons, bndPressure);
534 m_varConv->GetRhoFromPT(bndPressure, bndTotEngy, bndRho);
535 m_varConv->GetEFromRhoP(bndRho, bndPressure, bndIntEndy);
536 m_varConv->GetDynamicEnergy(bndCons, bndTotEngy);
538 Vmath::Vvtvp(nBndEdgePts, bndIntEndy, 1, bndCons[ndens], 1,
539 bndTotEngy, 1, bndTotEngy, 1);
542 tmp = consvar[nengy] + id2, 1);
558 size_t nConvectiveFields = inarray.size();
559 size_t nPts = inaverg[nConvectiveFields - 1].size();
561 for (
int i = 0; i < nConvectiveFields - 1; ++i)
563 nonZeroIndex[i] = i + 1;
569 m_varConv->GetTemperature(inaverg, temperature);
572 std::vector<NekDouble> inAvgTmp(nConvectiveFields);
573 std::vector<NekDouble> inTmp(nConvectiveFields);
574 std::vector<NekDouble> outTmp(nConvectiveFields);
575 for (
int d = 0; d < nDim; ++d)
577 for (
int nderiv = 0; nderiv < nDim; ++nderiv)
579 for (
size_t p = 0;
p < nPts; ++
p)
582 for (
int f = 0; f < nConvectiveFields; ++f)
584 inAvgTmp[f] = inaverg[f][
p];
585 inTmp[f] = inarray[f][
p];
589 inAvgTmp.data(), inTmp.data(),
590 mu[
p], outTmp.data());
592 for (
int f = 0; f < nConvectiveFields; ++f)
594 outarray[d][f][
p] += normal[nderiv][
p] * outTmp[f];
609 unsigned int nTracePts = uFwd[0].size();
612 unsigned int nVariables = uFwd.size();
615 uBwd[nVariables - 1], 1, tAve, 1);
624 for (
int i = 0; i < nVariables; ++i)
627 Vmath::Vsub(nTracePts, uFwd[i], 1, uBwd[i], 1, penaltyCoeff[i], 1);
629 if (i < nVariables - 1)
631 Vmath::Vmul(nTracePts, muAve, 1, penaltyCoeff[i], 1,
636 Vmath::Vmul(nTracePts, tcAve, 1, penaltyCoeff[i], 1,
650 auto nPts = temperature.size();
652 for (
size_t p = 0;
p < nPts; ++
p)
662 auto nTracePts =
m_fields[0]->GetTrace()->GetTotPoints();
663 if (nPts != nTracePts)
698 auto nDim = fields[0]->GetCoordim(0);
699 auto nVar = cnsVar.size();
700 auto nPts = cnsVar[0].size();
701 auto nPtsTrc = cnsVarFwd[0].size();
705 primVarBwd(nVar - 1);
707 for (
unsigned short d = 0; d < nVar - 2; ++d)
713 auto ergLoc = nVar - 2;
721 for (
unsigned short d = 0; d < nVar - 2; ++d)
724 for (
size_t p = 0;
p < nPts; ++
p)
726 primVar[d][
p] = cnsVar[d + 1][
p] / cnsVar[0][
p];
729 for (
size_t p = 0;
p < nPtsTrc; ++
p)
731 primVarFwd[d][
p] = cnsVarFwd[d + 1][
p] / cnsVarFwd[0][
p];
732 primVarBwd[d][
p] = cnsVarBwd[d + 1][
p] / cnsVarBwd[0][
p];
738 for (
unsigned short j = 0; j < nDim; ++j)
741 for (
unsigned short i = 0; i < nVar - 1; ++i)
748 m_diffusion->DiffuseCalcDerivative(fields, primVar, primVarDer, primVarFwd,
763 auto nDim = pVarDer.size();
764 auto nPts = div.size();
767 for (
size_t p = 0;
p < nPts; ++
p)
770 for (
unsigned short j = 0; j < nDim; ++j)
772 divTmp += pVarDer[j][j][
p];
780 for (
size_t p = 0;
p < nPts; ++
p)
804 curl0sqr += curl1sqr + curl2sqr;
806 curlSquare[
p] = curl0sqr;
811 for (
size_t p = 0;
p < nPts; ++
p)
820 curlSquare[
p] = curl * curl;
831 std::vector<std::string> &variables)
834 m_session->MatchSolverInfo(
"OutputExtraFields",
"True", extraFields,
true);
837 const int nPhys =
m_fields[0]->GetNpoints();
838 const int nCoeffs =
m_fields[0]->GetNcoeffs();
841 for (
int i = 0; i <
m_fields.size(); ++i)
859 m_varConv->GetVelocityVector(cnsVar, velocity);
861 m_varConv->GetTemperature(cnsVar, temperature);
863 m_varConv->GetSoundSpeed(cnsVar, soundspeed);
864 m_varConv->GetMach(cnsVar, soundspeed, mach);
867 m_session->LoadParameter(
"SensorOffset", sensorOffset, 1);
876 string velNames[3] = {
"u",
"v",
"w"};
879 m_fields[0]->FwdTransLocalElmt(velocity[i], velFwd[i]);
880 variables.push_back(velNames[i]);
881 fieldcoeffs.push_back(velFwd[i]);
885 m_fields[0]->FwdTransLocalElmt(temperature, TFwd);
886 m_fields[0]->FwdTransLocalElmt(entropy, sFwd);
887 m_fields[0]->FwdTransLocalElmt(soundspeed, aFwd);
888 m_fields[0]->FwdTransLocalElmt(mach, mFwd);
889 m_fields[0]->FwdTransLocalElmt(sensor, sensFwd);
891 variables.push_back(
"p");
892 variables.push_back(
"T");
893 variables.push_back(
"s");
894 variables.push_back(
"a");
895 variables.push_back(
"Mach");
896 variables.push_back(
"Sensor");
897 fieldcoeffs.push_back(pFwd);
898 fieldcoeffs.push_back(TFwd);
899 fieldcoeffs.push_back(sFwd);
900 fieldcoeffs.push_back(aFwd);
901 fieldcoeffs.push_back(mFwd);
902 fieldcoeffs.push_back(sensFwd);
911 variables.push_back(
"ArtificialVisc");
912 fieldcoeffs.push_back(sensorFwd);
921 for (
int i = 0; i <
m_fields.size(); ++i)
925 m_fields[i]->GetFwdBwdTracePhys(cnsVar[i], cnsVarFwd[i],
934 m_fields[0]->FwdTransLocalElmt(div, divFwd);
935 variables.push_back(
"div");
936 fieldcoeffs.push_back(divFwd);
939 m_fields[0]->FwdTransLocalElmt(curlSquare, curlFwd);
940 variables.push_back(
"curl^2");
941 fieldcoeffs.push_back(curlFwd);
947 variables.push_back(
"ArtificialVisc");
948 fieldcoeffs.push_back(muavFwd);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
NekDouble m_bndEvaluateTime
std::string m_shockCaptureType
SolverUtils::DiffusionSharedPtr m_diffusion
VariableConverterSharedPtr m_varConv
void SetBoundaryConditionsBwdWeight()
Set up a weight on physical boundaries for boundary condition applications.
ArtificialDiffusionSharedPtr m_artificialDiffusion
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
virtual ~NavierStokesCFE()
virtual void v_GetViscousFluxVector(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &derivatives, TensorOfArray3D< NekDouble > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
void GetDivCurlFromDvelT(const TensorOfArray3D< NekDouble > &pVarDer, Array< OneD, NekDouble > &div, Array< OneD, NekDouble > &curlSquare)
Get divergence and curl from velocity derivative tensor.
virtual void v_InitObject(bool DeclareField=true) override
Initialization object for CompressibleFlowSystem class.
void SpecialBndTreat(Array< OneD, Array< OneD, NekDouble >> &consvar)
For very special treatment. For general boundaries it does nothing But for WallViscous and WallAdiaba...
void GetDivCurlSquared(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &cnsVar, Array< OneD, NekDouble > &div, Array< OneD, NekDouble > &curlSquare, const Array< OneD, Array< OneD, NekDouble >> &cnsVarFwd, const Array< OneD, Array< OneD, NekDouble >> &cnsVarBwd)
Get divergence and curl squared.
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) override
Apply artificial diffusion (Laplacian operator)
void InitObject_Explicit()
NekDouble m_thermalConductivityRef
bool m_is_mu_variable
flag to switch between constant viscosity and Sutherland an enum could be added for more options
bool m_is_shockCaptPhys
flag for shock capturing switch on/off an enum could be added for more options
virtual void v_GetViscousFluxVectorDeAlias(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &derivatives, TensorOfArray3D< NekDouble > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
void GetViscosityAndThermalCondFromTemp(const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &thermalCond)
Update viscosity todo: add artificial viscosity here.
bool m_is_diffIP
flag to switch between IP and LDG an enum could be added for more options
void GetViscousSymmtrFluxConservVar(const int nSpaceDim, const Array< OneD, Array< OneD, NekDouble >> &inaverg, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &outarray, Array< OneD, int > &nonZeroIndex, const Array< OneD, Array< OneD, NekDouble >> &normals)
Calculate and return the Symmetric flux in IP method.
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables) override
void GetViscousFluxBilinearFormKernel(const unsigned short nDim, const unsigned short FluxDirection, const unsigned short DerivDirection, const T *inaverg, const T *injumpp, const T &mu, T *outarray)
Calculate diffusion flux using the Jacobian form.
virtual void v_GetFluxPenalty(const Array< OneD, const Array< OneD, NekDouble >> &uFwd, const Array< OneD, const Array< OneD, NekDouble >> &uBwd, Array< OneD, Array< OneD, NekDouble >> &penaltyCoeff)
Return the penalty vector for the LDGNS diffusion problem.
void GetViscosityAndThermalCondFromTempKernel(const T &temperature, T &mu, T &thermalCond)
std::string m_ViscosityType
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetNpoints()
Base class for unsteady solvers.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
DiffusionFactory & GetDiffusionFactory()
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
The above copyright notice and this permission notice shall be included.
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
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 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 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 Zero(int n, T *x, const int incx)
Zero vector.
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)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.