ExaDG
Loading...
Searching...
No Matches
kernels_and_operators.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef INCLUDE_EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_KERNELS_AND_OPERATORS_H_
23#define INCLUDE_EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_KERNELS_AND_OPERATORS_H_
24
25// C/C++
26#include <iostream>
27
28// deal.II
29#include <deal.II/lac/la_parallel_vector.h>
30
31// ExaDG
32#include <exadg/compressible_navier_stokes/user_interface/boundary_descriptor.h>
33#include <exadg/compressible_navier_stokes/user_interface/parameters.h>
34#include <exadg/functions_and_boundary_conditions/evaluate_functions.h>
35#include <exadg/matrix_free/integrators.h>
36#include <exadg/operators/interior_penalty_parameter.h>
37
38namespace ExaDG
39{
40namespace CompNS
41{
42template<int dim, typename Number>
43inline DEAL_II_ALWAYS_INLINE //
44 dealii::VectorizedArray<Number>
45 calculate_pressure(dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & rho_u,
46 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u,
47 dealii::VectorizedArray<Number> const & rho_E,
48 Number const & gamma)
49{
50 return (gamma - 1.0) * (rho_E - 0.5 * scalar_product(rho_u, u));
51}
52
53template<int dim, typename Number>
54inline DEAL_II_ALWAYS_INLINE //
55 dealii::VectorizedArray<Number>
56 calculate_pressure(dealii::VectorizedArray<Number> const & rho,
57 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u,
58 dealii::VectorizedArray<Number> const & E,
59 Number const & gamma)
60{
61 return (gamma - 1.0) * rho * (E - 0.5 * scalar_product(u, u));
62}
63
64template<typename Number>
65inline DEAL_II_ALWAYS_INLINE //
66 dealii::VectorizedArray<Number>
67 calculate_temperature(dealii::VectorizedArray<Number> const & p,
68 dealii::VectorizedArray<Number> const & rho,
69 Number const & R)
70{
71 return p / (rho * R);
72}
73
74template<int dim, typename Number>
75inline dealii::VectorizedArray<Number>
76calculate_energy(dealii::VectorizedArray<Number> const & T,
77 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u,
78 Number const & c_v)
79{
80 return c_v * T + 0.5 * scalar_product(u, u);
81}
82
83template<int dim, typename Number>
84inline DEAL_II_ALWAYS_INLINE //
85 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
86 calculate_grad_E(dealii::VectorizedArray<Number> const & rho_inverse,
87 dealii::VectorizedArray<Number> const & rho_E,
88 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & grad_rho,
89 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & grad_rho_E)
90{
91 dealii::VectorizedArray<Number> E = rho_inverse * rho_E;
92
93 return rho_inverse * (grad_rho_E - E * grad_rho);
94}
95
96template<int dim, typename Number>
97inline DEAL_II_ALWAYS_INLINE //
98 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
99 calculate_grad_u(dealii::VectorizedArray<Number> const & rho_inverse,
100 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & rho_u,
101 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & grad_rho,
102 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & grad_rho_u)
103{
104 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> out;
105 for(unsigned int d = 0; d < dim; ++d)
106 {
107 dealii::VectorizedArray<Number> ud = rho_inverse * rho_u[d];
108 for(unsigned int e = 0; e < dim; ++e)
109 out[d][e] = rho_inverse * (grad_rho_u[d][e] - ud * grad_rho[e]);
110 }
111
112 return out;
113}
114
115template<int dim, typename Number>
116inline DEAL_II_ALWAYS_INLINE //
117 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
118 calculate_grad_T(dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & grad_E,
119 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u,
120 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & grad_u,
121 Number const & gamma,
122 Number const & R)
123{
124 return (gamma - 1.0) / R * (grad_E - u * grad_u);
125}
126
127template<int dim, typename Number>
128inline DEAL_II_ALWAYS_INLINE //
129 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
130 calculate_stress_tensor(dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & grad_u,
131 Number const & viscosity)
132{
133 dealii::VectorizedArray<Number> const divu = (2. / 3.) * trace(grad_u);
134
135 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> out;
136 for(unsigned int d = 0; d < dim; ++d)
137 {
138 for(unsigned int e = 0; e < dim; ++e)
139 out[d][e] = viscosity * (grad_u[d][e] + grad_u[e][d]);
140 out[d][d] -= viscosity * divu;
141 }
142
143 return out;
144}
145
146/*
147 * Calculates exterior state "+" for a scalar/vectorial quantity depending on interior state "-" and
148 * boundary conditions.
149 */
150template<int dim, typename Number, int rank>
151inline DEAL_II_ALWAYS_INLINE //
152 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
153 calculate_exterior_value(
154 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>> const & value_m,
155 BoundaryType const boundary_type,
156 BoundaryDescriptorStd<dim> const & boundary_descriptor,
157 dealii::types::boundary_id const & boundary_id,
158 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_point,
159 Number const & time)
160{
161 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>> value_p;
162
163 if(boundary_type == BoundaryType::Dirichlet)
164 {
165 auto bc = boundary_descriptor.dirichlet_bc.find(boundary_id)->second;
166 auto g = FunctionEvaluator<rank, dim, Number>::value(*bc, q_point, time);
167
168 value_p = -value_m + dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>(2.0 * g);
169 }
170 else if(boundary_type == BoundaryType::Neumann)
171 {
172 value_p = value_m;
173 }
174 else
175 {
176 AssertThrow(false, dealii::ExcMessage("Not implemented."));
177 }
178
179 return value_p;
180}
181
182/*
183 * Calculates exterior state of normal gradient (Neumann type boundary conditions)
184 * depending on interior data and boundary conditions.
185 */
186template<int dim, typename Number, int rank>
187inline DEAL_II_ALWAYS_INLINE //
188 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>
189 calculate_exterior_normal_grad(
190 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>> const & grad_M_normal,
191 BoundaryType const & boundary_type,
192 BoundaryDescriptorStd<dim> const & boundary_descriptor,
193 dealii::types::boundary_id const & boundary_id,
194 dealii::Point<dim, dealii::VectorizedArray<Number>> const & q_point,
195 Number const & time)
196{
197 dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>> grad_P_normal;
198
199 if(boundary_type == BoundaryType::Dirichlet)
200 {
201 // do nothing
202 grad_P_normal = grad_M_normal;
203 }
204 else if(boundary_type == BoundaryType::Neumann)
205 {
206 auto bc = boundary_descriptor.neumann_bc.find(boundary_id)->second;
207 auto h = FunctionEvaluator<rank, dim, Number>::value(*bc, q_point, time);
208
209 grad_P_normal =
210 -grad_M_normal + dealii::Tensor<rank, dim, dealii::VectorizedArray<Number>>(2.0 * h);
211 }
212 else
213 {
214 AssertThrow(false, dealii::ExcMessage("Not implemented."));
215 }
216
217 return grad_P_normal;
218}
219
220/*
221 * This function calculates the Lax-Friedrichs flux for the momentum equation
222 */
223template<int dim, typename Number>
224inline DEAL_II_ALWAYS_INLINE //
225 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
226 calculate_flux(dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & momentum_flux_M,
227 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & momentum_flux_P,
228 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & rho_u_M,
229 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & rho_u_P,
230 dealii::VectorizedArray<Number> const & lambda,
231 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normal)
232{
233 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> out;
234 for(unsigned int d = 0; d < dim; ++d)
235 {
236 dealii::VectorizedArray<Number> sum = dealii::VectorizedArray<Number>();
237 for(unsigned int e = 0; e < dim; ++e)
238 sum += (momentum_flux_M[d][e] + momentum_flux_P[d][e]) * normal[e];
239 out[d] = 0.5 * (sum + lambda * (rho_u_M[d] - rho_u_P[d]));
240 }
241
242 return out;
243}
244
245/*
246 * This function calculates the Lax-Friedrichs flux for scalar quantities (density/energy)
247 */
248template<int dim, typename Number>
249inline DEAL_II_ALWAYS_INLINE //
250 dealii::VectorizedArray<Number>
251 calculate_flux(dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & flux_M,
252 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & flux_P,
253 dealii::VectorizedArray<Number> const & value_M,
254 dealii::VectorizedArray<Number> const & value_P,
255 dealii::VectorizedArray<Number> const & lambda,
256 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normal)
257{
258 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> average_flux = 0.5 * (flux_M + flux_P);
259
260 return average_flux * normal + 0.5 * lambda * (value_M - value_P);
261}
262
263/*
264 * Calculation of lambda for Lax-Friedrichs flux according to Hesthaven:
265 * lambda = max( |u_M| + sqrt(|gamma * p_M / rho_M|) , |u_P| + sqrt(|gamma * p_P / rho_P|) )
266 */
267template<int dim, typename Number>
268inline DEAL_II_ALWAYS_INLINE //
269 dealii::VectorizedArray<Number>
270 calculate_lambda(dealii::VectorizedArray<Number> const & rho_m,
271 dealii::VectorizedArray<Number> const & rho_p,
272 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u_m,
273 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u_p,
274 dealii::VectorizedArray<Number> const & p_m,
275 dealii::VectorizedArray<Number> const & p_p,
276 Number const & gamma)
277{
278 dealii::VectorizedArray<Number> lambda_m = u_m.norm() + std::sqrt(std::abs(gamma * p_m / rho_m));
279 dealii::VectorizedArray<Number> lambda_p = u_p.norm() + std::sqrt(std::abs(gamma * p_p / rho_p));
280
281 return std::max(lambda_m, lambda_p);
282}
283
284template<int dim>
285struct BodyForceOperatorData
286{
287 BodyForceOperatorData() : dof_index(0), quad_index(0)
288 {
289 }
290
291 unsigned int dof_index;
292 unsigned int quad_index;
293
294 std::shared_ptr<dealii::Function<dim>> rhs_rho;
295 std::shared_ptr<dealii::Function<dim>> rhs_u;
296 std::shared_ptr<dealii::Function<dim>> rhs_E;
297};
298
299template<int dim, typename Number>
300class BodyForceOperator
301{
302public:
303 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
304
305 typedef BodyForceOperator<dim, Number> This;
306
307 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
308 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
309
310 typedef dealii::VectorizedArray<Number> scalar;
311 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
312 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
313 typedef dealii::Point<dim, dealii::VectorizedArray<Number>> point;
314
315 BodyForceOperator() : matrix_free(nullptr), eval_time(0.0)
316 {
317 }
318
319 void
320 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
321 BodyForceOperatorData<dim> const & data_in)
322 {
323 this->matrix_free = &matrix_free_in;
324 this->data = data_in;
325 }
326
327 void
328 evaluate(VectorType & dst, VectorType const & src, double const evaluation_time) const
329 {
330 dst = 0;
331 evaluate_add(dst, src, evaluation_time);
332 }
333
334 void
335 evaluate_add(VectorType & dst, VectorType const & src, double const evaluation_time) const
336 {
337 this->eval_time = evaluation_time;
338
339 matrix_free->cell_loop(&This::cell_loop, this, dst, src);
340 }
341
342 inline DEAL_II_ALWAYS_INLINE //
343 std::tuple<scalar, vector, scalar>
344 get_volume_flux(CellIntegratorScalar & density,
345 CellIntegratorVector & momentum,
346 unsigned int const q) const
347 {
348 point q_points = density.quadrature_point(q);
349 scalar rho = density.get_value(q);
350 vector u = momentum.get_value(q) / rho;
351
352 scalar rhs_density =
353 FunctionEvaluator<0, dim, Number>::value(*(data.rhs_rho), q_points, eval_time);
354 vector rhs_momentum =
355 FunctionEvaluator<1, dim, Number>::value(*(data.rhs_u), q_points, eval_time);
356 scalar rhs_energy =
357 FunctionEvaluator<0, dim, Number>::value(*(data.rhs_E), q_points, eval_time);
358
359 return std::make_tuple(rhs_density, rhs_momentum, rhs_momentum * u + rhs_energy);
360 }
361
362private:
363 void
364 cell_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
365 VectorType & dst,
366 VectorType const & src,
367 std::pair<unsigned int, unsigned int> const & cell_range) const
368 {
369 CellIntegratorScalar density(matrix_free, data.dof_index, data.quad_index, 0);
370 CellIntegratorVector momentum(matrix_free, data.dof_index, data.quad_index, 1);
371 CellIntegratorScalar energy(matrix_free, data.dof_index, data.quad_index, 1 + dim);
372
373 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
374 {
375 density.reinit(cell);
376 density.gather_evaluate(src, dealii::EvaluationFlags::values);
377
378 momentum.reinit(cell);
379 momentum.gather_evaluate(src, dealii::EvaluationFlags::values);
380
381 energy.reinit(cell);
382
383 for(unsigned int q = 0; q < density.n_q_points; ++q)
384 {
385 std::tuple<scalar, vector, scalar> flux = get_volume_flux(density, momentum, q);
386
387 density.submit_value(std::get<0>(flux), q);
388 momentum.submit_value(std::get<1>(flux), q);
389 energy.submit_value(std::get<2>(flux), q);
390 }
391
392 density.integrate_scatter(dealii::EvaluationFlags::values, dst);
393 momentum.integrate_scatter(dealii::EvaluationFlags::values, dst);
394 energy.integrate_scatter(dealii::EvaluationFlags::values, dst);
395 }
396 }
397
398 dealii::MatrixFree<dim, Number> const * matrix_free;
399
400 BodyForceOperatorData<dim> data;
401
402 double mutable eval_time;
403};
404
405struct MassOperatorData
406{
407 MassOperatorData() : dof_index(0), quad_index(0)
408 {
409 }
410
411 unsigned int dof_index;
412 unsigned int quad_index;
413};
414
415template<int dim, typename Number>
416class MassOperator
417{
418public:
419 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
420
421 typedef MassOperator<dim, Number> This;
422
423 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
424 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
425
426 MassOperator() : matrix_free(nullptr)
427 {
428 }
429
430 void
431 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
432 MassOperatorData const & data_in)
433 {
434 this->matrix_free = &matrix_free_in;
435 this->data = data_in;
436 }
437
438 // apply matrix vector multiplication
439 void
440 apply(VectorType & dst, VectorType const & src) const
441 {
442 dst = 0;
443 apply_add(dst, src);
444 }
445
446 void
447 apply_add(VectorType & dst, VectorType const & src) const
448 {
449 matrix_free->cell_loop(&This::cell_loop, this, dst, src);
450 }
451
452private:
453 void
454 cell_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
455 VectorType & dst,
456 VectorType const & src,
457 std::pair<unsigned int, unsigned int> const & cell_range) const
458 {
459 CellIntegratorScalar density(matrix_free, data.dof_index, data.quad_index, 0);
460 CellIntegratorVector momentum(matrix_free, data.dof_index, data.quad_index, 1);
461 CellIntegratorScalar energy(matrix_free, data.dof_index, data.quad_index, 1 + dim);
462
463 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
464 {
465 density.reinit(cell);
466 density.gather_evaluate(src, dealii::EvaluationFlags::values);
467
468 momentum.reinit(cell);
469 momentum.gather_evaluate(src, dealii::EvaluationFlags::values);
470
471 energy.reinit(cell);
472 energy.gather_evaluate(src, dealii::EvaluationFlags::values);
473
474 for(unsigned int q = 0; q < density.n_q_points; ++q)
475 {
476 density.submit_value(density.get_value(q), q);
477 momentum.submit_value(momentum.get_value(q), q);
478 energy.submit_value(energy.get_value(q), q);
479 }
480
481 density.integrate_scatter(dealii::EvaluationFlags::values, dst);
482 momentum.integrate_scatter(dealii::EvaluationFlags::values, dst);
483 energy.integrate_scatter(dealii::EvaluationFlags::values, dst);
484 }
485 }
486
487 dealii::MatrixFree<dim, Number> const * matrix_free;
488 MassOperatorData data;
489};
490
491template<int dim>
492struct ConvectiveOperatorData
493{
494 ConvectiveOperatorData()
495 : dof_index(0), quad_index(0), heat_capacity_ratio(1.4), specific_gas_constant(287.0)
496 {
497 }
498
499 unsigned int dof_index;
500 unsigned int quad_index;
501
502 std::shared_ptr<BoundaryDescriptor<dim> const> bc;
503
504 double heat_capacity_ratio;
505 double specific_gas_constant;
506};
507
508template<int dim, typename Number>
509class ConvectiveOperator
510{
511public:
512 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
513
514 typedef ConvectiveOperator<dim, Number> This;
515
516 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
517 typedef FaceIntegrator<dim, 1, Number> FaceIntegratorScalar;
518 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
519 typedef FaceIntegrator<dim, dim, Number> FaceIntegratorVector;
520
521 typedef dealii::VectorizedArray<Number> scalar;
522 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
523 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
524 typedef dealii::Point<dim, dealii::VectorizedArray<Number>> point;
525
526 ConvectiveOperator() : matrix_free(nullptr)
527 {
528 }
529
530 void
531 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
532 ConvectiveOperatorData<dim> const & data_in)
533 {
534 this->matrix_free = &matrix_free_in;
535 this->data = data_in;
536
537 gamma = data.heat_capacity_ratio;
538 R = data.specific_gas_constant;
539 c_v = R / (gamma - 1.0);
540 }
541
542 void
543 evaluate(VectorType & dst, VectorType const & src, Number const evaluation_time) const
544 {
545 dst = 0;
546 evaluate_add(dst, src, evaluation_time);
547 }
548
549 void
550 evaluate_add(VectorType & dst, VectorType const & src, Number const evaluation_time) const
551 {
552 this->eval_time = evaluation_time;
553
554 matrix_free->loop(
555 &This::cell_loop, &This::face_loop, &This::boundary_face_loop, this, dst, src);
556 }
557
558 void
559 set_evaluation_time(double const & evaluation_time) const
560 {
561 eval_time = evaluation_time;
562 }
563
564 inline DEAL_II_ALWAYS_INLINE //
565 std::tuple<vector, tensor, vector>
566 get_volume_flux(CellIntegratorScalar & density,
567 CellIntegratorVector & momentum,
568 CellIntegratorScalar & energy,
569 unsigned int const q) const
570 {
571 scalar rho_inv = 1.0 / density.get_value(q);
572 vector rho_u = momentum.get_value(q);
573 scalar rho_E = energy.get_value(q);
574 vector u = rho_inv * rho_u;
575 scalar p = calculate_pressure(rho_u, u, rho_E, gamma);
576
577 tensor momentum_flux = outer_product(rho_u, u);
578 for(unsigned int d = 0; d < dim; ++d)
579 momentum_flux[d][d] += p;
580
581 vector energy_flux = (rho_E + p) * u;
582
583 return std::make_tuple(rho_u, momentum_flux, energy_flux);
584 }
585
586 inline DEAL_II_ALWAYS_INLINE //
587 std::tuple<scalar, vector, scalar>
588 get_flux(FaceIntegratorScalar & density_m,
589 FaceIntegratorScalar & density_p,
590 FaceIntegratorVector & momentum_m,
591 FaceIntegratorVector & momentum_p,
592 FaceIntegratorScalar & energy_m,
593 FaceIntegratorScalar & energy_p,
594 unsigned int const q) const
595 {
596 vector normal = momentum_m.get_normal_vector(q);
597
598 // get values
599 scalar rho_M = density_m.get_value(q);
600 scalar rho_P = density_p.get_value(q);
601 vector rho_u_M = momentum_m.get_value(q);
602 vector rho_u_P = momentum_p.get_value(q);
603 vector u_M = rho_u_M / rho_M;
604 vector u_P = rho_u_P / rho_P;
605 scalar rho_E_M = energy_m.get_value(q);
606 scalar rho_E_P = energy_p.get_value(q);
607
608 // calculate pressure
609 scalar p_M = calculate_pressure(rho_u_M, u_M, rho_E_M, gamma);
610 scalar p_P = calculate_pressure(rho_u_P, u_P, rho_E_P, gamma);
611
612 // calculate lambda
613 scalar lambda = calculate_lambda(rho_M, rho_P, u_M, u_P, p_M, p_P, gamma);
614
615 // flux density
616 scalar flux_density = calculate_flux(rho_u_M, rho_u_P, rho_M, rho_P, lambda, normal);
617
618 // flux momentum
619 tensor momentum_flux_M = outer_product(rho_u_M, u_M);
620 for(unsigned int d = 0; d < dim; ++d)
621 momentum_flux_M[d][d] += p_M;
622
623 tensor momentum_flux_P = outer_product(rho_u_P, u_P);
624 for(unsigned int d = 0; d < dim; ++d)
625 momentum_flux_P[d][d] += p_P;
626
627 vector flux_momentum =
628 calculate_flux(momentum_flux_M, momentum_flux_P, rho_u_M, rho_u_P, lambda, normal);
629
630 // flux energy
631 vector energy_flux_M = (rho_E_M + p_M) * u_M;
632 vector energy_flux_P = (rho_E_P + p_P) * u_P;
633
634 scalar flux_energy =
635 calculate_flux(energy_flux_M, energy_flux_P, rho_E_M, rho_E_P, lambda, normal);
636
637 return std::make_tuple(flux_density, flux_momentum, flux_energy);
638 }
639
640 inline DEAL_II_ALWAYS_INLINE //
641 std::tuple<scalar, vector, scalar>
642 get_flux_boundary(FaceIntegratorScalar & density,
643 FaceIntegratorVector & momentum,
644 FaceIntegratorScalar & energy,
645 BoundaryType const & boundary_type_density,
646 BoundaryType const & boundary_type_velocity,
647 BoundaryType const & boundary_type_pressure,
648 BoundaryType const & boundary_type_energy,
649 EnergyBoundaryVariable const & boundary_variable,
650 dealii::types::boundary_id const & boundary_id,
651 unsigned int const q) const
652 {
653 vector normal = momentum.get_normal_vector(q);
654
655 // element e⁻
656 scalar rho_M = density.get_value(q);
657 vector rho_u_M = momentum.get_value(q);
658 vector u_M = rho_u_M / rho_M;
659 scalar rho_E_M = energy.get_value(q);
660 scalar E_M = rho_E_M / rho_M;
661 scalar p_M = calculate_pressure(rho_M, u_M, E_M, gamma);
662
663 // element e⁺
664
665 // calculate rho_P
666 scalar rho_P = calculate_exterior_value<dim, Number, 0>(rho_M,
667 boundary_type_density,
668 data.bc->density,
669 boundary_id,
670 density.quadrature_point(q),
671 this->eval_time);
672
673 // calculate u_P
674 vector u_P = calculate_exterior_value<dim, Number, 1>(u_M,
675 boundary_type_velocity,
676 data.bc->velocity,
677 boundary_id,
678 momentum.quadrature_point(q),
679 this->eval_time);
680
681 vector rho_u_P = rho_P * u_P;
682
683 // calculate p_P
684 scalar p_P = calculate_exterior_value<dim, Number, 0>(p_M,
685 boundary_type_pressure,
686 data.bc->pressure,
687 boundary_id,
688 density.quadrature_point(q),
689 this->eval_time);
690
691 // calculate E_P
692 scalar E_P = dealii::make_vectorized_array<Number>(0.0);
693 if(boundary_variable == EnergyBoundaryVariable::Energy)
694 {
695 E_P = calculate_exterior_value<dim, Number, 0>(E_M,
696 boundary_type_energy,
697 data.bc->energy,
698 boundary_id,
699 energy.quadrature_point(q),
700 this->eval_time);
701 }
702 else if(boundary_variable == EnergyBoundaryVariable::Temperature)
703 {
704 scalar T_M = calculate_temperature(p_M, rho_M, R);
705 scalar T_P = calculate_exterior_value<dim, Number, 0>(T_M,
706 boundary_type_energy,
707 data.bc->energy,
708 boundary_id,
709 energy.quadrature_point(q),
710 this->eval_time);
711
712 E_P = calculate_energy(T_P, u_P, c_v);
713 }
714 scalar rho_E_P = rho_P * E_P;
715
716 // calculate lambda
717 scalar lambda = calculate_lambda(rho_M, rho_P, u_M, u_P, p_M, p_P, gamma);
718
719 // flux density
720 scalar flux_density = calculate_flux(rho_u_M, rho_u_P, rho_M, rho_P, lambda, normal);
721
722 // flux momentum
723 tensor momentum_flux_M = outer_product(rho_u_M, u_M);
724 for(unsigned int d = 0; d < dim; ++d)
725 momentum_flux_M[d][d] += p_M;
726
727 tensor momentum_flux_P = outer_product(rho_u_P, u_P);
728 for(unsigned int d = 0; d < dim; ++d)
729 momentum_flux_P[d][d] += p_P;
730
731 vector flux_momentum =
732 calculate_flux(momentum_flux_M, momentum_flux_P, rho_u_M, rho_u_P, lambda, normal);
733
734 // flux energy
735 vector energy_flux_M = (rho_E_M + p_M) * u_M;
736 vector energy_flux_P = (rho_E_P + p_P) * u_P;
737 scalar flux_energy =
738 calculate_flux(energy_flux_M, energy_flux_P, rho_E_M, rho_E_P, lambda, normal);
739
740 return std::make_tuple(flux_density, flux_momentum, flux_energy);
741 }
742
743private:
744 void
745 cell_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
746 VectorType & dst,
747 VectorType const & src,
748 std::pair<unsigned int, unsigned int> const & cell_range) const
749 {
750 CellIntegratorScalar density(matrix_free, data.dof_index, data.quad_index, 0);
751 CellIntegratorVector momentum(matrix_free, data.dof_index, data.quad_index, 1);
752 CellIntegratorScalar energy(matrix_free, data.dof_index, data.quad_index, 1 + dim);
753
754 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
755 {
756 density.reinit(cell);
757 density.gather_evaluate(src, dealii::EvaluationFlags::values);
758
759 momentum.reinit(cell);
760 momentum.gather_evaluate(src, dealii::EvaluationFlags::values);
761
762 energy.reinit(cell);
763 energy.gather_evaluate(src, dealii::EvaluationFlags::values);
764
765 for(unsigned int q = 0; q < momentum.n_q_points; ++q)
766 {
767 std::tuple<vector, tensor, vector> flux = get_volume_flux(density, momentum, energy, q);
768
769 density.submit_gradient(-std::get<0>(flux), q);
770 momentum.submit_gradient(-std::get<1>(flux), q);
771 energy.submit_gradient(-std::get<2>(flux), q);
772 }
773
774 density.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
775 momentum.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
776 energy.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
777 }
778 }
779
780 void
781 face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
782 VectorType & dst,
783 VectorType const & src,
784 std::pair<unsigned int, unsigned int> const & face_range) const
785 {
786 FaceIntegratorScalar density_m(matrix_free, true, data.dof_index, data.quad_index, 0);
787 FaceIntegratorScalar density_p(matrix_free, false, data.dof_index, data.quad_index, 0);
788
789 FaceIntegratorVector momentum_m(matrix_free, true, data.dof_index, data.quad_index, 1);
790 FaceIntegratorVector momentum_p(matrix_free, false, data.dof_index, data.quad_index, 1);
791
792 FaceIntegratorScalar energy_m(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
793 FaceIntegratorScalar energy_p(matrix_free, false, data.dof_index, data.quad_index, 1 + dim);
794
795 for(unsigned int face = face_range.first; face < face_range.second; face++)
796 {
797 // density
798 density_m.reinit(face);
799 density_m.gather_evaluate(src, dealii::EvaluationFlags::values);
800
801 density_p.reinit(face);
802 density_p.gather_evaluate(src, dealii::EvaluationFlags::values);
803
804 // velocity
805 momentum_m.reinit(face);
806 momentum_m.gather_evaluate(src, dealii::EvaluationFlags::values);
807
808 momentum_p.reinit(face);
809 momentum_p.gather_evaluate(src, dealii::EvaluationFlags::values);
810
811 // energy
812 energy_m.reinit(face);
813 energy_m.gather_evaluate(src, dealii::EvaluationFlags::values);
814
815 energy_p.reinit(face);
816 energy_p.gather_evaluate(src, dealii::EvaluationFlags::values);
817
818 for(unsigned int q = 0; q < momentum_m.n_q_points; ++q)
819 {
820 std::tuple<scalar, vector, scalar> flux =
821 get_flux(density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, q);
822
823 density_m.submit_value(std::get<0>(flux), q);
824 // - sign since n⁺ = -n⁻
825 density_p.submit_value(-std::get<0>(flux), q);
826
827 momentum_m.submit_value(std::get<1>(flux), q);
828 // - sign since n⁺ = -n⁻
829 momentum_p.submit_value(-std::get<1>(flux), q);
830
831 energy_m.submit_value(std::get<2>(flux), q);
832 // - sign since n⁺ = -n⁻
833 energy_p.submit_value(-std::get<2>(flux), q);
834 }
835
836 density_m.integrate_scatter(dealii::EvaluationFlags::values, dst);
837 density_p.integrate_scatter(dealii::EvaluationFlags::values, dst);
838
839 momentum_m.integrate_scatter(dealii::EvaluationFlags::values, dst);
840 momentum_p.integrate_scatter(dealii::EvaluationFlags::values, dst);
841
842 energy_m.integrate_scatter(dealii::EvaluationFlags::values, dst);
843 energy_p.integrate_scatter(dealii::EvaluationFlags::values, dst);
844 }
845 }
846
847 void
848 boundary_face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
849 VectorType & dst,
850 VectorType const & src,
851 std::pair<unsigned int, unsigned int> const & face_range) const
852 {
853 FaceIntegratorScalar density(matrix_free, true, data.dof_index, data.quad_index, 0);
854 FaceIntegratorVector momentum(matrix_free, true, data.dof_index, data.quad_index, 1);
855 FaceIntegratorScalar energy(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
856
857 for(unsigned int face = face_range.first; face < face_range.second; face++)
858 {
859 density.reinit(face);
860 density.gather_evaluate(src, dealii::EvaluationFlags::values);
861
862 momentum.reinit(face);
863 momentum.gather_evaluate(src, dealii::EvaluationFlags::values);
864
865 energy.reinit(face);
866 energy.gather_evaluate(src, dealii::EvaluationFlags::values);
867
868 dealii::types::boundary_id boundary_id = matrix_free.get_boundary_id(face);
869
870 BoundaryType boundary_type_density = data.bc->density.get_boundary_type(boundary_id);
871 BoundaryType boundary_type_velocity = data.bc->velocity.get_boundary_type(boundary_id);
872 BoundaryType boundary_type_pressure = data.bc->pressure.get_boundary_type(boundary_id);
873 BoundaryType boundary_type_energy = data.bc->energy.get_boundary_type(boundary_id);
874
875 EnergyBoundaryVariable boundary_variable = data.bc->energy.get_boundary_variable(boundary_id);
876
877 for(unsigned int q = 0; q < density.n_q_points; ++q)
878 {
879 std::tuple<scalar, vector, scalar> flux = get_flux_boundary(density,
880 momentum,
881 energy,
882 boundary_type_density,
883 boundary_type_velocity,
884 boundary_type_pressure,
885 boundary_type_energy,
886 boundary_variable,
887 boundary_id,
888 q);
889
890 density.submit_value(std::get<0>(flux), q);
891 momentum.submit_value(std::get<1>(flux), q);
892 energy.submit_value(std::get<2>(flux), q);
893 }
894
895 density.integrate_scatter(dealii::EvaluationFlags::values, dst);
896 momentum.integrate_scatter(dealii::EvaluationFlags::values, dst);
897 energy.integrate_scatter(dealii::EvaluationFlags::values, dst);
898 }
899 }
900
901 dealii::MatrixFree<dim, Number> const * matrix_free;
902
903 ConvectiveOperatorData<dim> data;
904
905 // heat capacity ratio
906 Number gamma;
907
908 // specific gas constant
909 Number R;
910
911 // specific heat at constant volume
912 Number c_v;
913
914 mutable Number eval_time;
915};
916
917
918template<int dim>
919struct ViscousOperatorData
920{
921 ViscousOperatorData()
922 : dof_index(0),
923 quad_index(0),
924 IP_factor(1.0),
925 dynamic_viscosity(1.0),
926 reference_density(1.0),
927 thermal_conductivity(0.0262),
928 heat_capacity_ratio(1.4),
929 specific_gas_constant(287.058)
930 {
931 }
932
933 unsigned int dof_index;
934 unsigned int quad_index;
935
936 double IP_factor;
937
938 std::shared_ptr<BoundaryDescriptor<dim> const> bc;
939
940 double dynamic_viscosity;
941 double reference_density;
942 double thermal_conductivity;
943 double heat_capacity_ratio;
944 double specific_gas_constant;
945};
946
947template<int dim, typename Number>
948class ViscousOperator
949{
950public:
951 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
952
953 typedef ViscousOperator<dim, Number> This;
954
955 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
956 typedef FaceIntegrator<dim, 1, Number> FaceIntegratorScalar;
957 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
958 typedef FaceIntegrator<dim, dim, Number> FaceIntegratorVector;
959
960 typedef dealii::VectorizedArray<Number> scalar;
961 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
962 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
963 typedef dealii::Point<dim, dealii::VectorizedArray<Number>> point;
964
965 ViscousOperator() : matrix_free(nullptr), degree(1)
966 {
967 }
968
969 void
970 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
971 ViscousOperatorData<dim> const & data_in)
972 {
973 this->matrix_free = &matrix_free_in;
974 this->data = data_in;
975
976 dealii::FiniteElement<dim> const & fe = matrix_free->get_dof_handler(data.dof_index).get_fe();
977 degree = fe.degree;
978
979 gamma = data.heat_capacity_ratio;
980 R = data.specific_gas_constant;
981 c_v = R / (gamma - 1.0);
982 mu = data.dynamic_viscosity;
983 nu = mu / data.reference_density;
984 lambda = data.thermal_conductivity;
985
986 IP::calculate_penalty_parameter<dim, Number>(array_penalty_parameter,
987 *matrix_free,
988 data.dof_index);
989 }
990
991 void
992 evaluate(VectorType & dst, VectorType const & src, Number const evaluation_time) const
993 {
994 dst = 0;
995 evaluate_add(dst, src, evaluation_time);
996 }
997
998 void
999 evaluate_add(VectorType & dst, VectorType const & src, Number const evaluation_time) const
1000 {
1001 this->eval_time = evaluation_time;
1002
1003 matrix_free->loop(
1004 &This::cell_loop, &This::face_loop, &This::boundary_face_loop, this, dst, src);
1005 }
1006
1007 void
1008 set_evaluation_time(double const & evaluation_time) const
1009 {
1010 eval_time = evaluation_time;
1011 }
1012
1013 inline DEAL_II_ALWAYS_INLINE //
1014 scalar
1015 get_penalty_parameter(FaceIntegratorScalar & fe_eval_m, FaceIntegratorScalar & fe_eval_p) const
1016 {
1017 scalar tau =
1018 std::max(fe_eval_m.read_cell_data(array_penalty_parameter),
1019 fe_eval_p.read_cell_data(array_penalty_parameter)) *
1020 IP::get_penalty_factor<dim, Number>(
1021 degree,
1022 get_element_type(
1023 fe_eval_m.get_matrix_free().get_dof_handler(data.dof_index).get_triangulation()),
1024 data.IP_factor) *
1025 nu;
1026
1027 return tau;
1028 }
1029
1030 inline DEAL_II_ALWAYS_INLINE //
1031 scalar
1032 get_penalty_parameter(FaceIntegratorScalar & fe_eval) const
1033 {
1034 scalar tau = fe_eval.read_cell_data(array_penalty_parameter) *
1035 IP::get_penalty_factor<dim, Number>(
1036 degree,
1037 get_element_type(
1038 fe_eval.get_matrix_free().get_dof_handler(data.dof_index).get_triangulation()),
1039 data.IP_factor) *
1040 nu;
1041
1042 return tau;
1043 }
1044
1045 inline DEAL_II_ALWAYS_INLINE //
1046 std::tuple<vector, tensor, vector>
1047 get_volume_flux(CellIntegratorScalar & density,
1048 CellIntegratorVector & momentum,
1049 CellIntegratorScalar & energy,
1050 unsigned int const q) const
1051 {
1052 scalar rho_inv = 1.0 / density.get_value(q);
1053 vector grad_rho = density.get_gradient(q);
1054
1055 vector rho_u = momentum.get_value(q);
1056 vector u = rho_inv * rho_u;
1057 tensor grad_rho_u = momentum.get_gradient(q);
1058
1059 scalar rho_E = energy.get_value(q);
1060 vector grad_rho_E = energy.get_gradient(q);
1061
1062 // calculate flux momentum
1063 tensor grad_u = calculate_grad_u(rho_inv, rho_u, grad_rho, grad_rho_u);
1064 tensor tau = calculate_stress_tensor(grad_u, mu);
1065
1066 // calculate flux energy
1067 vector grad_E = calculate_grad_E(rho_inv, rho_E, grad_rho, grad_rho_E);
1068 vector grad_T = calculate_grad_T(grad_E, u, grad_u, gamma, R);
1069 vector energy_flux = tau * u + lambda * grad_T;
1070
1071 return std::make_tuple(vector() /* dummy */, tau, energy_flux);
1072 }
1073
1074 inline DEAL_II_ALWAYS_INLINE //
1075 std::tuple<scalar, vector, scalar>
1076 get_gradient_flux(FaceIntegratorScalar & density_m,
1077 FaceIntegratorScalar & density_p,
1078 FaceIntegratorVector & momentum_m,
1079 FaceIntegratorVector & momentum_p,
1080 FaceIntegratorScalar & energy_m,
1081 FaceIntegratorScalar & energy_p,
1082 scalar const & tau_IP,
1083 unsigned int const q) const
1084 {
1085 vector normal = momentum_m.get_normal_vector(q);
1086
1087 // density
1088 scalar rho_M = density_m.get_value(q);
1089 vector grad_rho_M = density_m.get_gradient(q);
1090
1091 scalar rho_P = density_p.get_value(q);
1092 vector grad_rho_P = density_p.get_gradient(q);
1093
1094 // velocity
1095 vector rho_u_M = momentum_m.get_value(q);
1096 tensor grad_rho_u_M = momentum_m.get_gradient(q);
1097
1098 vector rho_u_P = momentum_p.get_value(q);
1099 tensor grad_rho_u_P = momentum_p.get_gradient(q);
1100
1101 // energy
1102 scalar rho_E_M = energy_m.get_value(q);
1103 vector grad_rho_E_M = energy_m.get_gradient(q);
1104
1105 scalar rho_E_P = energy_p.get_value(q);
1106 vector grad_rho_E_P = energy_p.get_gradient(q);
1107
1108 // flux density
1109 scalar jump_density = rho_M - rho_P;
1110 scalar gradient_flux_density = -tau_IP * jump_density;
1111
1112 // flux momentum
1113 scalar rho_inv_M = 1.0 / rho_M;
1114 tensor grad_u_M = calculate_grad_u(rho_inv_M, rho_u_M, grad_rho_M, grad_rho_u_M);
1115 tensor tau_M = calculate_stress_tensor(grad_u_M, mu);
1116
1117 scalar rho_inv_P = 1.0 / rho_P;
1118 tensor grad_u_P = calculate_grad_u(rho_inv_P, rho_u_P, grad_rho_P, grad_rho_u_P);
1119 tensor tau_P = calculate_stress_tensor(grad_u_P, mu);
1120
1121 vector jump_momentum = rho_u_M - rho_u_P;
1122 vector gradient_flux_momentum = 0.5 * (tau_M + tau_P) * normal - tau_IP * jump_momentum;
1123
1124 // flux energy
1125 vector u_M = rho_inv_M * rho_u_M;
1126 vector grad_E_M = calculate_grad_E(rho_inv_M, rho_E_M, grad_rho_M, grad_rho_E_M);
1127 vector grad_T_M = calculate_grad_T(grad_E_M, u_M, grad_u_M, gamma, R);
1128
1129 vector u_P = rho_inv_P * rho_u_P;
1130 vector grad_E_P = calculate_grad_E(rho_inv_P, rho_E_P, grad_rho_P, grad_rho_E_P);
1131 vector grad_T_P = calculate_grad_T(grad_E_P, u_P, grad_u_P, gamma, R);
1132
1133 vector flux_energy_average = 0.5 * (tau_M * u_M + tau_P * u_P + lambda * (grad_T_M + grad_T_P));
1134
1135 scalar jump_energy = rho_E_M - rho_E_P;
1136 scalar gradient_flux_energy = flux_energy_average * normal - tau_IP * jump_energy;
1137
1138 return std::make_tuple(gradient_flux_density, gradient_flux_momentum, gradient_flux_energy);
1139 }
1140
1141
1142 inline DEAL_II_ALWAYS_INLINE //
1143 std::tuple<scalar, vector, scalar>
1144 get_gradient_flux_boundary(FaceIntegratorScalar & density,
1145 FaceIntegratorVector & momentum,
1146 FaceIntegratorScalar & energy,
1147 scalar const & tau_IP,
1148 BoundaryType const & boundary_type_density,
1149 BoundaryType const & boundary_type_velocity,
1150 BoundaryType const & boundary_type_energy,
1151 EnergyBoundaryVariable const & boundary_variable,
1152 dealii::types::boundary_id const & boundary_id,
1153 unsigned int const q) const
1154 {
1155 vector normal = momentum.get_normal_vector(q);
1156
1157 // density
1158 scalar rho_M = density.get_value(q);
1159 vector grad_rho_M = density.get_gradient(q);
1160
1161 scalar rho_P = calculate_exterior_value<dim, Number, 0>(rho_M,
1162 boundary_type_density,
1163 data.bc->density,
1164 boundary_id,
1165 density.quadrature_point(q),
1166 this->eval_time);
1167
1168 scalar jump_density = rho_M - rho_P;
1169 scalar gradient_flux_density = -tau_IP * jump_density;
1170
1171 // velocity
1172 vector rho_u_M = momentum.get_value(q);
1173 tensor grad_rho_u_M = momentum.get_gradient(q);
1174
1175 scalar rho_inv_M = 1.0 / rho_M;
1176 vector u_M = rho_inv_M * rho_u_M;
1177
1178 vector u_P = calculate_exterior_value<dim, Number, 1>(u_M,
1179 boundary_type_velocity,
1180 data.bc->velocity,
1181 boundary_id,
1182 momentum.quadrature_point(q),
1183 this->eval_time);
1184
1185 vector rho_u_P = rho_P * u_P;
1186
1187 tensor grad_u_M = calculate_grad_u(rho_inv_M, rho_u_M, grad_rho_M, grad_rho_u_M);
1188 tensor tau_M = calculate_stress_tensor(grad_u_M, mu);
1189
1190 vector tau_P_normal =
1191 calculate_exterior_normal_grad<dim, Number, 1>(tau_M * normal,
1192 boundary_type_velocity,
1193 data.bc->velocity,
1194 boundary_id,
1195 momentum.quadrature_point(q),
1196 this->eval_time);
1197
1198 vector jump_momentum = rho_u_M - rho_u_P;
1199 vector gradient_flux_momentum = 0.5 * (tau_M * normal + tau_P_normal) - tau_IP * jump_momentum;
1200
1201 // energy
1202 scalar rho_E_M = energy.get_value(q);
1203 vector grad_rho_E_M = energy.get_gradient(q);
1204
1205 scalar E_M = rho_inv_M * rho_E_M;
1206 scalar E_P = dealii::make_vectorized_array<Number>(0.0);
1207 if(boundary_variable == EnergyBoundaryVariable::Energy)
1208 {
1209 E_P = calculate_exterior_value<dim, Number, 0>(E_M,
1210 boundary_type_energy,
1211 data.bc->energy,
1212 boundary_id,
1213 energy.quadrature_point(q),
1214 this->eval_time);
1215 }
1216 else if(boundary_variable == EnergyBoundaryVariable::Temperature)
1217 {
1218 scalar p_M = calculate_pressure(rho_M, u_M, E_M, gamma);
1219 scalar T_M = calculate_temperature(p_M, rho_M, R);
1220 scalar T_P = calculate_exterior_value<dim, Number, 0>(T_M,
1221 boundary_type_energy,
1222 data.bc->energy,
1223 boundary_id,
1224 energy.quadrature_point(q),
1225 this->eval_time);
1226
1227 E_P = calculate_energy(T_P, u_P, c_v);
1228 }
1229
1230 scalar rho_E_P = rho_P * E_P;
1231
1232 vector grad_E_M = calculate_grad_E(rho_inv_M, rho_E_M, grad_rho_M, grad_rho_E_M);
1233 vector grad_T_M = calculate_grad_T(grad_E_M, u_M, grad_u_M, gamma, R);
1234
1235 scalar grad_T_M_normal = grad_T_M * normal;
1236 scalar grad_T_P_normal =
1237 calculate_exterior_normal_grad<dim, Number, 0>(grad_T_M_normal,
1238 boundary_type_energy,
1239 data.bc->energy,
1240 boundary_id,
1241 energy.quadrature_point(q),
1242 this->eval_time);
1243
1244 scalar jump_energy = rho_E_M - rho_E_P;
1245 scalar gradient_flux_energy = 0.5 * (u_M * tau_M * normal + u_P * tau_P_normal +
1246 lambda * (grad_T_M * normal + grad_T_P_normal)) -
1247 tau_IP * jump_energy;
1248
1249 return std::make_tuple(gradient_flux_density, gradient_flux_momentum, gradient_flux_energy);
1250 }
1251
1252 inline DEAL_II_ALWAYS_INLINE //
1253 std::tuple<vector /*dummy_M*/,
1254 tensor /*value_flux_momentum_M*/,
1255 vector /*value_flux_energy_M*/,
1256 vector /*dummy_P*/,
1257 tensor /*value_flux_momentum_P*/,
1258 vector /*value_flux_energy_P*/>
1259 get_value_flux(FaceIntegratorScalar & density_m,
1260 FaceIntegratorScalar & density_p,
1261 FaceIntegratorVector & momentum_m,
1262 FaceIntegratorVector & momentum_p,
1263 FaceIntegratorScalar & energy_m,
1264 FaceIntegratorScalar & energy_p,
1265 unsigned int const q) const
1266 {
1267 vector normal = momentum_m.get_normal_vector(q);
1268
1269 // density
1270 scalar rho_M = density_m.get_value(q);
1271 scalar rho_P = density_p.get_value(q);
1272
1273 // velocity
1274 vector rho_u_M = momentum_m.get_value(q);
1275 vector rho_u_P = momentum_p.get_value(q);
1276
1277 // energy
1278 scalar rho_E_M = energy_m.get_value(q);
1279 scalar rho_E_P = energy_p.get_value(q);
1280
1281 vector jump_rho = (rho_M - rho_P) * normal;
1282 tensor jump_rho_u = outer_product(rho_u_M - rho_u_P, normal);
1283 vector jump_rho_E = (rho_E_M - rho_E_P) * normal;
1284
1285 scalar rho_inv_M = 1.0 / rho_M;
1286 scalar rho_inv_P = 1.0 / rho_P;
1287
1288 vector u_M = rho_inv_M * rho_u_M;
1289 vector u_P = rho_inv_P * rho_u_P;
1290
1291 // value flux momentum
1292 tensor grad_u_using_jumps_M = calculate_grad_u(rho_inv_M,
1293 rho_u_M,
1294 jump_rho /*instead of grad_rho*/,
1295 jump_rho_u /*instead of grad_rho_u*/);
1296
1297 tensor tau_using_jumps_M = calculate_stress_tensor(grad_u_using_jumps_M, mu);
1298 tensor value_flux_momentum_M = -0.5 * tau_using_jumps_M;
1299
1300 tensor grad_u_using_jumps_P = calculate_grad_u(rho_inv_P,
1301 rho_u_P,
1302 jump_rho /*instead of grad_rho*/,
1303 jump_rho_u /*instead of grad_rho_u*/);
1304
1305 tensor tau_using_jumps_P = calculate_stress_tensor(grad_u_using_jumps_P, mu);
1306 tensor value_flux_momentum_P = -0.5 * tau_using_jumps_P;
1307
1308 // value flux energy
1309 vector grad_E_using_jumps_M = calculate_grad_E(rho_inv_M,
1310 rho_E_M,
1311 jump_rho /*instead of grad_rho*/,
1312 jump_rho_E /*instead of grad_rho_E*/);
1313
1314 vector grad_T_using_jumps_M =
1315 calculate_grad_T(grad_E_using_jumps_M, u_M, grad_u_using_jumps_M, gamma, R);
1316 vector value_flux_energy_M = -0.5 * (tau_using_jumps_M * u_M + lambda * grad_T_using_jumps_M);
1317
1318 vector grad_E_using_jumps_P = calculate_grad_E(rho_inv_P,
1319 rho_E_P,
1320 jump_rho /*instead of grad_rho*/,
1321 jump_rho_E /*instead of grad_rho_E*/);
1322
1323 vector grad_T_using_jumps_P =
1324 calculate_grad_T(grad_E_using_jumps_P, u_P, grad_u_using_jumps_P, gamma, R);
1325 vector value_flux_energy_P = -0.5 * (tau_using_jumps_P * u_P + lambda * grad_T_using_jumps_P);
1326
1327 return std::make_tuple(vector() /*dummy*/,
1328 value_flux_momentum_M,
1329 value_flux_energy_M,
1330 vector() /*dummy*/,
1331 value_flux_momentum_P,
1332 value_flux_energy_P);
1333 }
1334
1335 inline DEAL_II_ALWAYS_INLINE //
1336 std::tuple<vector /*dummy_M*/, tensor /*value_flux_momentum_M*/, vector /*value_flux_energy_M*/>
1337 get_value_flux_boundary(FaceIntegratorScalar & density,
1338 FaceIntegratorVector & momentum,
1339 FaceIntegratorScalar & energy,
1340 BoundaryType const & boundary_type_density,
1341 BoundaryType const & boundary_type_velocity,
1342 BoundaryType const & boundary_type_energy,
1343 EnergyBoundaryVariable const & boundary_variable,
1344 dealii::types::boundary_id const & boundary_id,
1345 unsigned int const q) const
1346 {
1347 vector normal = momentum.get_normal_vector(q);
1348
1349 // density
1350 scalar rho_M = density.get_value(q);
1351 scalar rho_P = calculate_exterior_value<dim, Number, 0>(rho_M,
1352 boundary_type_density,
1353 data.bc->density,
1354 boundary_id,
1355 density.quadrature_point(q),
1356 this->eval_time);
1357
1358 scalar rho_inv_M = 1.0 / rho_M;
1359
1360 // velocity
1361 vector rho_u_M = momentum.get_value(q);
1362 vector u_M = rho_inv_M * rho_u_M;
1363
1364 vector u_P = calculate_exterior_value<dim, Number, 1>(u_M,
1365 boundary_type_velocity,
1366 data.bc->velocity,
1367 boundary_id,
1368 momentum.quadrature_point(q),
1369 this->eval_time);
1370
1371 vector rho_u_P = rho_P * u_P;
1372
1373 // energy
1374 scalar rho_E_M = energy.get_value(q);
1375 scalar E_M = rho_inv_M * rho_E_M;
1376
1377 scalar E_P = dealii::make_vectorized_array<Number>(0.0);
1378 if(boundary_variable == EnergyBoundaryVariable::Energy)
1379 {
1380 E_P = calculate_exterior_value<dim, Number, 0>(E_M,
1381 boundary_type_energy,
1382 data.bc->energy,
1383 boundary_id,
1384 energy.quadrature_point(q),
1385 this->eval_time);
1386 }
1387 else if(boundary_variable == EnergyBoundaryVariable::Temperature)
1388 {
1389 scalar p_M = calculate_pressure(rho_M, u_M, E_M, gamma);
1390 scalar T_M = calculate_temperature(p_M, rho_M, R);
1391 scalar T_P = calculate_exterior_value<dim, Number, 0>(T_M,
1392 boundary_type_energy,
1393 data.bc->energy,
1394 boundary_id,
1395 energy.quadrature_point(q),
1396 this->eval_time);
1397
1398 E_P = calculate_energy(T_P, u_P, c_v);
1399 }
1400 scalar rho_E_P = rho_P * E_P;
1401
1402 vector jump_rho = (rho_M - rho_P) * normal;
1403 tensor jump_rho_u = outer_product(rho_u_M - rho_u_P, normal);
1404 vector jump_rho_E = (rho_E_M - rho_E_P) * normal;
1405
1406 // value flux momentum
1407 tensor grad_u_using_jumps_M = calculate_grad_u(rho_inv_M,
1408 rho_u_M,
1409 jump_rho /*instead of grad_rho*/,
1410 jump_rho_u /*instead of grad_rho_u*/);
1411
1412 tensor tau_using_jumps_M = calculate_stress_tensor(grad_u_using_jumps_M, mu);
1413 tensor value_flux_momentum_M = -0.5 * tau_using_jumps_M;
1414
1415 // value flux energy
1416 vector grad_E_using_jumps_M = calculate_grad_E(rho_inv_M,
1417 rho_E_M,
1418 jump_rho /*instead of grad_rho*/,
1419 jump_rho_E /*instead of grad_rho_E*/);
1420
1421 vector grad_T_using_jumps_M =
1422 calculate_grad_T(grad_E_using_jumps_M, u_M, grad_u_using_jumps_M, gamma, R);
1423 vector value_flux_energy_M = -0.5 * (tau_using_jumps_M * u_M + lambda * grad_T_using_jumps_M);
1424
1425 return std::make_tuple(vector() /*dummy*/, value_flux_momentum_M, value_flux_energy_M);
1426 }
1427
1428private:
1429 void
1430 cell_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1431 VectorType & dst,
1432 VectorType const & src,
1433 std::pair<unsigned int, unsigned int> const & cell_range) const
1434 {
1435 CellIntegratorScalar density(matrix_free, data.dof_index, data.quad_index, 0);
1436 CellIntegratorVector momentum(matrix_free, data.dof_index, data.quad_index, 1);
1437 CellIntegratorScalar energy(matrix_free, data.dof_index, data.quad_index, 1 + dim);
1438
1439 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1440 {
1441 density.reinit(cell);
1442 density.gather_evaluate(src,
1443 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1444
1445 momentum.reinit(cell);
1446 momentum.gather_evaluate(src,
1447 dealii::EvaluationFlags::values |
1448 dealii::EvaluationFlags::gradients);
1449
1450 energy.reinit(cell);
1451 energy.gather_evaluate(src,
1452 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1453
1454 for(unsigned int q = 0; q < momentum.n_q_points; ++q)
1455 {
1456 std::tuple<vector, tensor, vector> flux = get_volume_flux(density, momentum, energy, q);
1457
1458 momentum.submit_gradient(std::get<1>(flux), q);
1459 energy.submit_gradient(std::get<2>(flux), q);
1460 }
1461
1462 momentum.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1463 energy.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1464 }
1465 }
1466
1467 void
1468 face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1469 VectorType & dst,
1470 VectorType const & src,
1471 std::pair<unsigned int, unsigned int> const & face_range) const
1472 {
1473 FaceIntegratorScalar density_m(matrix_free, true, data.dof_index, data.quad_index, 0);
1474 FaceIntegratorScalar density_p(matrix_free, false, data.dof_index, data.quad_index, 0);
1475
1476 FaceIntegratorVector momentum_m(matrix_free, true, data.dof_index, data.quad_index, 1);
1477 FaceIntegratorVector momentum_p(matrix_free, false, data.dof_index, data.quad_index, 1);
1478
1479 FaceIntegratorScalar energy_m(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
1480 FaceIntegratorScalar energy_p(matrix_free, false, data.dof_index, data.quad_index, 1 + dim);
1481
1482 for(unsigned int face = face_range.first; face < face_range.second; face++)
1483 {
1484 // density
1485 density_m.reinit(face);
1486 density_m.gather_evaluate(src,
1487 dealii::EvaluationFlags::values |
1488 dealii::EvaluationFlags::gradients);
1489
1490 density_p.reinit(face);
1491 density_p.gather_evaluate(src,
1492 dealii::EvaluationFlags::values |
1493 dealii::EvaluationFlags::gradients);
1494
1495 // momentum
1496 momentum_m.reinit(face);
1497 momentum_m.gather_evaluate(src,
1498 dealii::EvaluationFlags::values |
1499 dealii::EvaluationFlags::gradients);
1500
1501 momentum_p.reinit(face);
1502 momentum_p.gather_evaluate(src,
1503 dealii::EvaluationFlags::values |
1504 dealii::EvaluationFlags::gradients);
1505
1506 // energy
1507 energy_m.reinit(face);
1508 energy_m.gather_evaluate(src,
1509 dealii::EvaluationFlags::values |
1510 dealii::EvaluationFlags::gradients);
1511
1512 energy_p.reinit(face);
1513 energy_p.gather_evaluate(src,
1514 dealii::EvaluationFlags::values |
1515 dealii::EvaluationFlags::gradients);
1516
1517 scalar tau_IP = get_penalty_parameter(density_m, density_p);
1518
1519 for(unsigned int q = 0; q < density_m.n_q_points; ++q)
1520 {
1521 std::tuple<scalar, vector, scalar> gradient_flux = get_gradient_flux(
1522 density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, tau_IP, q);
1523
1524 std::tuple<vector, tensor, vector, vector, tensor, vector> value_flux =
1525 get_value_flux(density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, q);
1526
1527 density_m.submit_value(-std::get<0>(gradient_flux), q);
1528 // + sign since n⁺ = -n⁻
1529 density_p.submit_value(std::get<0>(gradient_flux), q);
1530
1531 momentum_m.submit_gradient(std::get<1>(value_flux), q);
1532 // note that value_flux_momentum is not conservative
1533 momentum_p.submit_gradient(std::get<4>(value_flux), q);
1534
1535 momentum_m.submit_value(-std::get<1>(gradient_flux), q);
1536 // + sign since n⁺ = -n⁻
1537 momentum_p.submit_value(std::get<1>(gradient_flux), q);
1538
1539 energy_m.submit_gradient(std::get<2>(value_flux), q);
1540 // note that value_flux_energy is not conservative
1541 energy_p.submit_gradient(std::get<5>(value_flux), q);
1542
1543 energy_m.submit_value(-std::get<2>(gradient_flux), q);
1544 // + sign since n⁺ = -n⁻
1545 energy_p.submit_value(std::get<2>(gradient_flux), q);
1546 }
1547
1548 density_m.integrate_scatter(dealii::EvaluationFlags::values, dst);
1549 density_p.integrate_scatter(dealii::EvaluationFlags::values, dst);
1550
1551 momentum_m.integrate_scatter(dealii::EvaluationFlags::values |
1552 dealii::EvaluationFlags::gradients,
1553 dst);
1554 momentum_p.integrate_scatter(dealii::EvaluationFlags::values |
1555 dealii::EvaluationFlags::gradients,
1556 dst);
1557
1558 energy_m.integrate_scatter(dealii::EvaluationFlags::values |
1559 dealii::EvaluationFlags::gradients,
1560 dst);
1561 energy_p.integrate_scatter(dealii::EvaluationFlags::values |
1562 dealii::EvaluationFlags::gradients,
1563 dst);
1564 }
1565 }
1566
1567 void
1568 boundary_face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1569 VectorType & dst,
1570 VectorType const & src,
1571 std::pair<unsigned int, unsigned int> const & face_range) const
1572 {
1573 FaceIntegratorScalar density(matrix_free, true, data.dof_index, data.quad_index, 0);
1574 FaceIntegratorVector momentum(matrix_free, true, data.dof_index, data.quad_index, 1);
1575 FaceIntegratorScalar energy(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
1576
1577 for(unsigned int face = face_range.first; face < face_range.second; face++)
1578 {
1579 density.reinit(face);
1580 density.gather_evaluate(src,
1581 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1582
1583 momentum.reinit(face);
1584 momentum.gather_evaluate(src,
1585 dealii::EvaluationFlags::values |
1586 dealii::EvaluationFlags::gradients);
1587
1588 energy.reinit(face);
1589 energy.gather_evaluate(src,
1590 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1591
1592 scalar tau_IP = get_penalty_parameter(density);
1593
1594 dealii::types::boundary_id boundary_id = matrix_free.get_boundary_id(face);
1595
1596 BoundaryType boundary_type_density = data.bc->density.get_boundary_type(boundary_id);
1597 BoundaryType boundary_type_velocity = data.bc->velocity.get_boundary_type(boundary_id);
1598 BoundaryType boundary_type_energy = data.bc->energy.get_boundary_type(boundary_id);
1599
1600 EnergyBoundaryVariable boundary_variable = data.bc->energy.get_boundary_variable(boundary_id);
1601
1602 for(unsigned int q = 0; q < density.n_q_points; ++q)
1603 {
1604 std::tuple<scalar, vector, scalar> gradient_flux =
1605 get_gradient_flux_boundary(density,
1606 momentum,
1607 energy,
1608 tau_IP,
1609 boundary_type_density,
1610 boundary_type_velocity,
1611 boundary_type_energy,
1612 boundary_variable,
1613 boundary_id,
1614 q);
1615
1616 std::tuple<vector, tensor, vector> value_flux =
1617 get_value_flux_boundary(density,
1618 momentum,
1619 energy,
1620 boundary_type_density,
1621 boundary_type_velocity,
1622 boundary_type_energy,
1623 boundary_variable,
1624 boundary_id,
1625 q);
1626
1627 density.submit_value(-std::get<0>(gradient_flux), q);
1628
1629 momentum.submit_gradient(std::get<1>(value_flux), q);
1630 momentum.submit_value(-std::get<1>(gradient_flux), q);
1631
1632 energy.submit_gradient(std::get<2>(value_flux), q);
1633 energy.submit_value(-std::get<2>(gradient_flux), q);
1634 }
1635
1636 density.integrate_scatter(dealii::EvaluationFlags::values, dst);
1637 momentum.integrate_scatter(dealii::EvaluationFlags::values |
1638 dealii::EvaluationFlags::gradients,
1639 dst);
1640 energy.integrate_scatter(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients,
1641 dst);
1642 }
1643 }
1644
1645 dealii::MatrixFree<dim, Number> const * matrix_free;
1646
1647 ViscousOperatorData<dim> data;
1648
1649 unsigned int degree;
1650
1651 // heat capacity ratio
1652 Number gamma;
1653
1654 // specific gas constant
1655 Number R;
1656
1657 // specific heat at constant volume
1658 Number c_v;
1659
1660 // dynamic viscosity
1661 Number mu;
1662
1663 // kinematic viscosity
1664 Number nu;
1665
1666 // thermal conductivity
1667 Number lambda;
1668
1669 dealii::AlignedVector<dealii::VectorizedArray<Number>> array_penalty_parameter;
1670
1671 mutable Number eval_time;
1672};
1673
1674template<int dim>
1675struct CombinedOperatorData
1676{
1677 CombinedOperatorData() : dof_index(0), quad_index(0)
1678 {
1679 }
1680
1681 unsigned int dof_index;
1682 unsigned int quad_index;
1683
1684 std::shared_ptr<BoundaryDescriptor<dim> const> bc;
1685};
1686
1687template<int dim, typename Number>
1688class CombinedOperator
1689{
1690public:
1691 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
1692
1693 typedef ConvectiveOperator<dim, Number> ConvectiveOp;
1694 typedef ViscousOperator<dim, Number> ViscousOp;
1695 typedef CombinedOperator<dim, Number> This;
1696
1697 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
1698 typedef FaceIntegrator<dim, 1, Number> FaceIntegratorScalar;
1699 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
1700 typedef FaceIntegrator<dim, dim, Number> FaceIntegratorVector;
1701
1702 typedef dealii::VectorizedArray<Number> scalar;
1703 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
1704 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
1705 typedef dealii::Point<dim, dealii::VectorizedArray<Number>> point;
1706
1707 CombinedOperator() : matrix_free(nullptr), convective_operator(nullptr), viscous_operator(nullptr)
1708 {
1709 }
1710
1711 void
1712 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
1713 CombinedOperatorData<dim> const & data_in,
1714 ConvectiveOp const & convective_operator_in,
1715 ViscousOp const & viscous_operator_in)
1716 {
1717 this->matrix_free = &matrix_free_in;
1718 this->data = data_in;
1719
1720 this->convective_operator = &convective_operator_in;
1721 this->viscous_operator = &viscous_operator_in;
1722 }
1723
1724 void
1725 evaluate(VectorType & dst, VectorType const & src, Number const evaluation_time) const
1726 {
1727 dst = 0;
1728 evaluate_add(dst, src, evaluation_time);
1729 }
1730
1731 void
1732 evaluate_add(VectorType & dst, VectorType const & src, Number const evaluation_time) const
1733 {
1734 convective_operator->set_evaluation_time(evaluation_time);
1735 viscous_operator->set_evaluation_time(evaluation_time);
1736
1737 matrix_free->loop(
1738 &This::cell_loop, &This::face_loop, &This::boundary_face_loop, this, dst, src);
1739
1740 // perform cell integrals only for performance measurements
1741 // matrix_free->cell_loop(&This::cell_loop, this, dst, src);
1742 }
1743
1744private:
1745 void
1746 cell_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1747 VectorType & dst,
1748 VectorType const & src,
1749 std::pair<unsigned int, unsigned int> const & cell_range) const
1750 {
1751 CellIntegratorScalar density(matrix_free, data.dof_index, data.quad_index, 0);
1752 CellIntegratorVector momentum(matrix_free, data.dof_index, data.quad_index, 1);
1753 CellIntegratorScalar energy(matrix_free, data.dof_index, data.quad_index, 1 + dim);
1754
1755 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1756 {
1757 density.reinit(cell);
1758 density.gather_evaluate(src,
1759 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1760
1761 momentum.reinit(cell);
1762 momentum.gather_evaluate(src,
1763 dealii::EvaluationFlags::values |
1764 dealii::EvaluationFlags::gradients);
1765
1766 energy.reinit(cell);
1767 energy.gather_evaluate(src,
1768 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1769
1770 for(unsigned int q = 0; q < momentum.n_q_points; ++q)
1771 {
1772 std::tuple<vector, tensor, vector> conv_flux =
1773 convective_operator->get_volume_flux(density, momentum, energy, q);
1774
1775 std::tuple<vector, tensor, vector> visc_flux =
1776 viscous_operator->get_volume_flux(density, momentum, energy, q);
1777
1778 density.submit_gradient(-std::get<0>(conv_flux), q);
1779 momentum.submit_gradient(-std::get<1>(conv_flux) + std::get<1>(visc_flux), q);
1780 energy.submit_gradient(-std::get<2>(conv_flux) + std::get<2>(visc_flux), q);
1781 }
1782
1783 density.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1784 momentum.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1785 energy.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1786 }
1787 }
1788
1789 void
1790 face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1791 VectorType & dst,
1792 VectorType const & src,
1793 std::pair<unsigned int, unsigned int> const & face_range) const
1794 {
1795 FaceIntegratorScalar density_m(matrix_free, true, data.dof_index, data.quad_index, 0);
1796 FaceIntegratorScalar density_p(matrix_free, false, data.dof_index, data.quad_index, 0);
1797 FaceIntegratorVector momentum_m(matrix_free, true, data.dof_index, data.quad_index, 1);
1798 FaceIntegratorVector momentum_p(matrix_free, false, data.dof_index, data.quad_index, 1);
1799 FaceIntegratorScalar energy_m(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
1800 FaceIntegratorScalar energy_p(matrix_free, false, data.dof_index, data.quad_index, 1 + dim);
1801
1802 for(unsigned int face = face_range.first; face < face_range.second; face++)
1803 {
1804 // density
1805 density_m.reinit(face);
1806 density_m.gather_evaluate(src,
1807 dealii::EvaluationFlags::values |
1808 dealii::EvaluationFlags::gradients);
1809
1810 density_p.reinit(face);
1811 density_p.gather_evaluate(src,
1812 dealii::EvaluationFlags::values |
1813 dealii::EvaluationFlags::gradients);
1814
1815 // momentum
1816 momentum_m.reinit(face);
1817 momentum_m.gather_evaluate(src,
1818 dealii::EvaluationFlags::values |
1819 dealii::EvaluationFlags::gradients);
1820
1821 momentum_p.reinit(face);
1822 momentum_p.gather_evaluate(src,
1823 dealii::EvaluationFlags::values |
1824 dealii::EvaluationFlags::gradients);
1825
1826 // energy
1827 energy_m.reinit(face);
1828 energy_m.gather_evaluate(src,
1829 dealii::EvaluationFlags::values |
1830 dealii::EvaluationFlags::gradients);
1831
1832 energy_p.reinit(face);
1833 energy_p.gather_evaluate(src,
1834 dealii::EvaluationFlags::values |
1835 dealii::EvaluationFlags::gradients);
1836
1837 scalar tau_IP = viscous_operator->get_penalty_parameter(density_m, density_p);
1838
1839 for(unsigned int q = 0; q < density_m.n_q_points; ++q)
1840 {
1841 std::tuple<scalar, vector, scalar> conv_flux = convective_operator->get_flux(
1842 density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, q);
1843
1844 std::tuple<scalar, vector, scalar> visc_grad_flux = viscous_operator->get_gradient_flux(
1845 density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, tau_IP, q);
1846
1847 std::tuple<vector, tensor, vector, vector, tensor, vector> visc_value_flux =
1848 viscous_operator->get_value_flux(
1849 density_m, density_p, momentum_m, momentum_p, energy_m, energy_p, q);
1850
1851 density_m.submit_value(std::get<0>(conv_flux) - std::get<0>(visc_grad_flux), q);
1852 // - sign since n⁺ = -n⁻
1853 density_p.submit_value(-std::get<0>(conv_flux) + std::get<0>(visc_grad_flux), q);
1854
1855 momentum_m.submit_value(std::get<1>(conv_flux) - std::get<1>(visc_grad_flux), q);
1856 // - sign since n⁺ = -n⁻
1857 momentum_p.submit_value(-std::get<1>(conv_flux) + std::get<1>(visc_grad_flux), q);
1858
1859 momentum_m.submit_gradient(std::get<1>(visc_value_flux), q);
1860 // note that value_flux_momentum is not conservative
1861 momentum_p.submit_gradient(std::get<4>(visc_value_flux), q);
1862
1863 energy_m.submit_value(std::get<2>(conv_flux) - std::get<2>(visc_grad_flux), q);
1864 // - sign since n⁺ = -n⁻
1865 energy_p.submit_value(-std::get<2>(conv_flux) + std::get<2>(visc_grad_flux), q);
1866
1867 energy_m.submit_gradient(std::get<2>(visc_value_flux), q);
1868 // note that value_flux_energy is not conservative
1869 energy_p.submit_gradient(std::get<5>(visc_value_flux), q);
1870 }
1871
1872 density_m.integrate_scatter(dealii::EvaluationFlags::values, dst);
1873 density_p.integrate_scatter(dealii::EvaluationFlags::values, dst);
1874
1875 momentum_m.integrate_scatter(dealii::EvaluationFlags::values |
1876 dealii::EvaluationFlags::gradients,
1877 dst);
1878 momentum_p.integrate_scatter(dealii::EvaluationFlags::values |
1879 dealii::EvaluationFlags::gradients,
1880 dst);
1881
1882 energy_m.integrate_scatter(dealii::EvaluationFlags::values |
1883 dealii::EvaluationFlags::gradients,
1884 dst);
1885 energy_p.integrate_scatter(dealii::EvaluationFlags::values |
1886 dealii::EvaluationFlags::gradients,
1887 dst);
1888 }
1889 }
1890
1891 void
1892 boundary_face_loop(dealii::MatrixFree<dim, Number> const & matrix_free,
1893 VectorType & dst,
1894 VectorType const & src,
1895 std::pair<unsigned int, unsigned int> const & face_range) const
1896 {
1897 FaceIntegratorScalar density(matrix_free, true, data.dof_index, data.quad_index, 0);
1898 FaceIntegratorVector momentum(matrix_free, true, data.dof_index, data.quad_index, 1);
1899 FaceIntegratorScalar energy(matrix_free, true, data.dof_index, data.quad_index, 1 + dim);
1900
1901 for(unsigned int face = face_range.first; face < face_range.second; face++)
1902 {
1903 density.reinit(face);
1904 density.gather_evaluate(src,
1905 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1906
1907 momentum.reinit(face);
1908 momentum.gather_evaluate(src,
1909 dealii::EvaluationFlags::values |
1910 dealii::EvaluationFlags::gradients);
1911
1912 energy.reinit(face);
1913 energy.gather_evaluate(src,
1914 dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients);
1915
1916 scalar tau_IP = viscous_operator->get_penalty_parameter(density);
1917
1918 dealii::types::boundary_id boundary_id = matrix_free.get_boundary_id(face);
1919
1920 BoundaryType boundary_type_density = data.bc->density.get_boundary_type(boundary_id);
1921 BoundaryType boundary_type_velocity = data.bc->velocity.get_boundary_type(boundary_id);
1922 BoundaryType boundary_type_pressure = data.bc->pressure.get_boundary_type(boundary_id);
1923 BoundaryType boundary_type_energy = data.bc->energy.get_boundary_type(boundary_id);
1924
1925 EnergyBoundaryVariable boundary_variable = data.bc->energy.get_boundary_variable(boundary_id);
1926
1927 for(unsigned int q = 0; q < density.n_q_points; ++q)
1928 {
1929 std::tuple<scalar, vector, scalar> conv_flux =
1930 convective_operator->get_flux_boundary(density,
1931 momentum,
1932 energy,
1933 boundary_type_density,
1934 boundary_type_velocity,
1935 boundary_type_pressure,
1936 boundary_type_energy,
1937 boundary_variable,
1938 boundary_id,
1939 q);
1940
1941 std::tuple<scalar, vector, scalar> visc_grad_flux =
1942 viscous_operator->get_gradient_flux_boundary(density,
1943 momentum,
1944 energy,
1945 tau_IP,
1946 boundary_type_density,
1947 boundary_type_velocity,
1948 boundary_type_energy,
1949 boundary_variable,
1950 boundary_id,
1951 q);
1952
1953 std::tuple<vector, tensor, vector> visc_value_flux =
1954 viscous_operator->get_value_flux_boundary(density,
1955 momentum,
1956 energy,
1957 boundary_type_density,
1958 boundary_type_velocity,
1959 boundary_type_energy,
1960 boundary_variable,
1961 boundary_id,
1962 q);
1963
1964 density.submit_value(std::get<0>(conv_flux) - std::get<0>(visc_grad_flux), q);
1965
1966 momentum.submit_value(std::get<1>(conv_flux) - std::get<1>(visc_grad_flux), q);
1967 momentum.submit_gradient(std::get<1>(visc_value_flux), q);
1968
1969 energy.submit_value(std::get<2>(conv_flux) - std::get<2>(visc_grad_flux), q);
1970 energy.submit_gradient(std::get<2>(visc_value_flux), q);
1971 }
1972
1973 density.integrate_scatter(dealii::EvaluationFlags::values, dst);
1974 momentum.integrate_scatter(dealii::EvaluationFlags::values |
1975 dealii::EvaluationFlags::gradients,
1976 dst);
1977 energy.integrate_scatter(dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients,
1978 dst);
1979 }
1980 }
1981
1982 dealii::MatrixFree<dim, Number> const * matrix_free;
1983
1984 CombinedOperatorData<dim> data;
1985
1986 ConvectiveOperator<dim, Number> const * convective_operator;
1987 ViscousOperator<dim, Number> const * viscous_operator;
1988};
1989
1990} // namespace CompNS
1991} // namespace ExaDG
1992
1993#endif /* INCLUDE_EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_KERNELS_AND_OPERATORS_H_ \
1994 */
ExaDG::CompNS::ConvectiveOperator
Definition kernels_and_operators.h:510
ExaDG::CompNS::ViscousOperator
Definition kernels_and_operators.h:949
ExaDG
Definition driver.cpp:33
ExaDG::get_element_type
ElementType get_element_type(dealii::Triangulation< dim > const &tria)
Definition grid_data.h:70
ExaDG::CompNS::BodyForceOperatorData
Definition kernels_and_operators.h:286
ExaDG::CompNS::BoundaryDescriptorStd
Definition boundary_descriptor.h:48
ExaDG::CompNS::CombinedOperatorData
Definition kernels_and_operators.h:1676
ExaDG::CompNS::ConvectiveOperatorData
Definition kernels_and_operators.h:493
ExaDG::CompNS::MassOperatorData
Definition kernels_and_operators.h:406
ExaDG::CompNS::ViscousOperatorData
Definition kernels_and_operators.h:920
Generated by doxygen 1.13.2