www.mooseframework.org
CNSFVHLLCSlipBoundaryFlux.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 object that computes the slip boundary flux using the HLLC "
17  "approximate Riemann solver.");
18 
19  params.addRequiredParam<UserObjectName>("bc_uo", "Name for boundary condition user object");
20 
21  params.addRequiredParam<UserObjectName>("fluid_properties",
22  "Name for fluid properties user object");
23 
24  return params;
25 }
26 
27 CNSFVHLLCSlipBoundaryFlux::CNSFVHLLCSlipBoundaryFlux(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 {
43  Real eps = 1e-6;
44 
46 
47  Real rho1 = uvec1[0];
48  Real rhou1 = uvec1[1];
49  Real rhov1 = uvec1[2];
50  Real rhow1 = uvec1[3];
51  Real rhoe1 = uvec1[4];
52 
53  Real nx = dwave(0);
54  Real ny = dwave(1);
55  Real nz = dwave(2);
56 
58 
59  flux.resize(5);
60 
62 
63  Real uadv1 = rhou1 / rho1;
64  Real vadv1 = rhov1 / rho1;
65  Real wadv1 = rhow1 / rho1;
66  Real vdov1 = uadv1 * uadv1 + vadv1 * vadv1 + wadv1 * wadv1;
67  Real v1 = 1. / rho1;
68  Real e1 = rhoe1 / rho1 - 0.5 * vdov1;
69  Real pres1 = _fp.pressure(v1, e1);
70  Real csou1 = _fp.c(v1, e1);
71  Real gamma = _fp.gamma(v1, e1);
72  Real enth1 = (rhoe1 + pres1) / rho1;
73 
75 
76  std::vector<Real> U2(5, 0.);
77 
78  U2 = _bc_uo.getGhostCellValue(iside, ielem, uvec1, dwave);
79 
80  Real rho2 = U2[0];
81  Real rhou2 = U2[1];
82  Real rhov2 = U2[2];
83  Real rhow2 = U2[3];
84  Real rhoe2 = U2[4];
85 
86  Real uadv2 = rhou2 / rho2;
87  Real vadv2 = rhov2 / rho2;
88  Real wadv2 = rhow2 / rho2;
89  Real vdov2 = uadv2 * uadv2 + vadv2 * vadv2 + wadv2 * wadv2;
90  Real v2 = 1. / rho2;
91  Real e2 = rhoe2 / rho2 - 0.5 * vdov2;
92  Real pres2 = _fp.pressure(v2, e2);
93 
94  Real csou2 = _fp.c(v2, e2);
95  Real enth2 = (rhoe2 + pres2) / rho2;
96 
98 
99  Real rhsca = std::sqrt(rho2 / rho1);
100  Real rmden = 1. / (rhsca + 1.);
101 
102  // Real rhoav = rhsca * rho1;
103  Real uaver = (rhsca * uadv2 + uadv1) * rmden;
104  Real vaver = (rhsca * vadv2 + vadv1) * rmden;
105  Real waver = (rhsca * wadv2 + wadv1) * rmden;
106  Real entav = (rhsca * enth2 + enth1) * rmden;
107 
109 
110  Real qave5 = 0.5 * (uaver * uaver + vaver * vaver + waver * waver);
111  Real cssa2 = std::max(eps, (gamma - 1.) * (entav - qave5));
112  Real cssav = std::sqrt(cssa2);
113 
115 
116  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
117  Real vdon2 = uadv2 * nx + vadv2 * ny + wadv2 * nz;
118  Real vnave = uaver * nx + vaver * ny + waver * nz;
119 
121 
122  Real s1 = std::min(vdon1 - csou1, vnave - cssav);
123  Real s2 = std::max(vdon2 + csou2, vnave + cssav);
124 
125  Real dsv1 = s1 - vdon1;
126  Real dsv2 = s2 - vdon2;
127 
129 
130  Real sm =
131  (rho2 * vdon2 * dsv2 - rho1 * vdon1 * dsv1 + pres1 - pres2) / (rho2 * dsv2 - rho1 * dsv1);
132 
134 
135  Real omeg1 = 1. / (s1 - sm);
136  Real omeg2 = 1. / (s2 - sm);
137  Real prsta = rho1 * dsv1 * (sm - vdon1) + pres1;
138 
139  Real prst1 = prsta - pres1;
140  Real prst2 = prsta - pres2;
141 
143 
144  Real rhol = omeg1 * dsv1 * rho1;
145  Real rhoul = omeg1 * (dsv1 * rhou1 + prst1 * nx);
146  Real rhovl = omeg1 * (dsv1 * rhov1 + prst1 * ny);
147  Real rhowl = omeg1 * (dsv1 * rhow1 + prst1 * nz);
148  Real rhoel = omeg1 * (dsv1 * rhoe1 - pres1 * vdon1 + prsta * sm);
149 
150  Real rhor = omeg2 * dsv2 * rho2;
151  Real rhour = omeg2 * (dsv2 * rhou2 + prst2 * nx);
152  Real rhovr = omeg2 * (dsv2 * rhov2 + prst2 * ny);
153  Real rhowr = omeg2 * (dsv2 * rhow2 + prst2 * nz);
154  Real rhoer = omeg2 * (dsv2 * rhoe2 - pres2 * vdon2 + prsta * sm);
155 
157 
158  if (s1 > 0.)
159  {
160  flux[0] = vdon1 * rho1;
161  flux[1] = vdon1 * rhou1 + pres1 * nx;
162  flux[2] = vdon1 * rhov1 + pres1 * ny;
163  flux[3] = vdon1 * rhow1 + pres1 * nz;
164  flux[4] = vdon1 * (rhoe1 + pres1);
165  }
166  else if (s1 <= 0. && sm > 0.)
167  {
168  flux[0] = sm * rhol;
169  flux[1] = sm * rhoul + prsta * nx;
170  flux[2] = sm * rhovl + prsta * ny;
171  flux[3] = sm * rhowl + prsta * nz;
172  flux[4] = sm * (rhoel + prsta);
173  }
174  else if (sm <= 0. && s2 >= 0.)
175  {
176  flux[0] = sm * rhor;
177  flux[1] = sm * rhour + prsta * nx;
178  flux[2] = sm * rhovr + prsta * ny;
179  flux[3] = sm * rhowr + prsta * nz;
180  flux[4] = sm * (rhoer + prsta);
181  }
182  else if (s2 < 0.)
183  {
184  flux[0] = vdon2 * rho2;
185  flux[1] = vdon2 * rhou2 + pres2 * nx;
186  flux[2] = vdon2 * rhov2 + pres2 * ny;
187  flux[3] = vdon2 * rhow2 + pres2 * nz;
188  flux[4] = vdon2 * (rhoe2 + pres2);
189  }
190  else
191  {
192  mooseError("Weird wave speed occured in ",
193  name(),
194  ": ",
195  __FUNCTION__,
196  "\n",
197  "iside = ",
198  iside,
199  "\n",
200  "ielem = ",
201  ielem,
202  "\n",
203  "rho1 = ",
204  rho1,
205  "\n",
206  "rhou1 = ",
207  rhou1,
208  "\n",
209  "rhov1 = ",
210  rhov1,
211  "\n",
212  "rhow1 = ",
213  rhow1,
214  "\n",
215  "rhoe1 = ",
216  rhoe1,
217  "\n",
218  "pres1 = ",
219  pres1,
220  "\n",
221  "enth1 = ",
222  enth1,
223  "\n",
224  "csou1 = ",
225  csou1,
226  "\n",
227  "rho2 = ",
228  rho2,
229  "\n",
230  "rhou2 = ",
231  rhou2,
232  "\n",
233  "rhov2 = ",
234  rhov2,
235  "\n",
236  "rhoe2 = ",
237  rhoe2,
238  "\n",
239  "pres2 = ",
240  pres2,
241  "\n",
242  "enth2 = ",
243  enth2,
244  "\n",
245  "csou2 = ",
246  csou2,
247  "\n",
248  "vdon1 = ",
249  vdon1,
250  "\n",
251  "vdon2 = ",
252  vdon2,
253  "\n",
254  "vnave = ",
255  vnave,
256  "\n",
257  "cssav = ",
258  cssav,
259  "\n",
260  "s1 = ",
261  s1,
262  "\n",
263  "s2 = ",
264  s2,
265  "\n",
266  "sm = ",
267  sm,
268  "\n",
269  "omeg1 = ",
270  omeg1,
271  "\n",
272  "omeg2 = ",
273  omeg2,
274  "\n",
275  "prsta = ",
276  prsta,
277  "\n",
278  "Please check before continuing!\n");
279  }
280 }
281 
282 void
284  dof_id_type ielem,
285  const std::vector<Real> & uvec1,
286  const RealVectorValue & dwave,
287  DenseMatrix<Real> & jac1) const
288 {
289  Real eps = 1e-6;
290 
292 
293  Real rho1 = uvec1[0];
294  Real rhou1 = uvec1[1];
295  Real rhov1 = uvec1[2];
296  Real rhow1 = uvec1[3];
297  Real rhoe1 = uvec1[4];
298 
299  Real nx = dwave(0);
300  Real ny = dwave(1);
301  Real nz = dwave(2);
302 
304 
305  jac1.resize(5, 5);
306 
308 
309  Real uadv1 = rhou1 / rho1;
310  Real vadv1 = rhov1 / rho1;
311  Real wadv1 = rhow1 / rho1;
312  Real vdov1 = uadv1 * uadv1 + vadv1 * vadv1 + wadv1 * wadv1;
313  Real v1 = 1. / rho1;
314  Real e1 = rhoe1 / rho1 - 0.5 * vdov1;
315  Real pres1 = _fp.pressure(v1, e1);
316  Real csou1 = _fp.c(v1, e1);
317  Real gamma = _fp.gamma(v1, e1);
318  Real gamm1 = gamma - 1.;
319  Real gamm2 = 2. - gamma;
320  Real enth1 = (rhoe1 + pres1) / rho1;
321  Real rq051 = 0.5 * gamm1 * vdov1;
322 
324 
325  std::vector<Real> U2(5, 0.);
326 
327  U2 = _bc_uo.getGhostCellValue(iside, ielem, uvec1, dwave);
328 
329  Real rho2 = U2[0];
330  Real rhou2 = U2[1];
331  Real rhov2 = U2[2];
332  Real rhow2 = U2[3];
333  Real rhoe2 = U2[4];
334 
335  Real uadv2 = rhou2 / rho2;
336  Real vadv2 = rhov2 / rho2;
337  Real wadv2 = rhow2 / rho2;
338  Real vdov2 = uadv2 * uadv2 + vadv2 * vadv2 + wadv2 * wadv2;
339  Real v2 = 1. / rho2;
340  Real e2 = rhoe2 / rho2 - 0.5 * vdov2;
341  Real pres2 = _fp.pressure(v2, e2);
342 
343  Real csou2 = _fp.c(v2, e2);
344  Real enth2 = (rhoe2 + pres2) / rho2;
345  Real rq052 = 0.5 * gamm1 * vdov2;
346 
348 
349  Real rhsca = std::sqrt(rho2 / rho1);
350  Real rmden = 1. / (rhsca + 1.);
351 
352  // Real rhoav = rhsca * rho1;
353  Real uaver = (rhsca * uadv2 + uadv1) * rmden;
354  Real vaver = (rhsca * vadv2 + vadv1) * rmden;
355  Real waver = (rhsca * wadv2 + wadv1) * rmden;
356  Real entav = (rhsca * enth2 + enth1) * rmden;
357 
359 
360  Real qave5 = 0.5 * (uaver * uaver + vaver * vaver + waver * waver);
361  Real cssa2 = std::max(eps, (gamma - 1.) * (entav - qave5));
362  Real cssav = std::sqrt(cssa2);
363 
365 
366  Real vdon1 = uadv1 * nx + vadv1 * ny + wadv1 * nz;
367  Real vdon2 = uadv2 * nx + vadv2 * ny + wadv2 * nz;
368  Real vnave = uaver * nx + vaver * ny + waver * nz;
369 
371 
372  Real s1 = std::min(vdon1 - csou1, vnave - cssav);
373  Real s2 = std::max(vdon2 + csou2, vnave + cssav);
374 
375  Real dsv1 = s1 - vdon1;
376  Real dsv2 = s2 - vdon2;
377 
379 
380  Real sm =
381  (rho2 * vdon2 * dsv2 - rho1 * vdon1 * dsv1 + pres1 - pres2) / (rho2 * dsv2 - rho1 * dsv1);
382 
384 
385  if (s1 > 0.)
386  {
388 
389  jac1(0, 0) = 0.;
390  jac1(0, 1) = nx;
391  jac1(0, 2) = ny;
392  jac1(0, 3) = nz;
393  jac1(0, 4) = 0.;
394 
395  jac1(1, 0) = rq051 * nx - uadv1 * vdon1;
396  jac1(1, 1) = gamm2 * nx * uadv1 + vdon1;
397  jac1(1, 2) = ny * uadv1 - vadv1 * gamm1 * nx;
398  jac1(1, 3) = nz * uadv1 - wadv1 * gamm1 * nx;
399  jac1(1, 4) = gamm1 * nx;
400 
401  jac1(2, 0) = rq051 * ny - vadv1 * vdon1;
402  jac1(2, 1) = nx * vadv1 - uadv1 * gamm1 * ny;
403  jac1(2, 2) = gamm2 * ny * vadv1 + vdon1;
404  jac1(2, 3) = nz * vadv1 - wadv1 * gamm1 * ny;
405  jac1(2, 4) = gamm1 * ny;
406 
407  jac1(3, 0) = rq051 * nz - wadv1 * vdon1;
408  jac1(3, 1) = nx * wadv1 - uadv1 * gamm1 * nz;
409  jac1(3, 2) = ny * wadv1 - vadv1 * gamm1 * nz;
410  jac1(3, 3) = gamm2 * nz * wadv1 + vdon1;
411  jac1(3, 4) = gamm1 * nz;
412 
413  jac1(4, 0) = (rq051 - enth1) * vdon1;
414  jac1(4, 1) = nx * enth1 - gamm1 * uadv1 * vdon1;
415  jac1(4, 2) = ny * enth1 - gamm1 * vadv1 * vdon1;
416  jac1(4, 3) = nz * enth1 - gamm1 * wadv1 * vdon1;
417  jac1(4, 4) = gamma * vdon1;
418  }
419  else if (s1 <= 0. && sm > 0.)
420  {
422 
423  Real omeg1 = 1. / (s1 - sm);
424  Real prsta = rho1 * dsv1 * (sm - vdon1) + pres1;
425 
426  Real prst1 = prsta - pres1;
427 
429 
430  Real rhol = omeg1 * dsv1 * rho1;
431  Real rhoul = omeg1 * (dsv1 * rhou1 + prst1 * nx);
432  Real rhovl = omeg1 * (dsv1 * rhov1 + prst1 * ny);
433  Real rhowl = omeg1 * (dsv1 * rhow1 + prst1 * nz);
434  Real rhoel = omeg1 * (dsv1 * rhoe1 - pres1 * vdon1 + prsta * sm);
435 
436  Real rhoepl = rhoel + prsta;
437 
439 
440  Real rhotd = std::max(eps, rho2 * (s2 - vdon2) - rho1 * (s1 - vdon1));
441  Real rhotm = 1. / rhotd;
442 
444 
445  Real sm_rho1 = rhotm * (-vdon1 * vdon1 + rq051 + sm * s1);
446  Real sm_rhou1 = rhotm * (nx * (2. * vdon1 - s1 - sm) - gamm1 * uadv1);
447  Real sm_rhov1 = rhotm * (ny * (2. * vdon1 - s1 - sm) - gamm1 * vadv1);
448  Real sm_rhow1 = rhotm * (nz * (2. * vdon1 - s1 - sm) - gamm1 * wadv1);
449  Real sm_rhoe1 = rhotm * gamm1;
450 
452 
453  Real sm_rho2 = rhotm * (vdon2 * vdon2 - rq052 - sm * s2);
454  Real sm_rhou2 = rhotm * (nx * (-2. * vdon2 + s2 + sm) + gamm1 * uadv2);
455  Real sm_rhov2 = rhotm * (ny * (-2. * vdon2 + s2 + sm) + gamm1 * vadv2);
456  Real sm_rhow2 = rhotm * (nz * (-2. * vdon2 + s2 + sm) + gamm1 * wadv2);
457  Real sm_rhoe2 = -rhotm * gamm1;
458 
460 
461  Real ps_rho1 = rho2 * dsv2 * sm_rho1;
462  Real ps_rhou1 = rho2 * dsv2 * sm_rhou1;
463  Real ps_rhov1 = rho2 * dsv2 * sm_rhov1;
464  Real ps_rhow1 = rho2 * dsv2 * sm_rhow1;
465  Real ps_rhoe1 = rho2 * dsv2 * sm_rhoe1;
466 
468 
469  Real ps_rho2 = rho1 * dsv1 * sm_rho2;
470  Real ps_rhou2 = rho1 * dsv1 * sm_rhou2;
471  Real ps_rhov2 = rho1 * dsv1 * sm_rhov2;
472  Real ps_rhow2 = rho1 * dsv1 * sm_rhow2;
473  Real ps_rhoe2 = rho1 * dsv1 * sm_rhoe2;
474 
476 
477  Real rhol_rho1 = omeg1 * (s1 + rhol * sm_rho1);
478  Real rhol_rhou1 = omeg1 * (-nx + rhol * sm_rhou1);
479  Real rhol_rhov1 = omeg1 * (-ny + rhol * sm_rhov1);
480  Real rhol_rhow1 = omeg1 * (-nz + rhol * sm_rhow1);
481  Real rhol_rhoe1 = omeg1 * (rhol * sm_rhoe1);
482 
484 
485  Real rhol_rho2 = omeg1 * rhol * sm_rho2;
486  Real rhol_rhou2 = omeg1 * rhol * sm_rhou2;
487  Real rhol_rhov2 = omeg1 * rhol * sm_rhov2;
488  Real rhol_rhow2 = omeg1 * rhol * sm_rhow2;
489  Real rhol_rhoe2 = omeg1 * rhol * sm_rhoe2;
490 
492 
493  Real rhoul_rho1 = omeg1 * (uadv1 * vdon1 - nx * rq051 + nx * ps_rho1 + rhoul * sm_rho1);
494  Real rhoul_rhou1 = omeg1 * (dsv1 - gamm2 * nx * uadv1 + nx * ps_rhou1 + rhoul * sm_rhou1);
495  Real rhoul_rhov1 =
496  omeg1 * (-uadv1 * ny + gamm1 * nx * vadv1 + nx * ps_rhov1 + rhoul * sm_rhov1);
497  Real rhoul_rhow1 =
498  omeg1 * (-uadv1 * nz + gamm1 * nx * wadv1 + nx * ps_rhow1 + rhoul * sm_rhow1);
499  Real rhoul_rhoe1 = omeg1 * (-gamm1 * nx + nx * ps_rhoe1 + rhoul * sm_rhoe1);
500 
502 
503  Real rhoul_rho2 = omeg1 * (nx * ps_rho2 + rhoul * sm_rho2);
504  Real rhoul_rhou2 = omeg1 * (nx * ps_rhou2 + rhoul * sm_rhou2);
505  Real rhoul_rhov2 = omeg1 * (nx * ps_rhov2 + rhoul * sm_rhov2);
506  Real rhoul_rhow2 = omeg1 * (nx * ps_rhow2 + rhoul * sm_rhow2);
507  Real rhoul_rhoe2 = omeg1 * (nx * ps_rhoe2 + rhoul * sm_rhoe2);
508 
510 
511  Real rhovl_rho1 = omeg1 * (vadv1 * vdon1 - ny * rq051 + ny * ps_rho1 + rhovl * sm_rho1);
512  Real rhovl_rhou1 =
513  omeg1 * (-vadv1 * nx + gamm1 * ny * uadv1 + ny * ps_rhou1 + rhovl * sm_rhou1);
514  Real rhovl_rhov1 = omeg1 * (dsv1 - gamm2 * ny * vadv1 + ny * ps_rhov1 + rhovl * sm_rhov1);
515  Real rhovl_rhow1 =
516  omeg1 * (-vadv1 * nz + gamm1 * ny * wadv1 + ny * ps_rhow1 + rhovl * sm_rhow1);
517  Real rhovl_rhoe1 = omeg1 * (-gamm1 * ny + ny * ps_rhoe1 + rhovl * sm_rhoe1);
518 
520 
521  Real rhovl_rho2 = omeg1 * (ny * ps_rho2 + rhovl * sm_rho2);
522  Real rhovl_rhou2 = omeg1 * (ny * ps_rhou2 + rhovl * sm_rhou2);
523  Real rhovl_rhov2 = omeg1 * (ny * ps_rhov2 + rhovl * sm_rhov2);
524  Real rhovl_rhow2 = omeg1 * (ny * ps_rhow2 + rhovl * sm_rhow2);
525  Real rhovl_rhoe2 = omeg1 * (ny * ps_rhoe2 + rhovl * sm_rhoe2);
526 
528 
529  Real rhowl_rho1 = omeg1 * (wadv1 * vdon1 - nz * rq051 + nz * ps_rho1 + rhowl * sm_rho1);
530  Real rhowl_rhou1 =
531  omeg1 * (-wadv1 * nx + gamm1 * nz * uadv1 + nz * ps_rhou1 + rhowl * sm_rhou1);
532  Real rhowl_rhov1 =
533  omeg1 * (-wadv1 * ny + gamm1 * nz * vadv1 + nz * ps_rhov1 + rhowl * sm_rhov1);
534  Real rhowl_rhow1 = omeg1 * (dsv1 - gamm2 * nz * wadv1 + nz * ps_rhow1 + rhowl * sm_rhow1);
535  Real rhowl_rhoe1 = omeg1 * (-gamm1 * nz + nz * ps_rhoe1 + rhowl * sm_rhoe1);
536 
538 
539  Real rhowl_rho2 = omeg1 * (nz * ps_rho2 + rhowl * sm_rho2);
540  Real rhowl_rhou2 = omeg1 * (nz * ps_rhou2 + rhowl * sm_rhou2);
541  Real rhowl_rhov2 = omeg1 * (nz * ps_rhov2 + rhowl * sm_rhov2);
542  Real rhowl_rhow2 = omeg1 * (nz * ps_rhow2 + rhowl * sm_rhow2);
543  Real rhowl_rhoe2 = omeg1 * (nz * ps_rhoe2 + rhowl * sm_rhoe2);
544 
548 
549  Real rhoel_rho1 = omeg1 * (vdon1 * enth1 - vdon1 * rq051 + sm * ps_rho1 + rhoepl * sm_rho1);
550  Real rhoel_rhou1 =
551  omeg1 * (-nx * enth1 + gamm1 * vdon1 * uadv1 + sm * ps_rhou1 + rhoepl * sm_rhou1);
552  Real rhoel_rhov1 =
553  omeg1 * (-ny * enth1 + gamm1 * vdon1 * vadv1 + sm * ps_rhov1 + rhoepl * sm_rhov1);
554  Real rhoel_rhow1 =
555  omeg1 * (-nz * enth1 + gamm1 * vdon1 * wadv1 + sm * ps_rhow1 + rhoepl * sm_rhow1);
556  Real rhoel_rhoe1 = omeg1 * (s1 - vdon1 * gamma + sm * ps_rhoe1 + rhoepl * sm_rhoe1);
557 
559 
560  Real rhoel_rho2 = omeg1 * (sm * ps_rho2 + rhoepl * sm_rho2);
561  Real rhoel_rhou2 = omeg1 * (sm * ps_rhou2 + rhoepl * sm_rhou2);
562  Real rhoel_rhov2 = omeg1 * (sm * ps_rhov2 + rhoepl * sm_rhov2);
563  Real rhoel_rhow2 = omeg1 * (sm * ps_rhow2 + rhoepl * sm_rhow2);
564  Real rhoel_rhoe2 = omeg1 * (sm * ps_rhoe2 + rhoepl * sm_rhoe2);
565 
567 
568  jac1(0, 0) = sm * rhol_rho1 + rhol * sm_rho1;
569  jac1(0, 1) = sm * rhol_rhou1 + rhol * sm_rhou1;
570  jac1(0, 2) = sm * rhol_rhov1 + rhol * sm_rhov1;
571  jac1(0, 3) = sm * rhol_rhow1 + rhol * sm_rhow1;
572  jac1(0, 4) = sm * rhol_rhoe1 + rhol * sm_rhoe1;
573 
574  jac1(1, 0) = sm * rhoul_rho1 + rhoul * sm_rho1 + nx * ps_rho1;
575  jac1(1, 1) = sm * rhoul_rhou1 + rhoul * sm_rhou1 + nx * ps_rhou1;
576  jac1(1, 2) = sm * rhoul_rhov1 + rhoul * sm_rhov1 + nx * ps_rhov1;
577  jac1(1, 3) = sm * rhoul_rhow1 + rhoul * sm_rhow1 + nx * ps_rhow1;
578  jac1(1, 4) = sm * rhoul_rhoe1 + rhoul * sm_rhoe1 + nx * ps_rhoe1;
579 
580  jac1(2, 0) = sm * rhovl_rho1 + rhovl * sm_rho1 + ny * ps_rho1;
581  jac1(2, 1) = sm * rhovl_rhou1 + rhovl * sm_rhou1 + ny * ps_rhou1;
582  jac1(2, 2) = sm * rhovl_rhov1 + rhovl * sm_rhov1 + ny * ps_rhov1;
583  jac1(2, 3) = sm * rhovl_rhow1 + rhovl * sm_rhow1 + ny * ps_rhow1;
584  jac1(2, 4) = sm * rhovl_rhoe1 + rhovl * sm_rhoe1 + ny * ps_rhoe1;
585 
586  jac1(3, 0) = sm * rhowl_rho1 + rhowl * sm_rho1 + nz * ps_rho1;
587  jac1(3, 1) = sm * rhowl_rhou1 + rhowl * sm_rhou1 + nz * ps_rhou1;
588  jac1(3, 2) = sm * rhowl_rhov1 + rhowl * sm_rhov1 + nz * ps_rhov1;
589  jac1(3, 3) = sm * rhowl_rhow1 + rhowl * sm_rhow1 + nz * ps_rhow1;
590  jac1(3, 4) = sm * rhowl_rhoe1 + rhowl * sm_rhoe1 + nz * ps_rhoe1;
591 
592  jac1(4, 0) = sm * (rhoel_rho1 + ps_rho1) + rhoepl * sm_rho1;
593  jac1(4, 1) = sm * (rhoel_rhou1 + ps_rhou1) + rhoepl * sm_rhou1;
594  jac1(4, 2) = sm * (rhoel_rhov1 + ps_rhov1) + rhoepl * sm_rhov1;
595  jac1(4, 3) = sm * (rhoel_rhow1 + ps_rhow1) + rhoepl * sm_rhow1;
596  jac1(4, 4) = sm * (rhoel_rhoe1 + ps_rhoe1) + rhoepl * sm_rhoe1;
597 
599 
600  DenseMatrix<Real> jac2(5, 5);
601 
602  jac2(0, 0) = sm * rhol_rho2 + rhol * sm_rho2;
603  jac2(0, 1) = sm * rhol_rhou2 + rhol * sm_rhou2;
604  jac2(0, 2) = sm * rhol_rhov2 + rhol * sm_rhov2;
605  jac2(0, 3) = sm * rhol_rhow2 + rhol * sm_rhow2;
606  jac2(0, 4) = sm * rhol_rhoe2 + rhol * sm_rhoe2;
607 
608  jac2(1, 0) = sm * rhoul_rho2 + rhoul * sm_rho2 + nx * ps_rho2;
609  jac2(1, 1) = sm * rhoul_rhou2 + rhoul * sm_rhou2 + nx * ps_rhou2;
610  jac2(1, 2) = sm * rhoul_rhov2 + rhoul * sm_rhov2 + nx * ps_rhov2;
611  jac2(1, 3) = sm * rhoul_rhow2 + rhoul * sm_rhow2 + nx * ps_rhow2;
612  jac2(1, 4) = sm * rhoul_rhoe2 + rhoul * sm_rhoe2 + nx * ps_rhoe2;
613 
614  jac2(2, 0) = sm * rhovl_rho2 + rhovl * sm_rho2 + ny * ps_rho2;
615  jac2(2, 1) = sm * rhovl_rhou2 + rhovl * sm_rhou2 + ny * ps_rhou2;
616  jac2(2, 2) = sm * rhovl_rhov2 + rhovl * sm_rhov2 + ny * ps_rhov2;
617  jac2(2, 3) = sm * rhovl_rhow2 + rhovl * sm_rhow2 + ny * ps_rhow2;
618  jac2(2, 4) = sm * rhovl_rhoe2 + rhovl * sm_rhoe2 + ny * ps_rhoe2;
619 
620  jac2(3, 0) = sm * rhowl_rho2 + rhowl * sm_rho2 + nz * ps_rho2;
621  jac2(3, 1) = sm * rhowl_rhou2 + rhowl * sm_rhou2 + nz * ps_rhou2;
622  jac2(3, 2) = sm * rhowl_rhov2 + rhowl * sm_rhov2 + nz * ps_rhov2;
623  jac2(3, 3) = sm * rhowl_rhow2 + rhowl * sm_rhow2 + nz * ps_rhow2;
624  jac2(3, 4) = sm * rhowl_rhoe2 + rhowl * sm_rhoe2 + nz * ps_rhoe2;
625 
626  jac2(4, 0) = sm * (rhoel_rho2 + ps_rho2) + rhoepl * sm_rho2;
627  jac2(4, 1) = sm * (rhoel_rhou2 + ps_rhou2) + rhoepl * sm_rhou2;
628  jac2(4, 2) = sm * (rhoel_rhov2 + ps_rhov2) + rhoepl * sm_rhov2;
629  jac2(4, 3) = sm * (rhoel_rhow2 + ps_rhow2) + rhoepl * sm_rhow2;
630  jac2(4, 4) = sm * (rhoel_rhoe2 + ps_rhoe2) + rhoepl * sm_rhoe2;
631 
633 
634  Real uu11 = 1. - 2. * nx * nx;
635  Real uu22 = 1. - 2. * ny * ny;
636  Real uu33 = 1. - 2. * nz * nz;
637  Real uu12 = -2. * nx * ny;
638  Real uu13 = -2. * nx * nz;
639  Real uu23 = -2. * ny * nz;
640 
642 
643  DenseMatrix<Real> dhdu(5, 5);
644 
645  dhdu(0, 0) = jac2(0, 0);
646  dhdu(0, 1) = jac2(0, 1) * uu11 + jac2(0, 2) * uu12 + jac2(0, 3) * uu13;
647  dhdu(0, 2) = jac2(0, 1) * uu12 + jac2(0, 2) * uu22 + jac2(0, 3) * uu23;
648  dhdu(0, 3) = jac2(0, 1) * uu13 + jac2(0, 2) * uu23 + jac2(0, 3) * uu33;
649  dhdu(0, 4) = jac2(0, 4);
650 
651  dhdu(1, 0) = jac2(1, 0);
652  dhdu(1, 1) = jac2(1, 1) * uu11 + jac2(1, 2) * uu12 + jac2(1, 3) * uu13;
653  dhdu(1, 2) = jac2(1, 1) * uu12 + jac2(1, 2) * uu22 + jac2(1, 3) * uu23;
654  dhdu(1, 3) = jac2(1, 1) * uu13 + jac2(1, 2) * uu23 + jac2(1, 3) * uu33;
655  dhdu(1, 4) = jac2(1, 4);
656 
657  dhdu(2, 0) = jac2(2, 0);
658  dhdu(2, 1) = jac2(2, 1) * uu11 + jac2(2, 2) * uu12 + jac2(2, 3) * uu13;
659  dhdu(2, 2) = jac2(2, 1) * uu12 + jac2(2, 2) * uu22 + jac2(2, 3) * uu23;
660  dhdu(2, 3) = jac2(2, 1) * uu13 + jac2(2, 2) * uu23 + jac2(2, 3) * uu33;
661  dhdu(2, 4) = jac2(2, 4);
662 
663  dhdu(3, 0) = jac2(3, 0);
664  dhdu(3, 1) = jac2(3, 1) * uu11 + jac2(3, 2) * uu12 + jac2(3, 3) * uu13;
665  dhdu(3, 2) = jac2(3, 1) * uu12 + jac2(3, 2) * uu22 + jac2(3, 3) * uu23;
666  dhdu(3, 3) = jac2(3, 1) * uu13 + jac2(3, 2) * uu23 + jac2(3, 3) * uu33;
667  dhdu(3, 4) = jac2(3, 4);
668 
669  dhdu(4, 0) = jac2(4, 0);
670  dhdu(4, 1) = jac2(4, 1) * uu11 + jac2(4, 2) * uu12 + jac2(4, 3) * uu13;
671  dhdu(4, 2) = jac2(4, 1) * uu12 + jac2(4, 2) * uu22 + jac2(4, 3) * uu23;
672  dhdu(4, 3) = jac2(4, 1) * uu13 + jac2(4, 2) * uu23 + jac2(4, 3) * uu33;
673  dhdu(4, 4) = jac2(4, 4);
674 
675  jac1 += dhdu;
676  }
677  else if (sm <= 0. && s2 >= 0.)
678  {
680 
681  Real omeg2 = 1. / (s2 - sm);
682  Real prsta = rho2 * dsv2 * (sm - vdon2) + pres2;
683 
684  Real prst2 = prsta - pres2;
685 
687 
688  Real rhor = omeg2 * dsv2 * rho2;
689  Real rhour = omeg2 * (dsv2 * rhou2 + prst2 * nx);
690  Real rhovr = omeg2 * (dsv2 * rhov2 + prst2 * ny);
691  Real rhowr = omeg2 * (dsv2 * rhow2 + prst2 * nz);
692  Real rhoer = omeg2 * (dsv2 * rhoe2 - pres2 * vdon2 + prsta * sm);
693 
694  Real rhoepr = rhoer + prsta;
695 
697 
698  Real rhotd = std::max(eps, rho2 * (s2 - vdon2) - rho1 * (s1 - vdon1));
699  Real rhotm = 1. / rhotd;
700 
702 
703  Real sm_rho1 = rhotm * (-vdon1 * vdon1 + rq051 + sm * s1);
704  Real sm_rhou1 = rhotm * (nx * (2. * vdon1 - s1 - sm) - gamm1 * uadv1);
705  Real sm_rhov1 = rhotm * (ny * (2. * vdon1 - s1 - sm) - gamm1 * vadv1);
706  Real sm_rhow1 = rhotm * (nz * (2. * vdon1 - s1 - sm) - gamm1 * wadv1);
707  Real sm_rhoe1 = rhotm * gamm1;
708 
710 
711  Real sm_rho2 = rhotm * (vdon2 * vdon2 - rq052 - sm * s2);
712  Real sm_rhou2 = rhotm * (nx * (-2. * vdon2 + s2 + sm) + gamm1 * uadv2);
713  Real sm_rhov2 = rhotm * (ny * (-2. * vdon2 + s2 + sm) + gamm1 * vadv2);
714  Real sm_rhow2 = rhotm * (nz * (-2. * vdon2 + s2 + sm) + gamm1 * wadv2);
715  Real sm_rhoe2 = -rhotm * gamm1;
716 
718 
719  Real ps_rho1 = rho2 * dsv2 * sm_rho1;
720  Real ps_rhou1 = rho2 * dsv2 * sm_rhou1;
721  Real ps_rhov1 = rho2 * dsv2 * sm_rhov1;
722  Real ps_rhow1 = rho2 * dsv2 * sm_rhow1;
723  Real ps_rhoe1 = rho2 * dsv2 * sm_rhoe1;
724 
726 
727  Real ps_rho2 = rho1 * dsv1 * sm_rho2;
728  Real ps_rhou2 = rho1 * dsv1 * sm_rhou2;
729  Real ps_rhov2 = rho1 * dsv1 * sm_rhov2;
730  Real ps_rhow2 = rho1 * dsv1 * sm_rhow2;
731  Real ps_rhoe2 = rho1 * dsv1 * sm_rhoe2;
732 
735 
737 
738  Real rhor_rho2 = omeg2 * (s2 + rhor * sm_rho2);
739  Real rhor_rhou2 = omeg2 * (-nx + rhor * sm_rhou2);
740  Real rhor_rhov2 = omeg2 * (-ny + rhor * sm_rhov2);
741  Real rhor_rhow2 = omeg2 * (-nz + rhor * sm_rhow2);
742  Real rhor_rhoe2 = omeg2 * (rhor * sm_rhoe2);
743 
745 
746  Real rhor_rho1 = omeg2 * rhor * sm_rho1;
747  Real rhor_rhou1 = omeg2 * rhor * sm_rhou1;
748  Real rhor_rhov1 = omeg2 * rhor * sm_rhov1;
749  Real rhor_rhow1 = omeg2 * rhor * sm_rhow1;
750  Real rhor_rhoe1 = omeg2 * rhor * sm_rhoe1;
751 
753 
754  Real rhour_rho2 = omeg2 * (uadv2 * vdon2 - nx * rq052 + nx * ps_rho2 + rhour * sm_rho2);
755  Real rhour_rhou2 = omeg2 * (dsv2 - gamm2 * nx * uadv2 + nx * ps_rhou2 + rhour * sm_rhou2);
756  Real rhour_rhov2 =
757  omeg2 * (-uadv2 * ny + gamm1 * nx * vadv2 + nx * ps_rhov2 + rhour * sm_rhov2);
758  Real rhour_rhow2 =
759  omeg2 * (-uadv2 * nz + gamm1 * nx * wadv2 + nx * ps_rhow2 + rhour * sm_rhow2);
760  Real rhour_rhoe2 = omeg2 * (-gamm1 * nx + nx * ps_rhoe2 + rhour * sm_rhoe2);
761 
763 
764  Real rhour_rho1 = omeg2 * (nx * ps_rho1 + rhour * sm_rho1);
765  Real rhour_rhou1 = omeg2 * (nx * ps_rhou1 + rhour * sm_rhou1);
766  Real rhour_rhov1 = omeg2 * (nx * ps_rhov1 + rhour * sm_rhov1);
767  Real rhour_rhow1 = omeg2 * (nx * ps_rhow1 + rhour * sm_rhow1);
768  Real rhour_rhoe1 = omeg2 * (nx * ps_rhoe1 + rhour * sm_rhoe1);
769 
771 
772  Real rhovr_rho2 = omeg2 * (vadv2 * vdon2 - ny * rq052 + ny * ps_rho2 + rhovr * sm_rho2);
773  Real rhovr_rhou2 =
774  omeg2 * (-vadv2 * nx + gamm1 * ny * uadv2 + ny * ps_rhou2 + rhovr * sm_rhou2);
775  Real rhovr_rhov2 = omeg2 * (dsv2 - gamm2 * ny * vadv2 + ny * ps_rhov2 + rhovr * sm_rhov2);
776  Real rhovr_rhow2 =
777  omeg2 * (-vadv2 * nz + gamm1 * ny * wadv2 + ny * ps_rhow2 + rhovr * sm_rhow2);
778  Real rhovr_rhoe2 = omeg2 * (-gamm1 * ny + ny * ps_rhoe2 + rhovr * sm_rhoe2);
779 
781 
782  Real rhovr_rho1 = omeg2 * (ny * ps_rho1 + rhovr * sm_rho1);
783  Real rhovr_rhou1 = omeg2 * (ny * ps_rhou1 + rhovr * sm_rhou1);
784  Real rhovr_rhov1 = omeg2 * (ny * ps_rhov1 + rhovr * sm_rhov1);
785  Real rhovr_rhow1 = omeg2 * (ny * ps_rhow1 + rhovr * sm_rhow1);
786  Real rhovr_rhoe1 = omeg2 * (ny * ps_rhoe1 + rhovr * sm_rhoe1);
787 
789 
790  Real rhowr_rho2 = omeg2 * (wadv2 * vdon2 - nz * rq052 + nz * ps_rho2 + rhowr * sm_rho2);
791  Real rhowr_rhou2 =
792  omeg2 * (-wadv2 * nx + gamm1 * nz * uadv2 + nz * ps_rhou2 + rhowr * sm_rhou2);
793  Real rhowr_rhov2 =
794  omeg2 * (-wadv2 * ny + gamm1 * nz * vadv2 + nz * ps_rhov2 + rhowr * sm_rhov2);
795  Real rhowr_rhow2 = omeg2 * (dsv2 - gamm2 * nz * wadv2 + nz * ps_rhow2 + rhowr * sm_rhow2);
796  Real rhowr_rhoe2 = omeg2 * (-gamm1 * nz + nz * ps_rhoe2 + rhowr * sm_rhoe2);
797 
799 
800  Real rhowr_rho1 = omeg2 * (nz * ps_rho1 + rhowr * sm_rho1);
801  Real rhowr_rhou1 = omeg2 * (nz * ps_rhou1 + rhowr * sm_rhou1);
802  Real rhowr_rhov1 = omeg2 * (nz * ps_rhov1 + rhowr * sm_rhov1);
803  Real rhowr_rhow1 = omeg2 * (nz * ps_rhow1 + rhowr * sm_rhow1);
804  Real rhowr_rhoe1 = omeg2 * (nz * ps_rhoe1 + rhowr * sm_rhoe1);
805 
809 
810  Real rhoer_rho2 = omeg2 * (vdon2 * enth2 - vdon2 * rq052 + sm * ps_rho2 + rhoepr * sm_rho2);
811  Real rhoer_rhou2 =
812  omeg2 * (-nx * enth2 + gamm1 * vdon2 * uadv2 + sm * ps_rhou2 + rhoepr * sm_rhou2);
813  Real rhoer_rhov2 =
814  omeg2 * (-ny * enth2 + gamm1 * vdon2 * vadv2 + sm * ps_rhov2 + rhoepr * sm_rhov2);
815  Real rhoer_rhow2 =
816  omeg2 * (-nz * enth2 + gamm1 * vdon2 * wadv2 + sm * ps_rhow2 + rhoepr * sm_rhow2);
817  Real rhoer_rhoe2 = omeg2 * (s2 - vdon2 * gamma + sm * ps_rhoe2 + rhoepr * sm_rhoe2);
818 
820 
821  Real rhoer_rho1 = omeg2 * (sm * ps_rho1 + rhoepr * sm_rho1);
822  Real rhoer_rhou1 = omeg2 * (sm * ps_rhou1 + rhoepr * sm_rhou1);
823  Real rhoer_rhov1 = omeg2 * (sm * ps_rhov1 + rhoepr * sm_rhov1);
824  Real rhoer_rhow1 = omeg2 * (sm * ps_rhow1 + rhoepr * sm_rhow1);
825  Real rhoer_rhoe1 = omeg2 * (sm * ps_rhoe1 + rhoepr * sm_rhoe1);
826 
828 
829  DenseMatrix<Real> jac2(5, 5);
830 
831  jac2(0, 0) = sm * rhor_rho2 + rhor * sm_rho2;
832  jac2(0, 1) = sm * rhor_rhou2 + rhor * sm_rhou2;
833  jac2(0, 2) = sm * rhor_rhov2 + rhor * sm_rhov2;
834  jac2(0, 3) = sm * rhor_rhow2 + rhor * sm_rhow2;
835  jac2(0, 4) = sm * rhor_rhoe2 + rhor * sm_rhoe2;
836 
837  jac2(1, 0) = sm * rhour_rho2 + rhour * sm_rho2 + nx * ps_rho2;
838  jac2(1, 1) = sm * rhour_rhou2 + rhour * sm_rhou2 + nx * ps_rhou2;
839  jac2(1, 2) = sm * rhour_rhov2 + rhour * sm_rhov2 + nx * ps_rhov2;
840  jac2(1, 3) = sm * rhour_rhow2 + rhour * sm_rhow2 + nx * ps_rhow2;
841  jac2(1, 4) = sm * rhour_rhoe2 + rhour * sm_rhoe2 + nx * ps_rhoe2;
842 
843  jac2(2, 0) = sm * rhovr_rho2 + rhovr * sm_rho2 + ny * ps_rho2;
844  jac2(2, 1) = sm * rhovr_rhou2 + rhovr * sm_rhou2 + ny * ps_rhou2;
845  jac2(2, 2) = sm * rhovr_rhov2 + rhovr * sm_rhov2 + ny * ps_rhov2;
846  jac2(2, 3) = sm * rhovr_rhow2 + rhovr * sm_rhow2 + ny * ps_rhow2;
847  jac2(2, 4) = sm * rhovr_rhoe2 + rhovr * sm_rhoe2 + ny * ps_rhoe2;
848 
849  jac2(3, 0) = sm * rhowr_rho2 + rhowr * sm_rho2 + nz * ps_rho2;
850  jac2(3, 1) = sm * rhowr_rhou2 + rhowr * sm_rhou2 + nz * ps_rhou2;
851  jac2(3, 2) = sm * rhowr_rhov2 + rhowr * sm_rhov2 + nz * ps_rhov2;
852  jac2(3, 3) = sm * rhowr_rhow2 + rhowr * sm_rhow2 + nz * ps_rhow2;
853  jac2(3, 4) = sm * rhowr_rhoe2 + rhowr * sm_rhoe2 + nz * ps_rhoe2;
854 
855  jac2(4, 0) = sm * (rhoer_rho2 + ps_rho2) + rhoepr * sm_rho2;
856  jac2(4, 1) = sm * (rhoer_rhou2 + ps_rhou2) + rhoepr * sm_rhou2;
857  jac2(4, 2) = sm * (rhoer_rhov2 + ps_rhov2) + rhoepr * sm_rhov2;
858  jac2(4, 3) = sm * (rhoer_rhow2 + ps_rhow2) + rhoepr * sm_rhow2;
859  jac2(4, 4) = sm * (rhoer_rhoe2 + ps_rhoe2) + rhoepr * sm_rhoe2;
860 
862 
863  jac1(0, 0) = sm * rhor_rho1 + rhor * sm_rho1;
864  jac1(0, 1) = sm * rhor_rhou1 + rhor * sm_rhou1;
865  jac1(0, 2) = sm * rhor_rhov1 + rhor * sm_rhov1;
866  jac1(0, 3) = sm * rhor_rhow1 + rhor * sm_rhow1;
867  jac1(0, 4) = sm * rhor_rhoe1 + rhor * sm_rhoe1;
868 
869  jac1(1, 0) = sm * rhour_rho1 + rhour * sm_rho1 + nx * ps_rho1;
870  jac1(1, 1) = sm * rhour_rhou1 + rhour * sm_rhou1 + nx * ps_rhou1;
871  jac1(1, 2) = sm * rhour_rhov1 + rhour * sm_rhov1 + nx * ps_rhov1;
872  jac1(1, 3) = sm * rhour_rhow1 + rhour * sm_rhow1 + nx * ps_rhow1;
873  jac1(1, 4) = sm * rhour_rhoe1 + rhour * sm_rhoe1 + nx * ps_rhoe1;
874 
875  jac1(2, 0) = sm * rhovr_rho1 + rhovr * sm_rho1 + ny * ps_rho1;
876  jac1(2, 1) = sm * rhovr_rhou1 + rhovr * sm_rhou1 + ny * ps_rhou1;
877  jac1(2, 2) = sm * rhovr_rhov1 + rhovr * sm_rhov1 + ny * ps_rhov1;
878  jac1(2, 3) = sm * rhovr_rhow1 + rhovr * sm_rhow1 + ny * ps_rhow1;
879  jac1(2, 4) = sm * rhovr_rhoe1 + rhovr * sm_rhoe1 + ny * ps_rhoe1;
880 
881  jac1(3, 0) = sm * rhowr_rho1 + rhowr * sm_rho1 + nz * ps_rho1;
882  jac1(3, 1) = sm * rhowr_rhou1 + rhowr * sm_rhou1 + nz * ps_rhou1;
883  jac1(3, 2) = sm * rhowr_rhov1 + rhowr * sm_rhov1 + nz * ps_rhov1;
884  jac1(3, 3) = sm * rhowr_rhow1 + rhowr * sm_rhow1 + nz * ps_rhow1;
885  jac1(3, 4) = sm * rhowr_rhoe1 + rhowr * sm_rhoe1 + nz * ps_rhoe1;
886 
887  jac1(4, 0) = sm * (rhoer_rho1 + ps_rho1) + rhoepr * sm_rho1;
888  jac1(4, 1) = sm * (rhoer_rhou1 + ps_rhou1) + rhoepr * sm_rhou1;
889  jac1(4, 2) = sm * (rhoer_rhov1 + ps_rhov1) + rhoepr * sm_rhov1;
890  jac1(4, 3) = sm * (rhoer_rhow1 + ps_rhow1) + rhoepr * sm_rhow1;
891  jac1(4, 4) = sm * (rhoer_rhoe1 + ps_rhoe1) + rhoepr * sm_rhoe1;
892 
894 
895  Real uu11 = 1. - 2. * nx * nx;
896  Real uu22 = 1. - 2. * ny * ny;
897  Real uu33 = 1. - 2. * nz * nz;
898  Real uu12 = -2. * nx * ny;
899  Real uu13 = -2. * nx * nz;
900  Real uu23 = -2. * ny * nz;
901 
903 
904  DenseMatrix<Real> dhdu(5, 5);
905 
906  dhdu(0, 0) = jac2(0, 0);
907  dhdu(0, 1) = jac2(0, 1) * uu11 + jac2(0, 2) * uu12 + jac2(0, 3) * uu13;
908  dhdu(0, 2) = jac2(0, 1) * uu12 + jac2(0, 2) * uu22 + jac2(0, 3) * uu23;
909  dhdu(0, 3) = jac2(0, 1) * uu13 + jac2(0, 2) * uu23 + jac2(0, 3) * uu33;
910  dhdu(0, 4) = jac2(0, 4);
911 
912  dhdu(1, 0) = jac2(1, 0);
913  dhdu(1, 1) = jac2(1, 1) * uu11 + jac2(1, 2) * uu12 + jac2(1, 3) * uu13;
914  dhdu(1, 2) = jac2(1, 1) * uu12 + jac2(1, 2) * uu22 + jac2(1, 3) * uu23;
915  dhdu(1, 3) = jac2(1, 1) * uu13 + jac2(1, 2) * uu23 + jac2(1, 3) * uu33;
916  dhdu(1, 4) = jac2(1, 4);
917 
918  dhdu(2, 0) = jac2(2, 0);
919  dhdu(2, 1) = jac2(2, 1) * uu11 + jac2(2, 2) * uu12 + jac2(2, 3) * uu13;
920  dhdu(2, 2) = jac2(2, 1) * uu12 + jac2(2, 2) * uu22 + jac2(2, 3) * uu23;
921  dhdu(2, 3) = jac2(2, 1) * uu13 + jac2(2, 2) * uu23 + jac2(2, 3) * uu33;
922  dhdu(2, 4) = jac2(2, 4);
923 
924  dhdu(3, 0) = jac2(3, 0);
925  dhdu(3, 1) = jac2(3, 1) * uu11 + jac2(3, 2) * uu12 + jac2(3, 3) * uu13;
926  dhdu(3, 2) = jac2(3, 1) * uu12 + jac2(3, 2) * uu22 + jac2(3, 3) * uu23;
927  dhdu(3, 3) = jac2(3, 1) * uu13 + jac2(3, 2) * uu23 + jac2(3, 3) * uu33;
928  dhdu(3, 4) = jac2(3, 4);
929 
930  dhdu(4, 0) = jac2(4, 0);
931  dhdu(4, 1) = jac2(4, 1) * uu11 + jac2(4, 2) * uu12 + jac2(4, 3) * uu13;
932  dhdu(4, 2) = jac2(4, 1) * uu12 + jac2(4, 2) * uu22 + jac2(4, 3) * uu23;
933  dhdu(4, 3) = jac2(4, 1) * uu13 + jac2(4, 2) * uu23 + jac2(4, 3) * uu33;
934  dhdu(4, 4) = jac2(4, 4);
935 
936  jac1 += dhdu;
937  }
938  else if (s2 < 0.)
939  {
941 
942  DenseMatrix<Real> jac2(5, 5);
943 
944  jac2(0, 0) = 0.;
945  jac2(0, 1) = nx;
946  jac2(0, 2) = ny;
947  jac2(0, 3) = nz;
948  jac2(0, 4) = 0.;
949 
950  jac2(1, 0) = rq052 * nx - uadv2 * vdon2;
951  jac2(1, 1) = gamm2 * nx * uadv2 + vdon2;
952  jac2(1, 2) = ny * uadv2 - vadv2 * gamm1 * nx;
953  jac2(1, 3) = nz * uadv2 - wadv2 * gamm1 * nx;
954  jac2(1, 4) = gamm1 * nx;
955 
956  jac2(2, 0) = rq052 * ny - vadv2 * vdon2;
957  jac2(2, 1) = nx * vadv2 - uadv2 * gamm1 * ny;
958  jac2(2, 2) = gamm2 * ny * vadv2 + vdon2;
959  jac2(2, 3) = nz * vadv2 - wadv2 * gamm1 * ny;
960  jac2(2, 4) = gamm1 * ny;
961 
962  jac2(3, 0) = rq052 * nz - wadv2 * vdon2;
963  jac2(3, 1) = nx * wadv2 - uadv2 * gamm1 * ny;
964  jac2(3, 2) = ny * wadv2 - vadv2 * gamm1 * ny;
965  jac2(3, 3) = gamm2 * nz * wadv2 + vdon2;
966  jac2(3, 4) = gamm1 * ny;
967 
968  jac2(4, 0) = (rq052 - enth2) * vdon2;
969  jac2(4, 1) = nx * enth2 - gamm1 * uadv2 * vdon2;
970  jac2(4, 2) = ny * enth2 - gamm1 * vadv2 * vdon2;
971  jac2(4, 3) = nz * enth2 - gamm1 * wadv2 * vdon2;
972  jac2(4, 4) = gamma * vdon2;
973 
975 
976  Real uu11 = 1. - 2. * nx * nx;
977  Real uu22 = 1. - 2. * ny * ny;
978  Real uu33 = 1. - 2. * nz * nz;
979  Real uu12 = -2. * nx * ny;
980  Real uu13 = -2. * nx * nz;
981  Real uu23 = -2. * ny * nz;
982 
984 
985  DenseMatrix<Real> dhdu(5, 5);
986 
987  dhdu(0, 0) = jac2(0, 0);
988  dhdu(0, 1) = jac2(0, 1) * uu11 + jac2(0, 2) * uu12 + jac2(0, 3) * uu13;
989  dhdu(0, 2) = jac2(0, 1) * uu12 + jac2(0, 2) * uu22 + jac2(0, 3) * uu23;
990  dhdu(0, 3) = jac2(0, 1) * uu13 + jac2(0, 2) * uu23 + jac2(0, 3) * uu33;
991  dhdu(0, 4) = jac2(0, 4);
992 
993  dhdu(1, 0) = jac2(1, 0);
994  dhdu(1, 1) = jac2(1, 1) * uu11 + jac2(1, 2) * uu12 + jac2(1, 3) * uu13;
995  dhdu(1, 2) = jac2(1, 1) * uu12 + jac2(1, 2) * uu22 + jac2(1, 3) * uu23;
996  dhdu(1, 3) = jac2(1, 1) * uu13 + jac2(1, 2) * uu23 + jac2(1, 3) * uu33;
997  dhdu(1, 4) = jac2(1, 4);
998 
999  dhdu(2, 0) = jac2(2, 0);
1000  dhdu(2, 1) = jac2(2, 1) * uu11 + jac2(2, 2) * uu12 + jac2(2, 3) * uu13;
1001  dhdu(2, 2) = jac2(2, 1) * uu12 + jac2(2, 2) * uu22 + jac2(2, 3) * uu23;
1002  dhdu(2, 3) = jac2(2, 1) * uu13 + jac2(2, 2) * uu23 + jac2(2, 3) * uu33;
1003  dhdu(2, 4) = jac2(2, 4);
1004 
1005  dhdu(3, 0) = jac2(3, 0);
1006  dhdu(3, 1) = jac2(3, 1) * uu11 + jac2(3, 2) * uu12 + jac2(3, 3) * uu13;
1007  dhdu(3, 2) = jac2(3, 1) * uu12 + jac2(3, 2) * uu22 + jac2(3, 3) * uu23;
1008  dhdu(3, 3) = jac2(3, 1) * uu13 + jac2(3, 2) * uu23 + jac2(3, 3) * uu33;
1009  dhdu(3, 4) = jac2(3, 4);
1010 
1011  dhdu(4, 0) = jac2(4, 0);
1012  dhdu(4, 1) = jac2(4, 1) * uu11 + jac2(4, 2) * uu12 + jac2(4, 3) * uu13;
1013  dhdu(4, 2) = jac2(4, 1) * uu12 + jac2(4, 2) * uu22 + jac2(4, 3) * uu23;
1014  dhdu(4, 3) = jac2(4, 1) * uu13 + jac2(4, 2) * uu23 + jac2(4, 3) * uu33;
1015  dhdu(4, 4) = jac2(4, 4);
1016 
1017  jac1 = dhdu;
1018  }
1019  else
1020  {
1022 
1023  Real omeg1 = 1. / (s1 - sm);
1024  Real omeg2 = 1. / (s2 - sm);
1025  Real prsta = rho1 * dsv1 * (sm - vdon1) + pres1;
1026 
1027  mooseError("Weird wave speed occured in ",
1028  name(),
1029  ": ",
1030  __FUNCTION__,
1031  "\n",
1032  "iside = ",
1033  iside,
1034  "\n",
1035  "ielem = ",
1036  ielem,
1037  "\n",
1038  "rho1 = ",
1039  rho1,
1040  "\n",
1041  "rhou1 = ",
1042  rhou1,
1043  "\n",
1044  "rhov1 = ",
1045  rhov1,
1046  "\n",
1047  "rhow1 = ",
1048  rhow1,
1049  "\n",
1050  "rhoe1 = ",
1051  rhoe1,
1052  "\n",
1053  "pres1 = ",
1054  pres1,
1055  "\n",
1056  "enth1 = ",
1057  enth1,
1058  "\n",
1059  "csou1 = ",
1060  csou1,
1061  "\n",
1062  "rho2 = ",
1063  rho2,
1064  "\n",
1065  "rhou2 = ",
1066  rhou2,
1067  "\n",
1068  "rhov2 = ",
1069  rhov2,
1070  "\n",
1071  "rhoe2 = ",
1072  rhoe2,
1073  "\n",
1074  "pres2 = ",
1075  pres2,
1076  "\n",
1077  "enth2 = ",
1078  enth2,
1079  "\n",
1080  "csou2 = ",
1081  csou2,
1082  "\n",
1083  "vdon1 = ",
1084  vdon1,
1085  "\n",
1086  "vdon2 = ",
1087  vdon2,
1088  "\n",
1089  "vnave = ",
1090  vnave,
1091  "\n",
1092  "cssav = ",
1093  cssav,
1094  "\n",
1095  "s1 = ",
1096  s1,
1097  "\n",
1098  "s2 = ",
1099  s2,
1100  "\n",
1101  "sm = ",
1102  sm,
1103  "\n",
1104  "omeg1 = ",
1105  omeg1,
1106  "\n",
1107  "omeg2 = ",
1108  omeg2,
1109  "\n",
1110  "prsta = ",
1111  prsta,
1112  "\n",
1113  "Please check before continuing!\n");
1114  }
1115 }
InputParameters validParams< BoundaryFluxBase >()
virtual Real c(Real v, Real u) const =0
Sound speed.
InputParameters validParams< CNSFVHLLCSlipBoundaryFlux >()
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.
A base class for computing/caching fluxes at boundaries.
const SinglePhaseFluidProperties & _fp
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.
Common class for single phase fluid properties.
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
CNSFVHLLCSlipBoundaryFlux(const InputParameters &parameters)
virtual Real gamma(Real v, Real u) const
Compute the ratio of specific heats.
virtual Real pressure(Real v, Real u) const =0
Pressure as a function of specific internal energy and specific volume.