45 "Riemann invariant boundary condition.");
51 const int pSpaceDim,
const int bcRegion,
const int cnt)
52 :
CFSBndCond(pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt)
55 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
73 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
80 NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
86 for (i = 0; i < nDimensions; ++i)
88 Vmath::Vdiv(nTracePts, Fwd[i + 1], 1, Fwd[0], 1, Vel, 1);
95 m_varConv->GetAbsoluteVelocity(Fwd, absVel);
102 m_varConv->GetSoundSpeed(Fwd, soundSpeed);
106 Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
111 int e, id1, id2, nBCEdgePts, pnt;
112 NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
120 for (e = 0; e < eMax; ++e)
133 for (i = 0; i < nBCEdgePts; i++)
141 if (Mach[pnt] < 1.00)
145 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
149 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
155 rPlus =
m_VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
159 rMinus =
m_VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
163 VNBC = 0.5 * (rPlus + rMinus);
164 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
169 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
170 pBC = rhoBC * cBC * cBC * gammaInv;
176 for (j = 0; j < nDimensions; ++j)
179 rhoVelBC[j] = rhoBC * velBC[j];
180 EkBC += 0.5 * rhoBC * velBC[j] * velBC[j];
184 EBC = pBC * gammaMinusOneInv + EkBC;
189 ->UpdatePhys())[id1 + i] = rhoBC;
190 for (j = 0; j < nDimensions; ++j)
194 ->UpdatePhys())[id1 + i] = rhoVelBC[j];
198 ->UpdatePhys())[id1 + i] = EBC;
203 if (Mach[pnt] < 1.00)
207 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
211 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
217 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
221 rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
225 VNBC = 0.5 * (rPlus + rMinus);
226 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
227 VDBC = VNBC - Vn[pnt];
231 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
232 pBC = rhoBC * cBC * cBC * gammaInv;
238 for (j = 0; j < nDimensions; ++j)
240 velBC[j] = Fwd[j + 1][pnt] / Fwd[0][pnt] +
242 rhoVelBC[j] = rhoBC * velBC[j];
243 EkBC += 0.5 * rhoBC * velBC[j] * velBC[j];
247 EBC = pBC * gammaMinusOneInv + EkBC;
252 ->UpdatePhys())[id1 + i] = rhoBC;
253 for (j = 0; j < nDimensions; ++j)
257 ->UpdatePhys())[id1 + i] = rhoVelBC[j];
261 ->UpdatePhys())[id1 + i] = EBC;
Encapsulates the user-defined boundary conditions for compressible flow solver.
int m_spacedim
Space dimension.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Trace normals.
Array< OneD, NekDouble > m_velInf
NekDouble m_gamma
Parameters of the flow.
int m_bcRegion
Id of the boundary region.
VariableConverterSharedPtr m_varConv
Auxiliary object to convert variables.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void v_Apply(Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time) override
static CFSBndCondSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const int pSpaceDim, const int bcRegion, const int cnt)
Creates an instance of this class.
static std::string className
Name of the class.
RiemannInvariantBC(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const int pSpaceDim, const int bcRegion, const int cnt)
Array< OneD, NekDouble > m_VnInf
Reference normal velocity.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |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 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.
scalarT< T > sqrt(scalarT< T > in)