39#include <boost/algorithm/string/predicate.hpp>
61 std::size_t nDim = pFields[0]->GetCoordim(0);
62 std::size_t nTracePts = pFields[0]->GetTrace()->GetTotPoints();
65 if (pSession->DefinesSolverInfo(
"HOMOGENEOUS"))
83 m_session->LoadSolverInfo(
"EquationOfState", eosType,
"IdealGas");
92 std::size_t nElements = pFields[0]->GetExpSize();
94 for (std::size_t e = 0; e < nElements; e++)
97 std::size_t expDim = pFields[0]->GetShapeDimension();
104 for (std::size_t i = 0; i < exp3D->GetNtraces(); ++i)
107 h, exp3D->GetGeom3D()->GetEdge(i)->GetVertex(0)->dist(*(
108 exp3D->GetGeom3D()->GetEdge(i)->GetVertex(1))));
117 for (std::size_t i = 0; i < exp2D->GetNtraces(); ++i)
120 h, exp2D->GetGeom2D()->GetEdge(i)->GetVertex(0)->dist(*(
121 exp2D->GetGeom2D()->GetEdge(i)->GetVertex(1))));
130 h = std::min(h, exp1D->
GetGeom1D()->GetVertex(0)->dist(
131 *(exp1D->GetGeom1D()->GetVertex(1))));
136 ASSERTL0(
false,
"Dimension out of bound.")
144 std::size_t nPts = pFields[0]->GetTotPoints();
147 for (std::size_t e = 0; e < pFields[0]->GetExpSize(); e++)
149 std::size_t nElmtPoints = pFields[0]->GetExp(e)->GetTotPoints();
150 std::size_t physOffset = pFields[0]->GetPhys_Offset(e);
151 Vmath::Fill(nElmtPoints, hEle[e], tmp = hElePts + physOffset, 1);
156 pFields[0]->GetFwdBwdTracePhys(hElePts, Fwd, Bwd);
159 std::size_t nBndRegions = pFields[0]->GetBndCondExpansions().size();
161 for (std::size_t i = 0; i < nBndRegions; ++i)
164 std::size_t nBndEdges =
165 pFields[0]->GetBndCondExpansions()[i]->GetExpSize();
167 if (pFields[0]->GetBndConditions()[i]->GetBoundaryConditionType() ==
174 for (std::size_t e = 0; e < nBndEdges; ++e)
176 std::size_t nBndEdgePts = pFields[0]
177 ->GetBndCondExpansions()[i]
181 std::size_t id2 = pFields[0]->GetTrace()->GetPhys_Offset(
182 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
210 const std::size_t nConvectiveFields,
217 std::size_t nCoeffs = fields[0]->GetNcoeffs();
219 for (std::size_t i = 0; i < nConvectiveFields; ++i)
224 for (std::size_t i = 0; i < nConvectiveFields; ++i)
226 fields[i]->BwdTrans(tmp2[i], outarray[i]);
231 const std::size_t nConvectiveFields,
238 std::size_t nDim = fields[0]->GetCoordim(0);
239 std::size_t nPts = fields[0]->GetTotPoints();
240 std::size_t nCoeffs = fields[0]->GetNcoeffs();
241 std::size_t nScalars = inarray.size();
242 std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
250 for (std::size_t i = 0; i < nScalars; ++i)
265 for (std::size_t i = 0; i < nScalars + 1; ++i)
271 for (std::size_t i = 0; i < nConvectiveFields; ++i)
287 for (std::size_t i = 0; i < nConvectiveFields; ++i)
291 for (std::size_t j = 0; j < nDim; ++j)
295 fields[i]->IProductWRTDerivBase(tmpIn, tmpOut);
299 fields[i]->AddTraceIntegral(viscousFlux[i], tmpOut);
300 fields[i]->SetPhysState(
false);
301 fields[i]->MultiplyByElmtInvMass(tmpOut, outarray[i]);
312 std::size_t nDim = fields[0]->GetCoordim(0);
313 std::size_t nCoeffs = fields[0]->GetNcoeffs();
314 std::size_t nScalars = inarray.size();
315 std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
316 std::size_t nConvectiveFields = fields.size();
326 for (std::size_t i = 0; i < nScalars; ++i)
334 for (std::size_t j = 0; j < nDim; ++j)
336 for (std::size_t i = 0; i < nScalars; ++i)
338 fields[i]->IProductWRTDerivBase(j, inarray[i], tmp1);
340 fields[i]->AddTraceIntegral(numericalFluxO1[j][i], tmp1);
341 fields[i]->SetPhysState(
false);
342 fields[i]->MultiplyByElmtInvMass(tmp1, tmp1);
343 fields[i]->BwdTrans(tmp1, qfields[j][i]);
349 for (std::size_t i = 0; i < nScalars; ++i)
389 std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
390 std::size_t nScalars = inarray.size();
394 for (std::size_t i = 0; i < nScalars; ++i)
396 numflux[i] = {pFwd[i]};
400 if (fields[0]->GetBndCondExpansions().size())
402 ApplyBCsO1(fields, inarray, pFwd, pBwd, numflux);
408 for (std::size_t i = 0; i < nScalars; ++i)
411 numericalFluxO1[j][i], 1);
427 std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
428 std::size_t nScalars = inarray.size();
437 for (std::size_t i = 0; i < nScalars; ++i)
443 for (std::size_t i = 0; i < nScalars - 1; ++i)
448 std::size_t nBndRegions = fields[i + 1]->GetBndCondExpansions().size();
449 for (std::size_t j = 0; j < nBndRegions; ++j)
452 ->GetBndConditions()[j]
458 std::size_t nBndEdges =
459 fields[i + 1]->GetBndCondExpansions()[j]->GetExpSize();
460 for (std::size_t e = 0; e < nBndEdges; ++e)
462 std::size_t nBndEdgePts = fields[i + 1]
463 ->GetBndCondExpansions()[j]
468 fields[i + 1]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
470 std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
471 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
475 fields[i]->GetBndConditions()[j]->GetUserDefined(),
478 fields[i]->GetBndConditions()[j]->GetUserDefined(),
482 Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
484 else if (boost::iequals(
485 fields[i]->GetBndConditions()[j]->GetUserDefined(),
488 fields[i]->GetBndConditions()[j]->GetUserDefined(),
497 for (std::size_t k = 0; k < nScalars - 1; ++k)
501 ->GetBndCondExpansions()[j]
502 ->UpdatePhys())[id1],
505 ->GetBndCondExpansions()[j]
506 ->UpdatePhys())[id1],
507 1, &scalarVariables[k][id2], 1);
509 1, &scalarVariables[k][id2], 1,
510 &tmp1[0], 1, &tmp1[0], 1);
512 Vmath::Smul(nBndEdgePts, -1.0, &tmp1[0], 1, &tmp1[0],
516 for (std::size_t k = 0; k < nScalars - 1; ++k)
520 &scalarVariables[k][id2], 1,
521 &scalarVariables[k][id2], 1);
526 ->GetBndConditions()[j]
527 ->GetBoundaryConditionType() ==
533 ->GetBndCondExpansions()[j]
534 ->UpdatePhys())[id1],
537 ->GetBndCondExpansions()[j]
538 ->UpdatePhys())[id1],
539 1, &scalarVariables[i][id2], 1);
544 ->GetBndConditions()[j]
545 ->GetBoundaryConditionType() ==
553 else if ((fields[i]->GetBndConditions()[j])
554 ->GetBoundaryConditionType() ==
557 Vmath::Vcopy(nBndEdgePts, &pFwd[i][id2], 1, &fluxO1[i][id2],
562 Vmath::Vmul(nBndEdgePts, &scalarVariables[i][id2], 1,
563 &scalarVariables[i][id2], 1, &tmp1[id2], 1);
565 Vmath::Smul(nBndEdgePts, 0.5, &tmp1[id2], 1, &tmp1[id2], 1);
567 Vmath::Vadd(nBndEdgePts, &tmp2[id2], 1, &tmp1[id2], 1,
575 std::size_t nBndRegions = fields[nScalars]->GetBndCondExpansions().size();
576 for (std::size_t j = 0; j < nBndRegions; ++j)
579 ->GetBndConditions()[j]
585 std::size_t nBndEdges =
586 fields[nScalars]->GetBndCondExpansions()[j]->GetExpSize();
587 for (std::size_t e = 0; e < nBndEdges; ++e)
589 std::size_t nBndEdgePts = fields[nScalars]
590 ->GetBndCondExpansions()[j]
595 fields[nScalars]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
597 std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
598 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
602 fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
606 &scalarVariables[nScalars - 1][id2], 1);
610 else if (fields[nScalars]
611 ->GetBndConditions()[j]
612 ->GetBoundaryConditionType() ==
617 for (std::size_t n = 0; n < nBndEdgePts; ++n)
620 ->GetBndCondExpansions()[j]
621 ->GetPhys()[id1 + n];
622 ene = fields[nScalars]
623 ->GetBndCondExpansions()[j]
624 ->GetPhys()[id1 + n] /
627 scalarVariables[nScalars - 1][id2 + n] =
628 m_eos->GetTemperature(rho, ene);
634 ->GetBndConditions()[j]
635 ->GetBoundaryConditionType() ==
638 fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
641 Vmath::Vcopy(nBndEdgePts, &scalarVariables[nScalars - 1][id2],
642 1, &fluxO1[nScalars - 1][id2], 1);
646 else if (((fields[nScalars]->GetBndConditions()[j])
647 ->GetBoundaryConditionType() ==
649 boost::iequals(fields[nScalars]
650 ->GetBndConditions()[j]
655 &fluxO1[nScalars - 1][id2], 1);
672 std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
673 std::size_t nVariables = fields.size();
677 for (std::size_t i = 0; i < nVariables - 1; ++i)
693 std::size_t nDim = fields[0]->GetCoordim(0);
694 for (std::size_t i = 1; i < nVariables; ++i)
696 for (std::size_t j = 0; j < nDim; ++j)
699 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
710 qfluxtemp, 1, qfluxtemp, 1);
713 if (fields[0]->GetBndCondExpansions().size())
715 ApplyBCsO2(fields, i, j, qfield[j][i], qFwd, qBwd, qfluxtemp);
719 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
730 const std::size_t var,
const std::size_t dir,
737 std::size_t nBndRegions = fields[var]->GetBndCondExpansions().size();
739 for (std::size_t i = 0; i < nBndRegions; ++i)
742 std::size_t nBndEdges =
743 fields[var]->GetBndCondExpansions()[i]->GetExpSize();
745 if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
752 for (std::size_t e = 0; e < nBndEdges; ++e)
754 std::size_t nBndEdgePts = fields[var]
755 ->GetBndCondExpansions()[i]
759 std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
760 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
765 ->GetBndConditions()[i]
766 ->GetBoundaryConditionType() ==
769 fields[var]->GetBndConditions()[i]->GetUserDefined(),
773 &qFwd[id2], 1, &penaltyflux[id2], 1);
777 else if ((fields[var]->GetBndConditions()[i])
778 ->GetBoundaryConditionType() ==
781 ASSERTL0(
false,
"Neumann bcs not implemented for LDGNS");
783 else if (boost::iequals(
784 fields[var]->GetBndConditions()[i]->GetUserDefined(),
795 &qFwd[id2], 1, &penaltyflux[id2], 1);
#define ASSERTL0(condition, msg)
void v_DiffuseCoeffs(const std::size_t nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, 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
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
void v_Diffuse(const std::size_t nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, 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
Calculate weak DG Diffusion in the LDG form for the Navier-Stokes (NS) equations:
EquationOfStateSharedPtr m_eos
Equation of system for computing temperature.
void NumericalFluxO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &numericalFluxO1, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
Builds the numerical flux for the 1st order derivatives.
LibUtilities::SessionReaderSharedPtr m_session
Array< OneD, Array< OneD, NekDouble > > m_traceVel
void NumericalFluxO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, TensorOfArray3D< NekDouble > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
Build the numerical flux for the 2nd order derivatives.
NekDouble m_C11
Penalty coefficient for LDGNS.
void v_DiffuseTraceFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &TraceFlux, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd, Array< OneD, int > &nonZeroIndex) override
Diffusion term Trace Flux.
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
void v_DiffuseCalcDerivative(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
Diffusion Flux, calculate the physical derivatives.
TensorOfArray3D< NekDouble > m_viscTensor
void ApplyBCsO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const std::size_t var, const std::size_t dir, const Array< OneD, const NekDouble > &qfield, const Array< OneD, const NekDouble > &qFwd, const Array< OneD, const NekDouble > &qBwd, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
void v_DiffuseVolumeFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex) override
Diffusion Volume Flux.
Array< OneD, NekDouble > m_traceOneOverH
h scaling for penalty term
void ApplyBCsO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd, Array< OneD, Array< OneD, NekDouble > > &flux01)
Imposes appropriate bcs for the 1st order derivatives.
static DiffusionSharedPtr create(std::string diffType)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
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.
SpatialDomains::Geometry1DSharedPtr GetGeom1D() const
SOLVER_UTILS_EXPORT void DiffuseCalcDerivative(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
SOLVER_UTILS_EXPORT void DiffuseVolumeFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
Diffusion Volume FLux.
DiffusionFluxVecCBNS m_fluxVectorNS
SOLVER_UTILS_EXPORT void DiffuseTraceFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &TraceFlux, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
Diffusion term Trace Flux.
DiffusionFluxPenaltyNS m_fluxPenaltyNS
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
DiffusionFactory & GetDiffusionFactory()
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (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 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 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 Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
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.
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)