44 RegisterCreatorFunction(
"RiemannInvariant",
45 RiemannInvariantBC::create,
46 "Riemann invariant boundary condition.");
48 RiemannInvariantBC::RiemannInvariantBC(
55 :
CFSBndCond(pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt)
58 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
78 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
86 NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
92 for (i = 0; i < nDimensions; ++i)
94 Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
101 m_varConv->GetAbsoluteVelocity(Fwd, absVel);
108 m_varConv->GetSoundSpeed(Fwd, pressure, soundSpeed);
112 Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
117 int e, id1, id2, nBCEdgePts, pnt;
118 NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
126 for (e = 0; e < eMax; ++e)
129 GetExp(e)->GetTotPoints();
136 for (i = 0; i < nBCEdgePts; i++)
144 if (Mach[pnt] < 1.00)
147 cPlus = sqrt(
m_gamma * pressure[pnt] / Fwd[0][pnt]);
148 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
152 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
158 rPlus =
m_VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
162 rMinus =
m_VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
166 VNBC = 0.5 * (rPlus + rMinus);
167 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
172 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
173 pBC = rhoBC * cBC * cBC * gammaInv;
179 for ( j = 0; j < nDimensions; ++j)
182 rhoVelBC[j] = rhoBC * velBC[j];
183 EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
187 EBC = pBC * gammaMinusOneInv + EkBC;
191 UpdatePhys())[id1+i] = rhoBC;
192 for (j = 0; j < nDimensions; ++j)
195 UpdatePhys())[id1+i] = rhoVelBC[j];
198 UpdatePhys())[id1+i] = EBC;
204 if (Mach[pnt] < 1.00)
207 cPlus = sqrt(
m_gamma * pressure[pnt] / Fwd[0][pnt]);
208 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
212 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
217 cPlus = sqrt(
m_gamma * pressure[pnt] / Fwd[0][pnt]);
218 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
221 cMinus = sqrt(
m_gamma * pressure[pnt] / Fwd[0][pnt]);
222 rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
226 VNBC = 0.5 * (rPlus + rMinus);
227 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
228 VDBC = VNBC - Vn[pnt];
231 sBC = pressure[pnt] / (pow(Fwd[0][pnt],
m_gamma));
232 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
233 pBC = rhoBC * cBC * cBC * gammaInv;
239 for ( j = 0; j < nDimensions; ++j)
241 velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
243 rhoVelBC[j] = rhoBC * velBC[j];
244 EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
248 EBC = pBC * gammaMinusOneInv + EkBC;
252 UpdatePhys())[id1+i] = rhoBC;
253 for (j = 0; j < nDimensions; ++j)
256 UpdatePhys())[id1+i] = rhoVelBC[j];
259 UpdatePhys())[id1+i] = EBC;
int m_spacedim
Space dimension.
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 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.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
virtual void v_Apply(Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time)
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
Array< OneD, NekDouble > m_VnInf
Reference normal velocity.
Encapsulates the user-defined boundary conditions for compressible flow solver.
NekDouble m_gamma
Parameters of the flow.
VariableConverterSharedPtr m_varConv
Auxiliary object to convert variables.
Array< OneD, NekDouble > m_velInf
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Trace normals.
int m_bcRegion
Id of the boundary region.