[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
reynolds_averaged_navier_stokes.cpp
1 #include <cmath>
2 #include <vector>
3 #include <complex> // for the jacobian
4 
5 #include "ADTypes.hpp"
6 
7 #include "model.h"
8 #include "reynolds_averaged_navier_stokes.h"
9 
10 namespace PHiLiP {
11 namespace Physics {
12 
13 //================================================================
14 // Reynolds Averaged Navier Stokes (RANS) Base Class
15 //================================================================
16 template <int dim, int nstate, typename real>
18  const Parameters::AllParameters *const parameters_input,
19  const double ref_length,
20  const double gamma_gas,
21  const double mach_inf,
22  const double angle_of_attack,
23  const double side_slip_angle,
24  const double prandtl_number,
25  const double reynolds_number_inf,
26  const bool use_constant_viscosity,
27  const double constant_viscosity,
28  const double turbulent_prandtl_number,
29  const double temperature_inf,
30  const double isothermal_wall_temperature,
31  const thermal_boundary_condition_enum thermal_boundary_condition_type,
32  std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function,
33  const two_point_num_flux_enum two_point_num_flux_type)
34  : ModelBase<dim,nstate,real>(manufactured_solution_function)
35  , turbulent_prandtl_number(turbulent_prandtl_number)
36  , navier_stokes_physics(std::make_unique < NavierStokes<dim,nstate_navier_stokes,real> > (
37  parameters_input,
38  ref_length,
39  gamma_gas,
40  mach_inf,
41  angle_of_attack,
42  side_slip_angle,
43  prandtl_number,
44  reynolds_number_inf,
45  use_constant_viscosity,
46  constant_viscosity,
47  temperature_inf,
48  isothermal_wall_temperature,
49  thermal_boundary_condition_type,
50  manufactured_solution_function,
51  two_point_num_flux_type))
52 {
53  static_assert(nstate>=dim+3, "ModelBase::ReynoldsAveragedNavierStokesBase() should be created with nstate>=dim+3");
54 }
55 //----------------------------------------------------------------
56 template <int dim, int nstate, typename real>
57 template<typename real2>
60  const dealii::Tensor<1,3,real2> &vector) const
61 {
62  real2 vector_magnitude_sqr; // complex initializes it as 0+0i
63  if(std::is_same<real2,real>::value){
64  vector_magnitude_sqr = 0.0;
65  }
66  for (int i=0; i<3; ++i) {
67  vector_magnitude_sqr += vector[i]*vector[i];
68  }
69  return vector_magnitude_sqr;
70 }
71 //----------------------------------------------------------------
72 template <int dim, int nstate, typename real>
73 template<typename real2>
76  const dealii::Tensor<2,dim,real2> &tensor) const
77 {
78  real2 tensor_magnitude_sqr; // complex initializes it as 0+0i
79  if(std::is_same<real2,real>::value){
80  tensor_magnitude_sqr = 0.0;
81  }
82  for (int i=0; i<dim; ++i) {
83  for (int j=0; j<dim; ++j) {
84  tensor_magnitude_sqr += tensor[i][j]*tensor[i][j];
85  }
86  }
87  return tensor_magnitude_sqr;
88 }
89 //----------------------------------------------------------------
90 template <int dim, int nstate, typename real>
91 std::array<dealii::Tensor<1,dim,real>,nstate> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
93  const std::array<real,nstate> &conservative_soln) const
94 {
95  return convective_flux_templated<real>(conservative_soln);
96 }
97 //----------------------------------------------------------------
98 template <int dim, int nstate, typename real>
99 template <typename real2>
100 std::array<dealii::Tensor<1,dim,real2>,nstate> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
102  const std::array<real2,nstate> &conservative_soln) const
103 {
104  const std::array<real2,nstate_navier_stokes> conservative_soln_rans = extract_rans_conservative_solution(conservative_soln);
105  const dealii::Tensor<1,dim,real2> vel = this->navier_stokes_physics->template compute_velocities<real2>(conservative_soln_rans); // from Euler
106  std::array<dealii::Tensor<1,dim,real2>,nstate> conv_flux;
107 
108  for (int flux_dim=0; flux_dim<nstate_navier_stokes; ++flux_dim) {
109  conv_flux[flux_dim] = 0.0; // No additional convective terms for RANS
110  }
111  // convective flux of additional RANS turbulence model
112  for (int flux_dim=nstate_navier_stokes; flux_dim<nstate; ++flux_dim) {
113  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim) {
114  conv_flux[flux_dim][velocity_dim] = conservative_soln[flux_dim]*vel[velocity_dim]; // Convective terms for turbulence model
115  }
116  }
117  return conv_flux;
118 }
119 //----------------------------------------------------------------
120 template <int dim, int nstate, typename real>
121 std::array<dealii::Tensor<1,dim,real>,nstate> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
123  const std::array<real,nstate> &conservative_soln,
124  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
125  const dealii::types::global_dof_index cell_index) const
126 {
127  return dissipative_flux_templated<real>(conservative_soln,solution_gradient,cell_index);
128 }
129 //----------------------------------------------------------------
130 template <int dim, int nstate, typename real>
131 template <typename real2>
132 std::array<dealii::Tensor<1,dim,real2>,nstate> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
134  const std::array<real2,nstate> &conservative_soln,
135  const std::array<dealii::Tensor<1,dim,real2>,nstate> &solution_gradient,
136  const dealii::types::global_dof_index /*cell_index*/) const
137 {
138 
139  const std::array<real2,nstate_navier_stokes> conservative_soln_rans = extract_rans_conservative_solution(conservative_soln);
140  const std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> solution_gradient_rans = extract_rans_solution_gradient(solution_gradient);
141 
142  // Step 1,2: Primitive solution and Gradient of primitive solution
143  const std::array<real2,nstate_navier_stokes> primitive_soln_rans = this->navier_stokes_physics->convert_conservative_to_primitive(conservative_soln_rans); // from Euler
144  const std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> primitive_soln_gradient_rans = this->navier_stokes_physics->convert_conservative_gradient_to_primitive_gradient(conservative_soln_rans, solution_gradient_rans);
145  const std::array<real2,nstate_turbulence_model> primitive_soln_turbulence_model = this->convert_conservative_to_primitive_turbulence_model(conservative_soln);
146  const std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> primitive_soln_gradient_turbulence_model = this->convert_conservative_gradient_to_primitive_gradient_turbulence_model(conservative_soln, solution_gradient);
147 
148  // Step 3: Viscous stress tensor, Velocities, Heat flux
149  const dealii::Tensor<1,dim,real2> vel = this->navier_stokes_physics->extract_velocities_from_primitive(primitive_soln_rans); // from Euler
150  // Templated virtual member functions
151  dealii::Tensor<2,dim,real2> viscous_stress_tensor;
152  dealii::Tensor<1,dim,real2> heat_flux;
153  if constexpr(std::is_same<real2,real>::value){
154  viscous_stress_tensor = compute_Reynolds_stress_tensor(primitive_soln_rans, primitive_soln_gradient_rans,primitive_soln_turbulence_model);
155  heat_flux = compute_Reynolds_heat_flux(primitive_soln_rans, primitive_soln_gradient_rans,primitive_soln_turbulence_model);
156  }
157  else if constexpr(std::is_same<real2,FadType>::value){
158  viscous_stress_tensor = compute_Reynolds_stress_tensor_fad(primitive_soln_rans, primitive_soln_gradient_rans,primitive_soln_turbulence_model);
159  heat_flux = compute_Reynolds_heat_flux_fad(primitive_soln_rans, primitive_soln_gradient_rans,primitive_soln_turbulence_model);
160  }
161  else{
162  std::cout << "ERROR in physics/reynolds_averaged_navier_stokes.cpp --> dissipative_flux_templated(): real2!=real || real2!=FadType)" << std::endl;
163  std::abort();
164  }
165 
166  // Step 4: Construct viscous flux; Note: sign corresponds to LHS
167  std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> viscous_flux_rans
168  = this->navier_stokes_physics->dissipative_flux_given_velocities_viscous_stress_tensor_and_heat_flux(vel,viscous_stress_tensor,heat_flux);
169 
170  std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> viscous_flux_turbulence_model
171  = this->dissipative_flux_turbulence_model(primitive_soln_rans,primitive_soln_turbulence_model,primitive_soln_gradient_turbulence_model);
172 
173  std::array<dealii::Tensor<1,dim,real2>,nstate> viscous_flux;
174  for(int flux_dim=0; flux_dim<nstate_navier_stokes; ++flux_dim)
175  {
176  viscous_flux[flux_dim] = viscous_flux_rans[flux_dim];
177  }
178  for(int flux_dim=nstate_navier_stokes; flux_dim<nstate; ++flux_dim)
179  {
180  viscous_flux[flux_dim] = viscous_flux_turbulence_model[flux_dim-(nstate_navier_stokes)];
181  }
182 
183  return viscous_flux;
184 }
185 //----------------------------------------------------------------
186 template <int dim, int nstate, typename real>
187 template <typename real2>
190  const std::array<real2,nstate> &conservative_soln) const
191 {
192  std::array<real2,nstate_navier_stokes> conservative_soln_rans;
193  for(int i=0; i<nstate_navier_stokes; ++i){
194  conservative_soln_rans[i] = conservative_soln[i];
195  }
196 
197  return conservative_soln_rans;
198 }
199 //----------------------------------------------------------------
200 template <int dim, int nstate, typename real>
201 template <typename real2>
202 std::array<dealii::Tensor<1,dim,real2>,dim+2> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
204  const std::array<dealii::Tensor<1,dim,real2>,nstate> &solution_gradient) const
205 {
206  std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> solution_gradient_rans;
207  for(int i=0; i<nstate_navier_stokes; ++i){
208  solution_gradient_rans[i] = solution_gradient[i];
209  }
210 
211  return solution_gradient_rans;
212 }
213 //----------------------------------------------------------------
214 template <int dim, int nstate, typename real>
215 template <typename real2>
216 std::array<dealii::Tensor<1,dim,real2>,nstate-(dim+2)> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
218  const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
219  const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model,
220  const std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> &primitive_solution_gradient_turbulence_model) const
221 {
222  std::array<real2,nstate_turbulence_model> effective_viscosity_turbulence_model;
223 
224  if constexpr(std::is_same<real2,real>::value){
225  effective_viscosity_turbulence_model = compute_effective_viscosity_turbulence_model(primitive_soln_rans, primitive_soln_turbulence_model);
226  }
227  else if constexpr(std::is_same<real2,FadType>::value){
228  effective_viscosity_turbulence_model = compute_effective_viscosity_turbulence_model_fad(primitive_soln_rans, primitive_soln_turbulence_model);
229  }
230  else{
231  std::cout << "ERROR in physics/reynolds_averaged_navier_stokes.cpp --> dissipative_flux_turbulence_model(): real2!=real || real2!=FadType)" << std::endl;
232  std::abort();
233  }
234 
235  std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> viscous_flux_turbulence_model;
236 
237  for(int i=0; i<nstate_turbulence_model; ++i){
238  for(int j=0; j<dim; ++j){
239  viscous_flux_turbulence_model[i][j] = -effective_viscosity_turbulence_model[i]*primitive_solution_gradient_turbulence_model[i][j];
240  }
241  }
242 
243  return viscous_flux_turbulence_model;
244 }
245 //----------------------------------------------------------------
246 template <int dim, int nstate, typename real>
247 template <typename real2>
248 std::array<real2,nstate-(dim+2)> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
250  const std::array<real2,nstate> &conservative_soln) const
251 {
252  std::array<real2,nstate_turbulence_model> primitive_soln_turbulence_model;
253  for(int i=0; i<nstate_turbulence_model; ++i){
254  primitive_soln_turbulence_model[i] = conservative_soln[nstate_navier_stokes+i]/conservative_soln[0];
255  }
256 
257  return primitive_soln_turbulence_model;
258 }
259 //----------------------------------------------------------------
260 template <int dim, int nstate, typename real>
261 template <typename real2>
262 std::array<dealii::Tensor<1,dim,real2>,nstate-(dim+2)> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
264  const std::array<real2,nstate> &conservative_soln,
265  const std::array<dealii::Tensor<1,dim,real2>,nstate> &solution_gradient) const
266 {
267  const std::array<real2,nstate_turbulence_model> primitive_soln_turbulence_model = this->convert_conservative_to_primitive_turbulence_model(conservative_soln);
268  std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> primitive_soln_gradient_turbulence_model;
269 
270  for(int i=0; i<nstate_turbulence_model; ++i){
271  for(int j=0; j<dim; ++j){
272  primitive_soln_gradient_turbulence_model[i][j] = (solution_gradient[nstate_navier_stokes+i][j]-primitive_soln_turbulence_model[i]*solution_gradient[0][j])/conservative_soln[0];
273  }
274  }
275 
276  return primitive_soln_gradient_turbulence_model;
277 }
278 //----------------------------------------------------------------
279 template <int dim, int nstate, typename real>
280 std::array<real,nstate-(dim+2)> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
281 ::compute_mean_turbulence_property(const std::array<real,nstate> &conservative_soln1,
282  const std::array<real,nstate> &conservative_soln2) const
283 {
284  std::array<real,nstate_turbulence_model> mean_turbulence_property;
285 
286  const std::array<real,nstate_turbulence_model> primitive_soln1_turbulence_model = convert_conservative_to_primitive_turbulence_model(conservative_soln1);
287  const std::array<real,nstate_turbulence_model> primitive_soln2_turbulence_model = convert_conservative_to_primitive_turbulence_model(conservative_soln2);
288 
289  for (int i=0;i<nstate_turbulence_model;++i){
290  mean_turbulence_property[i] = (primitive_soln1_turbulence_model[i]+primitive_soln2_turbulence_model[i])/2.0;
291  }
292 
293  return mean_turbulence_property;
294 }
295 //----------------------------------------------------------------
296 template <int dim, int nstate, typename real>
299  const std::array<real,nstate> &conservative_soln,
300  const dealii::Tensor<1,dim,real> &normal) const
301 {
302  const std::array<real,nstate_navier_stokes> conservative_soln_rans = extract_rans_conservative_solution(conservative_soln);
303  std::array<real,nstate> eig;
304  const real vel_dot_n = this->navier_stokes_physics->convective_eigenvalues(conservative_soln_rans,normal)[0];
305  for (int i=0; i<nstate_navier_stokes; ++i) {
306  eig[i] = 0.0;
307  }
308  for (int i=nstate_navier_stokes; i<nstate; ++i) {
309  eig[i] = vel_dot_n;
310  }
311  return eig;
312 }
313 //----------------------------------------------------------------
314 template <int dim, int nstate, typename real>
316 ::max_convective_eigenvalue (const std::array<real,nstate> &conservative_soln) const
317 {
318  const std::array<real,nstate_navier_stokes> conservative_soln_rans = extract_rans_conservative_solution(conservative_soln);
319 
320  const dealii::Tensor<1,dim,real> vel = this->navier_stokes_physics->template compute_velocities<real>(conservative_soln_rans);
321 
322  const real vel2 = this->navier_stokes_physics->template compute_velocity_squared<real>(vel);
323 
324  const real max_eig = sqrt(vel2);
325 
326  return max_eig;
327 }
328 //----------------------------------------------------------------
329 template <int dim, int nstate, typename real>
332  const std::array<real,nstate> &conservative_soln,
333  const dealii::Tensor<1,dim,real> &normal) const
334 {
335  const std::array<real,nstate_navier_stokes> conservative_soln_rans = extract_rans_conservative_solution(conservative_soln);
336 
337  const dealii::Tensor<1,dim,real> vel = this->navier_stokes_physics->template compute_velocities<real>(conservative_soln_rans);
338 
339  real vel_dot_n = 0.0;
340  for (int d=0;d<dim;++d) { vel_dot_n += vel[d]*normal[d]; };
341  const real max_normal_eig = sqrt(vel_dot_n*vel_dot_n);
342 
343  return max_normal_eig;
344 }
345 //----------------------------------------------------------------
346 template <int dim, int nstate, typename real>
349  const dealii::Point<dim,real> &pos,
350  const std::array<real,nstate> &conservative_soln,
351  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
352  const dealii::types::global_dof_index /*cell_index*/) const
353 {
354  std::array<real,nstate> physical_source;
355  physical_source = this->compute_production_dissipation_cross_term(pos, conservative_soln, solution_gradient);
356 
357  return physical_source;
358 }
359 //----------------------------------------------------------------
360 template <int dim, int nstate, typename real>
363  const dealii::Point<dim,real> &pos,
364  const std::array<real,nstate> &/*solution*/,
365  const real /*current_time*/,
366  const dealii::types::global_dof_index cell_index) const
367 {
368  std::array<real,nstate> conv_source_term = convective_source_term_computed_from_manufactured_solution(pos);
369  std::array<real,nstate> diss_source_term = dissipative_source_term_computed_from_manufactured_solution(pos,cell_index);
370  std::array<real,nstate> phys_source_source_term = physical_source_term_computed_from_manufactured_solution(pos,cell_index);
371  std::array<real,nstate> source_term;
372  for (int s=0; s<nstate; ++s)
373  {
374  source_term[s] = conv_source_term[s] + diss_source_term[s] - phys_source_source_term[s];
375  }
376  return source_term;
377 }
378 
379 //----------------------------------------------------------------
380 template <int dim, int nstate, typename real>
383  const dealii::Point<dim,real> &pos,
384  const std::array<real,nstate> &/*solution*/,
385  const dealii::types::global_dof_index cell_index) const
386 {
387  std::array<real,nstate> conv_source_term = convective_source_term_computed_from_manufactured_solution(pos);
388  std::array<real,nstate> diss_source_term = dissipative_source_term_computed_from_manufactured_solution(pos,cell_index);
389  std::array<real,nstate> convective_dissipative_source_term;
390  for (int s=0; s<nstate; ++s)
391  {
392  convective_dissipative_source_term[s] = conv_source_term[s] + diss_source_term[s];
393  }
395 }
396 
397 // Returns the value from a CoDiPack or Sacado variable.
398 template<typename real>
399 double getValue(const real &x) {
400  if constexpr(std::is_same<real,double>::value) {
401  return x;
402  }
403  else if constexpr(std::is_same<real,FadType>::value) {
404  return x.val(); // sacado
405  }
406  else if constexpr(std::is_same<real,FadFadType>::value) {
407  return x.val().val(); // sacado
408  }
409  else if constexpr(std::is_same<real,RadType>::value) {
410  return x.value(); // CoDiPack
411  }
412  else if(std::is_same<real,RadFadType>::value) {
413  return x.value().value(); // CoDiPack
414  }
415 }
416 //----------------------------------------------------------------
417 template <int dim, int nstate, typename real>
418 dealii::Tensor<2,nstate,real> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
420  const std::array<real,nstate> &conservative_soln,
421  const dealii::Tensor<1,dim,real> &normal) const
422 {
423  using adtype = FadType;
424 
425  // Initialize AD objects
426  std::array<adtype,nstate> AD_conservative_soln;
427  for (int s=0; s<nstate; ++s) {
428  adtype ADvar(nstate, s, getValue<real>(conservative_soln[s])); // create AD variable
429  AD_conservative_soln[s] = ADvar;
430  }
431 
432  // Compute AD convective flux
433  std::array<dealii::Tensor<1,dim,adtype>,nstate> AD_convective_flux = convective_flux_templated<adtype>(AD_conservative_soln);
434 
435  // Assemble the directional Jacobian
436  dealii::Tensor<2,nstate,real> jacobian;
437  for (int sp=0; sp<nstate; ++sp) {
438  // for each perturbed state (sp) variable
439  for (int s=0; s<nstate; ++s) {
440  jacobian[s][sp] = 0.0;
441  for (int d=0;d<dim;++d) {
442  // Compute directional jacobian
443  jacobian[s][sp] += AD_convective_flux[s][d].dx(sp)*normal[d];
444  }
445  }
446  }
447  return jacobian;
448 }
449 //----------------------------------------------------------------
450 template <int dim, int nstate, typename real>
451 dealii::Tensor<2,nstate,real> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
453  const std::array<real,nstate> &conservative_soln,
454  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
455  const dealii::Tensor<1,dim,real> &normal,
456  const dealii::types::global_dof_index cell_index) const
457 {
458  using adtype = FadType;
459 
460  // Initialize AD objects
461  std::array<adtype,nstate> AD_conservative_soln;
462  std::array<dealii::Tensor<1,dim,adtype>,nstate> AD_solution_gradient;
463  for (int s=0; s<nstate; ++s) {
464  adtype ADvar(nstate, s, getValue<real>(conservative_soln[s])); // create AD variable
465  AD_conservative_soln[s] = ADvar;
466  for (int d=0;d<dim;++d) {
467  AD_solution_gradient[s][d] = getValue<real>(solution_gradient[s][d]);
468  }
469  }
470 
471  // Compute AD dissipative flux
472  std::array<dealii::Tensor<1,dim,adtype>,nstate> AD_dissipative_flux = dissipative_flux_templated<adtype>(AD_conservative_soln, AD_solution_gradient, cell_index);
473 
474  // Assemble the directional Jacobian
475  dealii::Tensor<2,nstate,real> jacobian;
476  for (int sp=0; sp<nstate; ++sp) {
477  // for each perturbed state (sp) variable
478  for (int s=0; s<nstate; ++s) {
479  jacobian[s][sp] = 0.0;
480  for (int d=0;d<dim;++d) {
481  // Compute directional jacobian
482  jacobian[s][sp] += AD_dissipative_flux[s][d].dx(sp)*normal[d];
483  }
484  }
485  }
486  return jacobian;
487 }
488 //----------------------------------------------------------------
489 template <int dim, int nstate, typename real>
490 dealii::Tensor<2,nstate,real> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
492  const std::array<real,nstate> &conservative_soln,
493  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
494  const dealii::Tensor<1,dim,real> &normal,
495  const int d_gradient,
496  const dealii::types::global_dof_index cell_index) const
497 {
498  using adtype = FadType;
499 
500  // Initialize AD objects
501  std::array<adtype,nstate> AD_conservative_soln;
502  std::array<dealii::Tensor<1,dim,adtype>,nstate> AD_solution_gradient;
503  for (int s=0; s<nstate; ++s) {
504  AD_conservative_soln[s] = getValue<real>(conservative_soln[s]);
505  for (int d=0;d<dim;++d) {
506  if(d == d_gradient){
507  adtype ADvar(nstate, s, getValue<real>(solution_gradient[s][d])); // create AD variable
508  AD_solution_gradient[s][d] = ADvar;
509  }
510  else {
511  AD_solution_gradient[s][d] = getValue<real>(solution_gradient[s][d]);
512  }
513  }
514  }
515 
516  // Compute AD dissipative flux
517  std::array<dealii::Tensor<1,dim,adtype>,nstate> AD_dissipative_flux = dissipative_flux_templated<adtype>(AD_conservative_soln, AD_solution_gradient, cell_index);
518 
519  // Assemble the directional Jacobian
520  dealii::Tensor<2,nstate,real> jacobian;
521  for (int sp=0; sp<nstate; ++sp) {
522  // for each perturbed state (sp) variable
523  for (int s=0; s<nstate; ++s) {
524  jacobian[s][sp] = 0.0;
525  for (int d=0;d<dim;++d) {
526  // Compute directional jacobian
527  jacobian[s][sp] += AD_dissipative_flux[s][d].dx(sp)*normal[d];
528  }
529  }
530  }
531  return jacobian;
532 }
533 //----------------------------------------------------------------
534 template <int dim, int nstate, typename real>
537  const dealii::Point<dim,real> &pos) const
538 {
539  std::array<real,nstate> manufactured_solution;
540  for (int s=0; s<nstate; ++s) {
541  manufactured_solution[s] = this->manufactured_solution_function->value (pos, s);
542  if (s==0) {
543  assert(manufactured_solution[s] > 0);
544  }
545  }
546  return manufactured_solution;
547 }
548 //----------------------------------------------------------------
549 template <int dim, int nstate, typename real>
550 std::array<dealii::Tensor<1,dim,real>,nstate> ReynoldsAveragedNavierStokesBase<dim,nstate,real>
552  const dealii::Point<dim,real> &pos) const
553 {
554  std::vector<dealii::Tensor<1,dim,real>> manufactured_solution_gradient_dealii(nstate);
555  this->manufactured_solution_function->vector_gradient(pos,manufactured_solution_gradient_dealii);
556  std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient;
557  for (int d=0;d<dim;++d) {
558  for (int s=0; s<nstate; ++s) {
559  manufactured_solution_gradient[s][d] = manufactured_solution_gradient_dealii[s][d];
560  }
561  }
562  return manufactured_solution_gradient;
563 }
564 //----------------------------------------------------------------
565 template <int dim, int nstate, typename real>
568  const dealii::Point<dim,real> &pos) const
569 {
570  // Get Manufactured Solution values
571  const std::array<real,nstate> manufactured_solution = get_manufactured_solution_value(pos);
572 
573  // Get Manufactured Solution gradient
574  const std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient = get_manufactured_solution_gradient(pos);
575 
576  dealii::Tensor<1,nstate,real> convective_flux_divergence;
577  for (int d=0;d<dim;++d) {
578  dealii::Tensor<1,dim,real> normal;
579  normal[d] = 1.0;
580  const dealii::Tensor<2,nstate,real> jacobian = convective_flux_directional_jacobian(manufactured_solution, normal);
581 
582  //convective_flux_divergence += jacobian*manufactured_solution_gradient[d]; <-- needs second term! (jac wrt gradient)
583  for (int sr = 0; sr < nstate; ++sr) {
584  real jac_grad_row = 0.0;
585  for (int sc = 0; sc < nstate; ++sc) {
586  jac_grad_row += jacobian[sr][sc]*manufactured_solution_gradient[sc][d];
587  }
588  convective_flux_divergence[sr] += jac_grad_row;
589  }
590  }
592  for (int s=0; s<nstate; ++s) {
593  convective_source_term_computed_from_manufactured_solution[s] = convective_flux_divergence[s];
594  }
595 
597 }
598 //----------------------------------------------------------------
599 template <int dim, int nstate, typename real>
602  const dealii::Point<dim,real> &pos,
603  const dealii::types::global_dof_index cell_index) const
604 {
609  // Get Manufactured Solution values
610  const std::array<real,nstate> manufactured_solution = get_manufactured_solution_value(pos); // from Euler
611 
612  // Get Manufactured Solution gradient
613  const std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient = get_manufactured_solution_gradient(pos); // from Euler
614 
615  // Get Manufactured Solution hessian
616  std::array<dealii::SymmetricTensor<2,dim,real>,nstate> manufactured_solution_hessian;
617  for (int s=0; s<nstate; ++s) {
618  dealii::SymmetricTensor<2,dim,real> hessian = this->manufactured_solution_function->hessian(pos,s);
619  for (int dr=0;dr<dim;++dr) {
620  for (int dc=0;dc<dim;++dc) {
621  manufactured_solution_hessian[s][dr][dc] = hessian[dr][dc];
622  }
623  }
624  }
625 
626  // First term -- wrt to the conservative variables
627  // This is similar, should simply provide this function a flux_directional_jacobian() -- could restructure later
628  dealii::Tensor<1,nstate,real> dissipative_flux_divergence;
629  for (int d=0;d<dim;++d) {
630  dealii::Tensor<1,dim,real> normal;
631  normal[d] = 1.0;
632  const dealii::Tensor<2,nstate,real> jacobian = dissipative_flux_directional_jacobian(manufactured_solution, manufactured_solution_gradient, normal, cell_index);
633 
634  // get the directional jacobian wrt gradient
635  std::array<dealii::Tensor<2,nstate,real>,dim> jacobian_wrt_gradient;
636  for (int d_gradient=0;d_gradient<dim;++d_gradient) {
637 
638  // get the directional jacobian wrt gradient component (x,y,z)
639  const dealii::Tensor<2,nstate,real> jacobian_wrt_gradient_component = dissipative_flux_directional_jacobian_wrt_gradient_component(manufactured_solution, manufactured_solution_gradient, normal, d_gradient, cell_index);
640 
641  // store each component in jacobian_wrt_gradient -- could do this in the function used above
642  for (int sr = 0; sr < nstate; ++sr) {
643  for (int sc = 0; sc < nstate; ++sc) {
644  jacobian_wrt_gradient[d_gradient][sr][sc] = jacobian_wrt_gradient_component[sr][sc];
645  }
646  }
647  }
648 
649  //dissipative_flux_divergence += jacobian*manufactured_solution_gradient[d]; <-- needs second term! (jac wrt gradient)
650  for (int sr = 0; sr < nstate; ++sr) {
651  real jac_grad_row = 0.0;
652  for (int sc = 0; sc < nstate; ++sc) {
653  jac_grad_row += jacobian[sr][sc]*manufactured_solution_gradient[sc][d]; // Euler is the same as this
654  // Second term -- wrt to the gradient of conservative variables
655  // -- add the contribution of each gradient component (e.g. x,y,z for dim==3)
656  for (int d_gradient=0;d_gradient<dim;++d_gradient) {
657  jac_grad_row += jacobian_wrt_gradient[d_gradient][sr][sc]*manufactured_solution_hessian[sc][d_gradient][d]; // symmetric so d indexing works both ways
658  }
659  }
660  dissipative_flux_divergence[sr] += jac_grad_row;
661  }
662  }
664  for (int s=0; s<nstate; ++s) {
665  dissipative_source_term_computed_from_manufactured_solution[s] = dissipative_flux_divergence[s];
666  }
667 
669 }
670 //----------------------------------------------------------------
671 template <int dim, int nstate, typename real>
674  const dealii::Point<dim,real> &pos,
675  const dealii::types::global_dof_index cell_index) const
676 {
677  // Get Manufactured Solution values
678  const std::array<real,nstate> manufactured_solution = get_manufactured_solution_value(pos); // from Euler
679 
680  // Get Manufactured Solution gradient
681  const std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient = get_manufactured_solution_gradient(pos); // from Euler
682 
683  std::array<real,nstate> physical_source_source_term_computed_from_manufactured_solution;
684  for (int i=0;i<nstate;++i){
685  physical_source_source_term_computed_from_manufactured_solution = physical_source_term(pos, manufactured_solution, manufactured_solution_gradient, cell_index);
686  }
687 
688  return physical_source_source_term_computed_from_manufactured_solution;
689 }
690 //----------------------------------------------------------------
691 template <int dim, int nstate, typename real>
694  const dealii::Point<dim, real> &pos,
695  const dealii::Tensor<1,dim,real> &/*normal_int*/,
696  const std::array<real,nstate> &/*soln_int*/,
697  const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
698  std::array<real,nstate> &soln_bc,
699  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
700 {
701  // Manufactured solution
702  std::array<real,nstate> boundary_values;
703  std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
704  for (int i=0; i<nstate; ++i) {
705  boundary_values[i] = this->manufactured_solution_function->value (pos, i);
706  boundary_gradients[i] = this->manufactured_solution_function->gradient (pos, i);
707  }
708  for (int istate=nstate_navier_stokes; istate<nstate; ++istate) {
709  soln_bc[istate] = boundary_values[istate];
710  soln_grad_bc[istate] = boundary_gradients[istate];
711  }
712 }
713 //----------------------------------------------------------------
714 //----------------------------------------------------------------
715 //----------------------------------------------------------------
716 // Instantiate explicitly
717 // -- ReynoldsAveragedNavierStokesBase
723 //-------------------------------------------------------------------------------------
724 // Templated members used by derived classes, defined in respective parent classes
725 //-------------------------------------------------------------------------------------
726 // -- get_tensor_magnitude_sqr()
727 template double ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::get_tensor_magnitude_sqr< double >(const dealii::Tensor<2,PHILIP_DIM,double > &tensor) const;
728 template FadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadType >::get_tensor_magnitude_sqr< FadType >(const dealii::Tensor<2,PHILIP_DIM,FadType > &tensor) const;
729 template RadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::get_tensor_magnitude_sqr< RadType >(const dealii::Tensor<2,PHILIP_DIM,RadType > &tensor) const;
730 template FadFadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::get_tensor_magnitude_sqr< FadFadType >(const dealii::Tensor<2,PHILIP_DIM,FadFadType> &tensor) const;
731 template RadFadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::get_tensor_magnitude_sqr< RadFadType >(const dealii::Tensor<2,PHILIP_DIM,RadFadType> &tensor) const;
732 // -- instantiate all the real types with real2 = FadType for automatic differentiation
733 template FadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::get_tensor_magnitude_sqr< FadType >(const dealii::Tensor<2,PHILIP_DIM,FadType> &tensor) const;
734 template FadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::get_tensor_magnitude_sqr< FadType >(const dealii::Tensor<2,PHILIP_DIM,FadType> &tensor) const;
735 template FadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::get_tensor_magnitude_sqr< FadType >(const dealii::Tensor<2,PHILIP_DIM,FadType> &tensor) const;
736 template FadType ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::get_tensor_magnitude_sqr< FadType >(const dealii::Tensor<2,PHILIP_DIM,FadType> &tensor) const;
737 
738 // -- get_vector_magnitude_sqr()
739 template double ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::get_vector_magnitude_sqr< double >(const dealii::Tensor<1,3,double > &vector) const;
744 // -- instantiate all the real types with real2 = FadType for automatic differentiation
749 
750 // -- extract_rans_conservative_solution()
751 template std::array<double, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::extract_rans_conservative_solution< double >(const std::array<double, PHILIP_DIM+3> &conservative_soln) const;
752 template std::array<FadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadType >::extract_rans_conservative_solution< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
753 template std::array<RadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::extract_rans_conservative_solution< RadType >(const std::array<RadType, PHILIP_DIM+3> &conservative_soln) const;
754 template std::array<FadFadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::extract_rans_conservative_solution< FadFadType >(const std::array<FadFadType, PHILIP_DIM+3> &conservative_soln) const;
755 template std::array<RadFadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::extract_rans_conservative_solution< RadFadType >(const std::array<RadFadType, PHILIP_DIM+3> &conservative_soln) const;
756 // -- instantiate all the real types with real2 = FadType for automatic differentiation
757 template std::array<FadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::extract_rans_conservative_solution< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
758 template std::array<FadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::extract_rans_conservative_solution< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
759 template std::array<FadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::extract_rans_conservative_solution< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
760 template std::array<FadType, PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::extract_rans_conservative_solution< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
761 
762 // -- extract_rans_solution_gradient()
763 template std::array<dealii::Tensor<1,PHILIP_DIM,double >,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::extract_rans_solution_gradient< double >(const std::array<dealii::Tensor<1,PHILIP_DIM,double >,PHILIP_DIM+3> &solution_gradient) const;
764 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType >,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadType >::extract_rans_solution_gradient< FadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadType >,PHILIP_DIM+3> &solution_gradient) const;
765 template std::array<dealii::Tensor<1,PHILIP_DIM,RadType >,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::extract_rans_solution_gradient< RadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,RadType >,PHILIP_DIM+3> &solution_gradient) const;
766 template std::array<dealii::Tensor<1,PHILIP_DIM,FadFadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::extract_rans_solution_gradient< FadFadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadFadType>,PHILIP_DIM+3> &solution_gradient) const;
767 template std::array<dealii::Tensor<1,PHILIP_DIM,RadFadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::extract_rans_solution_gradient< RadFadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,RadFadType>,PHILIP_DIM+3> &solution_gradient) const;
768 // -- instantiate all the real types with real2 = FadType for automatic differentiation
769 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::extract_rans_solution_gradient< FadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
770 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::extract_rans_solution_gradient< FadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
771 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::extract_rans_solution_gradient< FadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
772 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+2> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::extract_rans_solution_gradient< FadType >(const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
773 
774 // -- convert_conservative_to_primitive_turbulence_model()
775 template std::array<double, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::convert_conservative_to_primitive_turbulence_model< double >(const std::array<double, PHILIP_DIM+3> &conservative_soln) const;
776 template std::array<FadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadType >::convert_conservative_to_primitive_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
777 template std::array<RadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::convert_conservative_to_primitive_turbulence_model< RadType >(const std::array<RadType, PHILIP_DIM+3> &conservative_soln) const;
778 template std::array<FadFadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::convert_conservative_to_primitive_turbulence_model< FadFadType >(const std::array<FadFadType, PHILIP_DIM+3> &conservative_soln) const;
779 template std::array<RadFadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::convert_conservative_to_primitive_turbulence_model< RadFadType >(const std::array<RadFadType, PHILIP_DIM+3> &conservative_soln) const;
780 // -- instantiate all the real types with real2 = FadType for automatic differentiation
781 template std::array<FadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::convert_conservative_to_primitive_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
782 template std::array<FadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::convert_conservative_to_primitive_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
783 template std::array<FadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::convert_conservative_to_primitive_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
784 template std::array<FadType, 1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::convert_conservative_to_primitive_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln) const;
785 
786 // -- convert_conservative_gradient_to_primitive_gradient_turbulence_model()
787 template std::array<dealii::Tensor<1,PHILIP_DIM,double >,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< double >(const std::array<double, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,double >,PHILIP_DIM+3> &solution_gradient) const;
788 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType >,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadType >,PHILIP_DIM+3> &solution_gradient) const;
789 template std::array<dealii::Tensor<1,PHILIP_DIM,RadType >,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< RadType >(const std::array<RadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,RadType >,PHILIP_DIM+3> &solution_gradient) const;
790 template std::array<dealii::Tensor<1,PHILIP_DIM,FadFadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadFadType >(const std::array<FadFadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadFadType>,PHILIP_DIM+3> &solution_gradient) const;
791 template std::array<dealii::Tensor<1,PHILIP_DIM,RadFadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< RadFadType >(const std::array<RadFadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,RadFadType>,PHILIP_DIM+3> &solution_gradient) const;
792 // -- instantiate all the real types with real2 = FadType for automatic differentiation
793 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, double >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
794 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
795 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, FadFadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
796 template std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,1> ReynoldsAveragedNavierStokesBase < PHILIP_DIM, PHILIP_DIM+3, RadFadType >::convert_conservative_gradient_to_primitive_gradient_turbulence_model< FadType >(const std::array<FadType, PHILIP_DIM+3> &conservative_soln, const std::array<dealii::Tensor<1,PHILIP_DIM,FadType>,PHILIP_DIM+3> &solution_gradient) const;
797 
798 
799 } // Physics namespace
800 } // PHiLiP namespace
virtual std::array< real, nstate > compute_production_dissipation_cross_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &conservative_solution, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient) const =0
Physical source term (production, dissipation source terms and source term with cross derivatives) in...
std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function
Manufactured solution function.
Definition: model.h:29
std::array< real, nstate > convective_dissipative_source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &conservative_solution, const dealii::types::global_dof_index cell_index) const
Convective and dissipative source term for manufactured solution functions.
void boundary_manufactured_solution(const dealii::Point< dim, real > &pos, const dealii::Tensor< 1, dim, real > &normal_int, const std::array< real, nstate > &soln_int, const std::array< dealii::Tensor< 1, dim, real >, nstate > &soln_grad_int, std::array< real, nstate > &soln_bc, std::array< dealii::Tensor< 1, dim, real >, nstate > &soln_grad_bc) const override
Evaluate the manufactured solution boundary conditions.
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &conservative_solution, const real current_time, const dealii::types::global_dof_index cell_index) const
Source term for manufactured solution functions.
std::array< real2, nstate-(dim+2)> convert_conservative_to_primitive_turbulence_model(const std::array< real2, nstate > &conservative_soln) const
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
Definition: ADTypes.hpp:12
dealii::Tensor< 2, nstate, real > dissipative_flux_directional_jacobian(const std::array< real, nstate > &conservative_soln, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient, const dealii::Tensor< 1, dim, real > &normal, const dealii::types::global_dof_index cell_index) const
virtual std::array< real, nstate_turbulence_model > compute_effective_viscosity_turbulence_model(const std::array< real, nstate_navier_stokes > &primitive_soln_rans, const std::array< real, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized effective (total) viscosities for the turbulence model.
Sacado::Fad::DFad< double > FadType
Sacado AD type for first derivatives.
Definition: ADTypes.hpp:11
Manufactured solution used for grid studies to check convergence orders.
codi_JacobianComputationType RadType
CoDiPaco reverse-AD type for first derivatives.
Definition: ADTypes.hpp:27
ReynoldsAveragedNavierStokesBase(const Parameters::AllParameters *const parameters_input, const double ref_length, const double gamma_gas, const double mach_inf, const double angle_of_attack, const double side_slip_angle, const double prandtl_number, const double reynolds_number_inf, const bool use_constant_viscosity, const double constant_viscosity, const double turbulent_prandtl_number, const double temperature_inf=273.15, const double isothermal_wall_temperature=1.0, const thermal_boundary_condition_enum thermal_boundary_condition_type=thermal_boundary_condition_enum::adiabatic, std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function=nullptr, const two_point_num_flux_enum two_point_num_flux_type=two_point_num_flux_enum::KG)
Constructor.
std::array< dealii::Tensor< 1, dim, real2 >, nstate-(dim+2)> dissipative_flux_turbulence_model(const std::array< real2, nstate_navier_stokes > &primitive_soln_rans, const std::array< real2, nstate_turbulence_model > &primitive_soln_turbulence_model, const std::array< dealii::Tensor< 1, dim, real2 >, nstate_turbulence_model > &primitive_solution_gradient_turbulence_model) const
Templated Additional viscous flux of RANS + viscous flux of turbulence model.
std::array< real, nstate > dissipative_source_term_computed_from_manufactured_solution(const dealii::Point< dim, real > &pos, const dealii::types::global_dof_index cell_index) const
std::array< real, nstate > get_manufactured_solution_value(const dealii::Point< dim, real > &pos) const
Get manufactured solution value.
std::array< real, nstate > physical_source_term_computed_from_manufactured_solution(const dealii::Point< dim, real > &pos, const dealii::types::global_dof_index cell_index) const
Files for the baseline physics.
Definition: ADTypes.hpp:10
std::array< real, nstate-(dim+2)> compute_mean_turbulence_property(const std::array< real, nstate > &conservative_soln1, const std::array< real, nstate > &conservative_soln2) const
Mean turbulence properties given two sets of conservative solutions.
real max_convective_eigenvalue(const std::array< real, nstate > &soln) const
Maximum convective eigenvalue of the additional models&#39; PDEs.
virtual std::array< FadType, nstate_turbulence_model > compute_effective_viscosity_turbulence_model_fad(const std::array< FadType, nstate_navier_stokes > &primitive_soln_rans, const std::array< FadType, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized effective (total) viscosities for the turbulence model (Automatic Differentiation ...
virtual dealii::Tensor< 2, dim, FadType > compute_Reynolds_stress_tensor_fad(const std::array< FadType, nstate_navier_stokes > &primitive_soln_rans, const std::array< dealii::Tensor< 1, dim, FadType >, nstate_navier_stokes > &primitive_soln_gradient_rans, const std::array< FadType, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized Reynolds stress tensor, (tau^reynolds)* (Automatic Differentiation Type: FadType) ...
Physics model additional terms and equations to the baseline physics.
Definition: model.h:18
Main parameter class that contains the various other sub-parameter classes.
dealii::Tensor< 2, nstate, real > dissipative_flux_directional_jacobian_wrt_gradient_component(const std::array< real, nstate > &conservative_soln, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient, const dealii::Tensor< 1, dim, real > &normal, const int d_gradient, const dealii::types::global_dof_index cell_index) const
std::array< dealii::Tensor< 1, dim, real2 >, dim+2 > extract_rans_solution_gradient(const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &solution_gradient) const
Returns the conservative solutions gradient of Reynolds-averaged Navier-Stokes equations (without add...
std::array< dealii::Tensor< 1, dim, real >, nstate > convective_flux(const std::array< real, nstate > &conservative_soln) const
Additional convective flux of RANS + convective flux of turbulence model.
virtual dealii::Tensor< 1, dim, FadType > compute_Reynolds_heat_flux_fad(const std::array< FadType, nstate_navier_stokes > &primitive_soln_rans, const std::array< dealii::Tensor< 1, dim, FadType >, nstate_navier_stokes > &primitive_soln_gradient_rans, const std::array< FadType, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized Reynolds heat flux, (q^reynolds)* (Automatic Differentiation Type: FadType) ...
TwoPointNumericalFlux
Two point numerical flux type for split form.
std::array< dealii::Tensor< 1, dim, real2 >, nstate > convective_flux_templated(const std::array< real2, nstate > &conservative_soln) const
Templated additional convective flux.
std::array< dealii::Tensor< 1, dim, real2 >, nstate > dissipative_flux_templated(const std::array< real2, nstate > &conservative_soln, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &solution_gradient, const dealii::types::global_dof_index cell_index) const
Templated additional dissipative (i.e. viscous) flux.
virtual dealii::Tensor< 2, dim, real > compute_Reynolds_stress_tensor(const std::array< real, nstate_navier_stokes > &primitive_soln_rans, const std::array< dealii::Tensor< 1, dim, real >, nstate_navier_stokes > &primitive_soln_gradient_rans, const std::array< real, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized Reynolds stress tensor, (tau^reynolds)*.
std::unique_ptr< NavierStokes< dim, nstate_navier_stokes, real > > navier_stokes_physics
Pointer to Navier-Stokes physics object.
static const int nstate_navier_stokes
Number of PDEs for RANS equations.
Reynolds-Averaged Navier-Stokes (RANS) equations. Derived from Navier-Stokes for modifying the stress...
std::array< real, nstate > convective_eigenvalues(const std::array< real, nstate > &, const dealii::Tensor< 1, dim, real > &) const
Convective eigenvalues of the additional models&#39; PDEs.
real max_convective_normal_eigenvalue(const std::array< real, nstate > &soln, const dealii::Tensor< 1, dim, real > &normal) const
Maximum convective normal eigenvalue (used in Lax-Friedrichs) of the additional models&#39; PDEs...
dealii::Tensor< 2, nstate, real > convective_flux_directional_jacobian(const std::array< real, nstate > &conservative_soln, const dealii::Tensor< 1, dim, real > &normal) const
real2 get_vector_magnitude_sqr(const dealii::Tensor< 1, 3, real2 > &vector) const
Returns the square of the magnitude of the vector.
virtual dealii::Tensor< 1, dim, real > compute_Reynolds_heat_flux(const std::array< real, nstate_navier_stokes > &primitive_soln_rans, const std::array< dealii::Tensor< 1, dim, real >, nstate_navier_stokes > &primitive_soln_gradient_rans, const std::array< real, nstate_turbulence_model > &primitive_soln_turbulence_model) const =0
Nondimensionalized Reynolds heat flux, (q^reynolds)*.
std::array< real, nstate > convective_source_term_computed_from_manufactured_solution(const dealii::Point< dim, real > &pos) const
static const int nstate_turbulence_model
Number of PDEs for RANS turbulence model.
ThermalBoundaryCondition
Types of thermal boundary conditions available.
std::array< dealii::Tensor< 1, dim, real2 >, nstate-(dim+2)> convert_conservative_gradient_to_primitive_gradient_turbulence_model(const std::array< real2, nstate > &conservative_soln, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &solution_gradient) const
real2 get_tensor_magnitude_sqr(const dealii::Tensor< 2, dim, real2 > &tensor) const
Returns the square of the magnitude of the tensor (i.e. the double dot product of a tensor with itsel...
std::array< real2, dim+2 > extract_rans_conservative_solution(const std::array< real2, nstate > &conservative_soln) const
Returns the conservative solutions of Reynolds-averaged Navier-Stokes equations (without additional R...
codi_HessianComputationType RadFadType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
Definition: ADTypes.hpp:28
Navier-Stokes equations. Derived from Euler for the convective terms, which is derived from PhysicsBa...
Definition: navier_stokes.h:12
std::array< dealii::Tensor< 1, dim, real >, nstate > get_manufactured_solution_gradient(const dealii::Point< dim, real > &pos) const
Get manufactured solution value.
std::array< real, nstate > physical_source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &conservative_solution, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient, const dealii::types::global_dof_index cell_index) const override
Physical source term.
std::array< dealii::Tensor< 1, dim, real >, nstate > dissipative_flux(const std::array< real, nstate > &conservative_soln, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient, const dealii::types::global_dof_index cell_index) const
Additional viscous flux of RANS + viscous flux of turbulence model.