www.mooseframework.org
CNSFVRiemannInvariantBoundaryFlux.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<BoundaryFluxBase>();
15 
16  params.addClassDescription("A user objec that computes the Riemann-invariant boundary flux.");
17 
18  params.addRequiredParam<UserObjectName>("bc_uo", "Name for boundary condition user object");
19 
20  params.addRequiredParam<UserObjectName>("fluid_properties",
21  "Name for fluid properties user object");
22 
23  return params;
24 }
25 
27  const InputParameters & parameters)
28  : BoundaryFluxBase(parameters),
29  _bc_uo(getUserObject<BCUserObject>("bc_uo")),
30  _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties"))
31 {
32 }
33 
35 
36 void
38  dof_id_type ielem,
39  const std::vector<Real> & uvec1,
40  const RealVectorValue & dwave,
41  std::vector<Real> & flux) const
42 {
44 
45  Real rho1 = uvec1[0];
46  Real rhou1 = uvec1[1];
47  Real rhov1 = uvec1[2];
48  Real rhow1 = uvec1[3];
49  Real rhoe1 = uvec1[4];
50 
51  Real nx = dwave(0);
52  Real ny = dwave(1);
53  Real nz = dwave(2);
54 
56 
57  flux.resize(5);
58 
60 
61  Real rhom1 = 1. / rho1;
62  Real uadv1 = rhou1 * rhom1;
63  Real vadv1 = rhov1 * rhom1;
64  Real wadv1 = rhow1 * rhom1;
65  Real vdov1 = uadv1 * uadv1 + vadv1 * vadv1 + wadv1 * wadv1;
66  Real eint1 = rhoe1 * rhom1 - 0.5 * vdov1;
67  Real pres1 = _fp.pressure(rhom1, eint1);
68  Real csou1 = _fp.c(rhom1, eint1);
69  Real gamma = _fp.gamma(rhom1, eint1);
70  Real gamm1 = gamma - 1.;
71  Real mach1 = std::sqrt(vdov1) / csou1;
72 
74 
75  if (std::abs(mach1) < 1.)
76  {
78 
79  std::vector<Real> U2(5, 0.);
80 
81  U2 = _bc_uo.getGhostCellValue(iside, ielem, uvec1, dwave);
82 
83  Real rho2 = U2[0];
84  Real rhou2 = U2[1];
85  Real rhov2 = U2[2];
86  Real rhow2 = U2[3];
87  Real rhoe2 = U2[4];
88 
89  Real uadv2 = rhou2 / rho2;
90  Real vadv2 = rhov2 / rho2;
91  Real wadv2 = rhow2 / rho2;
92  Real vdov2 = uadv2 * uadv2 + vadv2 * vadv2 + wadv2 * wadv2;
93  Real v2 = 1. / rho2;
94  Real e2 = rhoe2 / rho2 - 0.5 * vdov2;
95  Real pres2 = _fp.pressure(v2, e2);
96 
97  Real eint2 = _fp.e(pres2, rho2);
98  Real csou2 = _fp.c(1. / rho2, eint2);
99 
100  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
101  Real vdon2 = uadv2 * nx + vadv2 * ny + wadv2 * nz;
102 
103  Real rplus = vdon1 + 2. * csou1 / gamm1;
104  Real rmins = vdon2 - 2. * csou2 / gamm1;
105 
106  Real velob = (rplus + rmins) * 0.5;
107  Real csoub = (rplus - rmins) * 0.25 * gamm1;
108 
109  Real vdiff, uadvb, vadvb, wadvb, entrb;
110 
111  if (mach1 <= 0.)
112  {
114 
115  vdiff = velob - vdon2;
116  uadvb = uadv2 + vdiff * nx;
117  vadvb = vadv2 + vdiff * ny;
118  wadvb = wadv2 + vdiff * nz;
119  entrb = csou2 * csou2 / gamma / std::pow(rho2, gamm1);
120  }
121  else
122  {
124 
125  vdiff = velob - vdon1;
126  uadvb = uadv1 + vdiff * nx;
127  vadvb = vadv1 + vdiff * ny;
128  wadvb = wadv1 + vdiff * nz;
129  entrb = csou1 * csou1 / gamma / std::pow(rho1, gamm1);
130  }
131 
132  Real rhob = std::pow(csoub * csoub / gamma / entrb, 1. / gamm1);
133 
134  Real presb = rhob * csoub * csoub / gamma;
135 
136  Real vdonb = uadvb * nx + vadvb * ny + wadvb * nz;
137 
138  Real vdovb = uadvb * uadvb + vadvb * vadvb + wadvb * wadvb;
139 
140  Real rhoeb = rhob * _fp.e(presb, rhob) + 0.5 * rhob * vdovb;
141 
142  flux[0] = vdonb * rhob;
143  flux[1] = vdonb * rhob * uadvb + presb * nx;
144  flux[2] = vdonb * rhob * vadvb + presb * ny;
145  flux[3] = vdonb * rhob * wadvb + presb * nz;
146  flux[4] = vdonb * (rhoeb + presb);
147  }
148  else if (mach1 <= -1.)
149  {
151 
152  std::vector<Real> U2(5, 0.);
153 
154  U2 = _bc_uo.getGhostCellValue(iside, ielem, uvec1, dwave);
155 
156  Real rho2 = U2[0];
157  Real rhou2 = U2[1];
158  Real rhov2 = U2[2];
159  Real rhow2 = U2[3];
160  Real rhoe2 = U2[4];
161 
162  Real uadv2 = rhou2 / rho2;
163  Real vadv2 = rhov2 / rho2;
164  Real wadv2 = rhow2 / rho2;
165  Real vdov2 = uadv2 * uadv2 + vadv2 * vadv2 + wadv2 * wadv2;
166  Real v2 = 1. / rho2;
167  Real e2 = rhoe2 / rho2 - 0.5 * vdov2;
168  Real pres2 = _fp.pressure(v2, e2);
169 
170  Real vdon2 = uadv2 * nx + vadv2 * ny + wadv2 * nz;
171 
172  flux[0] = vdon2 * rho2;
173  flux[1] = vdon2 * rho2 * uadv2 + pres2 * nx;
174  flux[2] = vdon2 * rho2 * vadv2 + pres2 * ny;
175  flux[3] = vdon2 * rho2 * wadv2 + pres2 * nz;
176  flux[4] = vdon2 * (rhoe2 + pres2);
177  }
178  else if (mach1 >= 1.)
179  {
181 
182  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
183 
184  flux[0] = vdon1 * rho1;
185  flux[1] = vdon1 * rhou1 + pres1 * nx;
186  flux[2] = vdon1 * rhov1 + pres1 * ny;
187  flux[3] = vdon1 * rhow1 + pres1 * nz;
188  flux[4] = vdon1 * (rhoe1 + pres1);
189  }
190  else
191  mooseError("Something is wrong in ",
192  name(),
193  ": ",
194  __FUNCTION__,
195  "\n",
196  "ielem = ",
197  ielem,
198  "\n",
199  "iside = ",
200  iside,
201  "\n",
202  "mach1 = ",
203  mach1,
204  "\n");
205 }
206 
207 void
209  dof_id_type ielem,
210  const std::vector<Real> & uvec1,
211  const RealVectorValue & dwave,
212  DenseMatrix<Real> & jac1) const
213 {
215 
216  Real rho1 = uvec1[0];
217  Real rhou1 = uvec1[1];
218  Real rhov1 = uvec1[2];
219  Real rhow1 = uvec1[3];
220  Real rhoe1 = uvec1[4];
221 
222  Real nx = dwave(0);
223  Real ny = dwave(1);
224  Real nz = dwave(2);
225 
227 
228  jac1.resize(5, 5);
229 
231 
232  Real rhom1 = 1. / rho1;
233  Real uadv1 = rhou1 * rhom1;
234  Real vadv1 = rhov1 * rhom1;
235  Real wadv1 = rhow1 * rhom1;
236  Real vdov1 = uadv1 * uadv1 + vadv1 * vadv1 + wadv1 * wadv1;
237  Real eint1 = rhoe1 * rhom1 - 0.5 * vdov1;
238  Real pres1 = _fp.pressure(rhom1, eint1);
239  Real csou1 = _fp.c(rhom1, eint1);
240  Real gamma = _fp.gamma(rhom1, eint1);
241  Real gamm1 = gamma - 1.;
242  Real gamm2 = 2. - gamma;
243  Real mach1 = std::sqrt(vdov1) / csou1;
244 
246 
247  if (std::abs(mach1) < 1.)
248  {
250 
251  std::vector<Real> U2(5, 0.);
252 
253  U2 = _bc_uo.getGhostCellValue(iside, ielem, uvec1, dwave);
254 
255  Real rho2 = U2[0];
256  Real rhou2 = U2[1];
257  Real rhov2 = U2[2];
258  Real rhow2 = U2[3];
259  Real rhoe2 = U2[4];
260 
261  Real uadv2 = rhou2 / rho2;
262  Real vadv2 = rhov2 / rho2;
263  Real wadv2 = rhow2 / rho2;
264  Real vdov2 = uadv2 * uadv2 + vadv2 * vadv2 + wadv2 * wadv2;
265  Real v2 = 1. / rho2;
266  Real e2 = rhoe2 / rho2 - 0.5 * vdov2;
267  Real pres2 = _fp.pressure(v2, e2);
268 
269  Real eint2 = _fp.e(pres2, rho2);
270  Real csou2 = _fp.c(1. / rho2, eint2);
271 
272  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
273  Real vdon2 = uadv2 * nx + vadv2 * ny + wadv2 * nz;
274 
275  Real rplus = vdon1 + 2. * csou1 / gamm1;
276  Real rmins = vdon2 - 2. * csou2 / gamm1;
277 
278  Real velob = (rplus + rmins) * 0.5;
279  Real csoub = (rplus - rmins) * 0.25 * gamm1;
280 
281  Real vdiff, uadvb, vadvb, wadvb, entrb;
282 
283  if (mach1 <= 0.)
284  {
286 
287  vdiff = velob - vdon2;
288  uadvb = uadv2 + vdiff * nx;
289  vadvb = vadv2 + vdiff * ny;
290  wadvb = wadv2 + vdiff * nz;
291  entrb = csou2 * csou2 / gamma / std::pow(rho2, gamm1);
292  }
293  else
294  {
296 
297  vdiff = velob - vdon1;
298  uadvb = uadv1 + vdiff * nx;
299  vadvb = vadv1 + vdiff * ny;
300  wadvb = wadv1 + vdiff * nz;
301  entrb = csou1 * csou1 / gamma / std::pow(rho1, gamm1);
302  }
303 
304  Real rhob = std::pow(csoub * csoub / gamma / entrb, 1. / gamm1);
305 
306  Real presb = rhob * csoub * csoub / gamma;
307 
308  Real vdonb = uadvb * nx + vadvb * ny + wadvb * nz;
309 
310  Real vdovb = uadvb * uadvb + vadvb * vadvb + wadvb * wadvb;
311 
312  Real rhoeb = rhob * _fp.e(presb, rhob) + 0.5 * rhob * vdovb;
313 
314  Real rhoub = rhob * uadvb;
315  Real rhovb = rhob * vadvb;
316  Real rhowb = rhob * wadvb;
317 
318  Real tenthb = rhoeb + presb;
319 
321 
322  Real fmass_rhob = vdonb;
323  Real fmass_uadvb = nx * rhob;
324  Real fmass_vadvb = ny * rhob;
325  Real fmass_wadvb = nz * rhob;
326  Real fmass_presb = 0.;
327 
328  Real fmomx_rhob = vdonb * uadvb;
329  Real fmomx_uadvb = nx * rhoub * 2.;
330  Real fmomx_vadvb = ny * rhoub;
331  Real fmomx_wadvb = nz * rhoub;
332  Real fmomx_presb = nx;
333 
334  Real fmomy_rhob = vdonb * vadvb;
335  Real fmomy_uadvb = nx * rhovb;
336  Real fmomy_vadvb = ny * rhovb * 2.;
337  Real fmomy_wadvb = nz * rhovb;
338  Real fmomy_presb = ny;
339 
340  Real fmomz_rhob = vdonb * wadvb;
341  Real fmomz_uadvb = nx * rhowb;
342  Real fmomz_vadvb = ny * rhowb;
343  Real fmomz_wadvb = nz * rhowb * 2.;
344  Real fmomz_presb = nz;
345 
346  Real fetot_rhob = 0.5 * vdonb * vdovb;
347  Real fetot_uadvb = nx * tenthb + vdonb * rhoub;
348  Real fetot_vadvb = ny * tenthb + vdonb * rhovb;
349  Real fetot_wadvb = nz * tenthb + vdonb * rhowb;
350  Real fetot_presb = gamma / gamm1 * vdonb;
351 
353 
354  Real rho1_rho1 = 1.;
355  Real rho1_rhou1 = 0.;
356  Real rho1_rhov1 = 0.;
357  Real rho1_rhow1 = 0.;
358  Real rho1_rhoe1 = 0.;
359 
360  Real uadv1_rho1 = -rhom1 * uadv1;
361  Real uadv1_rhou1 = rhom1;
362  Real uadv1_rhov1 = 0.;
363  Real uadv1_rhow1 = 0.;
364  Real uadv1_rhoe1 = 0.;
365 
366  Real vadv1_rho1 = -rhom1 * vadv1;
367  Real vadv1_rhou1 = 0.;
368  Real vadv1_rhov1 = rhom1;
369  Real vadv1_rhow1 = 0.;
370  Real vadv1_rhoe1 = 0.;
371 
372  Real wadv1_rho1 = -rhom1 * wadv1;
373  Real wadv1_rhou1 = 0.;
374  Real wadv1_rhov1 = 0.;
375  Real wadv1_rhow1 = rhom1;
376  Real wadv1_rhoe1 = 0.;
377 
378  Real pres1_rho1 = gamm1 * vdov1 * 0.5;
379  Real pres1_rhou1 = -gamm1 * uadv1;
380  Real pres1_rhov1 = -gamm1 * vadv1;
381  Real pres1_rhow1 = -gamm1 * wadv1;
382  Real pres1_rhoe1 = gamm1;
383 
385 
386  Real vdon1_rho1 = -rhom1 * vdon1;
387  Real vdon1_rhou1 = rhom1 * nx;
388  Real vdon1_rhov1 = rhom1 * ny;
389  Real vdon1_rhow1 = rhom1 * nz;
390  Real vdon1_rhoe1 = 0.;
391 
393 
394  Real cp051 = 0.5 * csou1 / pres1;
395  Real prho1 = pres1 * rhom1;
396 
397  Real csou1_rho1 = cp051 * (pres1_rho1 - prho1 * rho1_rho1);
398  Real csou1_rhou1 = cp051 * (pres1_rhou1 - prho1 * rho1_rhou1);
399  Real csou1_rhov1 = cp051 * (pres1_rhov1 - prho1 * rho1_rhov1);
400  Real csou1_rhow1 = cp051 * (pres1_rhow1 - prho1 * rho1_rhow1);
401  Real csou1_rhoe1 = cp051 * (pres1_rhoe1 - prho1 * rho1_rhoe1);
402 
404 
405  Real Rplus_rho1 = vdon1_rho1 + 2. / gamm1 * csou1_rho1;
406  Real Rplus_rhou1 = vdon1_rhou1 + 2. / gamm1 * csou1_rhou1;
407  Real Rplus_rhov1 = vdon1_rhov1 + 2. / gamm1 * csou1_rhov1;
408  Real Rplus_rhow1 = vdon1_rhow1 + 2. / gamm1 * csou1_rhow1;
409  Real Rplus_rhoe1 = vdon1_rhoe1 + 2. / gamm1 * csou1_rhoe1;
410 
412 
413  Real csoub_rho1 = 0.25 * gamm1 * Rplus_rho1;
414  Real csoub_rhou1 = 0.25 * gamm1 * Rplus_rhou1;
415  Real csoub_rhov1 = 0.25 * gamm1 * Rplus_rhov1;
416  Real csoub_rhow1 = 0.25 * gamm1 * Rplus_rhow1;
417  Real csoub_rhoe1 = 0.25 * gamm1 * Rplus_rhoe1;
418 
420 
421  Real rhob_rho1, rhob_rhou1, rhob_rhov1, rhob_rhow1, rhob_rhoe1;
422  Real uadvb_rho1, uadvb_rhou1, uadvb_rhov1, uadvb_rhow1, uadvb_rhoe1;
423  Real vadvb_rho1, vadvb_rhou1, vadvb_rhov1, vadvb_rhow1, vadvb_rhoe1;
424  Real wadvb_rho1, wadvb_rhou1, wadvb_rhov1, wadvb_rhow1, wadvb_rhoe1;
425  Real presb_rho1, presb_rhou1, presb_rhov1, presb_rhow1, presb_rhoe1;
426 
427  if (mach1 <= 0.)
428  {
430 
432 
433  Real coeff = rhob / gamm1 * 2. / csoub;
434 
435  rhob_rho1 = coeff * csoub_rho1;
436  rhob_rhou1 = coeff * csoub_rhou1;
437  rhob_rhov1 = coeff * csoub_rhov1;
438  rhob_rhow1 = coeff * csoub_rhow1;
439  rhob_rhoe1 = coeff * csoub_rhoe1;
440 
441  uadvb_rho1 = nx * 0.5 * Rplus_rho1;
442  uadvb_rhou1 = nx * 0.5 * Rplus_rhou1;
443  uadvb_rhov1 = nx * 0.5 * Rplus_rhov1;
444  uadvb_rhow1 = nx * 0.5 * Rplus_rhow1;
445  uadvb_rhoe1 = nx * 0.5 * Rplus_rhoe1;
446 
447  vadvb_rho1 = ny * 0.5 * Rplus_rho1;
448  vadvb_rhou1 = ny * 0.5 * Rplus_rhou1;
449  vadvb_rhov1 = ny * 0.5 * Rplus_rhov1;
450  vadvb_rhow1 = ny * 0.5 * Rplus_rhow1;
451  vadvb_rhoe1 = ny * 0.5 * Rplus_rhoe1;
452 
453  wadvb_rho1 = nz * 0.5 * Rplus_rho1;
454  wadvb_rhou1 = nz * 0.5 * Rplus_rhou1;
455  wadvb_rhov1 = nz * 0.5 * Rplus_rhov1;
456  wadvb_rhow1 = nz * 0.5 * Rplus_rhow1;
457  wadvb_rhoe1 = nz * 0.5 * Rplus_rhoe1;
458 
459  presb_rho1 = presb * (1. / rhob * rhob_rho1 + 2. / csoub * csoub_rho1);
460  presb_rhou1 = presb * (1. / rhob * rhob_rhou1 + 2. / csoub * csoub_rhou1);
461  presb_rhov1 = presb * (1. / rhob * rhob_rhov1 + 2. / csoub * csoub_rhov1);
462  presb_rhow1 = presb * (1. / rhob * rhob_rhow1 + 2. / csoub * csoub_rhow1);
463  presb_rhoe1 = presb * (1. / rhob * rhob_rhoe1 + 2. / csoub * csoub_rhoe1);
464  }
465  else
466  {
468 
470 
471  Real coeff = csou1 / gamma / std::pow(rho1, gamma);
472 
473  Real entrb_rho1 = coeff * (2. * rho1 * csou1_rho1 - gamm1 * csou1 * rho1_rho1);
474  Real entrb_rhou1 = coeff * (2. * rho1 * csou1_rhou1 - gamm1 * csou1 * rho1_rhou1);
475  Real entrb_rhov1 = coeff * (2. * rho1 * csou1_rhov1 - gamm1 * csou1 * rho1_rhov1);
476  Real entrb_rhow1 = coeff * (2. * rho1 * csou1_rhow1 - gamm1 * csou1 * rho1_rhow1);
477  Real entrb_rhoe1 = coeff * (2. * rho1 * csou1_rhoe1 - gamm1 * csou1 * rho1_rhoe1);
478 
480 
481  rhob_rho1 = rhob / gamm1 * (2. / csoub * csoub_rho1 - 1. / entrb * entrb_rho1);
482  rhob_rhou1 = rhob / gamm1 * (2. / csoub * csoub_rhou1 - 1. / entrb * entrb_rhou1);
483  rhob_rhov1 = rhob / gamm1 * (2. / csoub * csoub_rhov1 - 1. / entrb * entrb_rhov1);
484  rhob_rhow1 = rhob / gamm1 * (2. / csoub * csoub_rhow1 - 1. / entrb * entrb_rhow1);
485  rhob_rhoe1 = rhob / gamm1 * (2. / csoub * csoub_rhoe1 - 1. / entrb * entrb_rhoe1);
486 
487  uadvb_rho1 = nx * (0.5 * Rplus_rho1 - vdon1_rho1) + uadv1_rho1;
488  uadvb_rhou1 = nx * (0.5 * Rplus_rhou1 - vdon1_rhou1) + uadv1_rhou1;
489  uadvb_rhov1 = nx * (0.5 * Rplus_rhov1 - vdon1_rhov1) + uadv1_rhov1;
490  uadvb_rhow1 = nx * (0.5 * Rplus_rhow1 - vdon1_rhow1) + uadv1_rhow1;
491  uadvb_rhoe1 = nx * (0.5 * Rplus_rhoe1 - vdon1_rhoe1) + uadv1_rhoe1;
492 
493  vadvb_rho1 = ny * (0.5 * Rplus_rho1 - vdon1_rho1) + vadv1_rho1;
494  vadvb_rhou1 = ny * (0.5 * Rplus_rhou1 - vdon1_rhou1) + vadv1_rhou1;
495  vadvb_rhov1 = ny * (0.5 * Rplus_rhov1 - vdon1_rhov1) + vadv1_rhov1;
496  vadvb_rhow1 = ny * (0.5 * Rplus_rhow1 - vdon1_rhow1) + vadv1_rhow1;
497  vadvb_rhoe1 = ny * (0.5 * Rplus_rhoe1 - vdon1_rhoe1) + vadv1_rhoe1;
498 
499  wadvb_rho1 = nz * (0.5 * Rplus_rho1 - vdon1_rho1) + wadv1_rho1;
500  wadvb_rhou1 = nz * (0.5 * Rplus_rhou1 - vdon1_rhou1) + wadv1_rhou1;
501  wadvb_rhov1 = nz * (0.5 * Rplus_rhov1 - vdon1_rhov1) + wadv1_rhov1;
502  wadvb_rhow1 = nz * (0.5 * Rplus_rhow1 - vdon1_rhow1) + wadv1_rhow1;
503  wadvb_rhoe1 = nz * (0.5 * Rplus_rhoe1 - vdon1_rhoe1) + wadv1_rhoe1;
504 
505  presb_rho1 = presb * (1. / rhob * rhob_rho1 + 2. / csoub * csoub_rho1);
506  presb_rhou1 = presb * (1. / rhob * rhob_rhou1 + 2. / csoub * csoub_rhou1);
507  presb_rhov1 = presb * (1. / rhob * rhob_rhov1 + 2. / csoub * csoub_rhov1);
508  presb_rhow1 = presb * (1. / rhob * rhob_rhow1 + 2. / csoub * csoub_rhow1);
509  presb_rhoe1 = presb * (1. / rhob * rhob_rhoe1 + 2. / csoub * csoub_rhoe1);
510  }
511 
513 
514  jac1(0, 0) = fmass_rhob * rhob_rho1 + fmass_uadvb * uadvb_rho1 + fmass_vadvb * vadvb_rho1 +
515  fmass_wadvb * wadvb_rho1 + fmass_presb * presb_rho1;
516 
517  jac1(0, 1) = fmass_rhob * rhob_rhou1 + fmass_uadvb * uadvb_rhou1 + fmass_vadvb * vadvb_rhou1 +
518  fmass_wadvb * wadvb_rhou1 + fmass_presb * presb_rhou1;
519 
520  jac1(0, 2) = fmass_rhob * rhob_rhov1 + fmass_uadvb * uadvb_rhov1 + fmass_vadvb * vadvb_rhov1 +
521  fmass_wadvb * wadvb_rhov1 + fmass_presb * presb_rhov1;
522 
523  jac1(0, 3) = fmass_rhob * rhob_rhow1 + fmass_uadvb * uadvb_rhow1 + fmass_vadvb * vadvb_rhow1 +
524  fmass_wadvb * wadvb_rhow1 + fmass_presb * presb_rhow1;
525 
526  jac1(0, 4) = fmass_rhob * rhob_rhoe1 + fmass_uadvb * uadvb_rhoe1 + fmass_vadvb * vadvb_rhoe1 +
527  fmass_wadvb * wadvb_rhoe1 + fmass_presb * presb_rhoe1;
528 
530 
531  jac1(1, 0) = fmomx_rhob * rhob_rho1 + fmomx_uadvb * uadvb_rho1 + fmomx_vadvb * vadvb_rho1 +
532  fmomx_wadvb * wadvb_rho1 + fmomx_presb * presb_rho1;
533 
534  jac1(1, 1) = fmomx_rhob * rhob_rhou1 + fmomx_uadvb * uadvb_rhou1 + fmomx_vadvb * vadvb_rhou1 +
535  fmomx_wadvb * wadvb_rhou1 + fmomx_presb * presb_rhou1;
536 
537  jac1(1, 2) = fmomx_rhob * rhob_rhov1 + fmomx_uadvb * uadvb_rhov1 + fmomx_vadvb * vadvb_rhov1 +
538  fmomx_wadvb * wadvb_rhov1 + fmomx_presb * presb_rhov1;
539 
540  jac1(1, 3) = fmomx_rhob * rhob_rhow1 + fmomx_uadvb * uadvb_rhow1 + fmomx_vadvb * vadvb_rhow1 +
541  fmomx_wadvb * wadvb_rhow1 + fmomx_presb * presb_rhow1;
542 
543  jac1(1, 4) = fmomx_rhob * rhob_rhoe1 + fmomx_uadvb * uadvb_rhoe1 + fmomx_vadvb * vadvb_rhoe1 +
544  fmomx_wadvb * wadvb_rhoe1 + fmomx_presb * presb_rhoe1;
545 
547 
548  jac1(2, 0) = fmomy_rhob * rhob_rho1 + fmomy_uadvb * uadvb_rho1 + fmomy_vadvb * vadvb_rho1 +
549  fmomy_wadvb * wadvb_rho1 + fmomy_presb * presb_rho1;
550 
551  jac1(2, 1) = fmomy_rhob * rhob_rhou1 + fmomy_uadvb * uadvb_rhou1 + fmomy_vadvb * vadvb_rhou1 +
552  fmomy_wadvb * wadvb_rhou1 + fmomy_presb * presb_rhou1;
553 
554  jac1(2, 2) = fmomy_rhob * rhob_rhov1 + fmomy_uadvb * uadvb_rhov1 + fmomy_vadvb * vadvb_rhov1 +
555  fmomy_wadvb * wadvb_rhov1 + fmomy_presb * presb_rhov1;
556 
557  jac1(2, 3) = fmomy_rhob * rhob_rhow1 + fmomy_uadvb * uadvb_rhow1 + fmomy_vadvb * vadvb_rhow1 +
558  fmomy_wadvb * wadvb_rhow1 + fmomy_presb * presb_rhow1;
559 
560  jac1(2, 4) = fmomy_rhob * rhob_rhoe1 + fmomy_uadvb * uadvb_rhoe1 + fmomy_vadvb * vadvb_rhoe1 +
561  fmomy_wadvb * wadvb_rhoe1 + fmomy_presb * presb_rhoe1;
562 
564 
565  jac1(3, 0) = fmomz_rhob * rhob_rho1 + fmomz_uadvb * uadvb_rho1 + fmomz_vadvb * vadvb_rho1 +
566  fmomz_wadvb * wadvb_rho1 + fmomz_presb * presb_rho1;
567 
568  jac1(3, 1) = fmomz_rhob * rhob_rhou1 + fmomz_uadvb * uadvb_rhou1 + fmomz_vadvb * vadvb_rhou1 +
569  fmomz_wadvb * wadvb_rhou1 + fmomz_presb * presb_rhou1;
570 
571  jac1(3, 2) = fmomz_rhob * rhob_rhov1 + fmomz_uadvb * uadvb_rhov1 + fmomz_vadvb * vadvb_rhov1 +
572  fmomz_wadvb * wadvb_rhov1 + fmomz_presb * presb_rhov1;
573 
574  jac1(3, 3) = fmomz_rhob * rhob_rhow1 + fmomz_uadvb * uadvb_rhow1 + fmomz_vadvb * vadvb_rhow1 +
575  fmomz_wadvb * wadvb_rhow1 + fmomz_presb * presb_rhow1;
576 
577  jac1(3, 4) = fmomz_rhob * rhob_rhoe1 + fmomz_uadvb * uadvb_rhoe1 + fmomz_vadvb * vadvb_rhoe1 +
578  fmomz_wadvb * wadvb_rhoe1 + fmomz_presb * presb_rhoe1;
579 
581 
582  jac1(4, 0) = fetot_rhob * rhob_rho1 + fetot_uadvb * uadvb_rho1 + fetot_vadvb * vadvb_rho1 +
583  fetot_wadvb * wadvb_rho1 + fetot_presb * presb_rho1;
584 
585  jac1(4, 1) = fetot_rhob * rhob_rhou1 + fetot_uadvb * uadvb_rhou1 + fetot_vadvb * vadvb_rhou1 +
586  fetot_wadvb * wadvb_rhou1 + fetot_presb * presb_rhou1;
587 
588  jac1(4, 2) = fetot_rhob * rhob_rhov1 + fetot_uadvb * uadvb_rhov1 + fetot_vadvb * vadvb_rhov1 +
589  fetot_wadvb * wadvb_rhov1 + fetot_presb * presb_rhov1;
590 
591  jac1(4, 3) = fetot_rhob * rhob_rhow1 + fetot_uadvb * uadvb_rhow1 + fetot_vadvb * vadvb_rhow1 +
592  fetot_wadvb * wadvb_rhow1 + fetot_presb * presb_rhow1;
593 
594  jac1(4, 4) = fetot_rhob * rhob_rhoe1 + fetot_uadvb * uadvb_rhoe1 + fetot_vadvb * vadvb_rhoe1 +
595  fetot_wadvb * wadvb_rhoe1 + fetot_presb * presb_rhoe1;
596  }
597  else if (mach1 <= -1.)
598  {
600 
602  }
603  else if (mach1 >= 1.)
604  {
606 
607  Real rq051 = 0.5 * gamm1 * vdov1;
608  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
609  Real enth1 = (rhoe1 + pres1) * rhom1;
610 
611  jac1(0, 0) = 0.;
612  jac1(0, 1) = nx;
613  jac1(0, 2) = ny;
614  jac1(0, 3) = nz;
615  jac1(0, 4) = 0.;
616 
617  jac1(1, 0) = rq051 * nx - uadv1 * vdon1;
618  jac1(1, 1) = gamm2 * nx * uadv1 + vdon1;
619  jac1(1, 2) = ny * uadv1 - vadv1 * gamm1 * nx;
620  jac1(1, 3) = nz * uadv1 - wadv1 * gamm1 * nx;
621  jac1(1, 4) = gamm1 * nx;
622 
623  jac1(2, 0) = rq051 * ny - vadv1 * vdon1;
624  jac1(2, 1) = nx * vadv1 - uadv1 * gamm1 * ny;
625  jac1(2, 2) = gamm2 * ny * vadv1 + vdon1;
626  jac1(2, 3) = nz * vadv1 - wadv1 * gamm1 * ny;
627  jac1(2, 4) = gamm1 * ny;
628 
629  jac1(3, 0) = rq051 * nz - wadv1 * vdon1;
630  jac1(3, 1) = nx * wadv1 - uadv1 * gamm1 * nz;
631  jac1(3, 2) = ny * wadv1 - vadv1 * gamm1 * nz;
632  jac1(3, 3) = gamm2 * nz * wadv1 + vdon1;
633  jac1(3, 4) = gamm1 * nz;
634 
635  jac1(4, 0) = (rq051 - enth1) * vdon1;
636  jac1(4, 1) = nx * enth1 - gamm1 * uadv1 * vdon1;
637  jac1(4, 2) = ny * enth1 - gamm1 * vadv1 * vdon1;
638  jac1(4, 3) = nz * enth1 - gamm1 * wadv1 * vdon1;
639  jac1(4, 4) = gamma * vdon1;
640  }
641  else
642  mooseError("Something is wrong in ",
643  name(),
644  ": ",
645  __FUNCTION__,
646  "\n",
647  "ielem = ",
648  ielem,
649  "\n",
650  "iside = ",
651  iside,
652  "\n",
653  "mach1 = ",
654  mach1,
655  "\n");
656 }
virtual void calcFlux(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, std::vector< Real > &flux) const
Solve the Riemann problem on the boundary face.
InputParameters validParams< BoundaryFluxBase >()
virtual Real c(Real v, Real u) const =0
Sound speed.
CNSFVRiemannInvariantBoundaryFlux(const InputParameters &parameters)
A base class for computing/caching fluxes at boundaries.
Common class for single phase fluid properties.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual std::vector< Real > getGhostCellValue(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave) const =0
compute the ghost cell variable values
A base class of user object for calculating the variable values in ghost element according to specifi...
Definition: BCUserObject.h:42
InputParameters validParams< CNSFVRiemannInvariantBoundaryFlux >()
virtual Real gamma(Real v, Real u) const
Compute the ratio of specific heats.
virtual void calcJacobian(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, DenseMatrix< Real > &jac1) const
Compute the Jacobian matrix on the boundary face.
virtual Real pressure(Real v, Real u) const =0
Pressure as a function of specific internal energy and specific volume.
virtual Real e(Real pressure, Real rho) const =0
Computes internal energy from pressure and density.