Skip to content

Commit a26a345

Browse files
keiwanjamalykeiwanjamaly
committed
Move FV reconstructor dimension to class template
1 parent c3a841b commit a26a345

7 files changed

Lines changed: 361 additions & 84 deletions

File tree

‎.clangd‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
CompileFlags:
2+
CompilationDatabase: build
3+
---
4+
If:
5+
PathMatch:
6+
- DiFfRG/build/.*
7+
Index:
8+
Background: Skip

‎DiFfRG/include/DiFfRG/discretization/FV/assembler/KurganovTadmor.hh‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,7 @@ namespace DiFfRG
113113
const std::array<std::array<autodiff::Real<1, NumberType>, n_components>, 2 * dim> &u_n)
114114
{
115115
const auto u_grad_deriv =
116-
Reconstructor::template compute_gradient_derivative<dim, n_components>(center, u_center, x_n, u_n);
116+
Reconstructor::template compute_gradient_derivative<n_components>(center, u_center, x_n, u_n);
117117
std::array<NumberType, n_components> result;
118118
for (size_t c = 0; c < n_components; ++c) {
119119
result[c] = derivative(u_center[c]);
@@ -427,8 +427,8 @@ namespace DiFfRG
427427
model.apply_boundary_conditions(neighbors_AD.u, neighbors_AD.x, boundary_ids, face_centers, cell_data_AD.u,
428428
cell_data_AD.x);
429429

430-
using ReconstructorAD = def::TVDReconstructor<typename Reconstructor::LimiterType, AD>;
431-
auto u_grad_cell_AD = ReconstructorAD::template compute_gradient<dim, n_components>(
430+
using ReconstructorAD = def::TVDReconstructor<dim, typename Reconstructor::LimiterType, AD>;
431+
auto u_grad_cell_AD = ReconstructorAD::template compute_gradient<n_components>(
432432
cell_data_AD.x, cell_data_AD.u, neighbors_AD.x, neighbors_AD.u);
433433

434434
std::array<dealii::Tensor<1, dim, AD>, n_components> ghost_grad_AD{};
@@ -451,7 +451,8 @@ namespace DiFfRG
451451
} // namespace internal
452452

453453
template <typename Discretization_, typename Model_,
454-
def::HasReconstructor Reconstructor_ = def::TVDReconstructor<def::MinModLimiter, double>,
454+
def::HasReconstructor Reconstructor_ =
455+
def::TVDReconstructor<Discretization_::dim, def::MinModLimiter, double>,
455456
def::HasWaveSpeed WaveSpeedStrategy_ = MaxEigenvalueWaveSpeed>
456457
requires MeshIsRectangular<typename Discretization_::Mesh>
457458
class Assembler : public AbstractAssembler<typename Discretization_::VectorType,
@@ -485,6 +486,7 @@ namespace DiFfRG
485486

486487
using Components = typename Discretization::Components;
487488
static constexpr uint dim = Discretization::dim;
489+
static_assert(Reconstructor::dim == dim, "Reconstructor dimension must match the discretization dimension.");
488490
static constexpr uint n_components = Components::count_fe_functions(0);
489491
static constexpr uint n_faces = GeometryInfo<dim>::faces_per_cell;
490492
// using CacheData = internal::Cache_Data<NumberType, dim, n_components>;
@@ -759,9 +761,9 @@ namespace DiFfRG
759761
const auto cell_data = get_cell_value(cell, solution_global);
760762
const auto ncell_data = get_cell_value(ncell, solution_global);
761763

762-
const GradientType u_grad_cell = Reconstructor::template compute_gradient<dim, n_components>(
764+
const GradientType u_grad_cell = Reconstructor::template compute_gradient<n_components>(
763765
cell_data.x, cell_data.u, cell_neighbors.x, cell_neighbors.u);
764-
const GradientType u_grad_ncell = Reconstructor::template compute_gradient<dim, n_components>(
766+
const GradientType u_grad_ncell = Reconstructor::template compute_gradient<n_components>(
765767
ncell_data.x, ncell_data.u, ncell_neighbors.x, ncell_neighbors.u);
766768

767769
const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
@@ -814,7 +816,7 @@ namespace DiFfRG
814816
const auto cell_neighbors = get_neighboring_cell_data(cell, solution_global, model);
815817
const auto cell_data = get_cell_value(cell, solution_global);
816818

817-
const GradientType u_grad_cell = Reconstructor::template compute_gradient<dim, n_components>(
819+
const GradientType u_grad_cell = Reconstructor::template compute_gradient<n_components>(
818820
cell_data.x, cell_data.u, cell_neighbors.x, cell_neighbors.u);
819821

820822
const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
@@ -1004,9 +1006,9 @@ namespace DiFfRG
10041006
copy_data_face.reinit(fe_iv, cell, ncell);
10051007

10061008
// Compute gradients and reconstructed interface values u_minus / u_plus
1007-
const GradientType u_grad_cell = Reconstructor::template compute_gradient<dim, n_components>(
1009+
const GradientType u_grad_cell = Reconstructor::template compute_gradient<n_components>(
10081010
cell_data.x, cell_data.u, cell_neighbors.x, cell_neighbors.u);
1009-
const GradientType u_grad_ncell = Reconstructor::template compute_gradient<dim, n_components>(
1011+
const GradientType u_grad_ncell = Reconstructor::template compute_gradient<n_components>(
10101012
ncell_data.x, ncell_data.u, ncell_neighbors.x, ncell_neighbors.u);
10111013

10121014
const std::array<NumberType, n_components> u_minus =
@@ -1104,7 +1106,7 @@ namespace DiFfRG
11041106
}
11051107

11061108
// Compute gradients and reconstructed interface values
1107-
const GradientType u_grad_cell = Reconstructor::template compute_gradient<dim, n_components>(
1109+
const GradientType u_grad_cell = Reconstructor::template compute_gradient<n_components>(
11081110
cell_data.x, cell_data.u, cell_neighbors.x, cell_neighbors.u);
11091111

11101112
// Ghost cell value and position from apply_boundary_conditions

‎DiFfRG/include/DiFfRG/discretization/FV/limiter/minmod_limiter.hh‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ namespace DiFfRG
2525
* @c Limiter template parameter of the TVDReconstructor:
2626
* @code
2727
* using Assembler = FV::KurganovTadmor::Assembler<Discretization, Model,
28-
* def::TVDReconstructor<def::MinModLimiter, double>>;
28+
* def::TVDReconstructor<Discretization::dim, def::MinModLimiter, double>>;
2929
* @endcode
3030
*/
3131
class MinModLimiter
@@ -45,4 +45,4 @@ namespace DiFfRG
4545
static_assert(HasSlopeLimiter<MinModLimiter>);
4646

4747
} // namespace def
48-
} // namespace DiFfRG
48+
} // namespace DiFfRG

‎DiFfRG/include/DiFfRG/discretization/FV/reconstructor/abstract_reconstructor.hh‎

Lines changed: 47 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,8 @@ namespace DiFfRG
1616
{
1717
namespace def
1818
{
19+
template <int dim> inline constexpr std::size_t n_faces = 2 * dim;
20+
1921
/**
2022
* @brief Per-component gradient type: one Tensor<1,dim> per solution component.
2123
*/
@@ -26,35 +28,58 @@ namespace DiFfRG
2628
*/
2729
template <typename T>
2830
concept HasReconstructor =
29-
// ---- compute_gradient (dim=1) ----
30-
requires(const dealii::Point<1> &center, const std::array<double, 1> &u_center,
31-
const std::array<dealii::Point<1>, 2> &x_n, const std::array<std::array<double, 1>, 2> &u_n) {
31+
requires {
32+
{
33+
T::dim
34+
} -> std::convertible_to<int>;
35+
} &&
36+
// ---- compute_gradient (n_components=1) ----
37+
requires(const dealii::Point<T::dim> &center, const std::array<double, 1> &u_center,
38+
const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
39+
const std::array<std::array<double, 1>, n_faces<T::dim>> &u_n) {
40+
{
41+
T::template compute_gradient<1>(center, u_center, x_n, u_n)
42+
} -> std::same_as<GradientType<T::dim, double, 1>>;
43+
} &&
44+
// ---- compute_gradient (n_components=2) ----
45+
requires(const dealii::Point<T::dim> &center, const std::array<double, 2> &u_center,
46+
const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
47+
const std::array<std::array<double, 2>, n_faces<T::dim>> &u_n) {
48+
{
49+
T::template compute_gradient<2>(center, u_center, x_n, u_n)
50+
} -> std::same_as<GradientType<T::dim, double, 2>>;
51+
} &&
52+
// ---- compute_gradient_at_point (n_components=1) ----
53+
requires(const dealii::Point<T::dim> &center, const dealii::Point<T::dim> &x,
54+
const std::array<double, 1> &u_center, const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
55+
const std::array<std::array<double, 1>, n_faces<T::dim>> &u_n) {
3256
{
33-
T::template compute_gradient<1, 1>(center, u_center, x_n, u_n)
34-
} -> std::same_as<GradientType<1, double, 1>>;
57+
T::template compute_gradient_at_point<1>(center, x, u_center, x_n, u_n)
58+
} -> std::same_as<GradientType<T::dim, double, 1>>;
3559
} &&
36-
// ---- compute_gradient (dim=2) ----
37-
requires(const dealii::Point<2> &center, const std::array<double, 2> &u_center,
38-
const std::array<dealii::Point<2>, 4> &x_n, const std::array<std::array<double, 2>, 4> &u_n) {
60+
// ---- compute_gradient_at_point (n_components=2) ----
61+
requires(const dealii::Point<T::dim> &center, const dealii::Point<T::dim> &x,
62+
const std::array<double, 2> &u_center, const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
63+
const std::array<std::array<double, 2>, n_faces<T::dim>> &u_n) {
3964
{
40-
T::template compute_gradient<2, 2>(center, u_center, x_n, u_n)
41-
} -> std::same_as<GradientType<2, double, 2>>;
65+
T::template compute_gradient_at_point<2>(center, x, u_center, x_n, u_n)
66+
} -> std::same_as<GradientType<T::dim, double, 2>>;
4267
} &&
43-
// ---- compute_gradient_derivative (dim=1) ----
44-
requires(const dealii::Point<1> &center, const std::array<autodiff::Real<1, double>, 1> &u_center,
45-
const std::array<dealii::Point<1>, 2> &x_n,
46-
const std::array<std::array<autodiff::Real<1, double>, 1>, 2> &u_n) {
68+
// ---- compute_gradient_derivative (n_components=1) ----
69+
requires(const dealii::Point<T::dim> &center, const std::array<autodiff::Real<1, double>, 1> &u_center,
70+
const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
71+
const std::array<std::array<autodiff::Real<1, double>, 1>, n_faces<T::dim>> &u_n) {
4772
{
48-
T::template compute_gradient_derivative<1, 1>(center, u_center, x_n, u_n)
49-
} -> std::same_as<GradientType<1, double, 1>>;
73+
T::template compute_gradient_derivative<1>(center, u_center, x_n, u_n)
74+
} -> std::same_as<GradientType<T::dim, double, 1>>;
5075
} &&
51-
// ---- compute_gradient_derivative (dim=2) ----
52-
requires(const dealii::Point<2> &center, const std::array<autodiff::Real<1, double>, 2> &u_center,
53-
const std::array<dealii::Point<2>, 4> &x_n,
54-
const std::array<std::array<autodiff::Real<1, double>, 2>, 4> &u_n) {
76+
// ---- compute_gradient_derivative (n_components=2) ----
77+
requires(const dealii::Point<T::dim> &center, const std::array<autodiff::Real<1, double>, 2> &u_center,
78+
const std::array<dealii::Point<T::dim>, n_faces<T::dim>> &x_n,
79+
const std::array<std::array<autodiff::Real<1, double>, 2>, n_faces<T::dim>> &u_n) {
5580
{
56-
T::template compute_gradient_derivative<2, 2>(center, u_center, x_n, u_n)
57-
} -> std::same_as<GradientType<2, double, 2>>;
81+
T::template compute_gradient_derivative<2>(center, u_center, x_n, u_n)
82+
} -> std::same_as<GradientType<T::dim, double, 2>>;
5883
};
5984

6085
} // namespace def

‎DiFfRG/include/DiFfRG/discretization/FV/reconstructor/tvd_reconstructor.hh‎

Lines changed: 46 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ namespace DiFfRG
3030
* When the limiter satisfies the TVD property (e.g. @c MinModLimiter),
3131
* the resulting reconstruction is Total Variation Diminishing.
3232
*
33+
* @tparam dim The spatial dimension of the reconstruction.
3334
* @tparam Limiter A type satisfying the @c HasSlopeLimiter concept
3435
* (e.g. @c MinModLimiter).
3536
* @tparam NumberType The numeric type used for computations
@@ -38,14 +39,17 @@ namespace DiFfRG
3839
* Usage with the KT assembler:
3940
* @code
4041
* using Assembler = FV::KurganovTadmor::Assembler<
41-
* Discretization, Model, def::TVDReconstructor<def::MinModLimiter, double>>;
42+
* Discretization, Model, def::TVDReconstructor<Discretization::dim, def::MinModLimiter, double>>;
4243
* @endcode
4344
*/
44-
template <HasSlopeLimiter Limiter, typename NumberType> class TVDReconstructor
45+
template <int dim_, HasSlopeLimiter Limiter, typename NumberType> class TVDReconstructor
4546
{
4647
using ADNumberType = autodiff::Real<1, NumberType>;
4748

4849
public:
50+
static constexpr int dim = dim_;
51+
static constexpr int n_faces = 2 * dim;
52+
4953
// Expose the Limiter type for use by other templates
5054
using LimiterType = Limiter;
5155
/**
@@ -55,22 +59,21 @@ namespace DiFfRG
5559
* two one-sided slopes from the cell centre to its axis-aligned
5660
* neighbours and returns the limiter-processed gradient.
5761
*
58-
* @tparam dim the spatial dimension
5962
* @tparam n_components the number of components in u
6063
*
6164
* @param center_pos position of the cell centre
6265
* @param u_center solution values at the cell centre
63-
* @param x_n positions of the 2*dim neighbouring cell centres,
66+
* @param x_n positions of the def::n_faces<dim> neighbouring cell centres,
6467
* ordered as pairs (left, right) for dim 0, then dim 1, …
6568
* @param u_n solution values at those neighbours
6669
*
6770
* @return per-component gradient vector
6871
*/
69-
template <int dim, int n_components>
72+
template <int n_components>
7073
static GradientType<dim, NumberType, n_components>
7174
compute_gradient(const dealii::Point<dim> &center_pos, const std::array<NumberType, n_components> &u_center,
72-
const std::array<dealii::Point<dim>, 2 * dim> &x_n,
73-
const std::array<std::array<NumberType, n_components>, 2 * dim> &u_n)
75+
const std::array<dealii::Point<dim>, n_faces> &x_n,
76+
const std::array<std::array<NumberType, n_components>, n_faces> &u_n)
7477
{
7578
GradientType<dim, NumberType, n_components> u_grad{};
7679
for (size_t c = 0; c < static_cast<size_t>(n_components); c++) {
@@ -92,6 +95,35 @@ namespace DiFfRG
9295
return u_grad;
9396
}
9497

98+
template <int n_components>
99+
static GradientType<dim, NumberType, n_components>
100+
compute_gradient_at_point(const dealii::Point<dim> &center_pos, const dealii::Point<dim> &x,
101+
const std::array<NumberType, n_components> &u_center,
102+
const std::array<dealii::Point<dim>, n_faces> &x_n,
103+
const std::array<std::array<NumberType, n_components>, n_faces> &u_n)
104+
{
105+
GradientType<dim, NumberType, n_components> u_grad{};
106+
107+
for (size_t c = 0; c < static_cast<size_t>(n_components); ++c) {
108+
const NumberType &u_val = u_center[c];
109+
for (int d = 0, i_n_1 = 0, i_n_2 = 1; d < dim; ++d, i_n_1 += 2, i_n_2 += 2) {
110+
const auto dx_1 = x_n[i_n_1] - center_pos;
111+
const NumberType du_1 = (u_n[i_n_1][c] - u_val) / dx_1[d];
112+
const auto dx_2 = x_n[i_n_2] - center_pos;
113+
const NumberType du_2 = (u_n[i_n_2][c] - u_val) / dx_2[d];
114+
115+
if (x[d] < center_pos[d])
116+
u_grad[c][d] = du_1;
117+
else if (x[d] > center_pos[d])
118+
u_grad[c][d] = du_2;
119+
else
120+
u_grad[c][d] = Limiter::slope_limit(du_1, du_2);
121+
}
122+
}
123+
124+
return u_grad;
125+
}
126+
95127
/**
96128
* @brief Compute the derivative of the limited gradient w.r.t. a single stencil DOF.
97129
*
@@ -103,25 +135,24 @@ namespace DiFfRG
103135
* reconstruction at a face point — it returns the per-component
104136
* derivative of the gradient itself.
105137
*
106-
* @tparam dim the spatial dimension
107138
* @tparam n_components the number of solution components
108139
*
109140
* @param center_pos position of the cell centre
110141
* @param u_center cell-centre values as seeded AD types
111-
* @param x_n positions of the 2*dim neighbouring cell centres
142+
* @param x_n positions of the def::n_faces<dim> neighbouring cell centres
112143
* @param u_n neighbour values as seeded AD types
113144
*
114145
* @return per-component gradient-derivative tensor d(grad u_c) / d(u_target)
115146
*/
116-
template <int dim, int n_components>
147+
template <int n_components>
117148
static GradientType<dim, NumberType, n_components>
118149
compute_gradient_derivative(const dealii::Point<dim> &center_pos,
119150
const std::array<ADNumberType, n_components> &u_center,
120-
const std::array<dealii::Point<dim>, 2 * dim> &x_n,
121-
const std::array<std::array<ADNumberType, n_components>, 2 * dim> &u_n)
151+
const std::array<dealii::Point<dim>, n_faces> &x_n,
152+
const std::array<std::array<ADNumberType, n_components>, n_faces> &u_n)
122153
{
123154
// Evaluate gradient reconstruction with AD types
124-
const auto u_grad_AD = TVDReconstructor<Limiter, ADNumberType>::template compute_gradient<dim, n_components>(
155+
const auto u_grad_AD = TVDReconstructor<dim, Limiter, ADNumberType>::template compute_gradient<n_components>(
125156
center_pos, u_center, x_n, u_n);
126157

127158
// Extract per-component gradient derivatives
@@ -135,7 +166,8 @@ namespace DiFfRG
135166
};
136167

137168
// Verify the default instantiation satisfies the concept.
138-
static_assert(HasReconstructor<TVDReconstructor<MinModLimiter, double>>);
169+
static_assert(HasReconstructor<TVDReconstructor<1, MinModLimiter, double>>);
170+
static_assert(HasReconstructor<TVDReconstructor<2, MinModLimiter, double>>);
139171

140172
} // namespace def
141173
} // namespace DiFfRG

‎DiFfRG/tests/discretization/FV/KurganovTadmor_tests.cc‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ namespace KT = DiFfRG::FV::KurganovTadmor;
1818
using WaveSpeedStrategy = DiFfRG::FV::KurganovTadmor::MaxEigenvalueWaveSpeed;
1919
using KT::internal::compute_numerical_flux;
2020
using KT::internal::reconstruct_u;
21-
using Reconstructor = DiFfRG::def::TVDReconstructor<DiFfRG::def::MinModLimiter, NumberType>;
21+
2222
struct CopyData {
2323
};
2424

@@ -67,6 +67,7 @@ TEST_CASE("reconstruct_u_derivative 1D single component", "[FV][KT]")
6767
constexpr int dim = 1;
6868
constexpr size_t n_components = 1;
6969
using AD = autodiff::Real<1, NumberType>;
70+
using Reconstructor = DiFfRG::def::TVDReconstructor<dim, DiFfRG::def::MinModLimiter, NumberType>;
7071

7172
const Point<dim> center(0.0);
7273
const Point<dim> x_q(0.5);
@@ -123,6 +124,7 @@ TEST_CASE("reconstruct_u_derivative at local extremum (vanishing limiter)", "[FV
123124
constexpr int dim = 1;
124125
constexpr size_t n_components = 1;
125126
using AD = autodiff::Real<1, NumberType>;
127+
using Reconstructor = DiFfRG::def::TVDReconstructor<dim, DiFfRG::def::MinModLimiter, NumberType>;
126128

127129
const Point<dim> center(0.0);
128130
const Point<dim> x_q(0.5);
@@ -162,6 +164,7 @@ TEST_CASE("reconstruct_u_derivative 2D single component", "[FV][KT]")
162164
constexpr int dim = 2;
163165
constexpr size_t n_components = 1;
164166
using AD = autodiff::Real<1, NumberType>;
167+
using Reconstructor = DiFfRG::def::TVDReconstructor<dim, DiFfRG::def::MinModLimiter, NumberType>;
165168

166169
const Point<dim> center(0.0, 0.0);
167170
const Point<dim> x_q(0.5, 0.3);

0 commit comments

Comments
 (0)