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
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_gravity_mode_no_gravity",
204 [](TConfig &self) {
205 self.gravity_config.gravity_mode = NoGravity;
206 })
207 .def(
208 "set_gravity_mode_cg",
209 [](TConfig &self) {
210 self.gravity_config.gravity_mode = CG;
211 })
212 .def(
213 "set_gravity_mode_pcg",
214 [](TConfig &self) {
215 self.gravity_config.gravity_mode = PCG;
216 })
217 .def(
218 "set_gravity_mode_bicgstab",
219 [](TConfig &self) {
220 self.gravity_config.gravity_mode = BICGSTAB;
221 })
222 .def("set_npscal_gas", [](TConfig &self, u32 npscal_gas) {
223 self.npscal_gas_config.npscal_gas = npscal_gas;
224 });
225
226 std::string sod_tube_analysis_name = name_model + "_AnalysisSodTube";
227 py::class_<TAnalysisSodTube>(m, sod_tube_analysis_name.c_str())
228 .def("compute_L2_dist", [](TAnalysisSodTube &self) -> std::tuple<Tscal, Tvec, Tscal> {
229 auto ret = self.compute_L2_dist();
230 return {ret.rho, ret.v, ret.P};
231 });
232
233 py::class_<T>(m, name_model.c_str())
234 .def("init", &T::init)
235 .def("init_scheduler", &T::init_scheduler)
236 .def("make_base_grid", &T::make_base_grid)
237 .def("dump_vtk", &T::dump_vtk)
238 .def("dump", &T::dump)
239 .def("load_from_dump", &T::load_from_dump)
240 .def("evolve_once_override_time", &T::evolve_once_time_expl)
241 .def("evolve_once", &T::evolve_once)
242 .def(
243 "evolve_until",
244 [](T &self, f64 target_time, i32 niter_max) {
245 return self.evolve_until(target_time, niter_max);
246 },
247 py::arg("target_time"),
248 py::kw_only(),
249 py::arg("niter_max") = -1)
250 .def("timestep", &T::timestep)
251 // .def("set_field_value_lambda_f64", &T::template set_field_value_lambda<f64>)
252 // .def("set_field_value_lambda_f64_3", &T::template set_field_value_lambda<f64_3>)
253 .def(
254 "set_field_value_lambda_f64",
255 [](T &self,
256 std::string field_name,
257 const std::function<f64(Tvec, Tvec)> pos_to_val,
258 const i32 offset) {
259 return self.template set_field_value_lambda<f64>(
260 field_name, pos_to_val, offset);
261 },
262 py::arg("field_name"),
263 py::arg("pos_to_val"),
264 py::arg("offset") = 0)
265 .def(
266 "set_field_value_lambda_f64_3",
267 [](T &self,
268 std::string field_name,
269 const std::function<f64_3(Tvec, Tvec)> pos_to_val,
270 const i32 offset) {
271 return self.template set_field_value_lambda<f64_3>(
272 field_name, pos_to_val, offset);
273 },
274 py::arg("field_name"),
275 py::arg("pos_to_val"),
276 py::arg("offset") = 0)
277 .def(
278 "gen_default_config",
279 [](T &self) -> TConfig {
280 return TConfig();
281 })
282 .def(
283 "set_solver_config",
284 [](T &self, TConfig cfg) {
285 if (self.ctx.is_scheduler_initialized()) {
286 shambase::throw_with_loc<std::runtime_error>(
287 "Cannot change solver config after scheduler is initialized");
288 }
289 cfg.check_config();
290 self.solver.solver_config = cfg;
291 })
292 .def(
293 "get_cell_coords",
294 [](T &self, std::pair<TgridVec, TgridVec> block_coord, u32 cell_local_id) {
295 return self.get_cell_coords(block_coord, cell_local_id);
296 })
297 .def(
298 "make_analysis_sodtube",
299 [](T &self,
300 shamphys::SodTube sod,
301 Tvec direction,
302 Tscal time_val,
303 Tscal x_ref,
304 Tscal x_min,
305 Tscal x_max) {
306 return std::make_unique<TAnalysisSodTube>(
307 self.ctx,
308 self.solver.solver_config,
309 self.solver.storage,
310 sod,
311 direction,
312 time_val,
313 x_ref,
314 x_min,
315 x_max);
316 })
317 .def(
318 "get_solver_tex",
319 [](T &self) {
320 return shambase::get_check_ref(self.solver.storage.solver_sequence).get_tex();
321 })
322 .def(
323 "get_solver_dot_graph",
324 [](T &self) {
325 return shambase::get_check_ref(self.solver.storage.solver_sequence)
326 .get_dot_graph();
327 })
328 .def(
329 "render_slice",
330 [](T &self, std::string name, std::string field_type, std::vector<Tvec> positions)
331 -> std::variant<std::vector<f64>, std::vector<f64_3>> {
332 if (field_type == "f64") {
333 ramses::modules::GridRender<Tvec, TgridVec, f64> render(
334 self.ctx, self.solver.solver_config, self.solver.storage);
335 return render.compute_slice(name, positions).copy_to_stdvec();
336 }
337
338 if (field_type == "f64_3") {
339 ramses::modules::GridRender<Tvec, TgridVec, f64_3> render(
340 self.ctx, self.solver.solver_config, self.solver.storage);
341 return render.compute_slice(name, positions).copy_to_stdvec();
342 }
343
344 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
345 })
346 .def(
347 "get_time",
348 [](T &self) {
349 return self.solver.solver_config.get_time();
350 })
351 .def(
352 "get_dt",
353 [](T &self) {
354 return self.solver.solver_config.get_dt();
355 })
356 .def(
357 "set_time",
358 [](T &self, Tscal t) {
359 return self.solver.solver_config.set_time(t);
360 })
361 .def("set_next_dt", [](T &self, Tscal dt) {
362 return self.solver.solver_config.set_next_dt(dt);
363 });
364 }
365} // namespace shammodels::basegodunov
366
367Register_pymod(pybasegodunovmodel) {
368
369 py::module mramses = m.def_submodule("model_ramses", "Shamrock Ramses solver");
370
371 std::string base_name = "RamsesModel";
372 using namespace shammodels::basegodunov;
373
374 add_instance<f64_3, i64_3>(
375 mramses, base_name + "_f64_3_i64_3_SolverConfig", base_name + "_f64_3_i64_3_Model");
376
377 using VariantAMRGodunovBind = std::variant<std::unique_ptr<Model<f64_3, i64_3>>>;
378
379 m.def(
380 "get_Model_Ramses",
381 [](ShamrockCtx &ctx,
382 std::string vector_type,
383 std::string grid_repr) -> VariantAMRGodunovBind {
384 VariantAMRGodunovBind ret;
385
386 if (vector_type == "f64_3" && grid_repr == "i64_3") {
387 ret = std::make_unique<Model<f64_3, i64_3>>(ctx);
388 } else {
389 throw shambase::make_except_with_loc<std::invalid_argument>(
390 "unknown combination of representation and grid_repr");
391 }
392
393 return ret;
394 },
395 py::kw_only(),
396 py::arg("context"),
397 py::arg("vector_type"),
398 py::arg("grid_repr"));
399}
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
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.
Register_pymod
#define Register_pymod(placeholdername)
Register a python module init function using static initialisation.
Definition pybindaliases.hpp:56
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 Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8