ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
indicator.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "hyperbolic_system.h"
9
10#include <compile_time_options.h>
12#include <simd.h>
13
14#include <deal.II/base/vectorization.h>
15
16
17namespace ryujin
18{
19 namespace ShallowWater
20 {
26 template <int dim, typename Number = double>
28 {
29 public:
33 static constexpr unsigned int problem_dimension =
34 HyperbolicSystem::problem_dimension<dim>;
35
40
44 static constexpr unsigned int n_precomputed_values =
45 HyperbolicSystem::n_precomputed_values<dim>;
46
51
56
61
63
81
85 Indicator(const HyperbolicSystem &hyperbolic_system,
87 &precomputed_values)
88 : hyperbolic_system(hyperbolic_system)
89 , precomputed_values(precomputed_values)
90 {
91 }
92
97 void reset(const unsigned int i, const state_type &U_i);
98
103 void add(const unsigned int *js,
104 const state_type &U_j,
105 const dealii::Tensor<1, dim, Number> &c_ij);
106
110 Number alpha(const Number h_i) const;
111
113
114 private:
119
120 const HyperbolicSystem &hyperbolic_system;
121
123 &precomputed_values;
124
125 Number h_i = 0.;
126 Number eta_i = 0.;
127 flux_type f_i;
128 state_type d_eta_i;
129 Number pressure_i = 0.;
130
131 Number left = 0.;
132 state_type right;
133
135 };
136
137
138 /* Inline definitions */
139
140
141 template <int dim, typename Number>
142 DEAL_II_ALWAYS_INLINE inline void
143 Indicator<dim, Number>::reset(const unsigned int i, const state_type &U_i)
144 {
145 h_i = hyperbolic_system.water_depth(U_i);
146 /* entropy viscosity commutator: */
147
148 const auto &[mathematical_entropy] =
149 precomputed_values.template get_tensor<Number, precomputed_type>(i);
150
151 eta_i = mathematical_entropy;
152
153 d_eta_i = hyperbolic_system.mathematical_entropy_derivative(U_i);
154 f_i = hyperbolic_system.f(U_i);
155 pressure_i = hyperbolic_system.pressure(U_i);
156
157 left = 0.;
158 right = 0.;
159 }
160
161
162 template <int dim, typename Number>
163 DEAL_II_ALWAYS_INLINE inline void
164 Indicator<dim, Number>::add(const unsigned int *js,
165 const state_type &U_j,
166 const dealii::Tensor<1, dim, Number> &c_ij)
167 {
168 /* entropy viscosity commutator: */
169
170 const auto &[eta_j] =
171 precomputed_values.template get_tensor<Number, precomputed_type>(js);
172
173 const auto velocity_j = hyperbolic_system.momentum(U_j) *
174 hyperbolic_system.inverse_water_depth_sharp(U_j);
175 const auto f_j = hyperbolic_system.f(U_j);
176 const auto pressure_j = hyperbolic_system.pressure(U_j);
177
178 left += (eta_j + pressure_j) * (velocity_j * c_ij);
179
180 for (unsigned int k = 0; k < problem_dimension; ++k)
181 right[k] += (f_j[k] - f_i[k]) * c_ij;
182 }
183
184
185 template <int dim, typename Number>
186 DEAL_II_ALWAYS_INLINE inline Number
187 Indicator<dim, Number>::alpha(const Number hd_i) const
188 {
190
191 Number my_sum = 0.;
192 for (unsigned int k = 0; k < problem_dimension; ++k) {
193 my_sum += d_eta_i[k] * right[k];
194 }
195
196 Number numerator = std::abs(left - my_sum);
197 Number denominator = std::abs(left) + std::abs(my_sum);
198
199 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
200 const auto quotient = std::abs(numerator + eps) /
201 (denominator + hd_i * std::abs(eta_i) + eps);
202
203 /* FIXME: this can be refactoring into a runtime parameter... */
204 const ScalarNumber evc_alpha_0_ = ScalarNumber(1.);
205
206 return std::min(Number(1.), evc_alpha_0_ * quotient);
207 }
208
209 } // namespace ShallowWater
210} // namespace ryujin
dealii::Tensor< 1, problem_dimension< dim >, dealii::Tensor< 1, dim, Number > > flux_type
std::array< Number, n_precomputed_values< dim > > precomputed_type
dealii::Tensor< 1, problem_dimension< dim >, Number > state_type
static constexpr unsigned int problem_dimension
Definition: indicator.h:33
static constexpr unsigned int n_precomputed_values
Definition: indicator.h:44
Number alpha(const Number h_i) const
Definition: indicator.h:187
HyperbolicSystem::flux_type< dim, Number > flux_type
Definition: indicator.h:55
void reset(const unsigned int i, const state_type &U_i)
Definition: indicator.h:143
HyperbolicSystem::precomputed_type< dim, Number > precomputed_type
Definition: indicator.h:50
Indicator(const HyperbolicSystem &hyperbolic_system, const MultiComponentVector< ScalarNumber, n_precomputed_values > &precomputed_values)
Definition: indicator.h:85
HyperbolicSystem::state_type< dim, Number > state_type
Definition: indicator.h:39
typename get_value_type< Number >::type ScalarNumber
Definition: indicator.h:60
void add(const unsigned int *js, const state_type &U_j, const dealii::Tensor< 1, dim, Number > &c_ij)
Definition: indicator.h:164