ExaDG
Loading...
Searching...
No Matches
calculators.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 EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_CALCULATORS_H_
23#define EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_CALCULATORS_H_
24
25// deal.II
26#include <deal.II/lac/la_parallel_vector.h>
27#include <deal.II/matrix_free/matrix_free.h>
28
29// ExaDG
30#include <exadg/matrix_free/integrators.h>
31
32namespace ExaDG
33{
34namespace CompNS
35{
36/*
37 * The DoF vector contains the vector of conserved quantities (rho, rho u, rho E).
38 * This class allows to transfer these quantities into the derived variables (p, u, T)
39 * by using L2-projections.
40 */
41template<int dim, typename Number>
42class p_u_T_Calculator
43{
44public:
45 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
46
47 typedef p_u_T_Calculator<dim, Number> This;
48
49 typedef CellIntegrator<dim, 1, Number> CellIntegratorScalar;
50 typedef CellIntegrator<dim, dim, Number> CellIntegratorVector;
51
52 typedef dealii::VectorizedArray<Number> scalar;
53 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
54
55 p_u_T_Calculator()
56 : matrix_free(nullptr),
57 dof_index_all(0),
58 dof_index_vector(1),
59 dof_index_scalar(2),
60 quad_index(0),
61 heat_capacity_ratio(-1.0),
62 specific_gas_constant(-1.0)
63 {
64 }
65
66 void
67 initialize(dealii::MatrixFree<dim, Number> const & matrix_free_in,
68 unsigned int const dof_index_all_in,
69 unsigned int const dof_index_vector_in,
70 unsigned int const dof_index_scalar_in,
71 unsigned int const quad_index_in,
72 double const heat_capacity_ratio_in,
73 double const specific_gas_constant_in)
74 {
75 matrix_free = &matrix_free_in;
76 dof_index_all = dof_index_all_in;
77 dof_index_vector = dof_index_vector_in;
78 dof_index_scalar = dof_index_scalar_in;
79 quad_index = quad_index_in;
80 heat_capacity_ratio = heat_capacity_ratio_in;
81 specific_gas_constant = specific_gas_constant_in;
82 }
83
84 void
85 compute_pressure(VectorType & pressure, VectorType const & solution_conserved) const
86 {
87 AssertThrow(heat_capacity_ratio > 0.0,
88 dealii::ExcMessage("heat capacity ratio has not been set!"));
89 AssertThrow(specific_gas_constant > 0.0,
90 dealii::ExcMessage("specific gas constant has not been set!"));
91
92 matrix_free->cell_loop(&This::local_apply_pressure, this, pressure, solution_conserved);
93 }
94
95 void
96 compute_velocity(VectorType & velocity, VectorType const & solution_conserved) const
97 {
98 matrix_free->cell_loop(&This::local_apply_velocity, this, velocity, solution_conserved);
99 }
100
101 void
102 compute_temperature(VectorType & temperature, VectorType const & solution_conserved) const
103 {
104 AssertThrow(heat_capacity_ratio > 0.0,
105 dealii::ExcMessage("heat capacity ratio has not been set!"));
106 AssertThrow(specific_gas_constant > 0.0,
107 dealii::ExcMessage("specific gas constant has not been set!"));
108
109 matrix_free->cell_loop(&This::local_apply_temperature, this, temperature, solution_conserved);
110 }
111
112private:
113 void
114 local_apply_pressure(dealii::MatrixFree<dim, Number> const & matrix_free,
115 VectorType & dst,
116 VectorType const & src,
117 std::pair<unsigned int, unsigned int> const & cell_range) const
118 {
119 // src-vector
120 CellIntegratorScalar density(matrix_free, dof_index_all, quad_index, 0);
121 CellIntegratorVector momentum(matrix_free, dof_index_all, quad_index, 1);
122 CellIntegratorScalar energy(matrix_free, dof_index_all, quad_index, 1 + dim);
123
124 // dst-vector
125 CellIntegratorScalar pressure(matrix_free, dof_index_scalar, quad_index, 0);
126
127 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
128 {
129 // src-vector
130 density.reinit(cell);
131 density.read_dof_values(src);
132 density.evaluate(dealii::EvaluationFlags::values);
133 momentum.reinit(cell);
134 momentum.read_dof_values(src);
135 momentum.evaluate(dealii::EvaluationFlags::values);
136 energy.reinit(cell);
137 energy.read_dof_values(src);
138 energy.evaluate(dealii::EvaluationFlags::values);
139
140 // dst-vector
141 pressure.reinit(cell);
142
143 for(unsigned int q = 0; q < energy.n_q_points; ++q)
144 {
145 // conserved variables
146 scalar rho = density.get_value(q);
147 vector rho_u = momentum.get_value(q);
148 scalar rho_E = energy.get_value(q);
149
150 // compute derived quantities in quadrature point
151 vector u = rho_u / rho;
152
153 scalar p = dealii::make_vectorized_array<Number>(heat_capacity_ratio - 1.0) *
154 (rho_E - 0.5 * rho * scalar_product(u, u));
155
156 pressure.submit_value(p, q);
157 }
158
159 pressure.integrate_scatter(dealii::EvaluationFlags::values, dst);
160 }
161 }
162
163 void
164 local_apply_velocity(dealii::MatrixFree<dim, Number> const & matrix_free,
165 VectorType & dst,
166 VectorType const & src,
167 std::pair<unsigned int, unsigned int> const & cell_range) const
168 {
169 // src-vector
170 CellIntegratorScalar density(matrix_free, dof_index_all, quad_index, 0);
171 CellIntegratorVector momentum(matrix_free, dof_index_all, quad_index, 1);
172
173 // dst-vector
174 CellIntegratorVector velocity(matrix_free, dof_index_vector, quad_index, 0);
175
176 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
177 {
178 // src-vector
179 density.reinit(cell);
180 density.read_dof_values(src);
181 density.evaluate(dealii::EvaluationFlags::values);
182 momentum.reinit(cell);
183 momentum.read_dof_values(src);
184 momentum.evaluate(dealii::EvaluationFlags::values);
185
186 // dst-vector
187 velocity.reinit(cell);
188
189 for(unsigned int q = 0; q < momentum.n_q_points; ++q)
190 {
191 // conserved variables
192 scalar rho = density.get_value(q);
193 vector rho_u = momentum.get_value(q);
194
195 // compute derived quantities in quadrature point
196 vector u = rho_u / rho;
197
198 velocity.submit_value(u, q);
199 }
200
201 velocity.integrate_scatter(dealii::EvaluationFlags::values, dst);
202 }
203 }
204
205 void
206 local_apply_temperature(dealii::MatrixFree<dim, Number> const & matrix_free,
207 VectorType & dst,
208 VectorType const & src,
209 std::pair<unsigned int, unsigned int> const & cell_range) const
210 {
211 // src-vector
212 CellIntegratorScalar density(matrix_free, dof_index_all, quad_index, 0);
213 CellIntegratorVector momentum(matrix_free, dof_index_all, quad_index, 1);
214 CellIntegratorScalar energy(matrix_free, dof_index_all, quad_index, 1 + dim);
215
216 // dst-vector
217 CellIntegratorScalar temperature(matrix_free, dof_index_scalar, quad_index, 0);
218
219 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
220 {
221 // src-vector
222 density.reinit(cell);
223 density.read_dof_values(src);
224 density.evaluate(dealii::EvaluationFlags::values);
225 momentum.reinit(cell);
226 momentum.read_dof_values(src);
227 momentum.evaluate(dealii::EvaluationFlags::values);
228 energy.reinit(cell);
229 energy.read_dof_values(src);
230 energy.evaluate(dealii::EvaluationFlags::values);
231
232 // dst-vector
233 temperature.reinit(cell);
234
235 for(unsigned int q = 0; q < energy.n_q_points; ++q)
236 {
237 // conserved variables
238 scalar rho = density.get_value(q);
239 vector rho_u = momentum.get_value(q);
240 scalar rho_E = energy.get_value(q);
241
242 // compute derived quantities in quadrature point
243 vector u = rho_u / rho;
244 scalar E = rho_E / rho;
245
246 scalar T = dealii::make_vectorized_array<Number>((heat_capacity_ratio - 1.0) /
247 specific_gas_constant) *
248 (E - 0.5 * scalar_product(u, u));
249
250 temperature.submit_value(T, q);
251 }
252
253 temperature.integrate_scatter(dealii::EvaluationFlags::values, dst);
254 }
255 }
256
257 dealii::MatrixFree<dim, Number> const * matrix_free;
258
259 unsigned int dof_index_all;
260 unsigned int dof_index_vector;
261 unsigned int dof_index_scalar;
262 unsigned int quad_index;
263
264 double heat_capacity_ratio;
265 double specific_gas_constant;
266};
267
268} // namespace CompNS
269} // namespace ExaDG
270
271#endif /* EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_CALCULATORS_H_ */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2