Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
SlopeLimitedGradientUtilities.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
20#include "shamcomm/logs.hpp"
21#include "shammath/riemann.hpp"
22#include "shammath/slopeLimiter.hpp"
23#include "shammodels/common/amr/NeighGraph.hpp"
24#include "shammodels/ramses/SolverConfig.hpp"
25#include <type_traits>
26
27namespace {
28 using AMRGraphLinkiterator = shammodels::basegodunov::modules::AMRGraph::ro_access;
29 using Direction = shammodels::basegodunov::modules::Direction;
30
31 template<class T>
32 inline T slope_function_van_leer_f_form(T sL, T sR) {
33 T st = sL + sR;
34
35 auto vanleer = [](T f) {
36 return 4. * f * (1. - f);
37 };
38
39 auto slopelim = [&](T f) {
40 if constexpr (std::is_same_v<T, f64_3>) {
41 f.x() = (f.x() >= 0 && f.x() <= 1) ? f.x() : 0;
42 f.y() = (f.y() >= 0 && f.y() <= 1) ? f.y() : 0;
43 f.z() = (f.z() >= 0 && f.z() <= 1) ? f.z() : 0;
44 } else {
45 f = (f >= 0 && f <= 1) ? f : 0;
46 }
47 return vanleer(f);
48 };
49
50 return slopelim(sL / st) * st * 0.5;
51 }
52
53 template<class T>
54 inline T slope_function_van_leer_symetric(T sL, T sR) {
55
56 if constexpr (std::is_same_v<T, f64_3>) {
57 return {
58 shammath::van_leer_slope_symetric(sL[0], sR[0]),
59 shammath::van_leer_slope_symetric(sL[1], sR[1]),
60 shammath::van_leer_slope_symetric(sL[2], sR[2])};
61 } else {
62 return shammath::van_leer_slope_symetric(sL, sR);
63 }
64 }
65
66 template<class T>
67 inline T slope_function_van_leer_standard(T sL, T sR) {
68
69 if constexpr (std::is_same_v<T, f64_3>) {
70 return {
71 shammath::van_leer_slope(sL[0], sR[0]),
72 shammath::van_leer_slope(sL[1], sR[1]),
73 shammath::van_leer_slope(sL[2], sR[2])};
74 } else {
75 return shammath::van_leer_slope(sL, sR);
76 }
77 }
78
79 template<class T>
80 inline T slope_function_minmod(T sL, T sR) {
81
82 if constexpr (std::is_same_v<T, f64_3>) {
83 return {
84 shammath::minmod(sL[0], sR[0]),
85 shammath::minmod(sL[1], sR[1]),
86 shammath::minmod(sL[2], sR[2])};
87 } else {
88 return shammath::minmod(sL, sR);
89 }
90 }
91
92 using SlopeMode = shammodels::basegodunov::SlopeMode;
93
94 template<class T, SlopeMode mode>
95 inline T slope_function(T sL, T sR) {
96 if constexpr (mode == SlopeMode::None) {
97 return sham::VectorProperties<T>::get_zero();
98 }
99
100 if constexpr (mode == SlopeMode::VanLeer_f) {
101 return slope_function_van_leer_f_form(sL, sR);
102 }
103
104 if constexpr (mode == SlopeMode::VanLeer_std) {
105 return slope_function_van_leer_standard(sL, sR);
106 }
107
108 if constexpr (mode == SlopeMode::VanLeer_sym) {
109 return slope_function_van_leer_symetric(sL, sR);
110 }
111
112 if constexpr (mode == SlopeMode::Minmod) {
113 return slope_function_minmod(sL, sR);
114 }
115 }
116
135 template<class Tfield, class Tvec, SlopeMode mode, class ACCField>
136 inline std::array<Tfield, 3> get_3d_grad(
137 const u32 cell_global_id,
138 const shambase::VecComponent<Tvec> delta_cell,
139 const AMRGraphLinkiterator &graph_iter_xp,
140 const AMRGraphLinkiterator &graph_iter_xm,
141 const AMRGraphLinkiterator &graph_iter_yp,
142 const AMRGraphLinkiterator &graph_iter_ym,
143 const AMRGraphLinkiterator &graph_iter_zp,
144 const AMRGraphLinkiterator &graph_iter_zm,
145 ACCField &&field_access) {
146
147 auto get_avg_neigh = [&](auto &graph_links) -> Tfield {
148 Tfield acc = shambase::VectorProperties<Tfield>::get_zero();
149 u32 cnt = graph_links.for_each_object_link_cnt(cell_global_id, [&](u32 id_b) {
150 acc += field_access(id_b);
151 });
152 return (cnt > 0) ? acc / cnt : shambase::VectorProperties<Tfield>::get_zero();
153 };
154
155 Tfield W_i = field_access(cell_global_id);
156 Tfield W_xp = get_avg_neigh(graph_iter_xp);
157 Tfield W_xm = get_avg_neigh(graph_iter_xm);
158 Tfield W_yp = get_avg_neigh(graph_iter_yp);
159 Tfield W_ym = get_avg_neigh(graph_iter_ym);
160 Tfield W_zp = get_avg_neigh(graph_iter_zp);
161 Tfield W_zm = get_avg_neigh(graph_iter_zm);
162
163 Tfield delta_W_x_p = W_xp - W_i;
164 Tfield delta_W_y_p = W_yp - W_i;
165 Tfield delta_W_z_p = W_zp - W_i;
166
167 Tfield delta_W_x_m = W_i - W_xm;
168 Tfield delta_W_y_m = W_i - W_ym;
169 Tfield delta_W_z_m = W_i - W_zm;
170
171 Tfield fact = 1. / Tfield(delta_cell);
172
173 Tfield lim_slope_W_x = slope_function<Tfield, mode>(delta_W_x_m * fact, delta_W_x_p * fact);
174 Tfield lim_slope_W_y = slope_function<Tfield, mode>(delta_W_y_m * fact, delta_W_y_p * fact);
175 Tfield lim_slope_W_z = slope_function<Tfield, mode>(delta_W_z_m * fact, delta_W_z_p * fact);
176
177 return {lim_slope_W_x, lim_slope_W_y, lim_slope_W_z};
178 }
179
202 template<class Tvec, SlopeMode mode, class ACCField1, class ACCField2, class ACCField3>
203 inline std::array<shammath::ConsState<Tvec>, 3> get_3d_grad_cons(
204 const u32 cell_global_id,
205 const shambase::VecComponent<Tvec> delta_cell,
206 const AMRGraphLinkiterator &graph_iter_xp,
207 const AMRGraphLinkiterator &graph_iter_xm,
208 const AMRGraphLinkiterator &graph_iter_yp,
209 const AMRGraphLinkiterator &graph_iter_ym,
210 const AMRGraphLinkiterator &graph_iter_zp,
211 const AMRGraphLinkiterator &graph_iter_zm,
212 ACCField1 &&field_access_rho,
213 ACCField2 &&field_access_rho_vel,
214 ACCField3 &&field_access_rhoe) {
215
216 using Tscal = shambase::VecComponent<Tvec>;
217
218 auto get_avg_neigh = [&](auto &graph_links) -> shammath::ConsState<Tvec> {
219 Tscal acc_rho = shambase::VectorProperties<Tscal>::get_zero();
220 Tscal acc_rhoe = shambase::VectorProperties<Tscal>::get_zero();
221 Tvec acc_rho_vel = shambase::VectorProperties<Tvec>::get_zero();
222 u32 cnt = graph_links.for_each_object_link_cnt(cell_global_id, [&](u32 id_b) {
223 acc_rho += field_access_rho(id_b);
224 acc_rho_vel += field_access_rho_vel(id_b);
225 acc_rhoe += field_access_rhoe(id_b);
226 });
227
228 shammath::ConsState<Tvec> res
229 = {shambase::VectorProperties<Tscal>::get_zero(),
230 shambase::VectorProperties<Tscal>::get_zero(),
231
232 {shambase::VectorProperties<Tscal>::get_zero(),
233 shambase::VectorProperties<Tscal>::get_zero(),
234 shambase::VectorProperties<Tscal>::get_zero()}};
235
236 if (cnt > 0) {
237 res = {acc_rho, acc_rhoe, acc_rho_vel};
238 res *= (1. / cnt);
239 }
240
241 return res;
242 };
243
244 shammath::ConsState<Tvec> W_i
245 = {field_access_rho(cell_global_id),
246 field_access_rhoe(cell_global_id),
247 field_access_rho_vel(cell_global_id)};
248
249 shammath::ConsState<Tvec> W_xp = get_avg_neigh(graph_iter_xp);
250 shammath::ConsState<Tvec> W_xm = get_avg_neigh(graph_iter_xm);
251 shammath::ConsState<Tvec> W_yp = get_avg_neigh(graph_iter_yp);
252 shammath::ConsState<Tvec> W_ym = get_avg_neigh(graph_iter_ym);
253 shammath::ConsState<Tvec> W_zp = get_avg_neigh(graph_iter_zp);
254 shammath::ConsState<Tvec> W_zm = get_avg_neigh(graph_iter_zm);
255
256 shammath::ConsState<Tvec> delta_W_x_p = W_xp - W_i;
257 shammath::ConsState<Tvec> delta_W_y_p = W_yp - W_i;
258 shammath::ConsState<Tvec> delta_W_z_p = W_zp - W_i;
259
260 shammath::ConsState<Tvec> delta_W_x_m = W_i - W_xm;
261 shammath::ConsState<Tvec> delta_W_y_m = W_i - W_ym;
262 shammath::ConsState<Tvec> delta_W_z_m = W_i - W_zm;
263
264 Tscal fact = 1. / delta_cell;
265
266 shammath::ConsState<Tvec> lim_slope_W_x
267 = {slope_function<Tscal, mode>(delta_W_x_m.rho * fact, delta_W_x_p.rho * fact),
268 slope_function<Tscal, mode>(delta_W_x_m.rhoe * fact, delta_W_x_p.rhoe * fact),
269 slope_function<Tvec, mode>(delta_W_x_m.rhovel * fact, delta_W_x_p.rhovel * fact)};
270
271 shammath::ConsState<Tvec> lim_slope_W_y
272 = {slope_function<Tscal, mode>(delta_W_y_m.rho * fact, delta_W_y_p.rho * fact),
273 slope_function<Tscal, mode>(delta_W_y_m.rhoe * fact, delta_W_y_p.rhoe * fact),
274 slope_function<Tvec, mode>(delta_W_y_m.rhovel * fact, delta_W_y_p.rhovel * fact)};
275
276 shammath::ConsState<Tvec> lim_slope_W_z
277 = {slope_function<Tscal, mode>(delta_W_z_m.rho * fact, delta_W_z_p.rho * fact),
278 slope_function<Tscal, mode>(delta_W_z_m.rhoe * fact, delta_W_z_p.rhoe * fact),
279 slope_function<Tvec, mode>(delta_W_z_m.rhovel * fact, delta_W_z_p.rhovel * fact)};
280
281 return {lim_slope_W_x, lim_slope_W_y, lim_slope_W_z};
282 }
283} // namespace
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
shammath::van_leer_slope
T van_leer_slope(T f, T g)
Van leer slope limiter.
Definition slopeLimiter.hpp:40
shammodels::basegodunov::SlopeMode
SlopeMode
Slope limiter modes.
Definition enum_SlopeMode.hpp:25
riemann.hpp
From original version by Thomas Guillet (T.A.Guillet@exeter.ac.uk)
slopeLimiter.hpp
Math header to compute slope limiters.
shammodels::basegodunov::modules::NeighGraph::ro_access
Definition NeighGraph.hpp:62
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8