Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
pyRamsesModel.cpp
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
19
20#include "shambindings/pybindaliases.hpp"
21#include "shambindings/pytypealias.hpp"
22#include "shammodels/common/shamrock_json_to_py_json.hpp"
23#include "shammodels/ramses/Model.hpp"
24#include "shammodels/ramses/Solver.hpp"
25#include "shammodels/ramses/modules/AnalysisSodTube.hpp"
26#include "shammodels/ramses/modules/render/GridRender.hpp"
27#include <pybind11/functional.h>
28#include <pybind11/numpy.h>
29#include <memory>
30
31namespace shammodels::basegodunov {
32 template<class Tvec, class TgridVec>
33 void add_instance(py::module &m, std::string name_config, std::string name_model) {
34
35 using Tscal = shambase::VecComponent<Tvec>;
36 using Tgridscal = shambase::VecComponent<TgridVec>;
37
38 using T = Model<Tvec, TgridVec>;
39 using TConfig = typename T::Solver::Config;
40 using TAnalysisSodTube = shammodels::basegodunov::modules::AnalysisSodTube<Tvec, TgridVec>;
41
42 shamlog_debug_ln("[Py]", "registering class :", name_config, typeid(T).name());
43 shamlog_debug_ln("[Py]", "registering class :", name_model, typeid(T).name());
44
45 py::class_<TConfig> config_cls(m, name_config.c_str());
46
47 shammodels::common::add_json_defs<TConfig>(config_cls);
48
49 config_cls
50 .def(
51 "set_scale_factor",
52 [](TConfig &self, Tscal scale_factor) {
53 self.grid_coord_to_pos_fact = scale_factor;
54 })
55 .def(
56 "set_Csafe",
57 [](TConfig &self, Tscal Csafe) {
58 self.Csafe = Csafe;
59 })
60 .def(
61 "set_eos_gamma",
62 [](TConfig &self, Tscal eos_gamma) {
63 self.set_eos_gamma(eos_gamma);
64 })
65 .def(
66 "set_riemann_solver_hll",
67 [](TConfig &self) {
68 self.riemann_config = HLL;
69 })
70 .def(
71 "set_riemann_solver_hllc",
72 [](TConfig &self) {
73 self.riemann_config = HLLC;
74 })
75 .def(
76 "set_riemann_solver_rusanov",
77 [](TConfig &self) {
78 self.riemann_config = Rusanov;
79 })
80 .def(
81 "set_slope_lim_none",
82 [](TConfig &self) {
83 self.slope_config = None;
84 })
85 .def(
86 "set_slope_lim_vanleer_f",
87 [](TConfig &self) {
88 self.slope_config = VanLeer_f;
89 })
90 .def(
91 "set_slope_lim_vanleer_std",
92 [](TConfig &self) {
93 self.slope_config = VanLeer_std;
94 })
95 .def(
96 "set_slope_lim_vanleer_sym",
97 [](TConfig &self) {
98 self.slope_config = VanLeer_sym;
99 })
100 .def(
101 "set_slope_lim_minmod",
102 [](TConfig &self) {
103 self.slope_config = Minmod;
104 })
105 .def(
106 "set_scheduler_config",
107 [](TConfig &self, u64 split_crit, u64 merge_crit) {
108 self.scheduler_conf.split_load_value = split_crit;
109 self.scheduler_conf.merge_load_value = merge_crit;
110 },
111 py::kw_only(),
112 py::arg("split_load_value"),
113 py::arg("merge_load_value"))
114 .def(
115 "set_face_time_interpolation",
116 [](TConfig &self, bool face_time_interpolate) {
117 self.face_half_time_interpolation = face_time_interpolate;
118 })
119 .def(
120 "set_boundary_condition",
121 [](TConfig &self, const std::string &axis, const std::string &bc_type) {
122 BCConfig::GhostType ghost_type;
123 if (bc_type == "periodic") {
124 ghost_type = BCConfig::GhostType::Periodic;
125 } else if (bc_type == "reflective") {
126 ghost_type = BCConfig::GhostType::Reflective;
127 } else if (bc_type == "outflow") {
128 ghost_type = BCConfig::GhostType::Outflow;
129 } else {
130 throw std::invalid_argument(
131 "Unsupported boundary condition type: " + bc_type);
132 }
133
134 if (axis == "x") {
135 self.bc_config.set_x(ghost_type);
136 } else if (axis == "y") {
137 self.bc_config.set_y(ghost_type);
138 } else if (axis == "z") {
139 self.bc_config.set_z(ghost_type);
140 } else {
141 throw std::invalid_argument("Unsupported axis: " + axis);
142 }
143 },
144 py::arg("axis"),
145 py::arg("bc_type"))
146 .def(
147 "set_dust_mode_dhll",
148 [](TConfig &self, u32 ndust) {
149 self.dust_config = {DHLL, ndust};
150 })
151 .def(
152 "set_dust_mode_hb",
153 [](TConfig &self, u32 ndust) {
154 self.dust_config = {HB, ndust};
155 })
156 .def(
157 "set_dust_mode_none",
158 [](TConfig &self) {
159 self.dust_config = {NoDust, 0};
160 })
161 .def(
162 "set_alpha_values",
163 [](TConfig &self, f32 alpha_values) {
164 return self.set_alphas_static(alpha_values);
165 })
166 .def(
167 "set_drag_mode_no_drag",
168 [](TConfig &self) {
169 self.drag_config.drag_solver_config = NoDrag;
170 self.drag_config.enable_frictional_heating = false;
171 })
172 .def(
173 "set_drag_mode_irk1",
174 [](TConfig &self, bool frictional_status) {
175 self.drag_config.drag_solver_config = IRK1;
176 self.drag_config.enable_frictional_heating = frictional_status;
177 })
178 .def(
179 "set_drag_mode_irk2",
180 [](TConfig &self, bool frictional_status) {
181 self.drag_config.drag_solver_config = IRK2;
182 self.drag_config.enable_frictional_heating = frictional_status;
183 })
184 .def(
185 "set_drag_mode_expo",
186 [](TConfig &self, bool frictional_status) {
187 self.drag_config.drag_solver_config = EXPO;
188 self.drag_config.enable_frictional_heating = frictional_status;
189 })
190 .def(
191 "set_amr_mode_none",
192 [](TConfig &self) {
193 self.amr_mode.set_refine_none();
194 })
195 .def(
196 "set_amr_mode_density_based",
197 [](TConfig &self, Tscal crit_mass) {
198 self.amr_mode.set_refine_density_based(crit_mass);
199 },
200 py::kw_only(),
201 py::arg("crit_mass"))
202 .def(
203 "set_amr_mode_pseudo_gradient_based",
204 [](TConfig &self, Tscal error_min, Tscal error_max) {
205 self.amr_mode.set_refine_pseudo_gradient_based(error_min, error_max);
206 },
207 py::kw_only(),
208 py::arg("error_min"),
209 py::arg("error_max"))
210 .def(
211 "set_amr_mode_jeans_length_based",
212 [](TConfig &self, u32 N_jeans, Tscal T_init) {
213 self.amr_mode.set_refine_jeans_length_based(N_jeans, T_init);
214 },
215 py::kw_only(),
216 py::arg("N_jeans"),
217 py::arg("T_init"))
218 .def(
219 "set_amr_mode_shear_based",
220 [](TConfig &self, Tscal threshold) {
221 self.amr_mode.set_refine_shear_based(threshold);
222 },
223 py::kw_only(),
224 py::arg("Threshold"))
225 .def(
226 "set_amr_mode_old",
227 [](TConfig &self, bool use_old_amr) {
228 self.amr_mode.old_amr = use_old_amr;
229 })
230 .def(
231 "set_gravity_mode_no_gravity",
232 [](TConfig &self) {
233 self.gravity_config.gravity_mode = NoGravity;
234 })
235 .def(
236 "set_gravity_mode_cg",
237 [](TConfig &self) {
238 self.gravity_config.gravity_mode = CG;
239 })
240 .def(
241 "set_gravity_mode_pcg",
242 [](TConfig &self) {
243 self.gravity_config.gravity_mode = PCG;
244 })
245 .def(
246 "set_gravity_mode_bicgstab",
247 [](TConfig &self) {
248 self.gravity_config.gravity_mode = BICGSTAB;
249 })
250 .def("set_npscal_gas", [](TConfig &self, u32 npscal_gas) {
251 self.npscal_gas_config.npscal_gas = npscal_gas;
252 });
253
254 std::string sod_tube_analysis_name = name_model + "_AnalysisSodTube";
255 py::class_<TAnalysisSodTube>(m, sod_tube_analysis_name.c_str())
256 .def("compute_L2_dist", [](TAnalysisSodTube &self) -> std::tuple<Tscal, Tvec, Tscal> {
257 auto ret = self.compute_L2_dist();
258 return {ret.rho, ret.v, ret.P};
259 });
260
261 py::class_<T>(m, name_model.c_str())
262 .def("init", &T::init)
263 .def("init_scheduler", &T::init_scheduler)
264 .def("make_base_grid", &T::make_base_grid)
265 .def("dump_vtk", &T::dump_vtk)
266 .def("dump", &T::dump)
267 .def("load_from_dump", &T::load_from_dump)
268 .def("evolve_once_override_time", &T::evolve_once_time_expl)
269 .def("evolve_once", &T::evolve_once)
270 .def(
271 "evolve_until",
272 [](T &self, f64 target_time, i32 niter_max) {
273 return self.evolve_until(target_time, niter_max);
274 },
275 py::arg("target_time"),
276 py::kw_only(),
277 py::arg("niter_max") = -1)
278 .def("timestep", &T::timestep)
279 .def(
280 "set_field_value_lambda_f64",
281 [](T &self,
282 std::string field_name,
283 const std::function<f64(Tvec, Tvec)> pos_to_val,
284 const i32 offset) {
285 return self.template set_field_value_lambda<f64>(
286 field_name, pos_to_val, offset);
287 },
288 py::arg("field_name"),
289 py::arg("pos_to_val"),
290 py::arg("offset") = 0)
291 .def(
292 "set_field_value_lambda_f64_3",
293 [](T &self,
294 std::string field_name,
295 const std::function<f64_3(Tvec, Tvec)> pos_to_val,
296 const i32 offset) {
297 return self.template set_field_value_lambda<f64_3>(
298 field_name, pos_to_val, offset);
299 },
300 py::arg("field_name"),
301 py::arg("pos_to_val"),
302 py::arg("offset") = 0)
303 .def(
304 "gen_default_config",
305 [](T &self) -> TConfig {
306 return TConfig();
307 })
308 .def(
309 "set_solver_config",
310 [](T &self, TConfig cfg) {
311 if (self.ctx.is_scheduler_initialized()) {
312 shambase::throw_with_loc<std::runtime_error>(
313 "Cannot change solver config after scheduler is initialized");
314 }
315 cfg.check_config();
316 self.solver.solver_config = cfg;
317 })
318 .def(
319 "get_cell_coords",
320 [](T &self, std::pair<TgridVec, TgridVec> block_coord, u32 cell_local_id) {
321 return self.get_cell_coords(block_coord, cell_local_id);
322 })
323 .def(
324 "make_analysis_sodtube",
325 [](T &self,
326 shamphys::SodTube sod,
327 Tvec direction,
328 Tscal time_val,
329 Tscal x_ref,
330 Tscal x_min,
331 Tscal x_max) {
332 return std::make_unique<TAnalysisSodTube>(
333 self.ctx,
334 self.solver.solver_config,
335 self.solver.storage,
336 sod,
337 direction,
338 time_val,
339 x_ref,
340 x_min,
341 x_max);
342 })
343 .def(
344 "get_solver_tex",
345 [](T &self) {
346 return shambase::get_check_ref(self.solver.storage.solver_sequence).get_tex();
347 })
348 .def(
349 "get_solver_dot_graph",
350 [](T &self) {
351 return shambase::get_check_ref(self.solver.storage.solver_sequence)
352 .get_dot_graph();
353 })
354 .def(
355 "render_slice",
356 [](T &self, std::string name, std::string field_type, std::vector<Tvec> positions)
357 -> std::variant<std::vector<f64>, std::vector<f64_3>> {
358 if (field_type == "f64") {
359 ramses::modules::GridRender<Tvec, TgridVec, f64> render(
360 self.ctx, self.solver.solver_config, self.solver.storage);
361 return render.compute_slice(name, positions).copy_to_stdvec();
362 }
363
364 if (field_type == "f64_3") {
365 ramses::modules::GridRender<Tvec, TgridVec, f64_3> render(
366 self.ctx, self.solver.solver_config, self.solver.storage);
367 return render.compute_slice(name, positions).copy_to_stdvec();
368 }
369
370 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
371 })
372 .def(
373 "get_time",
374 [](T &self) {
375 return self.solver.get_time();
376 })
377 .def(
378 "get_dt",
379 [](T &self) {
380 return self.solver.get_dt();
381 })
382 .def(
383 "set_time",
384 [](T &self, Tscal t) {
385 return self.solver.set_time(t);
386 })
387 .def("set_next_dt", [](T &self, Tscal dt) {
388 return self.solver.set_next_dt(dt);
389 });
390 }
391} // namespace shammodels::basegodunov
392
393ON_PYTHON_INIT {
394 auto &m = root_module;
395
396 py::module mramses = m.def_submodule("model_ramses", "Shamrock Ramses solver");
397
398 std::string base_name = "RamsesModel";
399 using namespace shammodels::basegodunov;
400
401 add_instance<f64_3, i64_3>(
402 mramses, base_name + "_f64_3_i64_3_SolverConfig", base_name + "_f64_3_i64_3_Model");
403
404 using VariantAMRGodunovBind = std::variant<std::unique_ptr<Model<f64_3, i64_3>>>;
405
406 m.def(
407 "get_Model_Ramses",
408 [](ShamrockCtx &ctx,
409 std::string vector_type,
410 std::string grid_repr) -> VariantAMRGodunovBind {
411 VariantAMRGodunovBind ret;
412
413 if (vector_type == "f64_3" && grid_repr == "i64_3") {
414 ret = std::make_unique<Model<f64_3, i64_3>>(ctx);
415 } else {
416 throw shambase::make_except_with_loc<std::invalid_argument>(
417 "unknown combination of representation and grid_repr");
418 }
419
420 return ret;
421 },
422 py::kw_only(),
423 py::arg("context"),
424 py::arg("vector_type"),
425 py::arg("grid_repr"));
426}
f64
double f64
Alias for double.
Definition aliases_float.hpp:20
f32
float f32
Alias for float.
Definition aliases_float.hpp:19
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
u64
std::uint64_t u64
64 bit unsigned integer
Definition aliases_int.hpp:26
i32
std::int32_t i32
32 bit integer
Definition aliases_int.hpp:23
shambase::throw_with_loc
void throw_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Throw an exception and append the source location to it.
Definition exception.hpp:132
shambase::get_check_ref
T & get_check_ref(const std::unique_ptr< T > &ptr, SourceLocation loc=SourceLocation())
Takes a std::unique_ptr and returns a reference to the object it holds. It throws a std::runtime_erro...
Definition memory.hpp:110
shambase::make_except_with_loc
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
Definition exception.hpp:105
shammodels::basegodunov
namespace for the basegodunov model
Definition AMRCellInfos.hpp:20
shammodels::basegodunov::HB
@ HB
Huang and Bai. Pressureless Riemann solver by Huang and Bai (2022) in Athena++.
Definition enum_DustRiemannSolverMode.hpp:33
shammodels::basegodunov::DHLL
@ DHLL
Dust HLL. This is merely the HLL solver for dust. It's then a Rusanov like.
Definition enum_DustRiemannSolverMode.hpp:31
shammodels::basegodunov::NoDust
@ NoDust
No dust, so no Riemann solver is used.
Definition enum_DustRiemannSolverMode.hpp:29
pybindaliases.hpp
Pybind11 include and definitions.
ON_PYTHON_INIT
#define ON_PYTHON_INIT
Register a Python module init function using static initialization.
Definition pybindaliases.hpp:69
shamrock_json_to_py_json.hpp
Utilities to convert JSON objects to Python objects and vice versa. TODO: try to convert directly wit...
Generated on for Shamrock by doxygen 1.15.0