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 INCLUDE_EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_CALCULATORS_H_
23#define INCLUDE_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(dealii::EvaluationFlags::values);
160 pressure.set_dof_values(dst);
161 }
162 }
163
164 void
165 local_apply_velocity(dealii::MatrixFree<dim, Number> const & matrix_free,
166 VectorType & dst,
167 VectorType const & src,
168 std::pair<unsigned int, unsigned int> const & cell_range) const
169 {
170 // src-vector
171 CellIntegratorScalar density(matrix_free, dof_index_all, quad_index, 0);
172 CellIntegratorVector momentum(matrix_free, dof_index_all, quad_index, 1);
173
174 // dst-vector
175 CellIntegratorVector velocity(matrix_free, dof_index_vector, quad_index, 0);
176
177 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
178 {
179 // src-vector
180 density.reinit(cell);
181 density.read_dof_values(src);
182 density.evaluate(dealii::EvaluationFlags::values);
183 momentum.reinit(cell);
184 momentum.read_dof_values(src);
185 momentum.evaluate(dealii::EvaluationFlags::values);
186
187 // dst-vector
188 velocity.reinit(cell);
189
190 for(unsigned int q = 0; q < momentum.n_q_points; ++q)
191 {
192 // conserved variables
193 scalar rho = density.get_value(q);
194 vector rho_u = momentum.get_value(q);
195
196 // compute derived quantities in quadrature point
197 vector u = rho_u / rho;
198
199 velocity.submit_value(u, q);
200 }
201
202 velocity.integrate(dealii::EvaluationFlags::values);
203 velocity.set_dof_values(dst);
204 }
205 }
206
207 void
208 local_apply_temperature(dealii::MatrixFree<dim, Number> const & matrix_free,
209 VectorType & dst,
210 VectorType const & src,
211 std::pair<unsigned int, unsigned int> const & cell_range) const
212 {
213 // src-vector
214 CellIntegratorScalar density(matrix_free, dof_index_all, quad_index, 0);
215 CellIntegratorVector momentum(matrix_free, dof_index_all, quad_index, 1);
216 CellIntegratorScalar energy(matrix_free, dof_index_all, quad_index, 1 + dim);
217
218 // dst-vector
219 CellIntegratorScalar temperature(matrix_free, dof_index_scalar, quad_index, 0);
220
221 for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
222 {
223 // src-vector
224 density.reinit(cell);
225 density.read_dof_values(src);
226 density.evaluate(dealii::EvaluationFlags::values);
227 momentum.reinit(cell);
228 momentum.read_dof_values(src);
229 momentum.evaluate(dealii::EvaluationFlags::values);
230 energy.reinit(cell);
231 energy.read_dof_values(src);
232 energy.evaluate(dealii::EvaluationFlags::values);
233
234 // dst-vector
235 temperature.reinit(cell);
236
237 for(unsigned int q = 0; q < energy.n_q_points; ++q)
238 {
239 // conserved variables
240 scalar rho = density.get_value(q);
241 vector rho_u = momentum.get_value(q);
242 scalar rho_E = energy.get_value(q);
243
244 // compute derived quantities in quadrature point
245 vector u = rho_u / rho;
246 scalar E = rho_E / rho;
247
248 scalar T = dealii::make_vectorized_array<Number>((heat_capacity_ratio - 1.0) /
249 specific_gas_constant) *
250 (E - 0.5 * scalar_product(u, u));
251
252 temperature.submit_value(T, q);
253 }
254
255 temperature.integrate(dealii::EvaluationFlags::values);
256 temperature.set_dof_values(dst);
257 }
258 }
259
260 dealii::MatrixFree<dim, Number> const * matrix_free;
261
262 unsigned int dof_index_all;
263 unsigned int dof_index_vector;
264 unsigned int dof_index_scalar;
265 unsigned int quad_index;
266
267 double heat_capacity_ratio;
268 double specific_gas_constant;
269};
270
271} // namespace CompNS
272} // namespace ExaDG
273
274#endif /* INCLUDE_EXADG_COMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_CALCULATORS_H_ \
275 */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2