71{
72 int i, j;
73 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
75
76 const Array<OneD, const int> &traceBndMap =
m_fields[0]->GetTraceBndMap();
77
80 NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
81
82
83
84 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
85 Array<OneD, NekDouble> Vel(nTracePts, 0.0);
86 for (i = 0; i < nDimensions; ++i)
87 {
88 Vmath::Vdiv(nTracePts, Fwd[i + 1], 1, Fwd[0], 1, Vel, 1);
90 }
91
92
93
94 Array<OneD, NekDouble> absVel(nTracePts, 0.0);
95 m_varConv->GetAbsoluteVelocity(Fwd, absVel);
96
97
98 Array<OneD, NekDouble>
pressure(nTracePts);
99 Array<OneD, NekDouble> soundSpeed(nTracePts);
100
102 m_varConv->GetSoundSpeed(Fwd, soundSpeed);
103
104
105 Array<OneD, NekDouble> Mach(nTracePts, 0.0);
106 Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
108
109
110 int eMax;
111 int e, id1, id2, nBCEdgePts, pnt;
112 NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
113 Array<OneD, NekDouble> velBC(nDimensions, 0.0);
114 Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
116
118
119
120 for (e = 0; e < eMax; ++e)
121 {
124 ->GetExp(e)
125 ->GetTotPoints();
126
127 id1 =
129 id2 =
131
132
133 for (i = 0; i < nBCEdgePts; i++)
134 {
135 pnt = id2 + i;
136
137
138 if (Vn[pnt] <= 0.0)
139 {
140
141 if (Mach[pnt] < 1.00)
142 {
143
145 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
146
147
149 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
150 }
151 else
152 {
153
155 rPlus =
m_VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
156
157
159 rMinus =
m_VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
160 }
161
162
163 VNBC = 0.5 * (rPlus + rMinus);
164 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
166
167
169 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
170 pBC = rhoBC * cBC * cBC * gammaInv;
171
172
174
175
176 for (j = 0; j < nDimensions; ++j)
177 {
179 rhoVelBC[j] = rhoBC * velBC[j];
180 EkBC += 0.5 * rhoBC * velBC[j] * velBC[j];
181 }
182
183
184 EBC = pBC * gammaMinusOneInv + EkBC;
185
186
189 ->UpdatePhys())[id1 + i] = rhoBC;
190 for (j = 0; j < nDimensions; ++j)
191 {
194 ->UpdatePhys())[id1 + i] = rhoVelBC[j];
195 }
198 ->UpdatePhys())[id1 + i] = EBC;
199 }
200 else
201 {
202
203 if (Mach[pnt] < 1.00)
204 {
205
207 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
208
209
211 rMinus =
m_VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
212 }
213 else
214 {
215
217 rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
218
219
221 rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
222 }
223
224
225 VNBC = 0.5 * (rPlus + rMinus);
226 cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
227 VDBC = VNBC - Vn[pnt];
228
229
231 rhoBC = pow((cBC * cBC) / (
m_gamma * sBC), gammaMinusOneInv);
232 pBC = rhoBC * cBC * cBC * gammaInv;
233
234
236
237
238 for (j = 0; j < nDimensions; ++j)
239 {
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];
244 }
245
246
247 EBC = pBC * gammaMinusOneInv + EkBC;
248
249
252 ->UpdatePhys())[id1 + i] = rhoBC;
253 for (j = 0; j < nDimensions; ++j)
254 {
257 ->UpdatePhys())[id1 + i] = rhoVelBC[j];
258 }
261 ->UpdatePhys())[id1 + i] = EBC;
262 }
263 }
264 }
265}
NekDouble m_gamma
Parameters of the flow.
int m_bcRegion
Id of the boundary region.
VariableConverterSharedPtr m_varConv
Auxiliary object to convert variables.
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)