/
euler_module.h
227 lines (184 loc) · 6.25 KB
/
euler_module.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
//
// SPDX-License-Identifier: MIT
// Copyright (C) 2020 - 2021 by the ryujin authors
//
#pragma once
#include <compile_time_options.h>
#include "convenience_macros.h"
#include "simd.h"
#include "limiter.h"
#include "initial_values.h"
#include "offline_data.h"
#include "problem_description.h"
#include "sparse_matrix_simd.h"
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/sparse_matrix.templates.h>
#include <deal.II/lac/vector.h>
namespace ryujin
{
/**
* Explicit (strong stability preserving) time-stepping for the
* compressible Euler equations described in ProblemDescription.
*
* This module is described in detail in @cite KronbichlerMaier2021, Alg.
* 1.
*
* @todo Write out some more documentation
*
* @ingroup EulerModule
*/
template <int dim, typename Number = double>
class EulerModule final : public dealii::ParameterAcceptor
{
public:
/**
* @copydoc ProblemDescription::problem_dimension
*/
// clang-format off
static constexpr unsigned int problem_dimension = ProblemDescription::problem_dimension<dim>;
// clang-format on
/**
* @copydoc ProblemDescription::rank1_type
*/
using rank1_type =
ProblemDescription::rank1_type<dim, Number>;
/**
* @copydoc ProblemDescription::rank2_type
*/
using rank2_type = ProblemDescription::rank2_type<dim, Number>;
/**
* @copydoc OfflineData::scalar_type
*/
using scalar_type = typename OfflineData<dim, Number>::scalar_type;
/**
* @copydoc OfflineData::vector_type
*/
using vector_type = typename OfflineData<dim, Number>::vector_type;
/**
* Constructor.
*/
EulerModule(const MPI_Comm &mpi_communicator,
std::map<std::string, dealii::Timer> &computing_timer,
const ryujin::OfflineData<dim, Number> &offline_data,
const ryujin::ProblemDescription &problem_description,
const ryujin::InitialValues<dim, Number> &initial_values,
const std::string &subsection = "EulerModule");
/**
* Prepare time stepping. A call to @ref prepare() allocates temporary
* storage and is necessary before any of the following time-stepping
* functions can be called.
*/
void prepare();
/**
* @name Functons for performing explicit time steps
*/
//@{
/**
* Given a reference to a previous state vector U perform an explicit
* euler step (and store the result in U). The function returns the
* computed maximal time step size tau_max.
*
* The time step is performed with either tau_max (if tau == 0), or tau
* (if tau != 0). Here, tau_max is the computed maximal time step size
* and tau is the optional third parameter.
*/
Number euler_step(vector_type &U, Number t, Number tau = 0.);
/**
* Given a reference to a previous state vector U perform an explicit
* Heun 2nd order step (and store the result in U). The function
* returns the computed maximal time step size tau_max.
*
* The time step is performed with either tau_max (if tau == 0), or tau
* (if tau != 0). Here, tau_max is the computed maximal time step size
* and tau is the optional third parameter.
*
* See @cite Shu1988, Eq. 2.15.
*/
Number ssph2_step(vector_type &U, Number t, Number tau = 0.);
/**
* Given a reference to a previous state vector U perform an explicit
* SSP Runge Kutta 3rd order step (and store the result in U). The
* function returns the computed maximal time step size tau_max.
*
* The time step is performed with either tau_max (if tau == 0), or tau
* (if tau != 0). Here, tau_max is the computed maximal time step size
* and tau is the optional third parameter.
*
* See @cite Shu1988, Eq. 2.18.
*/
Number ssprk3_step(vector_type &U, Number t, Number tau = 0.);
/**
* Given a reference to a previous state vector U perform an explicit
* time step (and store the result in U). The function returns the
* chosen time step size tau.
*
* This function switches between euler_step(), ssph2_step(), or
* ssprk3_step() depending on selected approximation order.
*
* The time step is performed with either tau_max (if tau == 0), or tau
* (if tau != 0). Here, tau_max is the computed maximal time step size
* and tau is the optional third parameter.
*/
Number step(vector_type &U, Number t, Number tau = 0.);
//@}
/**
* @name Functons for postprocessing
*/
//@{
/**
* Computes an estimate for the local residual viscosity:
* \f{align}
* \alpha_i\,\big(-d_{ii}^{L,n}\big)\frac{1}{\#\mathcal{I}(i)}\beta_{ii}^{-1}
* \f}
*/
void compute_residual_mu();
private:
//@}
/**
* @name Internally used time-stepping primitives
*/
//@{
Number single_step(vector_type &U, Number tau);
void apply_boundary_conditions(vector_type &U, Number t);
//@}
/**
* @name Run time options
*/
//@{
Number cfl_update_;
Number cfl_max_;
unsigned int time_step_order_;
unsigned int limiter_iter_;
bool enforce_noslip_;
//@}
//@}
/**
* @name Internal data
*/
//@{
const MPI_Comm &mpi_communicator_;
std::map<std::string, dealii::Timer> &computing_timer_;
dealii::SmartPointer<const ryujin::OfflineData<dim, Number>> offline_data_;
dealii::SmartPointer<const ryujin::ProblemDescription> problem_description_;
dealii::SmartPointer<const ryujin::InitialValues<dim, Number>>
initial_values_;
unsigned int n_restarts_;
ACCESSOR_READ_ONLY(n_restarts)
scalar_type residual_mu_;
ACCESSOR_READ_ONLY(residual_mu)
scalar_type alpha_;
scalar_type second_variations_;
scalar_type specific_entropies_;
scalar_type evc_entropies_;
MultiComponentVector<Number, Limiter<dim, Number>::n_bounds> bounds_;
vector_type r_;
SparseMatrixSIMD<Number> dij_matrix_;
SparseMatrixSIMD<Number> lij_matrix_;
SparseMatrixSIMD<Number> lij_matrix_next_;
SparseMatrixSIMD<Number, problem_dimension> pij_matrix_;
vector_type temp_euler_;
vector_type temp_ssp_;
//@}
};
} /* namespace ryujin */