www.mooseframework.org
Public Member Functions | Protected Attributes | List of all members
CNSFVHLLCSlipBoundaryFlux Class Reference

A user object that computes the slip boundary flux using the HLLC approximate Riemann solver. More...

#include <CNSFVHLLCSlipBoundaryFlux.h>

Inheritance diagram for CNSFVHLLCSlipBoundaryFlux:
[legend]

Public Member Functions

 CNSFVHLLCSlipBoundaryFlux (const InputParameters &parameters)
 
virtual ~CNSFVHLLCSlipBoundaryFlux ()
 
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. More...
 
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. More...
 
virtual void execute ()
 
virtual void initialize ()
 
virtual void finalize ()
 
virtual const std::vector< Real > & getFlux (unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, THREAD_ID tid) const
 Get the boundary flux vector. More...
 
virtual const DenseMatrix< Real > & getJacobian (unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, THREAD_ID tid) const
 Get the boundary Jacobian matrix. More...
 

Protected Attributes

const BCUserObject_bc_uo
 
const SinglePhaseFluidProperties_fp
 
unsigned int _cached_side_id
 
dof_id_type _cached_elem_id
 
std::vector< std::vector< Real > > _flux
 Threaded storage for fluxes. More...
 
std::vector< DenseMatrix< Real > > _jac1
 Threaded storage for jacobians. More...
 

Detailed Description

A user object that computes the slip boundary flux using the HLLC approximate Riemann solver.

Definition at line 24 of file CNSFVHLLCSlipBoundaryFlux.h.

Constructor & Destructor Documentation

CNSFVHLLCSlipBoundaryFlux::CNSFVHLLCSlipBoundaryFlux ( const InputParameters &  parameters)

Definition at line 27 of file CNSFVHLLCSlipBoundaryFlux.C.

28  : BoundaryFluxBase(parameters),
29  _bc_uo(getUserObject<BCUserObject>("bc_uo")),
30  _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties"))
31 {
32 }
const SinglePhaseFluidProperties & _fp
BoundaryFluxBase(const InputParameters &parameters)
CNSFVHLLCSlipBoundaryFlux::~CNSFVHLLCSlipBoundaryFlux ( )
virtual

Definition at line 34 of file CNSFVHLLCSlipBoundaryFlux.C.

34 {}

Member Function Documentation

void CNSFVHLLCSlipBoundaryFlux::calcFlux ( unsigned int  iside,
dof_id_type  ielem,
const std::vector< Real > &  uvec1,
const RealVectorValue &  dwave,
std::vector< Real > &  flux 
) const
virtual

Solve the Riemann problem on the boundary face.

Parameters
[in]isidelocal index of current side
[in]ielemglobal index of the current element
[in]uvec1vector of variables on the host side
[in]dwavevector of unit normal
[out]fluxflux vector for conservation equations

pass the inputs to local

assign the proper size for the flux vector

derived variables on the left

status in the ghost cell

get the so-called Roe averaged variables

averaged speed of sound

compute the eigenvalues at the left, right, and Roe states

get the S_L and S_R defined by Eq. (13).

compute the S_M defined by Eq. (12).

compute the Omega_l, Omega_r, and p^*

compute the U_l^, U_r^

compute the fluxes according to the wave speed

Implements BoundaryFluxBase.

Definition at line 37 of file CNSFVHLLCSlipBoundaryFlux.C.

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 }
virtual Real c(Real v, Real u) const =0
Sound speed.
const SinglePhaseFluidProperties & _fp
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
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.
void CNSFVHLLCSlipBoundaryFlux::calcJacobian ( unsigned int  iside,
dof_id_type  ielem,
const std::vector< Real > &  uvec1,
const RealVectorValue &  dwave,
DenseMatrix< Real > &  jac1 
) const
virtual

Compute the Jacobian matrix on the boundary face.

Parameters
[in]isidelocal index of current side
[in]ielemglobal index of the current element
[in]uvec1vector of variables on the host side
[in]dwavevector of unit normal
[out]jac1Jacobian matrix contribution

pass the inputs to local

assign the proper size for the Jacobian matrix

derived variables on the left

status in the ghost cell

get the so-called Roe averaged variables

averaged speed of sound

compute the eigenvalues at the left, right, and Roe states

get the S_L and S_R defined by Eq. (13).

compute the S_M defined by Eq. (12).

compute the flux Jacobians according to the wave speed

get the Jacobian matrix on the left only

compute the Omega_l and p^* by Eq. (10). and (11).

compute the U_l^ by Eq. (7).

compute rho^

compute a(S_M)/a(U_l) by Eq. (42).

compute a(S_M)/a(U_r) by Eq. (43).

compute a(p^*)/a(U_l) by Eq. (44).

compute a(p^*)/a(U_r) by Eq. (45).

compute a(rho_l^)/a(U_l) by Eq. (46).

compute a(rho_l^)/a(U_r) by Eq. (46a).

compute a(rhou_l^)/a(U_l) by Eq. (47).

compute a(rhou_l^)/a(U_r) by Eq. (48).

compute a(rhov_l^)/a(U_l) by Eq. (49).

compute a(rhov_l^)/a(U_r) by Eq. (50).

compute a(rhow_l^)/a(U_l) by Eq. (51).

compute a(rhow_l^)/a(U_r) by Eq. (52).

compute a(rhoe_l^)/a(U_l) by Eq. (53). note that enth1 = (rhoe1 + pres1)*rhom1 and rhoepl = rhoel + prsta

compute a(rhoe_l^)/a(U_r) by Eq. (54).

compute the HLLC Jacobians a(F_l^)/a(U_l) by Eq. (40).

compute the HLLC Jacobians a(F_l^)/a(U_r) by Eq. (41).

compute d(U_r)/d(U_l) by slip wall BC

compute a(F^)/a(U_r) * d(U_r)/d(U_l) by Eq.(39)

compute the Omega_r, and p^* by Eq. (10). and (11).

compute the U_r^ by Eq. (8).

compute rho^

compute a(S_M)/a(U_l) by Eq. (42).

compute a(S_M)/a(U_r) by Eq. (43).

compute a(p^*)/a(U_l) by Eq. (44).

compute a(p^*)/a(U_r) by Eq. (45).

interchange subscripts l <–> r and L <–> R in Eq. (46) through (54) and get Eq. (46*) through (54*).

compute a(rho_r^)/a(U_r) by Eq. (46*).

compute a(rho_r^)/a(U_l) by Eq. (46a*).

compute a(rhou_r^)/a(U_r) by Eq. (47*).

compute a(rhou_r^)/a(U_l) by Eq. (48*).

compute a(rhov_r^)/a(U_r) by Eq. (49*).

compute a(rhov_r^)/a(U_l) by Eq. (50*).

compute a(rhow_r^)/a(U_r) by Eq. (51*).

compute a(rhow_r^)/a(U_l) by Eq. (52*).

compute a(rhoe_r^)/a(U_r) by Eq. (53*). note that enth2 = (rhoe2 + pres2)*rhom2 and rhoepr = rhoer + prsta

compute a(rhoe_r^)/a(U_l) by Eq. (54*).

compute the HLLC Jacobians a(F_r^)/a(U_r) by Eq. (40*).

compute the HLLC Jacobians a(F_r^)/a(U_l) by Eq. (41*).

compute d(U_r)/d(U_l) by slip wall BC

compute a(F^)/a(U_r) * d(U_r)/d(U_l) by Eq.(39)

get the jacobian matrix at the point on the right only

compute d(U_r)/d(U_l) by slip wall BC

compute a(F^)/a(U_r) * d(U_r)/d(U_l) by Eq.(39)

compute the Omega_l and p^* by Eq. (10). and (11).

Implements BoundaryFluxBase.

Definition at line 283 of file CNSFVHLLCSlipBoundaryFlux.C.

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 }
virtual Real c(Real v, Real u) const =0
Sound speed.
const SinglePhaseFluidProperties & _fp
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
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.
void BoundaryFluxBase::execute ( )
virtualinherited

Definition at line 36 of file BoundaryFluxBase.C.

37 {
38 }
void BoundaryFluxBase::finalize ( )
virtualinherited

Definition at line 41 of file BoundaryFluxBase.C.

42 {
43 }
const std::vector< Real > & BoundaryFluxBase::getFlux ( unsigned int  iside,
dof_id_type  ielem,
const std::vector< Real > &  uvec1,
const RealVectorValue &  dwave,
THREAD_ID  tid 
) const
virtualinherited

Get the boundary flux vector.

Parameters
[in]isidelocal index of current side
[in]ielemglobal index of the current element
[in]uvec1vector of variables on the host side
[in]dwavevector of unit normal

Definition at line 46 of file BoundaryFluxBase.C.

Referenced by CNSFVBC::computeQpResidual(), and AEFVBC::computeQpResidual().

51 {
52  Threads::spin_mutex::scoped_lock lock(_mutex);
53  if (_cached_elem_id != ielem || _cached_side_id != iside)
54  {
55  _cached_elem_id = ielem;
56  _cached_side_id = iside;
57 
58  calcFlux(iside, ielem, uvec1, dwave, _flux[tid]);
59  }
60  return _flux[tid];
61 }
unsigned int _cached_side_id
static Threads::spin_mutex _mutex
std::vector< std::vector< Real > > _flux
Threaded storage for fluxes.
dof_id_type _cached_elem_id
virtual void calcFlux(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, std::vector< Real > &flux) const =0
Solve the Riemann problem on the boundary face.
const DenseMatrix< Real > & BoundaryFluxBase::getJacobian ( unsigned int  iside,
dof_id_type  ielem,
const std::vector< Real > &  uvec1,
const RealVectorValue &  dwave,
THREAD_ID  tid 
) const
virtualinherited

Get the boundary Jacobian matrix.

Parameters
[in]isidelocal index of current side
[in]ielemglobal index of the current element
[in]uvec1vector of variables on the host side
[in]dwavevector of unit normal

Definition at line 64 of file BoundaryFluxBase.C.

Referenced by CNSFVBC::computeQpJacobian(), AEFVBC::computeQpJacobian(), and CNSFVBC::computeQpOffDiagJacobian().

69 {
70  Threads::spin_mutex::scoped_lock lock(_mutex);
71  if (_cached_elem_id != ielem || _cached_side_id != iside)
72  {
73  _cached_elem_id = ielem;
74  _cached_side_id = iside;
75 
76  calcJacobian(iside, ielem, uvec1, dwave, _jac1[tid]);
77  }
78  return _jac1[tid];
79 }
unsigned int _cached_side_id
static Threads::spin_mutex _mutex
std::vector< DenseMatrix< Real > > _jac1
Threaded storage for jacobians.
virtual void calcJacobian(unsigned int iside, dof_id_type ielem, const std::vector< Real > &uvec1, const RealVectorValue &dwave, DenseMatrix< Real > &jac1) const =0
Compute the Jacobian matrix on the boundary face.
dof_id_type _cached_elem_id
void BoundaryFluxBase::initialize ( )
virtualinherited

Definition at line 29 of file BoundaryFluxBase.C.

30 {
31  _cached_elem_id = 0;
32  _cached_side_id = libMesh::invalid_uint;
33 }
unsigned int _cached_side_id
dof_id_type _cached_elem_id

Member Data Documentation

const BCUserObject& CNSFVHLLCSlipBoundaryFlux::_bc_uo
protected

Definition at line 43 of file CNSFVHLLCSlipBoundaryFlux.h.

Referenced by calcFlux(), and calcJacobian().

dof_id_type BoundaryFluxBase::_cached_elem_id
mutableprotectedinherited
unsigned int BoundaryFluxBase::_cached_side_id
mutableprotectedinherited
std::vector<std::vector<Real> > BoundaryFluxBase::_flux
mutableprotectedinherited

Threaded storage for fluxes.

Definition at line 98 of file BoundaryFluxBase.h.

Referenced by BoundaryFluxBase::BoundaryFluxBase(), and BoundaryFluxBase::getFlux().

const SinglePhaseFluidProperties& CNSFVHLLCSlipBoundaryFlux::_fp
protected

Definition at line 44 of file CNSFVHLLCSlipBoundaryFlux.h.

Referenced by calcFlux(), and calcJacobian().

std::vector<DenseMatrix<Real> > BoundaryFluxBase::_jac1
mutableprotectedinherited

Threaded storage for jacobians.

Definition at line 101 of file BoundaryFluxBase.h.

Referenced by BoundaryFluxBase::BoundaryFluxBase(), and BoundaryFluxBase::getJacobian().


The documentation for this class was generated from the following files: