MAST
fluid_elem_base.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // C++ includes
21 #include <iomanip>
22 
23 // MAST includes
24 #include "fluid/fluid_elem_base.h"
25 #include "fluid/primitive_fluid_solution.h"
26 #include "fluid/small_disturbance_primitive_fluid_solution.h"
27 #include "fluid/flight_condition.h"
28 #include "mesh/fe_base.h"
29 
30 // Basic include files
31 #include "libmesh/mesh.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/quadrature.h"
34 #include "libmesh/string_to_enum.h"
35 #include "libmesh/parameters.h"
36 
37 
38 
39 MAST::FluidElemBase::FluidElemBase(const unsigned int d,
40  const MAST::FlightCondition& f):
41 _if_viscous(f.gas_property.if_viscous),
42 _include_pressure_switch(false),
43 flight_condition(&f),
44 dim(d),
45 _dissipation_scaling(1.) {
46 
47 
48  // prepare the variable vector
49  _active_primitive_vars.push_back(RHO_PRIM);
50  _active_primitive_vars.push_back(VEL1);
51  _active_conservative_vars.push_back(RHO_CONS);
52  _active_conservative_vars.push_back(RHOVEL1);
53 
54  if (dim > 1)
55  {
56  _active_primitive_vars.push_back(VEL2);
57  _active_conservative_vars.push_back(RHOVEL2);
58  }
59 
60  if (dim > 2)
61  {
62  _active_primitive_vars.push_back(VEL3);
63  _active_conservative_vars.push_back(RHOVEL3);
64  }
65  _active_primitive_vars.push_back(TEMP);
66  _active_conservative_vars.push_back(ETOT);
67 
68 }
69 
70 
71 
72 
73 MAST::FluidElemBase::~FluidElemBase() {
74 
75 
76 }
77 
78 
79 
80 
81 
82 void MAST::FluidElemBase::
83 get_infinity_vars( RealVectorX& vars_inf ) const {
84 
85  vars_inf(0) = flight_condition->rho();
86 
87  vars_inf(1) = flight_condition->rho_u1();
88 
89  if (dim > 1)
90  vars_inf(2) = flight_condition->rho_u2();
91 
92  if (dim > 2)
93  vars_inf(3) = flight_condition->rho_u3();
94 
95  vars_inf(dim+2-1) = flight_condition->rho_e();
96 }
97 
98 
99 
100 
101 void
102 MAST::FluidElemBase::
103 update_solution_at_quadrature_point(const unsigned int qp,
104  const MAST::FEBase& fe,
105  const RealVectorX& elem_solution,
106  RealVectorX& conservative_sol,
107  PrimitiveSolution& primitive_sol,
108  MAST::FEMOperatorMatrix& B_mat,
109  std::vector<MAST::FEMOperatorMatrix>& dB_mat) {
110 
111  conservative_sol.setZero();
112 
113 
114  const std::vector<std::vector<Real> >& phi = fe.get_phi();
115  const unsigned int n_phi = (unsigned int)phi.size();
116 
117  RealVectorX
118  phi_vals = RealVectorX::Zero(n_phi);
119 
120  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
121  phi_vals(i_nd) = phi[i_nd][qp];
122 
123  B_mat.reinit(dim+2, phi_vals); // initialize the operator matrix
124 
125  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi =
126  fe.get_dphi();
127 
128 
129  for ( unsigned int i_dim=0; i_dim<dim; i_dim++ )
130  {
131  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
132  phi_vals(i_nd) = dphi[i_nd][qp](i_dim);
133  dB_mat[i_dim].reinit(dim+2, phi_vals);
134  }
135 
136  B_mat.vector_mult( conservative_sol, elem_solution );
137 
138  primitive_sol.zero();
139  primitive_sol.init(dim,
140  conservative_sol,
141  flight_condition->gas_property.cp,
142  flight_condition->gas_property.cv,
143  _if_viscous);
144 }
145 
146 
147 
148 
149 void
150 MAST::FluidElemBase::calculate_advection_flux(const unsigned int calculate_dim,
151  const MAST::PrimitiveSolution& sol,
152  RealVectorX& flux) {
153 
154  const unsigned int n1 = 2 + dim;
155 
156  const Real rho = sol.rho,
157  u1 = sol.u1,
158  u2 = sol.u2,
159  u3 = sol.u3,
160  p = sol.p,
161  e_tot = sol.e_tot;
162 
163  flux.setZero();
164 
165  // calculate the flux using given flow parameters at this point
166  switch (calculate_dim)
167  {
168  case 0:
169  {
170  flux(0) = rho * u1;
171  flux(n1-1) = u1 * (rho * e_tot + p);
172  switch (dim)
173  {
174  case 3:
175  flux(3) = rho * u1 * u3;
176  case 2:
177  flux(2) = rho * u1 * u2;
178  case 1:
179  flux(1) = rho * u1 * u1 + p;
180  }
181  }
182  break;
183 
184  case 1:
185  {
186  flux(0) = rho * u2;
187  flux(n1-1) = u2 * (rho * e_tot + p);
188  switch (dim)
189  {
190  case 3:
191  flux(3) = rho * u2 * u3;
192  case 2:
193  flux(2) = rho * u2 * u2 + p;
194  case 1:
195  flux(1) = rho * u2 * u1;
196  }
197  }
198  break;
199 
200  case 2:
201  {
202  flux(0) = rho * u3;
203  flux(n1-1) = u3 * (rho * e_tot + p);
204  switch (dim)
205  {
206  case 3:
207  flux(3) = rho * u3 * u3 + p;
208  case 2:
209  flux(2) = rho * u3 * u2;
210  case 1:
211  flux(1) = rho * u3 * u1;
212  }
213  }
214  break;
215  }
216 }
217 
218 
219 
220 
221 void
222 MAST::FluidElemBase::calculate_diffusion_flux(const unsigned int calculate_dim,
223  const MAST::PrimitiveSolution& sol,
224  const RealMatrixX& stress_tensor,
225  const RealVectorX& temp_gradient,
226  RealVectorX& flux) {
227 
228  const unsigned int n1 = 2 + dim;
229 
230  RealVectorX
231  uvec = RealVectorX::Zero(3);
232 
233  uvec(0) = sol.u1;
234  uvec(1) = sol.u2;
235  uvec(2) = sol.u3;
236 
237  flux.setZero();
238 
239  // the momentum flux is obtained from the stress tensor
240  for (unsigned int i_dim=0; i_dim<dim; i_dim++) {
241 
242  flux(1+i_dim) += stress_tensor(calculate_dim, i_dim); // tau_ij
243 
244  flux(n1-1) += uvec(i_dim) * stress_tensor(calculate_dim, i_dim); // u_j tau_ij
245  }
246 
247  flux(n1-1) += sol.k_thermal * temp_gradient(calculate_dim);
248 }
249 
250 
251 
252 
253 void
254 MAST::FluidElemBase::
255 calculate_diffusion_tensors(const RealVectorX& elem_sol,
256  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
257  const RealMatrixX& dprim_dcons,
258  const MAST::PrimitiveSolution& psol,
259  RealMatrixX& stress_tensor,
260  RealVectorX& temp_gradient) {
261 
262  const unsigned int n1 = dim+2;
263 
264  RealVectorX
265  dprim_dx = RealVectorX::Zero(n1),
266  dcons_dx = RealVectorX::Zero(n1);
267 
268  stress_tensor.setZero();
269  temp_gradient.setZero();
270 
271  for (unsigned int i_dim=0; i_dim<dim; i_dim++) {
272 
273  dB_mat[i_dim].vector_mult(dcons_dx, elem_sol); // dUcons/dx_i
274  dprim_dx = dprim_dcons * dcons_dx; // dUprim/dx_i
275 
276  for (unsigned int j_dim=0; j_dim<dim; j_dim++) {
277 
278  stress_tensor(i_dim, j_dim) += psol.mu * dprim_dx(j_dim+1); // mu * duj/dxi
279  stress_tensor(j_dim, i_dim) += psol.mu * dprim_dx(j_dim+1); // mu * dui/dxj
280 
281  stress_tensor(j_dim, j_dim) += psol.lambda * dprim_dx(i_dim+1); // + delta_ij lambda duk/dxk
282  }
283 
284  temp_gradient(i_dim) = dprim_dx(n1-1); // dT/dx_i
285  }
286 }
287 
288 
289 
290 void
291 MAST::FluidElemBase::
292 calculate_diffusion_tensors_sensitivity(const RealVectorX& elem_sol,
293  const RealVectorX& elem_sol_sens,
294  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
295  const RealMatrixX& dprim_dcons,
296  const MAST::PrimitiveSolution& psol,
297  const MAST::SmallPerturbationPrimitiveSolution<Real>& dsol,
298  RealMatrixX& stress_tensor_sens,
299  RealVectorX& temp_gradient_sens) {
300 
301  const unsigned int n1 = dim+2;
302 
303  RealVectorX
304  dprim_dx = RealVectorX::Zero(n1),
305  dcons_dx = RealVectorX::Zero(n1),
306  dprim_sens_dx = RealVectorX::Zero(n1),
307  dcons_sens_dx = RealVectorX::Zero(n1);
308 
309  stress_tensor_sens.setZero();
310  temp_gradient_sens.setZero();
311 
312  for (unsigned int i_dim=0; i_dim<dim; i_dim++) {
313 
314  dB_mat[i_dim].vector_mult(dcons_dx, elem_sol); // dUcons/dx_i
315  dprim_dx = dprim_dcons * dcons_dx; // dUprim/dx_i
316 
317  dB_mat[i_dim].vector_mult(dcons_sens_dx, elem_sol_sens); // dUcons/dx_i
318  dprim_sens_dx = dprim_dcons * dcons_sens_dx; // dUprim/dx_i
319 
320  for (unsigned int j_dim=0; j_dim<dim; j_dim++) {
321 
322  stress_tensor_sens(i_dim, j_dim) += (dsol.dmu * dprim_dx(j_dim+1) +
323  psol.mu * dprim_sens_dx(j_dim+1)); // mu * duj/dxi
324  stress_tensor_sens(j_dim, i_dim) += (dsol.dmu * dprim_dx(j_dim+1) +
325  psol.mu * dprim_sens_dx(j_dim+1)); // mu * dui/dxj
326 
327  stress_tensor_sens(j_dim, j_dim) += (dsol.dlambda * dprim_dx(i_dim+1) +
328  psol.lambda * dprim_sens_dx(i_dim+1)); // + delta_ij lambda duk/dxk
329  }
330 
331  temp_gradient_sens(i_dim) = dprim_sens_dx(n1-1); // dT/dx_i
332  }
333 }
334 
335 
336 
337 
338 void
339 MAST::FluidElemBase::
340 calculate_conservative_variable_jacobian(const MAST::PrimitiveSolution& sol,
341  RealMatrixX& dcons_dprim,
342  RealMatrixX& dprim_dcons) {
343 
344 
345  // calculate Ai = d F_adv / d x_i, where F_adv is the Euler advection flux vector
346 
347  const unsigned int n1 = 2 + dim;
348 
349  dcons_dprim.setZero();
350  dprim_dcons.setZero();
351 
352  const Real u1 = sol.u1,
353  u2 = sol.u2,
354  u3 = sol.u3,
355  rho = sol.rho,
356  k = sol.k,
357  e_tot = sol.e_tot,
358  cv = flight_condition->gas_property.cv;
359 
360  switch (dim)
361  {
362  case 3:
363  {
364  dcons_dprim(3, 0) = u3;
365  dcons_dprim(3, 3) = rho;
366 
367  dcons_dprim(n1-1, 3) = rho*u3;
368  }
369 
370  case 2:
371  {
372  dcons_dprim(2, 0) = u2;
373  dcons_dprim(2, 2) = rho;
374 
375  dcons_dprim(n1-1, 2) = rho*u2;
376  }
377 
378  case 1:
379  {
380  dcons_dprim(0, 0) = 1.;
381 
382  dcons_dprim(1, 0) = u1;
383  dcons_dprim(1, 1) = rho;
384 
385  dcons_dprim(n1-1, 0) = e_tot;
386  dcons_dprim(n1-1, 1) = rho*u1;
387  dcons_dprim(n1-1, n1-1) = rho*cv;
388  }
389  }
390 
391  switch (dim)
392  {
393  case 3:
394  {
395  dprim_dcons(3, 0) = -u3/rho;
396  dprim_dcons(3, 3) = 1./rho;
397 
398  dprim_dcons(n1-1, 3) = -u3/cv/rho;
399  }
400 
401  case 2:
402  {
403  dprim_dcons(2, 0) = -u2/rho;
404  dprim_dcons(2, 2) = 1./rho;
405 
406  dprim_dcons(n1-1, 2) = -u2/cv/rho;
407  }
408 
409  case 1:
410  {
411  dprim_dcons(0, 0) = 1.;
412 
413  dprim_dcons(1, 0) = -u1/rho;
414  dprim_dcons(1, 1) = 1./rho;
415 
416  dprim_dcons(n1-1, 0) = (-e_tot+2*k)/cv/rho;
417  dprim_dcons(n1-1, 1) = -u1/cv/rho;
418  dprim_dcons(n1-1, n1-1) = 1./cv/rho;
419  }
420  }
421 
422 }
423 
424 
425 
426 void
427 MAST::FluidElemBase::
428 calculate_advection_flux_jacobian(const unsigned int calculate_dim,
429  const MAST::PrimitiveSolution& sol,
430  RealMatrixX& mat) {
431 
432 
433  // calculate Ai = d F_adv / d x_i, where F_adv is the Euler advection flux vector
434 
435  const unsigned int n1 = 2 + dim;
436 
437  mat.setZero();
438 
439  const Real u1 = sol.u1,
440  u2 = sol.u2,
441  u3 = sol.u3,
442  k = sol.k,
443  e_tot = sol.e_tot,
444  T = sol.T,
445  gamma = flight_condition->gas_property.gamma,
446  R = flight_condition->gas_property.R,
447  cv = flight_condition->gas_property.cv;
448 
449  switch (calculate_dim)
450  {
451  case 0:
452  {
453  switch (dim)
454  {
455  case 3:
456  {
457  mat(1, 3) = -u3*R/cv;
458 
459  mat(3, 0) = -u1*u3;
460  mat(3, 1) = u3;
461  mat(3, 3) = u1;
462 
463  mat(n1-1, 3) = -u1*u3*R/cv;
464  }
465 
466  case 2:
467  {
468  mat(1, 2) = -u2*R/cv;
469 
470  mat(2, 0) = -u1*u2;
471  mat(2, 1) = u2;
472  mat(2, 2) = u1;
473 
474  mat(n1-1, 2) = -u1*u2*R/cv;
475  }
476 
477  case 1:
478  {
479  mat(0, 1) = 1.0; // d U / d (rho u1)
480 
481  mat(1, 0) = -u1*u1+R*k/cv;
482  mat(1, 1) = u1*(2.0-R/cv);
483  mat(1, n1-1) = R/cv;
484 
485  mat(n1-1, 0) = u1*(R*(-e_tot+2.0*k)-e_tot*cv)/cv;
486  mat(n1-1, 1) = e_tot+R*(T-u1*u1/cv);
487  mat(n1-1, n1-1) = u1*gamma;
488  }
489  break;
490  }
491  }
492  break;
493 
494  case 1:
495  {
496  switch (dim)
497  {
498  case 3:
499  {
500  mat(2, 3) = -u3*R/cv;
501 
502  mat(3, 0) = -u2*u3;
503  mat(3, 2) = u3;
504  mat(3, 3) = u2;
505 
506  mat(n1-1, 3) = -u2*u3*R/cv;
507  }
508 
509  case 2:
510  {
511  mat(0, 2) = 1.0; // d U / d (rho u2)
512 
513  mat(1, 0) = -u1*u2;
514  mat(1, 1) = u2;
515  mat(1, 2) = u1;
516 
517  mat(2, 0) = -u2*u2+R*k/cv;
518  mat(2, 1) = -u1*R/cv;
519  mat(2, 2) = u2*(2.0-R/cv);
520  mat(2, n1-1) = R/cv;
521 
522  mat(n1-1, 0) = u2*(R*(-e_tot+2.0*k)-e_tot*cv)/cv;
523  mat(n1-1, 1) = -u1*u2*R/cv;
524  mat(n1-1, 2) = e_tot+R*(T-u2*u2/cv);
525  mat(n1-1, n1-1) = u2*gamma;
526  }
527  break;
528 
529  case 1:
530  // if second coordinate divergence is being asked for, then the element is atleast 2D
531  libmesh_assert_msg(false, "invalid dim");
532  break;
533  }
534  }
535  break;
536 
537  case 2:
538  {
539  mat(0, 3) = 1.0; // d U / d (rho u3)
540 
541  mat(1, 0) = -u1*u3;
542  mat(1, 1) = u3;
543  mat(1, 3) = u1;
544 
545  mat(2, 0) = -u2*u3;
546  mat(2, 2) = u3;
547  mat(2, 3) = u2;
548 
549  mat(3, 0) = -u3*u3+R*k/cv;
550  mat(3, 1) = -u1*R/cv;
551  mat(3, 2) = -u2*R/cv;
552  mat(3, 3) = u3*(2.0-R/cv);
553  mat(3, n1-1) = R/cv;
554 
555  mat(n1-1, 0) = u3*(R*(-e_tot+2.0*k)-e_tot*cv)/cv;
556  mat(n1-1, 1) = -u1*u3*R/cv;
557  mat(n1-1, 2) = -u2*u3*R/cv;
558  mat(n1-1, 3) = e_tot+R*(T-u3*u3/cv);
559  mat(n1-1, n1-1) = u3*gamma;
560  }
561  break;
562 
563  default:
564  libmesh_assert_msg(false, "invalid dim");
565  break;
566  }
567 }
568 
569 
570 
571 
572 
573 void
574 MAST::FluidElemBase::
575 calculate_diffusion_flux_jacobian_primitive_vars (const unsigned int flux_dim,
576  const unsigned int deriv_dim,
577  const RealVectorX& uvec,
578  const bool zero_kth,
579  const MAST::PrimitiveSolution& sol,
580  RealMatrixX& mat) {
581 
582  const unsigned int n1 = 2 + dim;
583 
584  mat.setZero();
585 
586  const Real
587  u1 = uvec(0),
588  u2 = dim>1?uvec(1):0.,
589  u3 = dim>2?uvec(2):0.,
590  mu = sol.mu,
591  lambda = sol.lambda,
592  kth = zero_kth?0.:sol.k_thermal;
593 
594  switch (flux_dim)
595  {
596  case 0:
597  {
598  switch (deriv_dim)
599  {
600  case 0: // K11
601  {
602  switch (dim)
603  {
604  case 3:
605  {
606  mat(3,3) = mu;
607 
608  mat(n1-1,3) = u3*mu;
609  }
610 
611  case 2:
612  {
613  mat(2,2) = mu;
614 
615  mat(n1-1,2) = u2*mu;
616  }
617 
618  case 1:
619  {
620  mat(1,1) = (lambda+2.*mu);
621 
622  mat(n1-1,1) = u1*(lambda+2.*mu);
623  mat(n1-1,n1-1) = kth;
624  }
625  }
626  }
627  break;
628 
629  case 1: // K12
630  {
631  mat(1,2) = lambda;
632 
633  mat(2,1) = mu;
634 
635  mat(n1-1,1) = u2*mu;
636  mat(n1-1,2) = u1*lambda;
637  }
638  break;
639 
640  case 2: // K13
641  {
642  mat(1,3) = lambda;
643 
644  mat(3,1) = mu;
645 
646  mat(n1-1,1) = u3*mu;
647  mat(n1-1,3) = u1*lambda;
648  }
649  break;
650  }
651  }
652  break;
653 
654  case 1:
655  {
656  switch (deriv_dim)
657  {
658  case 0: // K21
659  {
660  mat(1,2) = mu;
661 
662  mat(2,1) = lambda;
663 
664  mat(n1-1,1) = u2*lambda;
665  mat(n1-1,2) = u1*mu;
666  }
667  break;
668 
669  case 1: // K22
670  {
671  switch (dim)
672  {
673  case 3:
674  {
675  mat(3,3) = mu;
676 
677  mat(n1-1,3) = u3*mu;
678  }
679 
680  case 2:
681  case 1:
682  {
683  mat(1,1) = mu;
684 
685  mat(2,2) = (lambda+2.*mu);
686 
687  mat(n1-1,1) = u1*mu;
688  mat(n1-1,2) = u2*(lambda+2.*mu);
689  mat(n1-1,n1-1) = kth;
690  }
691  }
692  }
693  break;
694 
695  case 2: // K23
696  {
697  mat(2,3) = lambda;
698 
699  mat(3,2) = mu;
700 
701  mat(n1-1,2) = u3*mu;
702  mat(n1-1,3) = u2*lambda;
703  }
704  break;
705  }
706  }
707  break;
708 
709  case 2:
710  {
711  switch (deriv_dim)
712  {
713  case 0: // K31
714  {
715  mat(1,3) = mu;
716 
717  mat(3,1) = lambda;
718 
719  mat(n1-1,1) = u3*lambda;
720  mat(n1-1,3) = u1*mu;
721  }
722  break;
723 
724  case 1: // K32
725  {
726  mat(2,3) = mu;
727 
728  mat(3,2) = lambda;
729 
730  mat(n1-1,2) = u3*lambda;
731  mat(n1-1,3) = u2*mu;
732  }
733  break;
734 
735  case 2: // K33
736  {
737  mat(1,1) = mu;
738 
739  mat(2,2) = mu;
740 
741  mat(3,3) = (lambda+2.*mu);
742 
743  mat(n1-1,1) = u1*mu;
744  mat(n1-1,2) = u2*mu;
745  mat(n1-1,3) = u3*(lambda+2.*mu);
746  mat(n1-1,n1-1) = kth;
747  }
748  break;
749  }
750  }
751  break;
752  }
753 }
754 
755 
756 
757 void
758 MAST::FluidElemBase::
759 calculate_diffusion_flux_jacobian (const unsigned int flux_dim,
760  const unsigned int deriv_dim,
761  const MAST::PrimitiveSolution& sol,
762  RealMatrixX& mat) {
763 
764  const unsigned int n1 = 2 + dim;
765 
766  mat.setZero();
767 
768  const Real rho = sol.rho,
769  u1 = sol.u1,
770  u2 = sol.u2,
771  u3 = sol.u3,
772  k = sol.k,
773  e_tot = sol.e_tot,
774  mu = sol.mu,
775  lambda = sol.lambda,
776  kth = sol.k_thermal,
777  cv = flight_condition->gas_property.cv;
778 
779  switch (flux_dim)
780  {
781  case 0:
782  {
783  switch (deriv_dim)
784  {
785  case 0: // K11
786  {
787  switch (dim)
788  {
789  case 3:
790  {
791  mat(3,0) = -u3*mu/rho;
792  mat(3,3) = mu/rho;
793 
794  mat(n1-1,3) = u3*(-kth+cv*mu)/cv/rho;
795  }
796 
797  case 2:
798  {
799  mat(2,0) = -u2*mu/rho;
800  mat(2,2) = mu/rho;
801 
802  mat(n1-1,2) = u2*(-kth+cv*mu)/cv/rho;
803  }
804 
805  case 1:
806  {
807  mat(1,0) = -u1*(lambda+2.*mu)/rho;
808  mat(1,1) = (lambda+2.*mu)/rho;
809 
810  mat(n1-1,0) = (kth*(2.*k-e_tot)-cv*(2.*k*mu+u1*u1*(mu+lambda)))/cv/rho;
811  mat(n1-1,1) = u1*(-kth+cv*(lambda+2.*mu))/cv/rho;
812  mat(n1-1,n1-1) = kth/cv/rho;
813  }
814  }
815  }
816  break;
817 
818  case 1: // K12
819  {
820  mat(1,0) = -u2*lambda/rho;
821  mat(1,2) = lambda/rho;
822 
823  mat(2,0) = -u1*mu/rho;
824  mat(2,1) = mu/rho;
825 
826  mat(n1-1,0) = -u1*u2*(lambda+mu)/rho;
827  mat(n1-1,1) = u2*mu/rho;
828  mat(n1-1,2) = u1*lambda/rho;
829  }
830  break;
831 
832  case 2: // K13
833  {
834  mat(1,0) = -u3*lambda/rho;
835  mat(1,3) = lambda/rho;
836 
837  mat(3,0) = -u1*mu/rho;
838  mat(3,1) = mu/rho;
839 
840  mat(n1-1,0) = -u1*u3*(lambda+mu)/rho;
841  mat(n1-1,1) = u3*mu/rho;
842  mat(n1-1,3) = u1*lambda/rho;
843  }
844  break;
845  }
846  }
847  break;
848 
849  case 1:
850  {
851  switch (deriv_dim)
852  {
853  case 0: // K21
854  {
855  mat(1,0) = -u2*mu/rho;
856  mat(1,2) = mu/rho;
857 
858  mat(2,0) = -u1*lambda/rho;
859  mat(2,1) = lambda/rho;
860 
861  mat(n1-1,0) = -u1*u2*(lambda+mu)/rho;
862  mat(n1-1,1) = u2*lambda/rho;
863  mat(n1-1,2) = u1*mu/rho;
864  }
865  break;
866 
867  case 1: // K22
868  {
869  switch (dim)
870  {
871  case 3:
872  {
873  mat(3,0) = -u3*mu/rho;
874  mat(3,3) = mu/rho;
875 
876  mat(n1-1,3) = u3*(-kth+cv*mu)/cv/rho;
877  }
878 
879  case 2:
880  case 1:
881  {
882  mat(1,0) = -u1*mu/rho;
883  mat(1,1) = mu/rho;
884 
885  mat(2,0) = -u2*(lambda+2.*mu)/rho;
886  mat(2,2) = (lambda+2.*mu)/rho;
887 
888  mat(n1-1,0) = (kth*(2.*k-e_tot)-cv*(2.*k*mu+u2*u2*(mu+lambda)))/cv/rho;
889  mat(n1-1,1) = u1*(-kth+cv*mu)/cv/rho;
890  mat(n1-1,2) = u2*(-kth+cv*(lambda+2.*mu))/cv/rho;
891  mat(n1-1,n1-1) = kth/cv/rho;
892  }
893  }
894  }
895  break;
896 
897  case 2: // K23
898  {
899  mat(2,0) = -u3*lambda/rho;
900  mat(2,3) = lambda/rho;
901 
902  mat(3,0) = -u2*mu/rho;
903  mat(3,2) = mu/rho;
904 
905  mat(n1-1,0) = -u2*u3*(lambda+mu)/rho;
906  mat(n1-1,2) = u3*mu/rho;
907  mat(n1-1,3) = u2*lambda/rho;
908  }
909  break;
910  }
911  }
912  break;
913 
914  case 2:
915  {
916  switch (deriv_dim)
917  {
918  case 0: // K31
919  {
920  mat(1,0) = -u3*mu/rho;
921  mat(1,3) = mu/rho;
922 
923  mat(3,0) = -u1*lambda/rho;
924  mat(3,1) = lambda/rho;
925 
926  mat(n1-1,0) = -u1*u3*(lambda+mu)/rho;
927  mat(n1-1,1) = u3*lambda/rho;
928  mat(n1-1,3) = u1*mu/rho;
929  }
930  break;
931 
932  case 1: // K32
933  {
934  mat(2,0) = -u3*mu/rho;
935  mat(2,3) = mu/rho;
936 
937  mat(3,0) = -u2*lambda/rho;
938  mat(3,2) = lambda/rho;
939 
940  mat(n1-1,0) = -u2*u3*(lambda+mu)/rho;
941  mat(n1-1,2) = u3*lambda/rho;
942  mat(n1-1,3) = u2*mu/rho;
943  }
944  break;
945 
946  case 2: // K33
947  {
948  mat(1,0) = -u1*mu/rho;
949  mat(1,1) = mu/rho;
950 
951  mat(2,0) = -u2*mu/rho;
952  mat(2,2) = mu/rho;
953 
954  mat(3,0) = -u3*(lambda+2.*mu)/rho;
955  mat(3,3) = (lambda+2.*mu)/rho;
956 
957  mat(n1-1,0) = (kth*(2.*k-e_tot)-cv*(2.*k*mu+u3*u3*(mu+lambda)))/cv/rho;
958  mat(n1-1,1) = u1*(-kth+cv*mu)/cv/rho;
959  mat(n1-1,2) = u2*(-kth+cv*mu)/cv/rho;
960  mat(n1-1,3) = u3*(-kth+cv*(lambda+2.*mu))/cv/rho;
961  mat(n1-1,n1-1) = kth/cv/rho;
962  }
963  break;
964  }
965  }
966  break;
967  }
968 }
969 
970 
971 
972 
973 void
974 MAST::FluidElemBase::
975 calculate_advection_flux_jacobian_sensitivity_for_conservative_variable
976 (const unsigned int calculate_dim,
977  const MAST::PrimitiveSolution& sol,
978  std::vector<RealMatrixX >& jac) {
979 
980  const unsigned int n1 = 2 + dim;
981 
982  RealMatrixX
983  dprim_dcons = RealMatrixX::Zero(n1, n1),
984  mat = RealMatrixX::Zero(n1, n1);
985 
986  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
987  jac[i_cvar].setZero();
988 
989  this->calculate_conservative_variable_jacobian(sol, mat, dprim_dcons);
990 
991  // calculate based on chain rule of the primary variables
992  for (unsigned int i_pvar=0; i_pvar<n1; i_pvar++) // iterate over the primitive variables for chain rule
993  {
994  this->calculate_advection_flux_jacobian_sensitivity_for_primitive_variable
995  (calculate_dim, i_pvar, sol, mat);
996  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
997  {
998  if (fabs(dprim_dcons(i_pvar, i_cvar)) > 0.0)
999  jac[i_cvar] += dprim_dcons(i_pvar, i_cvar) * mat;
1000  }
1001  }
1002 }
1003 
1004 
1005 
1006 void
1007 MAST::FluidElemBase::
1008 calculate_advection_flux_jacobian_sensitivity_for_primitive_variable
1009 (const unsigned int calculate_dim,
1010  const unsigned int primitive_var,
1011  const MAST::PrimitiveSolution& sol,
1012  RealMatrixX& mat) {
1013 
1014  // calculate Ai = d F_adv / d x_i, where F_adv is the Euler advection flux vector
1015 
1016  const unsigned int n1 = 2 + dim;
1017 
1018  mat.setZero();
1019 
1020  const Real u1 = sol.u1,
1021  u2 = sol.u2,
1022  u3 = sol.u3,
1023  k = sol.k,
1024  e_tot = sol.e_tot,
1025  R = flight_condition->gas_property.R,
1026  cv = flight_condition->gas_property.cv;
1027 
1028  switch (calculate_dim)
1029  {
1030  case 0: // Ax
1031  {
1032  switch (_active_primitive_vars[primitive_var])
1033  {
1034  case RHO_PRIM:
1035  {
1036  // nothing to be done for density; matrix is zero
1037  }
1038  break;
1039 
1040  case VEL1:
1041  {
1042  switch (dim)
1043  {
1044  case 3:
1045  {
1046  mat(3, 0) = -u3;
1047  mat(3, 3) = 1.0;
1048 
1049  mat(n1-1, 3) = -u3*R/cv;
1050  }
1051 
1052  case 2:
1053  {
1054  mat(2, 0) = -u2;
1055  mat(2, 2) = 1.0;
1056 
1057  mat(n1-1, 2) = -u2*R/cv;
1058  }
1059 
1060  case 1:
1061  {
1062  mat(1, 0) = -u1*(2.0-R/cv);
1063  mat(1, 1) = (2.0-R/cv);
1064 
1065  mat(n1-1, 0) = (-e_tot*(cv+R)+2.0*R*k + u1*u1*(-cv+R))/cv;
1066  mat(n1-1, 1) = u1*(1.0-2.0*R/cv);
1067  mat(n1-1, n1-1) = (cv+R)/cv;
1068  }
1069  break;
1070  }
1071  }
1072  break;
1073 
1074  case VEL2:
1075  {
1076  mat(1, 0) = u2*R/cv;
1077  mat(1, 2) = -R/cv;
1078 
1079  mat(2, 0) = -u1;
1080  mat(2, 1) = 1.0;
1081 
1082  mat(n1-1, 0) = (-cv+R)*u1*u2/cv;
1083  mat(n1-1, 1) = u2;
1084  mat(n1-1, 2) = -u1*R/cv;
1085  }
1086  break;
1087 
1088  case VEL3:
1089  {
1090  mat(1, 0) = u3*R/cv;
1091  mat(1, 3) = -R/cv;
1092 
1093  mat(3, 0) = -u1;
1094  mat(3, 1) = 1.0;
1095 
1096  mat(n1-1, 0) = (-cv+R)*u1*u3/cv;
1097  mat(n1-1, 1) = u3;
1098  mat(n1-1, 3) = -u1*R/cv;
1099  }
1100  break;
1101 
1102  case TEMP:
1103  {
1104  // all dimensions have the same format
1105  mat(n1-1, 0) = -u1*(cv+R);
1106  mat(n1-1, 1) = cv+R;
1107  }
1108  break;
1109 
1110  default:
1111  libmesh_assert_msg(false, "Invalid primitive variable number");
1112  break;
1113  }
1114  }
1115  break;
1116 
1117  case 1: // Ay
1118  {
1119  switch (_active_primitive_vars[primitive_var])
1120  {
1121  case RHO_PRIM:
1122  {
1123  // nothing to be done for density; matrix is zero
1124  }
1125  break;
1126 
1127  case VEL1:
1128  {
1129  mat(1, 0) = -u2;
1130  mat(1, 2) = 1.0;
1131 
1132  mat(2, 0) = u1*R/cv;
1133  mat(2, 1) = -R/cv;
1134 
1135 
1136  mat(n1-1, 0) = (-cv+R)*u1*u2/cv;
1137  mat(n1-1, 1) = -u2*R/cv;
1138  mat(n1-1, 2) = u1;
1139  }
1140  break;
1141 
1142  case VEL2:
1143  {
1144  switch (dim)
1145  {
1146  case 3:
1147  {
1148  mat(3, 0) = -u3;
1149  mat(3, 3) = 1.0;
1150 
1151  mat(n1-1, 3) = -u3*R/cv;
1152  }
1153 
1154  case 1:
1155  case 2:
1156  {
1157  mat(1, 0) = -u1;
1158  mat(1, 1) = 1.0;
1159 
1160  mat(2, 0) = -u2*(2.0-R/cv);
1161  mat(2, 2) = (2.0-R/cv);
1162 
1163  mat(n1-1, 0) = (-e_tot*(cv+R)+2.0*R*k + u2*u2*(-cv+R))/cv;
1164  mat(n1-1, 1) = -u1*R/cv;
1165  mat(n1-1, 2) = u2*(1.0-2.0*R/cv);
1166  mat(n1-1, n1-1) = (cv+R)/cv;
1167  }
1168  break;
1169  }
1170  }
1171  break;
1172 
1173  case VEL3:
1174  {
1175  mat(2, 0) = u3*R/cv;
1176  mat(2, 3) = -R/cv;
1177 
1178  mat(3, 0) = -u2;
1179  mat(3, 2) = 1.0;
1180 
1181  mat(n1-1, 0) = (-cv+R)*u2*u3/cv;
1182  mat(n1-1, 2) = u3;
1183  mat(n1-1, 3) = -u2*R/cv;
1184  }
1185  break;
1186 
1187  case TEMP:
1188  {
1189  // all dimensions have the same format
1190  mat(n1-1, 0) = -u2*(cv+R);
1191  mat(n1-1, 2) = cv+R;
1192  }
1193  break;
1194 
1195  default:
1196  libmesh_assert_msg(false, "Invalid primitive variable number");
1197  break;
1198  }
1199  }
1200  break;
1201 
1202  case 2: // Az
1203  {
1204  switch (_active_primitive_vars[primitive_var])
1205  {
1206  case RHO_PRIM:
1207  {
1208  // nothing to be done for density; matrix is zero
1209  }
1210  break;
1211 
1212  case VEL1:
1213  {
1214  mat(1, 0) = -u3;
1215  mat(1, 3) = 1.0;
1216 
1217  mat(3, 0) = u1*R/cv;
1218  mat(3, 1) = -R/cv;
1219 
1220  mat(n1-1, 0) = (-cv+R)*u1*u3/cv;
1221  mat(n1-1, 1) = -u3*R/cv;
1222  mat(n1-1, 3) = u1;
1223  }
1224  break;
1225 
1226  case VEL2:
1227  {
1228  mat(2, 0) = -u3;
1229  mat(2, 3) = 1.0;
1230 
1231  mat(3, 0) = u2*R/cv;
1232  mat(3, 2) = -R/cv;
1233 
1234  mat(n1-1, 0) = (-cv+R)*u2*u3/cv;
1235  mat(n1-1, 2) = -u3*R/cv;
1236  mat(n1-1, 3) = u2;
1237  }
1238  break;
1239 
1240  case VEL3:
1241  {
1242  mat(1, 0) = -u1;
1243  mat(1, 1) = 1.0;
1244 
1245  mat(2, 0) = -u2;
1246  mat(2, 2) = 1.0;
1247 
1248  mat(3, 0) = -u3*(2.0-R/cv);
1249  mat(3, 3) = (2.0-R/cv);
1250 
1251  mat(n1-1, 0) = (-e_tot*(cv+R)+2.0*R*k + u3*u3*(-cv+R))/cv;
1252  mat(n1-1, 1) = -u1*R/cv;
1253  mat(n1-1, 2) = -u2*R/cv;
1254  mat(n1-1, 3) = u3*(1.0-2.0*R/cv);
1255  mat(n1-1, n1-1) = (cv+R)/cv;
1256  }
1257  break;
1258 
1259  case TEMP:
1260  {
1261  // all dimensions have the same format
1262  mat(n1-1, 0) = -u3*(cv+R);
1263  mat(n1-1, 3) = cv+R;
1264  }
1265  break;
1266 
1267  default:
1268  libmesh_assert_msg(false, "Invalid primitive variable number");
1269  break;
1270  }
1271  }
1272  break;
1273 
1274  default:
1275  libmesh_assert_msg(false, "invalid dim");
1276  break;
1277  }
1278 }
1279 
1280 
1281 
1282 void
1283 MAST::FluidElemBase::
1284 calculate_advection_left_eigenvector_and_inverse_for_normal
1285 (const MAST::PrimitiveSolution& sol,
1286  const libMesh::Point& normal,
1287  RealVectorX& eig_vals,
1288  RealMatrixX& l_eig_mat,
1289  RealMatrixX& l_eig_mat_inv_tr) {
1290 
1291 
1292  const unsigned int n1 = 2 + dim;
1293 
1294  eig_vals.setZero(); l_eig_mat.setZero(); l_eig_mat_inv_tr.setZero();
1295 
1296  Real nx=0., ny=0., nz=0., u=0.;
1297  unsigned int dim_for_eig_vec=100; // initializing with arbitrarily high value
1298 
1299  const Real u1 = sol.u1,
1300  u2 = sol.u2,
1301  u3 = sol.u3,
1302  k = sol.k,
1303  T = sol.T,
1304  a = sol.a,
1305  gamma = flight_condition->gas_property.gamma,
1306  cv = flight_condition->gas_property.cv;
1307 
1308  // initialize the values
1309  switch (dim)
1310  {
1311  case 3:
1312  {
1313  nz = normal(2);
1314  u += u3*nz;
1315  }
1316 
1317  case 2:
1318  {
1319  ny = normal(1);
1320  u += u2*ny;
1321  }
1322 
1323  case 1:
1324  {
1325  nx = normal(0);
1326  u += u1*nx;
1327  }
1328  }
1329 
1330  // select the largest value of surface normal component, so that the appropriate section can be chosen
1331  if ((fabs(nx)>=fabs(ny)) && (fabs(nx)>=fabs(nz)))
1332  dim_for_eig_vec = 1;
1333  else if ((fabs(ny)>=fabs(nx)) && (fabs(ny)>=fabs(nz)))
1334  dim_for_eig_vec = 2;
1335  else if ((fabs(nz)>=fabs(nx)) && (fabs(nz)>=fabs(ny)))
1336  dim_for_eig_vec = 3;
1337 
1338  // set eigenvalues
1339  switch (dim)
1340  {
1341  case 3:
1342  eig_vals(2) = u;
1343  case 2:
1344  eig_vals(1) = u;
1345  case 1:
1346  {
1347  eig_vals(0) = u;
1348  eig_vals(n1-2) = u-a;
1349  eig_vals(n1-1) = u+a;
1350  }
1351  }
1352 
1353  // set last two columns of the eigenvector matrices, and all columns of the eigenvector inverse.
1354  // Note that the dim_for_eig_vec column of the inverse matrix will be overwritten by the appropriate matrix
1355  switch (dim)
1356  {
1357  case 3:
1358  {
1359  // for u-a
1360  l_eig_mat(3, n1-2) = -u3-nz*a/(gamma-1.0);
1361 
1362  // for u+a
1363  l_eig_mat(3, n1-1) = -u3+nz*a/(gamma-1.0);
1364 
1365  // for u
1366  l_eig_mat_inv_tr(3, 0) = (gamma-1.0)*u1*u3/(a*a)-nx*nz;
1367 
1368  // for u
1369  l_eig_mat_inv_tr(3, 1) = (gamma-1.0)*u2*u3/(a*a)-ny*nz;
1370 
1371  // for u
1372  l_eig_mat_inv_tr(0, 2) = (gamma-1.0)*u3/(a*a);
1373  l_eig_mat_inv_tr(1, 2) = (gamma-1.0)*u3*u1/(a*a)-nz*nx;
1374  l_eig_mat_inv_tr(2, 2) = (gamma-1.0)*u3*u2/(a*a)-nz*ny;
1375  l_eig_mat_inv_tr(3, 2) = (gamma-1.0)*u3*u3/(a*a)+1.0-nz*nz;
1376  l_eig_mat_inv_tr(n1-1, 2) = (gamma-1.0)*k*u3/(a*a)+u3-nz*u;
1377 
1378  // overwrite for u
1379  l_eig_mat_inv_tr(3, dim_for_eig_vec-1) = u3;
1380 
1381  // for u-a
1382  l_eig_mat_inv_tr(3, n1-2) = 0.5*(gamma-1.0)*(-nz+u3/a)/a;
1383 
1384  // for u+a
1385  l_eig_mat_inv_tr(3, n1-1) = 0.5*(gamma-1.0)*(nz+u3/a)/a;
1386 
1387  }
1388 
1389  case 2:
1390  {
1391  // for u-a
1392  l_eig_mat(2, n1-2) = -u2-ny*a/(gamma-1.0);
1393 
1394  // for u+a
1395  l_eig_mat(2, n1-1) = -u2+ny*a/(gamma-1.0);
1396 
1397  // for u
1398  l_eig_mat_inv_tr(2, 0) = (gamma-1.0)*u1*u2/(a*a)-nx*ny;
1399 
1400  // for u
1401  l_eig_mat_inv_tr(0, 1) = (gamma-1.0)*u2/(a*a);
1402  l_eig_mat_inv_tr(1, 1) = (gamma-1.0)*u2*u1/(a*a)-ny*nx;
1403  l_eig_mat_inv_tr(2, 1) = (gamma-1.0)*u2*u2/(a*a)+1.0-ny*ny;
1404  l_eig_mat_inv_tr(n1-1, 1) = (gamma-1.0)*k*u2/(a*a)+u2-ny*u;
1405 
1406  // overwrite for u
1407  l_eig_mat_inv_tr(2, dim_for_eig_vec-1) = u2;
1408 
1409  // for u-a
1410  l_eig_mat_inv_tr(2, n1-2) = 0.5*(gamma-1.0)*(-ny+u2/a)/a;
1411 
1412  // for u+a
1413  l_eig_mat_inv_tr(2, n1-1) = 0.5*(gamma-1.0)*(ny+u2/a)/a;
1414  }
1415 
1416  case 1:
1417  {
1418  // for u-a
1419  l_eig_mat(0, n1-2) = k+u*a/(gamma-1.0);
1420  l_eig_mat(1, n1-2) = -u1-nx*a/(gamma-1.0);
1421  l_eig_mat(n1-1, n1-2) = 1.0;
1422 
1423  // for u+a
1424  l_eig_mat(0, n1-1) = k-u*a/(gamma-1.0);
1425  l_eig_mat(1, n1-1) = -u1+nx*a/(gamma-1.0);
1426  l_eig_mat(n1-1, n1-1) = 1.0;
1427 
1428  // for u
1429  l_eig_mat_inv_tr(0, 0) = (gamma-1.0)*u1/(a*a);
1430  l_eig_mat_inv_tr(1, 0) = (gamma-1.0)*u1*u1/(a*a)+1.0-nx*nx;
1431  l_eig_mat_inv_tr(n1-1, 0) = (gamma-1.0)*k*u1/(a*a)+u1-nx*u;
1432 
1433  // overwrite for u
1434  l_eig_mat_inv_tr(0, dim_for_eig_vec-1) = 1.0;
1435  l_eig_mat_inv_tr(1, dim_for_eig_vec-1) = u1;
1436  l_eig_mat_inv_tr(n1-1, dim_for_eig_vec-1) = k;
1437  l_eig_mat_inv_tr.col(dim_for_eig_vec-1) *= -1.0/(cv*T*gamma);
1438 
1439  // for u-a
1440  l_eig_mat_inv_tr(0, n1-2) = 1.0/(2.0*cv*T*gamma);
1441  l_eig_mat_inv_tr(1, n1-2) = 0.5*(gamma-1.0)*(-nx+u1/a)/a;
1442  l_eig_mat_inv_tr(n1-1, n1-2) = 0.5*(1.0+(gamma-1.0)*(-u+k/a)/a);
1443 
1444  // for u+a
1445  l_eig_mat_inv_tr(0, n1-1) = 1.0/(2.0*cv*T*gamma);
1446  l_eig_mat_inv_tr(1, n1-1) = 0.5*(gamma-1.0)*(nx+u1/a)/a;
1447  l_eig_mat_inv_tr(n1-1, n1-1) = 0.5*(1.0+(gamma-1.0)*(u+k/a)/a);
1448 
1449  }
1450  }
1451 
1452 
1453  switch (dim_for_eig_vec)
1454  {
1455  case 1:
1456  {
1457  // set values in the left eigenvector matrix and eigenvalue matrix
1458  switch (dim)
1459  {
1460  case 3:
1461  {
1462  // for u
1463  l_eig_mat(0, 2) = nz*u1/nx-u3;
1464  l_eig_mat(1, 2) = -nz/nx;
1465  l_eig_mat(3, 2) = 1.0;
1466  }
1467 
1468  case 2:
1469  {
1470  // for u
1471  l_eig_mat(0, 1) = ny*u1/nx-u2;
1472  l_eig_mat(1, 1) = -ny/nx;
1473  l_eig_mat(2, 1) = 1.0;
1474  }
1475 
1476  case 1:
1477  {
1478  // for u
1479  l_eig_mat(0, 0) = u*u1/nx-(cv*T*gamma+k);
1480  l_eig_mat(1, 0) = -u/nx;
1481  l_eig_mat(n1-1, 0) = 1.0;
1482  }
1483  }
1484  }
1485  break;
1486 
1487  case 2:
1488  {
1489  // set values in the left eigenvector matrix and eigenvalue matrix
1490  switch (dim)
1491  {
1492  case 3:
1493  {
1494  // for u
1495  l_eig_mat(0, 2) = nz*u2/ny-u3;
1496  l_eig_mat(2, 2) = -nz/ny;
1497  l_eig_mat(3, 2) = 1.0;
1498  }
1499 
1500  case 2:
1501  case 1:
1502  {
1503  // for u
1504  l_eig_mat(0, 0) = nx*u2/ny-u1;
1505  l_eig_mat(1, 0) = 1.0;
1506  l_eig_mat(2, 0) = -nx/ny;
1507 
1508  // for u
1509  l_eig_mat(0, 1) = u*u2/ny-(cv*T*gamma+k);
1510  l_eig_mat(2, 1) = -u/ny;
1511  l_eig_mat(n1-1, 1) = 1.0;
1512  }
1513  }
1514  }
1515  break;
1516 
1517  case 3:
1518  {
1519  // set values in the left eigenvector matrix and eigenvalue matrix
1520 
1521  // for u
1522  l_eig_mat(0, 0) = nx*u3/nz-u1;
1523  l_eig_mat(1, 0) = 1.0;
1524  l_eig_mat(3, 0) = -nx/nz;
1525 
1526  // for u
1527  l_eig_mat(0, 1) = ny*u3/nz-u2;
1528  l_eig_mat(2, 1) = 1.0;
1529  l_eig_mat(3, 1) = -ny/nz;
1530 
1531  // for u
1532  l_eig_mat(0, 2) = u*u3/nz-(cv*T*gamma+k);
1533  l_eig_mat(3, 2) = -u/nz;
1534  l_eig_mat(n1-1, 2) = 1.0;
1535  }
1536  break;
1537  }
1538 }
1539 
1540 
1541 
1542 void
1543 MAST::FluidElemBase::
1544 calculate_pressure_derivative_wrt_conservative_variables(const MAST::PrimitiveSolution& sol,
1545  RealVectorX& dpdX) {
1546 
1547  const unsigned int n1 = 2 + dim;
1548 
1549  dpdX.setZero();
1550  const Real rho = sol.rho,
1551  u1 = sol.u1,
1552  u2 = sol.u2,
1553  u3 = sol.u3,
1554  k = sol.k,
1555  R = flight_condition->gas_property.R,
1556  cv= flight_condition->gas_property.cv,
1557  e = sol.e_tot,
1558  p = sol.p;
1559 
1560  switch (dim)
1561  {
1562  case 3:
1563  dpdX(3) = -R/cv*u3;
1564 
1565  case 2:
1566  dpdX(2) = -R/cv*u2;
1567 
1568  case 1: {
1569 
1570  dpdX(0) = R/cv*k;
1571  dpdX(1) = -R/cv*u1;
1572  dpdX(n1-1) = R/cv;
1573  }
1574  }
1575 }
1576 
1577 
1578 
1579 
1580 void
1581 MAST::FluidElemBase::
1582 calculate_advection_flux_jacobian_for_moving_solid_wall_boundary
1583 (const MAST::PrimitiveSolution& sol,
1584  const Real ui_ni,
1585  const libMesh::Point& nvec,
1586  const RealVectorX& dnormal,
1587  RealMatrixX& mat) {
1588 
1589  // the solid wall boundary flux is defined as
1590  // fn = ui_ni (U_c + {0,0,0,0,p}^T) + {0, p n_hat, 0}
1591  // then,
1592  // dfn/dU = ui_ni (I + {0,0,0,0,dp/dU}^T) +
1593  // dui_ni/dU (U_c + {0,0,0,0,p}^T) +
1594  // {0, dp/dU n1, dp/dU n2, dp/dU n3, 0}
1595  //
1596  // and using ui_ni = wdot (ni + dni) - ui dni
1597  // dui_ni/dU = - dni (dui/dU)
1598  //
1599  // the first and third terms of dfn/dU are calculated first.
1600  // the last section calculates the contribution of the second term
1601  // dfn/dU
1602 
1603 
1604  const unsigned int n1 = 2 + dim;
1605 
1606  mat.setZero();
1607  const Real rho = sol.rho,
1608  u1 = sol.u1,
1609  u2 = sol.u2,
1610  u3 = sol.u3,
1611  k = sol.k,
1612  R = flight_condition->gas_property.R,
1613  cv= flight_condition->gas_property.cv,
1614  e = sol.e_tot,
1615  p = sol.p;
1616 
1617  switch (dim)
1618  {
1619  case 3:
1620  {
1621  mat(1, 3) = -R/cv*nvec(0)*u3; // dp/dU n1
1622 
1623  mat(2, 3) = -R/cv*nvec(1)*u3; // dp/dU n2
1624 
1625  mat(3, 0) = nvec(2)*R/cv*k; // dp/dU n3
1626  mat(3, 1) = -R/cv*nvec(2)*u1; // dp/dU n3
1627  mat(3, 2) = -R/cv*nvec(2)*u2; // dp/dU n3
1628  mat(3, 3) = ui_ni-R/cv*nvec(2)*u3; // ui_ni dU/dU + dp/dU n3
1629  mat(3, n1-1) = R/cv*nvec(2); // dp/dU n3
1630 
1631  mat(n1-1, 3) = -ui_ni*R/cv*u3; // ui_ni dp/dU
1632  }
1633 
1634  case 2:
1635  {
1636  mat(1, 2) = -R/cv*nvec(0)*u2; // dp/dU n1
1637 
1638  mat(2, 0) = nvec(1)*R/cv*k; // dp/dU n2
1639  mat(2, 1) = -R/cv*nvec(1)*u1; // dp/dU n2
1640  mat(2, 2) = ui_ni-R/cv*nvec(1)*u2; // ui_ni dU/dU + dp/dU n2
1641  mat(2, n1-1) = R/cv*nvec(1); // dp/dU n2
1642 
1643  mat(n1-1, 2) = -ui_ni*R/cv*u2; // ui_ni dp/dU
1644  }
1645 
1646  case 1:
1647  {
1648  mat(0, 0) = ui_ni; // ui_ni dU/dU
1649 
1650  mat(1, 0) = nvec(0)*R/cv*k; // dp/dU n1
1651  mat(1, 1) = ui_ni-R/cv*nvec(0)*u1; // ui_ni dU/dU + dp/dU n1
1652  mat(1, n1-1) = R/cv*nvec(0); // dp/dU n1
1653 
1654  mat(n1-1, 0) = ui_ni*R/cv*k; // ui_ni dp/dU
1655  mat(n1-1, 1) = -ui_ni*R/cv*u1; // ui_ni dp/dU
1656  mat(n1-1, n1-1) = ui_ni*(1.+R/cv); // ui_ni dp/dU
1657  }
1658  break;
1659  }
1660 
1661 
1662  //
1663  // this calculates the Jacobian contribution from
1664  // -dni (dui/dU) (U_c + {0,0,0,0,p}^T)
1665  //
1666  //
1667  // note that
1668  // ui * ni = wi_dot * (ni + dni) - (ui * dni)
1669  // hence,
1670  // d (ui * ni)/d U = - d ui / dU * dni
1671  //
1672  // from primitive variable Jacobian,
1673  // du1/dU = 1/rho {-u1, 1, 0, 0, 0}
1674  // du2/dU = 1/rho {-u2, 0, 1, 0, 0}
1675  // du3/dU = 1/rho {-u3, 0, 0, 1, 0}
1676  //
1677  // Hence,
1678  // dni dui/dU = 1/rho{-dn1 u1 - dn2 u2 - dn3 u3, dn1, dn2 dn3, 0}
1679  //
1680  // and
1681  // (U_c + {0,0,0,0,p}^T) dni dui/dU =
1682  // {rho, rho u1, rho u2, rho u3, rho E + p}^T . 1/rho{-dn1 u1 - dn2 u2 - dn3 u3, dn1, dn2 dn3, 0}
1683  // = {1, u1, u2, u3, E+p/rho} {-dn1 u1 - dn2 u2 - dn3 u3, dn1, dn2 dn3, 0}
1684  //
1685  if (dnormal.norm() > 0.)
1686  {
1687  RealVectorX
1688  duini_dU = RealVectorX::Zero(n1),
1689  tmp = RealVectorX::Zero(n1);
1690 
1691  // initialze the Jacobian of ui_ni wrt sol variables
1692  tmp(0) = 1.;
1693  tmp(n1-1) = e+p/rho;
1694  for (unsigned int i_dim=0; i_dim<dim; i_dim++) {
1695 
1696  switch (i_dim) {
1697  case 0: {
1698 
1699  duini_dU(0) -= u1 * dnormal(0);
1700  tmp(1) = u1;
1701  }
1702  break;
1703  case 1: {
1704 
1705  duini_dU(0) -= u2 * dnormal(1);
1706  tmp(2) = u2;
1707  }
1708  break;
1709  case 2: {
1710 
1711  duini_dU(0) -= u3 * dnormal(2);
1712  tmp(3) = u3;
1713  }
1714  break;
1715  }
1716  duini_dU(i_dim+1) = dnormal(i_dim);
1717  }
1718 
1719  // now add to the Jacobian matrix
1720  for (unsigned int i=0; i<n1; i++)
1721  for (unsigned int j=0; j<n1; j++)
1722  // look at the note above to see reason for negative sign
1723  mat(i,j) -= tmp(i)*duini_dU(j);
1724  }
1725 }
1726 
1727 
1728 
1729 
1730 void
1731 MAST::FluidElemBase::
1732 calculate_entropy_variable_jacobian(const MAST::PrimitiveSolution& sol,
1733  RealMatrixX& dUdV,
1734  RealMatrixX& dVdU) {
1735 
1736  // calculates dU/dV where V is the Entropy variable vector
1737 
1738  // calculate A0 = d U / d Y , where U = conservation variable vector,
1739  // and Y = unknown variable vector
1740  // note that for conservation variables as the unknown, this is an
1741  // identity matrix
1742 
1743  const unsigned int n1 = 2 + dim;
1744  const Real rho = sol.rho,
1745  u1 = sol.u1,
1746  u2 = sol.u2,
1747  u3 = sol.u3,
1748  k = sol.k,
1749  T = sol.T,
1750  e_tot = sol.e_tot,
1751  gamma = flight_condition->gas_property.gamma,
1752  cv = flight_condition->gas_property.cv;
1753 
1754  dUdV.setZero(); dVdU.setZero();
1755 
1756  // du/dv
1757  switch (dim)
1758  {
1759  case 3:
1760  {
1761  dUdV(0, 3) = u3;
1762 
1763  dUdV(1, 3) = u1*u3;
1764 
1765  dUdV(2, 3) = u2*u3;
1766 
1767  dUdV(3, 0) = dUdV(0, 3);
1768  dUdV(3, 1) = dUdV(1, 3);
1769  dUdV(3, 2) = dUdV(2, 3);
1770  dUdV(3, 3) = u3*u3+cv*T*(gamma-1.0);
1771  dUdV(3, n1-1) = u3*(cv*T*gamma+k);
1772 
1773  dUdV(n1-1, 3) = dUdV(3, n1-1);
1774  }
1775 
1776  case 2:
1777  {
1778  dUdV(0, 2) = u2;
1779 
1780  dUdV(1, 2) = u1*u2;
1781 
1782  dUdV(2, 0) = dUdV(0, 2);
1783  dUdV(2, 1) = dUdV(1, 2);
1784  dUdV(2, 2) = u2*u2+cv*T*(gamma-1.0);
1785  dUdV(2, n1-1) = u2*(cv*T*gamma+k);
1786 
1787  dUdV(n1-1, 2) = dUdV(2, n1-1);
1788  }
1789 
1790  case 1:
1791  {
1792  dUdV(0, 0) = 1.0;
1793  dUdV(0, 1) = u1;
1794  dUdV(0, n1-1) = e_tot;
1795 
1796  dUdV(1, 0) = dUdV(0, 1);
1797  dUdV(1, 1) = u1*u1+cv*T*(gamma-1.0);
1798  dUdV(1, n1-1) = u1*(cv*T*gamma+k);
1799 
1800  dUdV(n1-1, 0) = dUdV(0, n1-1);
1801  dUdV(n1-1, 1) = dUdV(1, n1-1);
1802  dUdV(n1-1, n1-1) = k*k+gamma*cv*T*(cv*T+2*k);
1803 
1804  }
1805  break;
1806 
1807  default:
1808  break;
1809  }
1810 
1811  dUdV *= rho/(gamma-1.0);
1812 
1813 
1814  // dv/du
1815  switch (dim)
1816  {
1817  case 3:
1818  {
1819  dVdU(0, 3) = -u3*k;
1820 
1821  dVdU(1, 3) = u1*u3;
1822 
1823  dVdU(2, 3) = u2*u3;
1824 
1825  dVdU(3, 0) = dVdU(0, 3);
1826  dVdU(3, 1) = dVdU(1, 3);
1827  dVdU(3, 2) = dVdU(2, 3);
1828  dVdU(3, 3) = u3*u3+cv*T;
1829  dVdU(3, n1-1) = -u3;
1830 
1831  dVdU(n1-1, 3) = dVdU(3, n1-1);
1832  }
1833 
1834  case 2:
1835  {
1836  dVdU(0, 2) = -u2*k;
1837 
1838  dVdU(1, 2) = u1*u2;
1839 
1840  dVdU(2, 0) = dVdU(0, 2);
1841  dVdU(2, 1) = dVdU(1, 2);
1842  dVdU(2, 2) = u2*u2+cv*T;
1843  dVdU(2, n1-1) = -u2;
1844 
1845  dVdU(n1-1, 2) = dVdU(2, n1-1);
1846  }
1847 
1848  case 1:
1849  {
1850  dVdU(0, 0) = k*k+cv*cv*T*T*gamma;
1851  dVdU(0, 1) = -u1*k;
1852  dVdU(0, n1-1) = -e_tot+2.0*k;
1853 
1854  dVdU(1, 0) = dVdU(0, 1);
1855  dVdU(1, 1) = u1*u1+cv*T;
1856  dVdU(1, n1-1) = -u1;
1857 
1858  dVdU(n1-1, 0) = dVdU(0, n1-1);
1859  dVdU(n1-1, 1) = dVdU(1, n1-1);
1860  dVdU(n1-1, n1-1) = 1.0;
1861  }
1862  break;
1863 
1864  default:
1865  break;
1866  }
1867 
1868  dVdU /= (rho*cv*cv*T*T);
1869 }
1870 
1871 
1872 
1873 
1874 bool
1875 MAST::FluidElemBase::
1876 calculate_barth_tau_matrix (const unsigned int qp,
1877  const MAST::FEBase& fe,
1878  const MAST::PrimitiveSolution& sol,
1879  RealMatrixX& tau,
1880  std::vector<RealMatrixX >& tau_sens) {
1881 
1882  const unsigned int n1 = 2 + dim;
1883 
1884  libMesh::Point nvec;
1885  RealVectorX
1886  eig_val = RealVectorX::Zero(n1);
1887 
1888  RealMatrixX
1889  l_eig_vec = RealMatrixX::Zero(n1, n1),
1890  l_eig_vec_inv_tr = RealMatrixX::Zero(n1, n1),
1891  tmp1 = RealMatrixX::Zero(n1, n1);
1892 
1893  Real nval = 0.;
1894 
1895  const std::vector<std::vector<libMesh::RealVectorValue> >&
1896  dphi = fe.get_dphi(); // assuming that all variables have the same interpolation
1897 
1898  for (unsigned int i_node=0; i_node<dphi.size(); i_node++)
1899  {
1900  nvec = dphi[i_node][qp];
1901  nval = nvec.norm();
1902  if (nval > 0.) {
1903 
1904  nvec /= nval;
1905  this->calculate_advection_left_eigenvector_and_inverse_for_normal
1906  (sol, nvec, eig_val, l_eig_vec, l_eig_vec_inv_tr);
1907 
1908  for (unsigned int i_var=0; i_var<n1; i_var++)
1909  l_eig_vec_inv_tr.col(i_var) *= fabs(eig_val(i_var)); // L^-T [omaga]
1910 
1911  l_eig_vec_inv_tr *= l_eig_vec.transpose(); // A = L^-T [omaga] L^T
1912 
1913  tmp1 += nval * l_eig_vec_inv_tr; // sum_inode | A_i |
1914  }
1915  }
1916 
1917  // add the viscous contribution
1918  if (if_viscous()) {
1919 
1920  for (unsigned int i_node=0; i_node<dphi.size(); i_node++) {
1921 
1922  nvec = dphi[i_node][qp];
1923 
1924  for (unsigned int i=0; i<dim; i++)
1925  for (unsigned int j=0; j<dim; j++) {
1926 
1927  calculate_diffusion_flux_jacobian(i, j, sol, l_eig_vec);
1928  tmp1 += l_eig_vec * nvec(i) * nvec(j);
1929  }
1930  }
1931  }
1932 
1933 
1934  // now invert the tmp matrix to get the tau matrix
1935  RealVectorX
1936  x = RealVectorX::Zero(n1),
1937  b = RealVectorX::Zero(n1);
1938 
1939  for (unsigned int i_var=0; i_var<n1; i_var++) {
1940 
1941  tau_sens[i_var].setZero(); // zero the sensitivity matrix for now
1942  b.setZero(); b(i_var) = 1.;
1943  x = tmp1.lu().solve(b);
1944  tau.col(i_var) = x;
1945  }
1946 
1947  return false;
1948 }
1949 
1950 
1951 
1952 bool
1953 MAST::FluidElemBase::
1954 calculate_aliabadi_tau_matrix (const unsigned int qp,
1955  const MAST::FEBase& fe,
1956  const MAST::PrimitiveSolution& sol,
1957  RealMatrixX& tau,
1958  std::vector<RealMatrixX >& tau_sens) {
1959 
1960  const unsigned int n1 = 2 + dim;
1961 
1962  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi =
1963  fe.get_dphi(); // assuming that all variables have the same interpolation
1964 
1965  RealVectorX
1966  u = RealVectorX::Zero(dim),
1967  dN = RealVectorX::Zero(dim);
1968 
1969  const Real u1 = sol.u1,
1970  u2 = sol.u2,
1971  u3 = sol.u3,
1972  a = sol.a,
1973  dt = 0., //c.get_deltat_value(),
1974  cp = flight_condition->gas_property.cp,
1975  R = flight_condition->gas_property.R,
1976  gamma = flight_condition->gas_property.gamma;
1977 
1978  tau.setZero();
1979 
1980  // calculate the gradients
1981  switch (dim)
1982  {
1983  case 3:
1984  u(2) = u3;
1985 
1986  case 2:
1987  u(1) = u2;
1988 
1989  case 1:
1990  u(0) = u1;
1991  break;
1992 
1993  default:
1994  break;
1995  }
1996 
1997  // calculate the dot product of velocity times gradient of shape function
1998  Real h = 0, u_val = u.norm(), tau_rho, tau_m, tau_e;
1999  u /= u_val;
2000 
2001  for (unsigned int i_nodes=0; i_nodes<dphi.size(); i_nodes++)
2002  {
2003  // set value of shape function gradient
2004  for (unsigned int i_dim=0; i_dim<dim; i_dim++)
2005  dN(i_dim) = dphi[i_nodes][qp](i_dim);
2006 
2007  h += fabs(dN.dot(u));
2008  }
2009 
2010  h = 2.0/h;
2011 
2012  // now set the value of streamwise dissipation
2013  tau_rho = 1.0/sqrt(pow(2.0/dt, 2)+ pow(2.0/h*(u_val+a), 2));
2014  if (!_if_viscous)
2015  {
2016  tau_m = 1.0/sqrt(pow(2.0/dt, 2)+ pow(2.0/h*(u_val+a), 2));
2017  tau_e = 1.0/sqrt(pow(2.0/dt, 2)+ pow(2.0/h*(u_val+a), 2));
2018  }
2019  else
2020  {
2021  tau_m = 1.0/sqrt(pow(2.0/dt, 2)+ pow(2.0/h*(u_val+a), 2) +
2022  pow(4.0*sol.mu/sol.rho/h/h, 2));
2023  tau_e = 1.0/sqrt(pow(2.0/dt, 2)+ pow(2.0/h*(u_val+a), 2)+
2024  pow(4.0*sol.k_thermal/sol.rho/cp/h/h, 2));
2025  }
2026 
2027  tau(0, 0) = tau_rho;
2028  for (unsigned int i_dim=0; i_dim<dim; i_dim++)
2029  tau(1+i_dim, 1+i_dim) = tau_m;
2030  tau(n1-1, n1-1) = tau_e;
2031 
2032  // calculation sensitivity of the tau matrix for each conservative variable
2033  std::vector<Real> primitive_sens(n1), cons_sens(n1);
2034 
2035  // sensitivity wrt primitive variables
2036  for (unsigned int i_pvar=0; i_pvar<n1; i_pvar++)
2037  {
2038  switch (_active_primitive_vars[i_pvar])
2039  {
2040  case RHO_PRIM:
2041  primitive_sens[i_pvar] = 0.;
2042  break;
2043 
2044  case VEL1:
2045  primitive_sens[i_pvar] = 2.0*u1/sqrt(2.0*sol.k)/h;
2046  break;
2047 
2048  case VEL2:
2049  primitive_sens[i_pvar] = 2.0*u2/sqrt(2.0*sol.k)/h;
2050  break;
2051 
2052  case VEL3:
2053  primitive_sens[i_pvar] = 2.0*u3/sqrt(2.0*sol.k)/h;
2054  break;
2055 
2056  case TEMP:
2057  primitive_sens[i_pvar] = sqrt(gamma*R/sol.T)/h;
2058  break;
2059 
2060  default:
2061  break;
2062  }
2063 
2064  RealMatrixX
2065  dprim_dcons = RealMatrixX::Zero(n1, n1),
2066  dcons_dprim = RealMatrixX::Zero(n1, n1);
2067 
2068  std::fill(cons_sens.begin(), cons_sens.end(), 0.);
2069 
2070  // apply chain rule
2071  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
2072  {
2073  for (unsigned int i_pvar=0; i_pvar<n1; i_pvar++)
2074  cons_sens[i_cvar] += primitive_sens[i_pvar] *
2075  dprim_dcons(i_pvar, i_cvar);
2076  }
2077 
2078  // set the values in the sensitivity matrices
2079  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
2080  {
2081  tau_sens[i_cvar].setZero();
2082  cons_sens[i_cvar] *= -pow(2.0/h*(u_val+a),-2.0);
2083  for (unsigned int i=0; i<n1; i++)
2084  tau_sens[i_cvar](i,i) = cons_sens[i_cvar];
2085  }
2086  }
2087 
2088  return true;
2089 }
2090 
2091 
2092 
2093 void
2094 MAST::FluidElemBase::
2095 calculate_dxidX (const unsigned int qp, const MAST::FEBase& fe,
2096  RealMatrixX& dxi_dX,
2097  RealMatrixX& dX_dxi) {
2098 
2099 
2100  // initialize dxi_dX and dX_dxi
2101  dxi_dX.setZero(); dX_dxi.setZero();
2102  Real val=0., val2=0.;
2103 
2104  for (unsigned int i_dim=0; i_dim<dim; i_dim++)
2105  for (unsigned int j_dim=0; j_dim<dim; j_dim++)
2106  {
2107  switch (i_dim)
2108  {
2109  case 0:
2110  {
2111  switch (j_dim)
2112  {
2113  case 0:
2114  val = fe.get_dxidx()[qp];
2115  val2 = fe.get_dxyzdxi()[qp](i_dim);
2116  break;
2117  case 1:
2118  val = fe.get_dxidy()[qp];
2119  val2 = fe.get_dxyzdeta()[qp](i_dim);
2120  break;
2121  case 2:
2122  val = fe.get_dxidz()[qp];
2123  val2 = fe.get_dxyzdzeta()[qp](i_dim);
2124  break;
2125  }
2126  }
2127  break;
2128  case 1:
2129  {
2130  switch (j_dim)
2131  {
2132  case 0:
2133  val = fe.get_detadx()[qp];
2134  val2 = fe.get_dxyzdxi()[qp](i_dim);
2135  break;
2136  case 1:
2137  val = fe.get_detady()[qp];
2138  val2 = fe.get_dxyzdeta()[qp](i_dim);
2139  break;
2140  case 2:
2141  val = fe.get_detadz()[qp];
2142  val2 = fe.get_dxyzdzeta()[qp](i_dim);
2143  break;
2144  }
2145  }
2146  break;
2147  case 2:
2148  {
2149  switch (j_dim)
2150  {
2151  case 0:
2152  val = fe.get_dzetadx()[qp];
2153  val2 = fe.get_dxyzdxi()[qp](i_dim);
2154  break;
2155  case 1:
2156  val = fe.get_dzetady()[qp];
2157  val2 = fe.get_dxyzdeta()[qp](i_dim);
2158  break;
2159  case 2:
2160  val = fe.get_dzetadz()[qp];
2161  val2 = fe.get_dxyzdzeta()[qp](i_dim);
2162  break;
2163  }
2164  }
2165  break;
2166  }
2167  dxi_dX(i_dim, j_dim) = val;
2168  dX_dxi(i_dim, j_dim) = val2;
2169  }
2170 }
2171 
2172 
2173 
2174 
2175 void MAST::FluidElemBase::
2176 calculate_hartmann_discontinuity_operator (const unsigned int qp,
2177  const MAST::FEBase& fe,
2178  const MAST::PrimitiveSolution& sol,
2179  const RealVectorX& elem_solution,
2180  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2181  const RealMatrixX& Ai_Bi_advection,
2182  RealVectorX& discontinuity_val) {
2183 
2184  discontinuity_val.setZero();
2185  const unsigned int n1 = 2 + dim;
2186 
2187  RealMatrixX
2188  tmpmat1_n1n1 = RealMatrixX::Zero(dim+2, dim+2),
2189  dxi_dX = RealMatrixX::Zero(dim, dim),
2190  dX_dxi = RealMatrixX::Zero(dim, dim),
2191  dpdc = RealMatrixX::Zero(n1, n1),
2192  dcdp = RealMatrixX::Zero(n1, n1);
2193  RealVectorX
2194  vec1 = RealVectorX::Zero(n1),
2195  vec2 = RealVectorX::Zero(n1),
2196  dpress_dp = RealVectorX::Zero(n1),
2197  dp = RealVectorX::Zero(dim),
2198  hk = RealVectorX::Zero(dim);
2199 
2200 
2201  // residual of the strong form of the equation: assuming steady flow
2202  vec2 = Ai_Bi_advection * elem_solution;// Ai dU/dxi
2203 
2204  // calculate dp/dU
2205  calculate_conservative_variable_jacobian(sol, dcdp, dpdc);
2206  dpress_dp(0) = (sol.cp - sol.cv)*sol.T; // R T
2207  dpress_dp(n1-1) = (sol.cp - sol.cv)*sol.rho; // R rho
2208  vec1 = dpdc.transpose() * dpress_dp; // dpress/dprimitive * dprimitive/dconservative
2209 
2210  Real rp = fabs(vec2.dot(vec1)) / sol.p; // sum_i=1..m (dp/dU_i)*R_i / p
2211 
2212  // calculate the approximate element size
2213  this->calculate_dxidX (qp, fe, dxi_dX, dX_dxi);
2214  for (unsigned int i=0; i<dim; i++) {
2215 
2216  for (unsigned int j=0; j<dim; j++)
2217  hk(i) = fmax(hk(i), fabs(dX_dxi(i, j))); // get the mqximum value out of each xi
2218 
2219  // calculate gradient of p = dp/dprim * dprim/dcons * dcons/dx_i
2220  dB_mat[i].vector_mult(vec1, elem_solution);
2221  vec2 = dpdc * vec1;
2222  dp(i) = vec2.dot(dpress_dp);
2223  }
2224 
2225  const unsigned int fe_order = fe.get_fe_type().order;
2226 
2227  // now set the value of the dissipation coefficient
2228  for (unsigned int i=0; i<dim; i++)
2229  discontinuity_val(i) =
2230  (dp.norm() * hk(i) / sol.p / (fe_order + 1)) * // pressure switch
2231  pow(hk(i), 2) * rp;
2232  discontinuity_val *= _dissipation_scaling;
2233 }
2234 
2235 
2236 
2237 
2238 void MAST::FluidElemBase::
2239 calculate_aliabadi_discontinuity_operator(const unsigned int qp,
2240  const MAST::FEBase& fe,
2241  const MAST::PrimitiveSolution& sol,
2242  const RealVectorX& elem_solution,
2243  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2244  const RealMatrixX& Ai_Bi_advection,
2245  RealVectorX& discontinuity_val) {
2246 
2247 
2248  discontinuity_val.setZero();
2249  const unsigned int n1 = 2 + dim;
2250 
2251  std::vector<RealVectorX >
2252  diff_vec(3);
2253  RealMatrixX
2254  A_inv_entropy = RealMatrixX::Zero(dim+2, dim+2),
2255  A_entropy = RealMatrixX::Zero(dim+2, dim+2),
2256  dxi_dX = RealMatrixX::Zero(dim, dim),
2257  dX_dxi = RealMatrixX::Zero(dim, dim),
2258  tmpmat1_n1n1 = RealMatrixX::Zero(n1, n1);
2259  RealVectorX
2260  vec1 = RealVectorX::Zero(n1),
2261  vec2 = RealVectorX::Zero(n1);
2262 
2263  for (unsigned int i=0; i<dim; i++) diff_vec[i].setZero(n1);
2264 
2265  Real dval;
2266 
2267  this->calculate_dxidX (qp, fe, dxi_dX, dX_dxi);
2268  this->calculate_entropy_variable_jacobian ( sol, A_entropy, A_inv_entropy );
2269 
2270 
2271  for (unsigned int i=0; i<dim; i++)
2272  dB_mat[i].vector_mult(diff_vec[i], elem_solution); // dU/dxi
2273  vec1 = Ai_Bi_advection * elem_solution; // Ai dU/dxi
2274 
2275  // TODO: divergence of diffusive flux
2276 
2277  // add the velocity and calculate the numerator of the discontinuity
2278  // capturing term coefficient
2279  //vec2 += c.elem_solution; // add velocity TODO: how to get the
2280  vec2 = A_inv_entropy * vec1;
2281  dval = vec1.dot(vec2); // this is the numerator term
2282 
2283  // now evaluate the dissipation factor for the discontinuity capturing term
2284  // this is the denominator term
2285 
2286  // add a small number to avoid division of zero by zero
2287  Real val1 = 1.0e-6;
2288  for (unsigned int i=0; i<dim; i++) {
2289  vec1.setZero();
2290 
2291  for (unsigned int j=0; j<dim; j++)
2292  vec1 += dxi_dX(i, j) * diff_vec[j];
2293 
2294  // calculate the value of denominator
2295  vec2 = A_inv_entropy * vec1;
2296  val1 += vec1.dot(vec2);
2297  }
2298 
2299  /* // now calculate the Ducros shock sensor
2300  // Real d_ducros = 0., div = 0.;
2301  // RealVectorX u, dudx, dudy, dudz, dN, curl;
2302  // u.resize(dim); dudx.resize(3); dudy.resize(3); dudz.resize(3); curl.resize(3);
2303  // sol.get_uvec(u);
2304  //
2305  // for (unsigned int i=0; i<dim; i++) {
2306  // dudx(i) = diff_vec[0](i+1); // du/dx
2307  // dudx(i) -= diff_vec[0](0)*u(i); // -drho/dx * u
2308  // dudx(i) /= sol.rho;
2309  // }
2310  // div += dudx(0);
2311  // if (dim > 1) {
2312  // for (unsigned int i=0; i<dim; i++) {
2313  // dudy(i) = diff_vec[1](i+2); // du/dy
2314  // dudy(i) -= diff_vec[1](0)*u(i); // -drho/dy * u
2315  // dudy(i) /= sol.rho;
2316  // }
2317  // div += dudy(1);
2318  // }
2319  // if (dim > 2) {
2320  // for (unsigned int i=0; i<dim; i++) {
2321  // dudz(i) = diff_vec[2](i+3); // du/dz
2322  // dudz(i) -= diff_vec[2](0)*u(i); // -drho/dz * u
2323  // dudz(i) /= sol.rho;
2324  // }
2325  // div += dudz(2);
2326  // }
2327  // curl(0) = dudy(2)-dudz(1);
2328  // curl(1) = -(dudx(2)-dudz(0));
2329  // curl(2) = dudx(1)-dudy(0);
2330  // d_ducros = (div*div) / (div*div + pow(curl.l2_norm(),2) + 1.0e-6); */
2331 
2332  dval = sqrt(dval/val1);
2333 
2334  if (_include_pressure_switch) {
2335  // also add a pressure switch q
2336  RealMatrixX
2337  dpdc = RealMatrixX::Zero(n1, n1),
2338  dcdp = RealMatrixX::Zero(n1, n1);
2339 
2340  RealVectorX
2341  dpress_dp = RealVectorX::Zero(n1, n1),
2342  dp = RealVectorX::Zero(dim);
2343 
2344  Real p_sensor = 0., hk = 0.;
2345  calculate_conservative_variable_jacobian(sol, dcdp, dpdc);
2346  dpress_dp(0) = (sol.cp - sol.cv)*sol.T; // R T
2347  dpress_dp(n1-1) = (sol.cp - sol.cv)*sol.rho; // R rho
2348  for (unsigned int i=0; i<dim; i++) {
2349  dB_mat[i].vector_mult(vec1, elem_solution);
2350  vec2 = dpdc * vec1;
2351  dp(i) = vec2.dot(dpress_dp);
2352  for (unsigned int j=0; j<dim; j++)
2353  hk = fmax(hk, fabs(dX_dxi(i, j)));
2354  }
2355  p_sensor = dp.norm() * hk / sol.p;
2356  dval *= (p_sensor * _dissipation_scaling);
2357  }
2358 
2359  dval *= _dissipation_scaling;
2360 
2361  // set value in all three dimensions to be the same
2362  const unsigned int fe_order = fe.get_fe_type().order;
2363  for (unsigned int i=0; i<dim; i++)
2364  discontinuity_val(i) = dval/fe_order;
2365 }
2366 
2367 
2368 
2369 
2370 template <typename ValType>
2371 void
2372 MAST::FluidElemBase::
2373 calculate_small_disturbance_aliabadi_discontinuity_operator
2374 (const unsigned int qp,
2375  const MAST::FEBase& fe,
2376  const MAST::PrimitiveSolution& sol,
2377  const SmallPerturbationPrimitiveSolution<ValType>& dsol,
2378  const RealVectorX& elem_solution,
2379  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2380  const RealMatrixX& Ai_Bi_advection,
2381  RealVectorX& discontinuity_val) {
2382 
2383 
2384  discontinuity_val.setZero();
2385  const unsigned int n1 = 2 + dim;
2386 
2387  std::vector<RealVectorX >
2388  diff_vec(3);
2389 
2390  RealMatrixX
2391  A_inv_entropy = RealMatrixX::Zero(dim+2, dim+2),
2392  A_entropy = RealMatrixX::Zero(dim+2, dim+2),
2393  dxi_dX = RealMatrixX::Zero(dim, dim),
2394  dX_dxi = RealMatrixX::Zero(dim, dim),
2395  tmpmat1_n1n1 = RealMatrixX::Zero(dim+2, dim+2);
2396  RealVectorX
2397  vec1 = RealVectorX::Zero(n1),
2398  vec2 = RealVectorX::Zero(n1);
2399 
2400  for (unsigned int i=0; i<dim; i++) diff_vec[i].setZero(n1);
2401 
2402  Real dval;
2403 
2404  this->calculate_dxidX (qp, fe, dxi_dX, dX_dxi);
2405  this->calculate_entropy_variable_jacobian ( sol, A_entropy, A_inv_entropy );
2406 
2407 
2408  for (unsigned int i=0; i<dim; i++)
2409  dB_mat[i].vector_mult(diff_vec[i], elem_solution); // dU/dxi
2410  vec1 = Ai_Bi_advection * elem_solution; // Ai dU/dxi
2411 
2412  // TODO: divergence of diffusive flux
2413 
2414  // add the velocity and calculate the numerator of the discontinuity
2415  // capturing term coefficient
2416  //vec2 += c.elem_solution; // add velocity TODO: how to get the
2417  vec2 = A_inv_entropy * vec1;
2418  dval = vec1.dot(vec2); // this is the numerator term
2419 
2420  // now evaluate the dissipation factor for the discontinuity capturing term
2421  // this is the denominator term
2422 
2423  // add a small number to avoid division of zero by zero
2424  Real val1 = 1.0e-6;
2425  for (unsigned int i=0; i<dim; i++) {
2426  vec1.setZero();
2427 
2428  for (unsigned int j=0; j<dim; j++)
2429  vec1 += dxi_dX(i, j) * diff_vec[j];
2430 
2431  // calculate the value of denominator
2432  vec2 = A_inv_entropy * vec1;
2433  val1 += vec1.dot(vec2);
2434  }
2435 
2436  dval = std::min(1., sqrt(dval/val1));
2437 
2438  if (_include_pressure_switch) {
2439  // also add a pressure switch q
2440  RealMatrixX
2441  dpdc = RealMatrixX::Zero(n1, n1),
2442  dcdp = RealMatrixX::Zero(n1, n1);
2443  RealVectorX
2444  dpress_dp = RealVectorX::Zero(n1),
2445  dp = RealVectorX::Zero(dim);
2446 
2447  Real p_sensor = 0., hk = 0.;
2448  calculate_conservative_variable_jacobian(sol, dcdp, dpdc);
2449  dpress_dp(0) = (sol.cp - sol.cv)*sol.T; // R T
2450  dpress_dp(n1-1) = (sol.cp - sol.cv)*sol.rho; // R rho
2451  for (unsigned int i=0; i<dim; i++) {
2452  dB_mat[i].vector_mult(vec1, elem_solution);
2453  vec2 = dpdc * vec1;
2454  dp(i) = vec2.dot(dpress_dp);
2455  for (unsigned int j=0; j<dim; j++)
2456  hk = fmax(hk, fabs(dX_dxi(i, j)));
2457  }
2458  p_sensor = std::min(1.,
2459  dp.norm() * hk /
2460  (sol.p + abs(dsol.dp)));
2461  dval *= (p_sensor * _dissipation_scaling);
2462  }
2463 
2464  dval *= _dissipation_scaling;
2465 
2466  // set value in all three dimensions to be the same
2467  const unsigned int fe_order = fe.get_fe_type().order;
2468  for (unsigned int i=0; i<dim; i++)
2469  discontinuity_val(i) = dval/fe_order;
2470 }
2471 
2472 
2473 
2474 void MAST::FluidElemBase::
2475 calculate_differential_operator_matrix (const unsigned int qp,
2476  const MAST::FEBase& fe,
2477  const RealVectorX& elem_solution,
2478  const MAST::PrimitiveSolution& sol,
2479  const MAST::FEMOperatorMatrix& B_mat,
2480  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2481  const std::vector<RealMatrixX >& Ai_advection,
2482  const RealMatrixX& Ai_Bi_advection,
2483  const std::vector<std::vector<RealMatrixX > >& Ai_sens,
2484  RealMatrixX& LS_operator,
2485  RealMatrixX& LS_sens) {
2486 
2487  const unsigned int n1 = 2 + dim, n2 = B_mat.n();
2488 
2489  RealMatrixX
2490  mat = RealMatrixX::Zero(n1, n2),
2491  mat2 = RealMatrixX::Zero(n1, n2),
2492  tau = RealMatrixX::Zero(n1, n1);
2493  RealVectorX
2494  vec1 = RealVectorX::Zero(n1),
2495  vec2 = RealVectorX::Zero(n1),
2496  vec3 = RealVectorX::Zero(n1),
2497  vec4_n2 = RealVectorX::Zero(n2);
2498 
2499 
2500  const std::vector<std::vector<Real> >& phi =
2501  fe.get_phi(); // assuming that all variables have the same interpolation
2502  const unsigned int n_phi = (unsigned int) phi.size();
2503  std::vector<RealMatrixX > tau_sens(n1);
2504  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
2505  tau_sens[i_cvar].setZero(n1, n1);
2506 
2507  // contribution of unsteady term
2508  LS_operator.setZero();
2509  LS_sens.setZero();
2510 
2511  bool if_diagonal_tau = false;
2512 
2513  vec2.setZero();
2514  vec2 = Ai_Bi_advection * elem_solution; // sum A_i dU/dx_i
2515 
2516  //if_diagonal_tau = this->calculate_aliabadi_tau_matrix
2517  //(qp, c, sol, tau, tau_sens);
2518  if_diagonal_tau = this->calculate_barth_tau_matrix
2519  (qp, fe, sol, tau, tau_sens);
2520 
2521  // contribution of advection flux term
2522  for (unsigned int i=0; i<dim; i++)
2523  {
2524  mat = Ai_advection[i].transpose();
2525  dB_mat[i].left_multiply(mat2, mat);
2526  LS_operator += mat2; // A_i^T dB/dx_i
2527 
2528 
2529  // sensitivity of the LS operator times strong form of residual
2530  // Bi^T dAi/dalpha tau
2531  vec1 = tau * vec2;
2532  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
2533  {
2534  vec3 = Ai_sens[i][i_cvar] * vec1;
2535  dB_mat[i].vector_mult_transpose(vec4_n2, vec3);
2536  for (unsigned int i_phi=0; i_phi<n_phi; i_phi++)
2537  LS_sens.col((n_phi*i_cvar)+i_phi) += phi[i_phi][qp] * vec4_n2;
2538  }
2539 
2540  // Bi^T Ai dtau/dalpha
2541  for (unsigned int i_cvar=0; i_cvar<n1; i_cvar++)
2542  {
2543  vec1 = tau_sens[i_cvar] * vec2;
2544  vec3 = Ai_advection[i] * vec1;
2545  dB_mat[i].vector_mult_transpose(vec4_n2, vec3);
2546  for (unsigned int i_phi=0; i_phi<n_phi; i_phi++)
2547  LS_sens.col((n_phi*i_cvar)+i_phi) += phi[i_phi][qp] * vec4_n2;
2548  }
2549  }
2550 
2551  // scale the LS matrix with the correct factor
2552  if (if_diagonal_tau) {
2553 
2554  for (unsigned int i=0; i<n1; i++)
2555  LS_operator.row(i) *= tau(i,i);
2556  }
2557  else
2558  LS_operator = tau.transpose() * LS_operator;
2559 }
2560 
2561 
2562 
2563 // template instantiations
2564 template void
2565 MAST::FluidElemBase::
2566 calculate_small_disturbance_aliabadi_discontinuity_operator<Real>
2567 (const unsigned int qp,
2568  const MAST::FEBase& fe,
2569  const MAST::PrimitiveSolution& sol,
2570  const SmallPerturbationPrimitiveSolution<Real>& dsol,
2571  const RealVectorX& elem_solution,
2572  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2573  const RealMatrixX& Ai_Bi_advection,
2574  RealVectorX& discontinuity_val);
2575 
2576 
2577 
2578 template void
2579 MAST::FluidElemBase::
2580 calculate_small_disturbance_aliabadi_discontinuity_operator<Complex>
2581 (const unsigned int qp,
2582  const MAST::FEBase& fe,
2583  const MAST::PrimitiveSolution& sol,
2584  const SmallPerturbationPrimitiveSolution<Complex>& dsol,
2585  const RealVectorX& elem_solution,
2586  const std::vector<MAST::FEMOperatorMatrix>& dB_mat,
2587  const RealMatrixX& Ai_Bi_advection,
2588  RealVectorX& discontinuity_val);
2589 
2590 
2591 
MAST::FEBase::get_detadz
virtual const std::vector< Real > & get_detadz() const
Definition: fe_base.cpp:311
MAST::FEBase::get_dzetadx
virtual const std::vector< Real > & get_dzetadx() const
Definition: fe_base.cpp:319
MAST::SmallPerturbationPrimitiveSolution
Class defines basic operations and calculation of the small disturbance primitive variables...
Definition: fluid_elem_base.h:42
MAST::FluidElemBase::calculate_advection_flux
void calculate_advection_flux(const unsigned int calculate_dim, const MAST::PrimitiveSolution &sol, RealVectorX &flux)
Definition: fluid_elem_base.cpp:150
small_disturbance_primitive_fluid_solution.h
MAST::FluidElemBase::if_viscous
bool if_viscous() const
Definition: fluid_elem_base.h:100
MAST::FluidElemBase::calculate_advection_flux_jacobian_for_moving_solid_wall_boundary
void calculate_advection_flux_jacobian_for_moving_solid_wall_boundary(const MAST::PrimitiveSolution &sol, const Real ui_ni, const libMesh::Point &nvec, const RealVectorX &dnvec, RealMatrixX &mat)
Definition: fluid_elem_base.cpp:1583
fe_base.h
MAST::FEBase::get_dxyzdzeta
virtual const std::vector< libMesh::RealVectorValue > & get_dxyzdzeta() const
Definition: fe_base.cpp:359
MAST::FluidElemBase::FluidElemBase
FluidElemBase(const unsigned int dimension, const MAST::FlightCondition &f)
Constructor.
Definition: fluid_elem_base.cpp:39
MAST::FEMOperatorMatrix::reinit
void reinit(unsigned int n_interpolated_vars, unsigned int n_discrete_vars, unsigned int n_discrete_dofs_per_var)
this initializes the operator for number of rows and variables, assuming that all variables has the s...
Definition: fem_operator_matrix.h:210
MAST::FEBase::get_dxidy
virtual const std::vector< Real > & get_dxidy() const
Definition: fe_base.cpp:279
MAST::FluidElemBase::calculate_diffusion_tensors_sensitivity
void calculate_diffusion_tensors_sensitivity(const RealVectorX &elem_sol, const RealVectorX &elem_sol_sens, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const RealMatrixX &dprim_dcons, const MAST::PrimitiveSolution &psol, const MAST::SmallPerturbationPrimitiveSolution< Real > &dsol, RealMatrixX &stress_tensor_sens, RealVectorX &temp_gradient_sens)
calculates and returns the stress tensor in stress_tensor.
Definition: fluid_elem_base.cpp:292
MAST::FEBase::get_dxidz
virtual const std::vector< Real > & get_dxidz() const
Definition: fe_base.cpp:287
MAST::PrimitiveSolution
Class defines the conversion and some basic operations on primitive fluid variables used in calculati...
Definition: primitive_fluid_solution.h:34
MAST::FluidElemBase::calculate_conservative_variable_jacobian
void calculate_conservative_variable_jacobian(const MAST::PrimitiveSolution &sol, RealMatrixX &dcons_dprim, RealMatrixX &dprim_dcons)
Definition: fluid_elem_base.cpp:340
MAST::FEBase::get_detadx
virtual const std::vector< Real > & get_detadx() const
Definition: fe_base.cpp:295
MAST::FluidElemBase::calculate_advection_flux_jacobian_sensitivity_for_primitive_variable
void calculate_advection_flux_jacobian_sensitivity_for_primitive_variable(const unsigned int calculate_dim, const unsigned int primitive_var, const MAST::PrimitiveSolution &sol, RealMatrixX &mat)
Definition: fluid_elem_base.cpp:1009
MAST::FluidElemBase::update_solution_at_quadrature_point
void update_solution_at_quadrature_point(const unsigned int qp, const MAST::FEBase &fe, const RealVectorX &elem_solution, RealVectorX &conservative_sol, MAST::PrimitiveSolution &primitive_sol, MAST::FEMOperatorMatrix &B_mat, std::vector< MAST::FEMOperatorMatrix > &dB_mat)
Definition: fluid_elem_base.cpp:103
MAST::FluidElemBase::calculate_aliabadi_tau_matrix
bool calculate_aliabadi_tau_matrix(const unsigned int qp, const MAST::FEBase &fe, const MAST::PrimitiveSolution &sol, RealMatrixX &tau, std::vector< RealMatrixX > &tau_sens)
Definition: fluid_elem_base.cpp:1954
MAST::FluidElemBase::calculate_advection_flux_jacobian_sensitivity_for_conservative_variable
void calculate_advection_flux_jacobian_sensitivity_for_conservative_variable(const unsigned int calculate_dim, const MAST::PrimitiveSolution &sol, std::vector< RealMatrixX > &mat)
Definition: fluid_elem_base.cpp:976
MAST::FEBase::get_dxyzdeta
virtual const std::vector< libMesh::RealVectorValue > & get_dxyzdeta() const
Definition: fe_base.cpp:351
Real
libMesh::Real Real
Definition: mast_data_types.h:32
MAST::FEBase::get_dzetady
virtual const std::vector< Real > & get_dzetady() const
Definition: fe_base.cpp:327
MAST::FluidElemBase::calculate_diffusion_flux
void calculate_diffusion_flux(const unsigned int calculate_dim, const MAST::PrimitiveSolution &sol, const RealMatrixX &stress_tensor, const RealVectorX &temp_gradient, RealVectorX &flux)
Definition: fluid_elem_base.cpp:222
MAST::FluidElemBase::calculate_diffusion_tensors
void calculate_diffusion_tensors(const RealVectorX &elem_sol, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const RealMatrixX &dprim_dcons, const MAST::PrimitiveSolution &psol, RealMatrixX &stress_tensor, RealVectorX &temp_gradient)
calculates and returns the stress tensor in stress_tensor.
Definition: fluid_elem_base.cpp:255
MAST::FEMOperatorMatrix::n
unsigned int n() const
Definition: fem_operator_matrix.h:52
MAST::FluidElemBase::calculate_diffusion_flux_jacobian
void calculate_diffusion_flux_jacobian(const unsigned int flux_dim, const unsigned int deriv_dim, const MAST::PrimitiveSolution &sol, RealMatrixX &mat)
Definition: fluid_elem_base.cpp:759
MAST::FluidElemBase::_active_conservative_vars
std::vector< FluidConservativeVars > _active_conservative_vars
Definition: fluid_elem_base.h:446
MAST::FluidElemBase::flight_condition
const MAST::FlightCondition * flight_condition
This defines the surface motion for use with the nonlinear fluid solver.
Definition: fluid_elem_base.h:96
MAST::FEBase::get_dzetadz
virtual const std::vector< Real > & get_dzetadz() const
Definition: fe_base.cpp:335
MAST::FluidElemBase::calculate_small_disturbance_aliabadi_discontinuity_operator
void calculate_small_disturbance_aliabadi_discontinuity_operator(const unsigned int qp, const MAST::FEBase &fe, const MAST::PrimitiveSolution &sol, const MAST::SmallPerturbationPrimitiveSolution< ValType > &dsol, const RealVectorX &elem_solution, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const RealMatrixX &Ai_Bi_advection, RealVectorX &discontinuity_val)
Definition: fluid_elem_base.cpp:2374
MAST::FluidElemBase::calculate_entropy_variable_jacobian
void calculate_entropy_variable_jacobian(const MAST::PrimitiveSolution &sol, RealMatrixX &dUdV, RealMatrixX &dVdU)
Definition: fluid_elem_base.cpp:1732
MAST::FluidElemBase::calculate_advection_flux_jacobian
void calculate_advection_flux_jacobian(const unsigned int calculate_dim, const MAST::PrimitiveSolution &sol, RealMatrixX &mat)
Definition: fluid_elem_base.cpp:428
RealMatrixX
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Definition: mast_data_types.h:44
primitive_fluid_solution.h
MAST::FluidElemBase::calculate_aliabadi_discontinuity_operator
void calculate_aliabadi_discontinuity_operator(const unsigned int qp, const MAST::FEBase &fe, const MAST::PrimitiveSolution &sol, const RealVectorX &elem_solution, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const RealMatrixX &Ai_Bi_advection, RealVectorX &discontinuity_val)
Definition: fluid_elem_base.cpp:2239
MAST::FEMOperatorMatrix::vector_mult
void vector_mult(T &res, const T &v) const
res = [this] * v
Definition: fem_operator_matrix.h:280
RealVectorX
Matrix< Real, Dynamic, 1 > RealVectorX
Definition: mast_data_types.h:35
MAST::FluidElemBase::calculate_differential_operator_matrix
void calculate_differential_operator_matrix(const unsigned int qp, const MAST::FEBase &fe, const RealVectorX &elem_solution, const MAST::PrimitiveSolution &sol, const MAST::FEMOperatorMatrix &B_mat, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const std::vector< RealMatrixX > &Ai_advection, const RealMatrixX &Ai_Bi_advection, const std::vector< std::vector< RealMatrixX > > &Ai_sens, RealMatrixX &LS_operator, RealMatrixX &LS_sens)
Definition: fluid_elem_base.cpp:2475
fluid_elem_base.h
MAST::FEBase::get_detady
virtual const std::vector< Real > & get_detady() const
Definition: fe_base.cpp:303
MAST::PrimitiveSolution::init
void init(const unsigned int dim, const RealVectorX &conservative_sol, const Real cp_val, const Real cv_val, bool if_viscous)
Definition: primitive_fluid_solution.cpp:62
MAST::FEBase::get_dxidx
virtual const std::vector< Real > & get_dxidx() const
Definition: fe_base.cpp:271
MAST::FluidElemBase::calculate_pressure_derivative_wrt_conservative_variables
void calculate_pressure_derivative_wrt_conservative_variables(const MAST::PrimitiveSolution &sol, RealVectorX &dpdX)
Definition: fluid_elem_base.cpp:1544
flight_condition.h
MAST::FEBase::get_dphi
virtual const std::vector< std::vector< libMesh::RealVectorValue > > & get_dphi() const
Definition: fe_base.cpp:255
MAST::FlightCondition::gas_property
GasProperty gas_property
Ambient air properties.
Definition: flight_condition.h:67
MAST::FluidElemBase::calculate_hartmann_discontinuity_operator
void calculate_hartmann_discontinuity_operator(const unsigned int qp, const MAST::FEBase &fe, const MAST::PrimitiveSolution &sol, const RealVectorX &elem_solution, const std::vector< MAST::FEMOperatorMatrix > &dB_mat, const RealMatrixX &Ai_Bi_advection, RealVectorX &discontinuity_val)
Definition: fluid_elem_base.cpp:2176
MAST::FEBase::get_dxyzdxi
virtual const std::vector< libMesh::RealVectorValue > & get_dxyzdxi() const
Definition: fe_base.cpp:343
MAST::FluidElemBase::calculate_dxidX
void calculate_dxidX(const unsigned int qp, const MAST::FEBase &fe, RealMatrixX &dxi_dX, RealMatrixX &dX_dxi)
Definition: fluid_elem_base.cpp:2095
MAST::FEBase::get_phi
virtual const std::vector< std::vector< Real > > & get_phi() const
Definition: fe_base.cpp:247
MAST::FluidElemBase::calculate_diffusion_flux_jacobian_primitive_vars
void calculate_diffusion_flux_jacobian_primitive_vars(const unsigned int flux_dim, const unsigned int deriv_dim, const RealVectorX &uvec, const bool zero_kth, const MAST::PrimitiveSolution &sol, RealMatrixX &mat)
Definition: fluid_elem_base.cpp:575
MAST::FluidElemBase::calculate_barth_tau_matrix
bool calculate_barth_tau_matrix(const unsigned int qp, const MAST::FEBase &fe, const MAST::PrimitiveSolution &sol, RealMatrixX &tau, std::vector< RealMatrixX > &tau_sens)
Definition: fluid_elem_base.cpp:1876
MAST::FEBase::get_fe_type
libMesh::FEType get_fe_type() const
Definition: fe_base.cpp:212
MAST::FluidElemBase::calculate_advection_left_eigenvector_and_inverse_for_normal
void calculate_advection_left_eigenvector_and_inverse_for_normal(const MAST::PrimitiveSolution &sol, const libMesh::Point &normal, RealVectorX &eig_vals, RealMatrixX &l_eig_mat, RealMatrixX &l_eig_mat_inv_tr)
Definition: fluid_elem_base.cpp:1285
MAST::FluidElemBase::_active_primitive_vars
std::vector< MAST::FluidPrimitiveVars > _active_primitive_vars
Definition: fluid_elem_base.h:444
MAST::FluidElemBase::get_infinity_vars
void get_infinity_vars(RealVectorX &vars_inf) const
Definition: fluid_elem_base.cpp:83