Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
pySPHModel.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
18#include "shambase/exception.hpp"
19#include "shambase/logs/loglevels.hpp"
20#include "shambase/memory.hpp"
21#include "shambindings/pybindaliases.hpp"
22#include "shambindings/pytypealias.hpp"
23#include "shamcomm/worldInfo.hpp"
24#include "shammath/sphkernels.hpp"
25#include "shammodels/common/shamrock_json_to_py_json.hpp"
26#include "shammodels/sph/Model.hpp"
27#include "shammodels/sph/io/PhantomDump.hpp"
28#include "shammodels/sph/modules/AnalysisAngularMomentum.hpp"
29#include "shammodels/sph/modules/AnalysisBarycenter.hpp"
30#include "shammodels/sph/modules/AnalysisDisc.hpp"
31#include "shammodels/sph/modules/AnalysisEnergyKinetic.hpp"
32#include "shammodels/sph/modules/AnalysisEnergyPotential.hpp"
33#include "shammodels/sph/modules/AnalysisSodTube.hpp"
34#include "shammodels/sph/modules/AnalysisTotalMomentum.hpp"
35#include "shammodels/sph/modules/render/CartesianRender.hpp"
36#include "shammodels/sph/modules/render/RenderFieldGetter.hpp"
37#include "shamphys/SodTube.hpp"
38#include "shamrock/scheduler/PatchScheduler.hpp"
39#include <pybind11/cast.h>
40#include <pybind11/numpy.h>
41#include <pybind11/pytypes.h>
42#include <memory>
43#include <optional>
44#include <random>
45#include <utility>
46
47template<class Tvec, template<class> class SPHKernel>
48void add_instance(py::module &m, std::string name_config, std::string name_model) {
49 using namespace shammodels::sph;
50
51 using Tscal = shambase::VecComponent<Tvec>;
52
53 using T = Model<Tvec, SPHKernel>;
54
55 using TAnalysisSodTube = shammodels::sph::modules::AnalysisSodTube<Tvec, SPHKernel>;
56 using TAnalysisDisc = shammodels::sph::modules::AnalysisDisc<Tvec, SPHKernel>;
57 using TSPHSetup = shammodels::sph::modules::SPHSetup<Tvec, SPHKernel>;
58 using TConfig = typename T::Solver::Config;
59
60 using custom_getter_t = std::function<pybind11::array_t<f64>(size_t, pybind11::dict &)>;
61
62 shamlog_debug_ln("[Py]", "registering class :", name_config, typeid(T).name());
63 shamlog_debug_ln("[Py]", "registering class :", name_model, typeid(T).name());
64
65 py::class_<TConfig> config_cls(m, name_config.c_str());
66
67 shammodels::common::add_json_defs<TConfig>(config_cls);
68
69 config_cls.def("print_status", &TConfig::print_status)
70 .def("set_particle_tracking", &TConfig::set_particle_tracking)
71 .def(
72 "set_scheduler_config",
73 [](TConfig &self, u64 split_crit, u64 merge_crit) {
74 self.scheduler_conf.split_load_value = split_crit;
75 self.scheduler_conf.merge_load_value = merge_crit;
76 },
77 py::kw_only(),
78 py::arg("split_load_value"),
79 py::arg("merge_load_value"))
80 .def("set_tree_reduction_level", &TConfig::set_tree_reduction_level)
81 .def("set_two_stage_search", &TConfig::set_two_stage_search)
82 .def("set_show_neigh_stats", &TConfig::set_show_neigh_stats)
83 .def(
84 "set_max_neigh_cache_size",
85 [](TConfig &self, const py::object &max_neigh_cache_size) {
86 ON_RANK_0(shamlog_warn_ln(
87 "SPH",
88 ".set_max_neigh_cache_size() is deprecated,\n"
89 " -> calling this is a no-op,\n"
90 " -> you can remove the call to that function"););
91 })
92 .def("set_smoothing_length_density_based", &TConfig::set_smoothing_length_density_based)
93 .def(
94 "set_smoothing_length_density_based_neigh_lim",
95 &TConfig::set_smoothing_length_density_based_neigh_lim)
96 .def("set_enable_particle_reordering", &TConfig::set_enable_particle_reordering)
97 .def("set_particle_reordering_step_freq", &TConfig::set_particle_reordering_step_freq)
98 .def("set_show_ghost_zone_graph", &TConfig::set_show_ghost_zone_graph)
99 .def("use_luminosity", &TConfig::use_luminosity)
100 .def("set_save_dt_to_fields", &TConfig::set_save_dt_to_fields)
101 .def("should_save_dt_to_fields", &TConfig::should_save_dt_to_fields)
102 .def("set_eos_isothermal", &TConfig::set_eos_isothermal)
103 .def("set_eos_adiabatic", &TConfig::set_eos_adiabatic)
104 .def("set_eos_polytropic", &TConfig::set_eos_polytropic)
105 .def("set_eos_locally_isothermal", &TConfig::set_eos_locally_isothermal)
106 .def(
107 "set_eos_locally_isothermalLP07",
108 [](TConfig &self, Tscal cs0, Tscal q, Tscal r0) {
109 self.set_eos_locally_isothermalLP07(cs0, q, r0);
110 },
111 py::kw_only(),
112 py::arg("cs0"),
113 py::arg("q"),
114 py::arg("r0"))
115 .def(
116 "set_eos_locally_isothermalFA2014",
117 [](TConfig &self, Tscal h_over_r) {
118 self.set_eos_locally_isothermalFA2014(h_over_r);
119 },
120 py::kw_only(),
121 py::arg("h_over_r"))
122 .def(
123 "set_eos_locally_isothermalFA2014_extended",
124 [](TConfig &self, Tscal cs0, Tscal q, Tscal r0, u32 n_sinks) {
125 self.set_eos_locally_isothermalFA2014_extended(cs0, q, r0, n_sinks);
126 },
127 py::kw_only(),
128 py::arg("cs0"),
129 py::arg("q"),
130 py::arg("r0"),
131 py::arg("n_sinks"))
132 .def(
133 "set_eos_fermi",
134 [](TConfig &self, Tscal mu_e) {
135 self.set_eos_fermi(mu_e);
136 },
137 py::kw_only(),
138 py::arg("mu_e"))
139 .def("set_artif_viscosity_None", &TConfig::set_artif_viscosity_None)
140 .def(
141 "set_artif_viscosity_Constant",
142 [](TConfig &self, Tscal alpha_u, Tscal alpha_AV, Tscal beta_AV) {
143 self.set_artif_viscosity_Constant({alpha_u, alpha_AV, beta_AV});
144 },
145 py::kw_only(),
146 py::arg("alpha_u"),
147 py::arg("alpha_AV"),
148 py::arg("beta_AV"))
149 .def(
150 "set_artif_viscosity_VaryingMM97",
151 [](TConfig &self,
152 Tscal alpha_min,
153 Tscal alpha_max,
154 Tscal sigma_decay,
155 Tscal alpha_u,
156 Tscal beta_AV) {
157 self.set_artif_viscosity_VaryingMM97(
158 {alpha_min, alpha_max, sigma_decay, alpha_u, beta_AV});
159 },
160 py::kw_only(),
161 py::arg("alpha_min"),
162 py::arg("alpha_max"),
163 py::arg("sigma_decay"),
164 py::arg("alpha_u"),
165 py::arg("beta_AV"))
166 .def(
167 "set_artif_viscosity_VaryingCD10",
168 [](TConfig &self,
169 Tscal alpha_min,
170 Tscal alpha_max,
171 Tscal sigma_decay,
172 Tscal alpha_u,
173 Tscal beta_AV) {
174 self.set_artif_viscosity_VaryingCD10(
175 {alpha_min, alpha_max, sigma_decay, alpha_u, beta_AV});
176 },
177 py::kw_only(),
178 py::arg("alpha_min"),
179 py::arg("alpha_max"),
180 py::arg("sigma_decay"),
181 py::arg("alpha_u"),
182 py::arg("beta_AV"))
183 .def(
184 "set_artif_viscosity_ConstantDisc",
185 [](TConfig &self, Tscal alpha_AV, Tscal alpha_u, Tscal beta_AV) {
186 self.set_artif_viscosity_ConstantDisc({alpha_AV, alpha_u, beta_AV});
187 },
188 py::kw_only(),
189 py::arg("alpha_AV"),
190 py::arg("alpha_u"),
191 py::arg("beta_AV"))
192 .def(
193 "set_IdealMHD",
194 [](TConfig &self, Tscal sigma_mhd, Tscal sigma_u) {
195 self.set_IdealMHD({sigma_mhd, sigma_u});
196 },
197 py::kw_only(),
198 py::arg("sigma_mhd"),
199 py::arg("sigma_u"))
200 .def(
201 "set_self_gravity_none",
202 [](TConfig &self) {
203 self.self_grav_config.set_none();
204 })
205 .def(
206 "set_self_gravity_direct",
207 [](TConfig &self, bool reference_mode = false) {
208 self.self_grav_config.set_direct(reference_mode);
209 },
210 py::kw_only(),
211 py::arg("reference_mode") = false)
212 .def(
213 "set_self_gravity_mm",
214 [](TConfig &self, u32 mm_order, f64 opening_angle, u32 reduction_level) {
215 self.self_grav_config.set_mm(mm_order, opening_angle, reduction_level);
216 },
217 py::kw_only(),
218 py::arg("order"),
219 py::arg("opening_angle"),
220 py::arg("reduction_level") = 3)
221 .def(
222 "set_self_gravity_fmm",
223 [](TConfig &self, u32 order, f64 opening_angle, u32 reduction_level) {
224 self.self_grav_config.set_fmm(order, opening_angle, reduction_level);
225 },
226 py::kw_only(),
227 py::arg("order"),
228 py::arg("opening_angle"),
229 py::arg("reduction_level") = 3)
230 .def(
231 "set_self_gravity_sfmm",
232 [](TConfig &self,
233 u32 sfmm_order,
234 f64 opening_angle,
235 bool leaf_lowering,
236 u32 reduction_level) {
237 self.self_grav_config.set_sfmm(
238 sfmm_order, opening_angle, leaf_lowering, reduction_level);
239 },
240 py::kw_only(),
241 py::arg("order"),
242 py::arg("opening_angle"),
243 py::arg("leaf_lowering") = true,
244 py::arg("reduction_level") = 3)
245 .def(
246 "set_softening_plummer",
247 [](TConfig &self, f64 epsilon) {
248 self.self_grav_config.set_softening_plummer(epsilon);
249 },
250 py::kw_only(),
251 py::arg("epsilon"))
252 .def(
253 "set_softening_none",
254 [](TConfig &self) {
255 self.self_grav_config.set_softening_none();
256 })
257 .def("set_boundary_free", &TConfig::set_boundary_free)
258 .def("set_boundary_periodic", &TConfig::set_boundary_periodic)
259 .def("set_boundary_shearing_periodic", &TConfig::set_boundary_shearing_periodic)
260 .def(
261 "set_dust_mode_none",
262 [](TConfig &self) {
263 self.dust_config.set_none();
264 })
265 .def(
266 "set_dust_mode_monofluid_tvi",
267 [](TConfig &self, u32 nvar) {
268 self.dust_config.set_monofluid_tvi(nvar);
269 })
270 .def(
271 "set_dust_mode_monofluid_complete",
272 [](TConfig &self, u32 ndust) {
273 self.dust_config.set_monofluid_complete(ndust);
274 },
275 py::kw_only(),
276 py::arg("ndust"))
277 .def("add_ext_force_point_mass", &TConfig::add_ext_force_point_mass)
278 .def(
279 "add_ext_force_lense_thirring",
280 [](TConfig &self, Tscal central_mass, Tscal Racc, Tscal a_spin, Tvec dir_spin) {
281 self.add_ext_force_lense_thirring(central_mass, Racc, a_spin, dir_spin);
282 },
283 py::kw_only(),
284 py::arg("central_mass"),
285 py::arg("Racc"),
286 py::arg("a_spin"),
287 py::arg("dir_spin"))
288 .def(
289 "add_ext_force_shearing_box",
290 [](TConfig &self, Tscal Omega_0, Tscal eta, Tscal q) {
291 self.add_ext_force_shearing_box(Omega_0, eta, q);
292 },
293 py::kw_only(),
294 py::arg("Omega_0"),
295 py::arg("eta"),
296 py::arg("q"))
297 .def(
298 "add_ext_force_velocity_dissipation",
299 [](TConfig &self, Tscal eta) {
300 self.ext_force_config.add_velocity_dissipation(eta);
301 },
302 py::kw_only(),
303 py::arg("eta"))
304 .def(
305 "add_ext_force_vertical_disc_potential",
306 [](TConfig &self, Tscal central_mass, Tscal R0) {
307 self.ext_force_config.add_vertical_disc_potential(central_mass, R0);
308 },
309 py::kw_only(),
310 py::arg("central_mass"),
311 py::arg("R0"))
312 .def("set_units", &TConfig::set_units)
313 .def(
314 "get_units",
315 [](TConfig &self) {
316 return self.unit_sys;
317 })
318 .def(
319 "set_cfl_cour",
320 [](TConfig &self, Tscal cfl_cour) {
321 self.cfl_config.cfl_cour = cfl_cour;
322 })
323 .def(
324 "set_cfl_force",
325 [](TConfig &self, Tscal cfl_force) {
326 self.cfl_config.cfl_force = cfl_force;
327 })
328 .def(
329 "set_eta_sink",
330 [](TConfig &self, Tscal eta_sink) {
331 self.cfl_config.eta_sink = eta_sink;
332 })
333 .def("set_cfl_multipler", &TConfig::set_cfl_multipler)
334 .def("set_cfl_mult_stiffness", &TConfig::set_cfl_mult_stiffness)
335 .def(
336 "set_particle_mass",
337 [](TConfig &self, Tscal gpart_mass) {
338 self.gpart_mass = gpart_mass;
339 })
340 .def(
341 "add_kill_sphere",
342 [](TConfig &self, const Tvec &center, Tscal radius) {
343 self.particle_killing.add_kill_sphere(center, radius);
344 },
345 py::kw_only(),
346 py::arg("center"),
347 py::arg("radius"));
348
349 std::string sod_tube_analysis_name = name_model + "_AnalysisSodTube";
350 py::class_<TAnalysisSodTube>(m, sod_tube_analysis_name.c_str())
351 .def("compute_L2_dist", [](TAnalysisSodTube &self) -> std::tuple<Tscal, Tvec, Tscal> {
352 auto ret = self.compute_L2_dist();
353 return {ret.rho, ret.v, ret.P};
354 });
355
356 std::string disc_analysis_name = name_model + "_AnalysisDisc";
357 py::class_<TAnalysisDisc>(m, disc_analysis_name.c_str())
358 .def(
359 "collect_data",
360 [](TAnalysisDisc &self, Tscal Rmin, Tscal Rmax, u32 Nbin, ShamrockCtx &ctx) {
361 auto anal = self.compute_analysis(Rmin, Rmax, Nbin, ctx);
362 py::dict dic_out;
363
364 auto radius = anal.radius.copy_to_stdvec();
365 auto counter = anal.counter.copy_to_stdvec();
366 auto Sigma = anal.Sigma.copy_to_stdvec();
367 auto lx = anal.lx.copy_to_stdvec();
368 auto ly = anal.ly.copy_to_stdvec();
369 auto lz = anal.lz.copy_to_stdvec();
370 auto tilt = anal.tilt.copy_to_stdvec();
371 auto twist = anal.twist.copy_to_stdvec();
372 auto psi = anal.psi.copy_to_stdvec();
373 auto Hsq = anal.Hsq.copy_to_stdvec();
374
375 dic_out["radius"] = radius;
376 dic_out["counter"] = counter;
377 dic_out["Sigma"] = Sigma;
378 dic_out["lx"] = lx;
379 dic_out["ly"] = ly;
380 dic_out["lz"] = lz;
381 dic_out["tilt"] = tilt;
382 dic_out["twist"] = twist;
383 dic_out["psi"] = psi;
384 dic_out["Hsq"] = Hsq;
385
386 return dic_out;
387 });
388
389 std::string setup_name = name_model + "_SPHSetup";
390 py::class_<TSPHSetup>(m, setup_name.c_str())
391 .def(
392 "make_generator_lattice_hcp",
393 [](TSPHSetup &self, Tscal dr, Tvec box_min, Tvec box_max, bool discontinuous) {
394 return self.make_generator_lattice_hcp(dr, {box_min, box_max}, discontinuous);
395 },
396 py::arg("dr"),
397 py::arg("box_min"),
398 py::arg("box_max"),
399 py::arg("discontinuous") = true)
400 .def(
401 "make_generator_lattice_cubic",
402 [](TSPHSetup &self, Tscal dr, Tvec box_min, Tvec box_max) {
403 return self.make_generator_lattice_cubic(dr, {box_min, box_max});
404 })
405 .def(
406 "make_generator_disc_mc",
407 [](TSPHSetup &self,
408 Tscal part_mass,
409 Tscal disc_mass,
410 Tscal r_in,
411 Tscal r_out,
412 std::function<Tscal(Tscal)> sigma_profile,
413 std::function<Tscal(Tscal)> H_profile,
414 std::function<Tscal(Tscal)> rot_profile,
415 std::function<Tscal(Tscal)> cs_profile,
416 u64 random_seed,
417 Tscal init_h_factor) {
418 return self.make_generator_disc_mc(
419 part_mass,
420 disc_mass,
421 r_in,
422 r_out,
423 sigma_profile,
424 H_profile,
425 rot_profile,
426 cs_profile,
427 std::mt19937_64(random_seed),
428 init_h_factor);
429 },
430 py::kw_only(),
431 py::arg("part_mass"),
432 py::arg("disc_mass"),
433 py::arg("r_in"),
434 py::arg("r_out"),
435 py::arg("sigma_profile"),
436 py::arg("H_profile"),
437 py::arg("rot_profile"),
438 py::arg("cs_profile"),
439 py::arg("random_seed"),
440 py::arg("init_h_factor") = 0.8)
441 .def(
442 "make_generator_from_context",
443 [](TSPHSetup &self, ShamrockCtx &context_other) {
444 return self.make_generator_from_context(context_other);
445 })
446 .def(
447 "make_combiner_add",
448 [](TSPHSetup &self,
449 shammodels::sph::modules::SetupNodePtr parent1,
450 shammodels::sph::modules::SetupNodePtr parent2) {
451 return self.make_combiner_add(parent1, parent2);
452 })
453 .def(
454 "make_modifier_warp_disc",
455 [](TSPHSetup &self,
456 shammodels::sph::modules::SetupNodePtr parent,
457 Tscal Rwarp,
458 Tscal Hwarp,
459 Tscal inclination,
460 Tscal posangle) {
461 return self.make_modifier_warp_disc(parent, Rwarp, Hwarp, inclination, posangle);
462 },
463 py::kw_only(),
464 py::arg("parent"),
465 py::arg("Rwarp"),
466 py::arg("Hwarp"),
467 py::arg("inclination"),
468 py::arg("posangle") = 0.)
469 .def(
470 "make_modifier_custom_warp",
471 [](TSPHSetup &self,
472 shammodels::sph::modules::SetupNodePtr parent,
473 std::function<Tscal(Tscal)> inc_profile,
474 std::function<Tscal(Tscal)> psi_profile,
475 std::function<Tvec(Tscal)> k_profile) {
476 return self.make_modifier_custom_warp(parent, inc_profile, psi_profile, k_profile);
477 },
478 py::kw_only(),
479 py::arg("parent"),
480 py::arg("inc_profile"),
481 py::arg("psi_profile"),
482 py::arg("k_profile"))
483 .def(
484 "make_modifier_offset",
485 [](TSPHSetup &self,
486 shammodels::sph::modules::SetupNodePtr parent,
487 Tvec offset_postion,
488 Tvec offset_velocity) {
489 return self.make_modifier_add_offset(parent, offset_postion, offset_velocity);
490 },
491 py::kw_only(),
492 py::arg("parent"),
493 py::arg("offset_position"),
494 py::arg("offset_velocity"))
495 .def(
496 "make_modifier_filter",
497 [](TSPHSetup &self,
498 shammodels::sph::modules::SetupNodePtr parent,
499 std::function<bool(Tvec)> filter) {
500 return self.make_modifier_filter(parent, filter);
501 },
502 py::kw_only(),
503 py::arg("parent"),
504 py::arg("filter"))
505 .def(
506 "make_modifier_split_part",
507 [](TSPHSetup &self,
508 shammodels::sph::modules::SetupNodePtr parent,
509 u64 n_split,
510 u64 seed,
511 Tscal h_scaling) {
512 return self.make_modifier_split_part(parent, n_split, seed, h_scaling);
513 },
514 py::kw_only(),
515 py::arg("parent"),
516 py::arg("n_split"),
517 py::arg("seed"),
518 py::arg("h_scaling") = 0.6)
519 .def(
520 "apply_setup",
521 [](TSPHSetup &self,
522 shammodels::sph::modules::SetupNodePtr setup,
523 bool part_reordering,
524 std::optional<u32> gen_step,
525 std::optional<u32> insert_step,
526 std::optional<u64> msg_count_limit,
527 std::optional<u64> msg_size_limit,
528 std::optional<u64> max_msg_size,
529 bool do_setup_log,
530 bool use_new_setup,
531 bool speculative_balancing) {
532 if (use_new_setup) {
533 return self.apply_setup_new(
534 setup,
535 part_reordering,
536 gen_step,
537 insert_step,
538 msg_count_limit,
539 msg_size_limit,
540 max_msg_size,
541 do_setup_log,
542 speculative_balancing);
543 } else {
544 if (bool(gen_step)) {
545 ON_RANK_0(
546 logger::warn_ln(
547 "SPHSetup", "gen_step is ignored when using old setup"));
548 }
549 if (bool(msg_count_limit)) {
550 ON_RANK_0(
551 logger::warn_ln(
552 "SPHSetup", "msg_count_limit is ignored when using old setup"));
553 }
554 if (bool(msg_size_limit)) {
555 ON_RANK_0(
556 logger::warn_ln(
557 "SPHSetup", "msg_size_limit is ignored when using old setup"));
558 }
559 if (bool(max_msg_size)) {
560 ON_RANK_0(
561 logger::warn_ln(
562 "SPHSetup", "max_msg_size is ignored when using old setup"));
563 }
564 if (bool(do_setup_log)) {
565 ON_RANK_0(
566 logger::warn_ln(
567 "SPHSetup", "do_setup_log is ignored when using old setup"));
568 }
569 return self.apply_setup(setup, part_reordering, insert_step);
570 }
571 },
572 py::arg("setup"),
573 py::kw_only(),
574 py::arg("part_reordering") = true,
575 py::arg("gen_step") = std::nullopt,
576 py::arg("insert_step") = std::nullopt,
577 py::arg("msg_count_limit") = std::nullopt,
578 py::arg("rank_comm_size_limit") = std::nullopt,
579 py::arg("max_msg_size") = std::nullopt,
580 py::arg("do_setup_log") = false,
581 py::arg("use_new_setup") = true,
582 py::arg("speculative_balancing") = false);
583
584 py::class_<T>(m, name_model.c_str())
585 .def(py::init([](ShamrockCtx &ctx) {
586 return std::make_unique<T>(ctx);
587 }))
588 .def("init", &T::init)
589 .def("init_scheduler", &T::init_scheduler)
590
591 .def(
592 "evolve_once_override_time",
593 &T::evolve_once_time_expl,
594 py::arg("t_curr"),
595 py::arg("dt_input"))
596 .def("evolve_once", &T::evolve_once)
597 .def(
598 "evolve_until",
599 [](T &self, f64 target_time, i32 niter_max) {
600 return self.evolve_until(target_time, niter_max);
601 },
602 py::arg("target_time"),
603 py::kw_only(),
604 py::arg("niter_max") = -1)
605 .def(
606 "set_dt",
607 [](T &self, f64 dt) {
608 self.solver.solver_config.set_next_dt(dt);
609 })
610 .def("timestep", &T::timestep)
611 .def("set_cfl_cour", &T::set_cfl_cour, py::arg("cfl_cour"))
612 .def("set_cfl_force", &T::set_cfl_force, py::arg("cfl_force"))
613 .def("set_eta_sink", &T::set_eta_sink, py::arg("eta_sink"))
614 .def("set_particle_mass", &T::set_particle_mass, py::arg("gpart_mass"))
615 .def("get_particle_mass", &T::get_particle_mass)
616 .def("rho_h", &T::rho_h)
617 .def("get_hfact", &T::get_hfact)
618 .def(
619 "get_solver_tex",
620 [](T &self) {
621 return shambase::get_check_ref(self.solver.storage.solver_sequence).get_tex();
622 })
623 .def(
624 "get_solver_dot_graph",
625 [](T &self) {
626 return shambase::get_check_ref(self.solver.storage.solver_sequence).get_dot_graph();
627 })
628 .def(
629 "get_box_dim_fcc_3d",
630 [](T &self, f64 dr, u32 xcnt, u32 ycnt, u32 zcnt) {
631 return self.get_box_dim_fcc_3d(dr, xcnt, ycnt, zcnt);
632 })
633 .def(
634 "get_ideal_fcc_box",
635 [](T &self, f64 dr, f64_3 box_min, f64_3 box_max) {
636 return self.get_ideal_fcc_box(dr, {box_min, box_max});
637 })
638 .def(
639 "get_ideal_hcp_box",
640 [](T &self, f64 dr, f64_3 box_min, f64_3 box_max) {
641 return self.get_ideal_hcp_box(dr, {box_min, box_max});
642 })
643 .def(
644 "resize_simulation_box",
645 [](T &self, f64_3 box_min, f64_3 box_max) {
646 return self.resize_simulation_box({box_min, box_max});
647 })
648 .def(
649 "push_particle",
650 [](T &self, std::vector<f64_3> pos, std::vector<f64> hpart, std::vector<f64> upart) {
651 return self.push_particle(pos, hpart, upart);
652 })
653 .def(
654 "push_particle_mhd",
655 [](T &self,
656 std::vector<f64_3> pos,
657 std::vector<f64> hpart,
658 std::vector<f64> upart,
659 std::vector<f64_3> B_on_rho,
660 std::vector<f64> psi_on_ch) {
661 return self.push_particle_mhd(pos, hpart, upart, B_on_rho, psi_on_ch);
662 })
663 .def(
664 "add_cube_fcc_3d",
665 [](T &self, f64 dr, f64_3 box_min, f64_3 box_max) {
666 return self.add_cube_fcc_3d(dr, {box_min, box_max});
667 })
668 .def(
669 "add_cube_hcp_3d",
670 [](T &self, f64 dr, f64_3 box_min, f64_3 box_max) {
671 return self.add_cube_hcp_3d(dr, {box_min, box_max});
672 })
673 .def(
674 "add_cube_hcp_3d_v2",
675 [](T &self, f64 dr, f64_3 box_min, f64_3 box_max) {
676 return self.add_cube_hcp_3d_v2(dr, {box_min, box_max});
677 })
678 .def(
679 "add_disc_3d_keplerian",
680 [](T &self,
681 Tvec center,
682 u32 Npart,
683 Tscal p,
684 Tscal rho_0,
685 Tscal m,
686 Tscal r_in,
687 Tscal r_out,
688 Tscal q,
689 Tscal cmass) {
690 return self.add_cube_disc_3d(center, Npart, p, rho_0, m, r_in, r_out, q, cmass);
691 })
692 .def(
693 "add_disc_3d",
694 [](T &self,
695 Tvec center,
696 Tscal central_mass,
697 u32 Npart,
698 Tscal r_in,
699 Tscal r_out,
700 Tscal disc_mass,
701 Tscal p,
702 Tscal H_r_in,
703 Tscal q) {
704 return self.add_disc_3d(
705 center, central_mass, Npart, r_in, r_out, disc_mass, p, H_r_in, q);
706 })
707 .def(
708 "add_big_disc_3d",
709 [](T &self,
710 Tvec center,
711 Tscal central_mass,
712 u32 Npart,
713 Tscal r_in,
714 Tscal r_out,
715 Tscal disc_mass,
716 Tscal p,
717 Tscal H_r_in,
718 Tscal q,
719 u16 seed) {
720 self.add_big_disc_3d(
721 center,
722 central_mass,
723 Npart,
724 r_in,
725 r_out,
726 disc_mass,
727 p,
728 H_r_in,
729 q,
730 std::mt19937{seed});
731 return disc_mass / Npart;
732 })
733 .def("get_total_part_count", &T::get_total_part_count)
734 .def("total_mass_to_part_mass", &T::total_mass_to_part_mass)
735 .def(
736 "set_value_in_a_box",
737 [](T &self,
738 const std::string &field_name,
739 const std::string &field_type,
740 const pybind11::object &value,
741 f64_3 box_min,
742 f64_3 box_max,
743 u32 ivar) {
744 if (field_type == "f64") {
745 f64 val = value.cast<f64>();
746 self.set_value_in_a_box(field_name, val, {box_min, box_max}, ivar);
747 } else if (field_type == "f64_3") {
748 f64_3 val = value.cast<f64_3>();
749 self.set_value_in_a_box(field_name, val, {box_min, box_max}, ivar);
750 } else {
751 throw shambase::make_except_with_loc<std::invalid_argument>(
752 "unknown field type");
753 }
754 },
755 py::arg("field_name"),
756 py::arg("field_type"),
757 py::arg("value"),
758 py::arg("box_min"),
759 py::arg("box_max"),
760 py::kw_only(),
761 py::arg("ivar") = 0)
762 .def(
763 "set_value_in_sphere",
764 [](T &self,
765 const std::string &field_name,
766 const std::string &field_type,
767 const pybind11::object &value,
768 f64_3 center,
769 f64 radius) {
770 if (field_type == "f64") {
771 f64 val = value.cast<f64>();
772 self.set_value_in_sphere(field_name, val, center, radius);
773 } else if (field_type == "f64_3") {
774 f64_3 val = value.cast<f64_3>();
775 self.set_value_in_sphere(field_name, val, center, radius);
776 } else {
777 throw shambase::make_except_with_loc<std::invalid_argument>(
778 "unknown field type");
779 }
780 })
781 .def(
782 "set_field_value_lambda_f64",
783 [](T &self,
784 std::string field_name,
785 const std::function<f64(Tvec)> pos_to_val,
786 const u32 offset) {
787 return self.template set_field_value_lambda<f64>(
788 std::move(field_name), pos_to_val, offset);
789 },
790 py::arg("field_name"),
791 py::arg("pos_to_val"),
792 py::arg("offset") = 0)
793 .def(
794 "set_field_value_lambda_f64_3",
795 [](T &self,
796 std::string field_name,
797 const std::function<f64_3(Tvec)> pos_to_val,
798 const u32 offset) {
799 return self.template set_field_value_lambda<f64_3>(
800 std::move(field_name), pos_to_val, offset);
801 },
802 py::arg("field_name"),
803 py::arg("pos_to_val"),
804 py::arg("offset") = 0)
805 .def("overwrite_field_value_f64", &T::template overwrite_field_value<f64>)
806 .def("overwrite_field_value_f64_3", &T::template overwrite_field_value<f64_3>)
807 .def("remap_positions", &T::remap_positions)
808 //.def("set_field_value_lambda_f64_3",[](T&self,std::string field_name, const
809 // std::function<f64_3 (Tscal, Tscal , Tscal)> pos_to_val){
810 // self.template set_field_value_lambda<f64_3>(field_name, [=](Tvec v){
811 // return pos_to_val(v.x(), v.y(),v.z());
812 // });
813 //})
814 .def(
815 "add_kernel_value",
816 [](T &self,
817 const std::string &field_name,
818 const std::string &field_type,
819 const pybind11::object &value,
820 f64_3 center,
821 f64 h_ker) {
822 if (field_type == "f64") {
823 f64 val = value.cast<f64>();
824 self.add_kernel_value(field_name, val, center, h_ker);
825 } else if (field_type == "f64_3") {
826 f64_3 val = value.cast<f64_3>();
827 self.add_kernel_value(field_name, val, center, h_ker);
828 } else {
829 throw shambase::make_except_with_loc<std::invalid_argument>(
830 "unknown field type");
831 }
832 })
833 .def(
834 "get_sum",
835 [](T &self, const std::string &field_name, const std::string &field_type) {
836 if (field_type == "f64") {
837 return py::cast(self.template get_sum<f64>(field_name));
838 } else if (field_type == "f64_3") {
839 return py::cast(self.template get_sum<f64_3>(field_name));
840 } else {
841 throw shambase::make_except_with_loc<std::invalid_argument>(
842 "unknown field type");
843 }
844 })
845 .def(
846 "get_closest_part_to",
847 [](T &self, f64_3 pos) -> f64_3 {
848 return self.get_closest_part_to(pos);
849 })
850 .def(
851 "gen_default_config",
852 [](T &self) {
853 return typename T::Solver::Config{};
854 })
855 .def(
856 "get_current_config",
857 [](T &self) {
858 return self.solver.solver_config;
859 })
860 .def("set_solver_config", &T::set_solver_config)
861 .def("add_sink", &T::add_sink)
862 .def(
863 "get_sinks",
864 [](T &self) {
865 py::list list_out;
866
867 if (!self.solver.storage.sinks.is_empty()) {
868 for (auto &sink : self.solver.storage.sinks.get()) {
869 py::dict sink_dic;
870 sink_dic["pos"] = sink.pos;
871 sink_dic["velocity"] = sink.velocity;
872 sink_dic["sph_acceleration"] = sink.sph_acceleration;
873 sink_dic["ext_acceleration"] = sink.ext_acceleration;
874 sink_dic["mass"] = sink.mass;
875 sink_dic["angular_momentum"] = sink.angular_momentum;
876 sink_dic["accretion_radius"] = sink.accretion_radius;
877 list_out.append(sink_dic);
878 }
879 }
880
881 return list_out;
882 })
883 .def(
884 "get_units",
885 [](T &self) {
886 return self.solver.solver_config.unit_sys;
887 })
888 .def(
889 "render_slice",
890 [](T &self,
891 const std::string &name,
892 const std::string &field_type,
893 const std::vector<Tvec> &positions,
894 const std::optional<custom_getter_t> &custom_getter)
895 -> std::variant<std::vector<f64>, std::vector<f64_3>> {
896 if (custom_getter.has_value()) {
897 if (!(name == "custom" && field_type == "f64")) {
898 throw shambase::make_except_with_loc<std::invalid_argument>(
899 "custom_getter only available for name=custom and field_type=f64");
900 }
901 }
902
903 if (field_type == "f64") {
904 modules::CartesianRender<Tvec, f64, SPHKernel> render(
905 self.ctx, self.solver.solver_config, self.solver.storage);
906 return render.compute_slice(name, positions, custom_getter).copy_to_stdvec();
907 }
908
909 if (field_type == "f64_3") {
910 modules::CartesianRender<Tvec, f64_3, SPHKernel> render(
911 self.ctx, self.solver.solver_config, self.solver.storage);
912 return render.compute_slice(name, positions, std::nullopt).copy_to_stdvec();
913 }
914
915 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
916 },
917 py::arg("name"),
918 py::arg("field_type"),
919 py::arg("positions"),
920 py::arg("custom_getter") = std::nullopt)
921 .def(
922 "render_column_integ",
923 [](T &self,
924 const std::string &name,
925 const std::string &field_type,
926 const std::vector<shammath::Ray<Tvec>> &rays,
927 const std::optional<custom_getter_t> &custom_getter)
928 -> std::variant<std::vector<f64>, std::vector<f64_3>> {
929 if (custom_getter.has_value()) {
930 if (!(name == "custom" && field_type == "f64")) {
931 throw shambase::make_except_with_loc<std::invalid_argument>(
932 "custom_getter only available for name=custom and field_type=f64");
933 }
934 }
935
936 if (field_type == "f64") {
937 modules::CartesianRender<Tvec, f64, SPHKernel> render(
938 self.ctx, self.solver.solver_config, self.solver.storage);
939 return render.compute_column_integ(name, rays, custom_getter).copy_to_stdvec();
940 }
941
942 if (field_type == "f64_3") {
943 modules::CartesianRender<Tvec, f64_3, SPHKernel> render(
944 self.ctx, self.solver.solver_config, self.solver.storage);
945 return render.compute_column_integ(name, rays, std::nullopt).copy_to_stdvec();
946 }
947
948 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
949 },
950 py::arg("name"),
951 py::arg("field_type"),
952 py::arg("rays"),
953 py::arg("custom_getter") = std::nullopt)
954 .def(
955 "compute_field",
956 [](T &self,
957 const std::string &name,
958 const std::string &field_type,
959 const std::optional<custom_getter_t> &custom_getter)
960 -> std::variant<
961 shamrock::solvergraph::Field<f64>,
962 shamrock::solvergraph::Field<f64_3>> {
963 if (custom_getter.has_value()) {
964 if (!(name == "custom" && field_type == "f64")) {
965 throw shambase::make_except_with_loc<std::invalid_argument>(
966 "custom_getter only available for name=custom and field_type=f64");
967 }
968 }
969
970 if (field_type == "f64") {
971 modules::RenderFieldGetter<Tvec, f64, SPHKernel> render_field_getter(
972 self.ctx, self.solver.solver_config, self.solver.storage);
973 return render_field_getter.build_field(name, custom_getter);
974 }
975
976 if (field_type == "f64_3") {
977 modules::RenderFieldGetter<Tvec, f64_3, SPHKernel> render_field_getter(
978 self.ctx, self.solver.solver_config, self.solver.storage);
979 return render_field_getter.build_field(name, custom_getter);
980 }
981
982 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
983 },
984 py::arg("name"),
985 py::arg("field_type"),
986 py::arg("custom_getter") = std::nullopt)
987 .def(
988 "render_azymuthal_integ",
989 [](T &self,
990 const std::string &name,
991 const std::string &field_type,
992 const std::vector<shammath::RingRay<Tvec>> &ring_rays,
993 const std::optional<custom_getter_t> &custom_getter)
994 -> std::variant<std::vector<f64>, std::vector<f64_3>> {
995 if (custom_getter.has_value()) {
996 if (!(name == "custom" && field_type == "f64")) {
997 throw shambase::make_except_with_loc<std::invalid_argument>(
998 "custom_getter only available for name=custom and field_type=f64");
999 }
1000 }
1001
1002 if (field_type == "f64") {
1003 modules::CartesianRender<Tvec, f64, SPHKernel> render(
1004 self.ctx, self.solver.solver_config, self.solver.storage);
1005 return render.compute_azymuthal_integ(name, ring_rays, custom_getter)
1006 .copy_to_stdvec();
1007 }
1008
1009 if (field_type == "f64_3") {
1010 modules::CartesianRender<Tvec, f64_3, SPHKernel> render(
1011 self.ctx, self.solver.solver_config, self.solver.storage);
1012 return render.compute_azymuthal_integ(name, ring_rays, std::nullopt)
1013 .copy_to_stdvec();
1014 }
1015
1016 throw shambase::make_except_with_loc<std::runtime_error>("unknown field type");
1017 },
1018 py::arg("name"),
1019 py::arg("field_type"),
1020 py::arg("ring_rays"),
1021 py::arg("custom_getter") = std::nullopt)
1022 .def(
1023 "render_cartesian_slice",
1024 [](T &self,
1025 const std::string &name,
1026 const std::string &field_type,
1027 Tvec center,
1028 Tvec delta_x,
1029 Tvec delta_y,
1030 u32 nx,
1031 u32 ny,
1032 const std::optional<custom_getter_t> &custom_getter)
1033 -> std::variant<py::array_t<Tscal>> {
1034 if (custom_getter.has_value()) {
1035 if (!(name == "custom" && field_type == "f64")) {
1036 throw shambase::make_except_with_loc<std::invalid_argument>(
1037 "custom_getter only available for name=custom and field_type=f64");
1038 }
1039 }
1040
1041 if (field_type == "f64") {
1042 py::array_t<Tscal> ret({ny, nx});
1043
1044 modules::CartesianRender<Tvec, f64, SPHKernel> render(
1045 self.ctx, self.solver.solver_config, self.solver.storage);
1046
1047 std::vector<f64> slice
1048 = render
1049 .compute_slice(name, center, delta_x, delta_y, nx, ny, custom_getter)
1050 .copy_to_stdvec();
1051
1052 for (u32 iy = 0; iy < ny; iy++) {
1053 for (u32 ix = 0; ix < nx; ix++) {
1054 ret.mutable_at(iy, ix) = slice[ix + nx * iy];
1055 }
1056 }
1057
1058 return ret;
1059 }
1060
1061 if (field_type == "f64_3") {
1062 py::array_t<Tscal> ret({ny, nx, 3_u32});
1063
1064 modules::CartesianRender<Tvec, f64_3, SPHKernel> render(
1065 self.ctx, self.solver.solver_config, self.solver.storage);
1066
1067 std::vector<f64_3> slice
1068 = render.compute_slice(name, center, delta_x, delta_y, nx, ny, std::nullopt)
1069 .copy_to_stdvec();
1070
1071 for (u32 iy = 0; iy < ny; iy++) {
1072 for (u32 ix = 0; ix < nx; ix++) {
1073 ret.mutable_at(iy, ix, 0) = slice[ix + nx * iy][0];
1074 ret.mutable_at(iy, ix, 1) = slice[ix + nx * iy][1];
1075 ret.mutable_at(iy, ix, 2) = slice[ix + nx * iy][2];
1076 }
1077 }
1078
1079 return ret;
1080 }
1081
1082 shambase::throw_with_loc<std::runtime_error>("unknown field type");
1083 return py::array_t<Tscal>({nx, ny});
1084 },
1085 py::arg("name"),
1086 py::arg("field_type"),
1087 py::arg("center"),
1088 py::arg("delta_x"),
1089 py::arg("delta_y"),
1090 py::arg("nx"),
1091 py::arg("ny"),
1092 py::arg("custom_getter") = std::nullopt)
1093 .def(
1094 "render_cartesian_column_integ",
1095 [](T &self,
1096 const std::string &name,
1097 const std::string &field_type,
1098 Tvec center,
1099 Tvec delta_x,
1100 Tvec delta_y,
1101 u32 nx,
1102 u32 ny,
1103 const std::optional<custom_getter_t> &custom_getter)
1104 -> std::variant<py::array_t<Tscal>> {
1105 if (custom_getter.has_value()) {
1106 if (!(name == "custom" && field_type == "f64")) {
1107 throw shambase::make_except_with_loc<std::invalid_argument>(
1108 "custom_getter only available for name=custom and field_type=f64");
1109 }
1110 }
1111
1112 if (field_type == "f64") {
1113 py::array_t<Tscal> ret({ny, nx});
1114
1115 modules::CartesianRender<Tvec, f64, SPHKernel> render(
1116 self.ctx, self.solver.solver_config, self.solver.storage);
1117
1118 std::vector<f64> slice
1119 = render
1120 .compute_column_integ(
1121 name, center, delta_x, delta_y, nx, ny, custom_getter)
1122 .copy_to_stdvec();
1123
1124 for (u32 iy = 0; iy < ny; iy++) {
1125 for (u32 ix = 0; ix < nx; ix++) {
1126 ret.mutable_at(iy, ix) = slice[ix + nx * iy];
1127 }
1128 }
1129
1130 return ret;
1131 }
1132
1133 if (field_type == "f64_3") {
1134 py::array_t<Tscal> ret({ny, nx, 3_u32});
1135
1136 modules::CartesianRender<Tvec, f64_3, SPHKernel> render(
1137 self.ctx, self.solver.solver_config, self.solver.storage);
1138
1139 std::vector<f64_3> slice
1140 = render
1141 .compute_column_integ(
1142 name, center, delta_x, delta_y, nx, ny, std::nullopt)
1143 .copy_to_stdvec();
1144
1145 for (u32 iy = 0; iy < ny; iy++) {
1146 for (u32 ix = 0; ix < nx; ix++) {
1147 ret.mutable_at(iy, ix, 0) = slice[ix + nx * iy][0];
1148 ret.mutable_at(iy, ix, 1) = slice[ix + nx * iy][1];
1149 ret.mutable_at(iy, ix, 2) = slice[ix + nx * iy][2];
1150 }
1151 }
1152
1153 return ret;
1154 }
1155
1156 shambase::throw_with_loc<std::runtime_error>("unknown field type");
1157 return py::array_t<Tscal>({nx, ny});
1158 },
1159 py::arg("name"),
1160 py::arg("field_type"),
1161 py::arg("center"),
1162 py::arg("delta_x"),
1163 py::arg("delta_y"),
1164 py::arg("nx"),
1165 py::arg("ny"),
1166 py::arg("custom_getter") = std::nullopt)
1167 .def(
1168 "gen_config_from_phantom_dump",
1169 [](T &self, PhantomDump &dump, bool bypass_error) {
1170 return self.gen_config_from_phantom_dump(dump, bypass_error);
1171 },
1172 py::arg("dump"),
1173 py::arg("bypass_error") = false,
1174 R"==(
1175 This function generate a shamrock sph solver config from a phantom dump
1176
1177 Parameters
1178 ----------
1179 PhantomDump dump
1180 bypass_error = false (default) bypass any error in the config
1181)==")
1182 .def(
1183 "init_from_phantom_dump",
1184 [](T &self, PhantomDump &dump, Tscal hpart_fact_load) {
1185 self.init_from_phantom_dump(dump, hpart_fact_load);
1186 },
1187 py::arg("dump"),
1188 py::arg("hpart_fact_load") = 1.0)
1189 .def(
1190 "make_phantom_dump",
1191 [](T &self) {
1192 return self.make_phantom_dump();
1193 })
1194 .def("do_vtk_dump", &T::do_vtk_dump)
1195 .def("set_debug_dump", &T::set_debug_dump)
1196 .def("solver_logs_last_rate", &T::solver_logs_last_rate)
1197 .def("solver_logs_last_obj_count", &T::solver_logs_last_obj_count)
1198 .def(
1199 "solver_logs_last_system_metrics",
1200 [&](T &self) {
1201 auto system_metrics = self.solver.solve_logs.get_last_system_metrics();
1202 py::dict ret;
1203 ret["duration"] = system_metrics.wall_time;
1204 if (system_metrics.rank_energy_consummed.has_value()) {
1205 ret["rank_energy_consummed"] = system_metrics.rank_energy_consummed.value();
1206 }
1207 if (system_metrics.gpu_energy_consummed.has_value()) {
1208 ret["gpu_energy_consummed"] = system_metrics.gpu_energy_consummed.value();
1209 }
1210 if (system_metrics.cpu_energy_consummed.has_value()) {
1211 ret["cpu_energy_consummed"] = system_metrics.cpu_energy_consummed.value();
1212 }
1213 if (system_metrics.dram_energy_consummed.has_value()) {
1214 ret["dram_energy_consummed"] = system_metrics.dram_energy_consummed.value();
1215 }
1216 return ret;
1217 })
1218 .def("solver_logs_cumulated_step_time", &T::solver_logs_cumulated_step_time)
1219 .def("solver_logs_reset_cumulated_step_time", &T::solver_logs_reset_cumulated_step_time)
1220 .def("solver_logs_step_count", &T::solver_logs_step_count)
1221 .def("solver_logs_reset_step_count", &T::solver_logs_reset_step_count)
1222 .def(
1223 "get_time",
1224 [](T &self) {
1225 return self.solver.solver_config.get_time();
1226 })
1227 .def(
1228 "get_dt",
1229 [](T &self) {
1230 return self.solver.solver_config.get_dt_sph();
1231 })
1232 .def(
1233 "set_time",
1234 [](T &self, Tscal t) {
1235 return self.solver.solver_config.set_time(t);
1236 })
1237 .def(
1238 "set_next_dt",
1239 [](T &self, Tscal dt) {
1240 return self.solver.solver_config.set_next_dt(dt);
1241 })
1242 .def(
1243 "set_cfl_multipler",
1244 [](T &self, Tscal lambda) {
1245 return self.solver.solver_config.set_cfl_multipler(lambda);
1246 },
1247 py::arg("lambda"))
1248 .def(
1249 "set_cfl_mult_stiffness",
1250 [](T &self, Tscal cstiff) {
1251 return self.solver.solver_config.set_cfl_mult_stiffness(cstiff);
1252 },
1253 py::arg("cstiff"))
1254 .def(
1255 "change_htolerance",
1256 [](T &self, Tscal in) {
1257 ON_RANK_0(shamlog_warn_ln(
1258 "SPH",
1259 ".change_htolerance(val) is deprecated,\n"
1260 " -> calling this is replaced internally by "
1261 ".change_htolerances(coarse=val, fine=min(val, 1.1))\n"
1262 " see: "
1263 "https://shamrock-code.github.io/Shamrock/mkdocs/models/sph/"
1264 "smoothing_length_tolerance"););
1265 self.change_htolerances(in, std::min(in, (Tscal) 1.1));
1266 })
1267 .def(
1268 "change_htolerances",
1269 [](T &self, Tscal coarse, Tscal fine) {
1270 self.change_htolerances(coarse, fine);
1271 },
1272 py::kw_only(),
1273 py::arg("coarse"),
1274 py::arg("fine"))
1275 .def(
1276 "make_analysis_sodtube",
1277 [](T &self,
1278 shamphys::SodTube sod,
1279 Tvec direction,
1280 Tscal time_val,
1281 Tscal x_ref,
1282 Tscal x_min,
1283 Tscal x_max) {
1284 return std::make_unique<TAnalysisSodTube>(
1285 self.ctx,
1286 self.solver.solver_config,
1287 self.solver.storage,
1288 sod,
1289 direction,
1290 time_val,
1291 x_ref,
1292 x_min,
1293 x_max);
1294 },
1295 py::arg("sod"),
1296 py::arg("direction"),
1297 py::arg("time_val"),
1298 py::arg("x_ref"),
1299 py::arg("x_min"),
1300 py::arg("x_max"))
1301 .def(
1302 "make_analysis_disc",
1303 [](T &self) {
1304 return std::make_unique<TAnalysisDisc>(
1305 self.ctx, self.solver.solver_config, self.solver.storage);
1306 })
1307 .def("load_from_dump", &T::load_from_dump)
1308 .def("dump", &T::dump)
1309 .def("get_setup", &T::get_setup)
1310 .def(
1311 "get_patch_transform",
1312 [](T &self) {
1313 PatchScheduler &sched = shambase::get_check_ref(self.ctx.sched);
1314 return sched.get_patch_transform<Tvec>();
1315 })
1316 .def("apply_momentum_offset", &T::apply_momentum_offset)
1317 .def("apply_position_offset", &T::apply_position_offset)
1318 .def(
1319 "add_timestep_callback",
1320 [](T &self,
1321 std::optional<std::function<void(void)>> step_begin_callback,
1322 std::optional<std::function<void(void)>> step_end_callback) {
1323 self.solver.timestep_callbacks.push_back(
1324 {std::move(step_begin_callback), std::move(step_end_callback)});
1325 },
1326 py::kw_only(),
1327 py::arg("step_begin") = std::nullopt,
1328 py::arg("step_end") = std::nullopt);
1329}
1330
1331template<class Tvec, template<class> class SPHKernel>
1332void add_analysisBarycenter_instance(py::module &m, const std::string &name_model) {
1333 using namespace shammodels::sph;
1334
1335 using Tscal = shambase::VecComponent<Tvec>;
1336
1337 using T = Model<Tvec, SPHKernel>;
1338
1339 py::class_<modules::AnalysisBarycenter<Tvec, SPHKernel>>(m, name_model.c_str())
1340 .def(py::init([](T &model) {
1341 return std::make_unique<modules::AnalysisBarycenter<Tvec, SPHKernel>>(model);
1342 }))
1343 .def("get_barycenter", [](modules::AnalysisBarycenter<Tvec, SPHKernel> &self) {
1344 auto result = self.get_barycenter();
1345 return py::make_tuple(result.barycenter, result.mass_disc);
1346 });
1347}
1348
1349template<class Tvec, template<class> class SPHKernel>
1350void add_analysisEnergyKinetic_instance(py::module &m, const std::string &name_model) {
1351 using namespace shammodels::sph;
1352
1353 using Tscal = shambase::VecComponent<Tvec>;
1354 using T = Model<Tvec, SPHKernel>;
1355
1356 py::class_<modules::AnalysisEnergyKinetic<Tvec, SPHKernel>>(m, name_model.c_str())
1357 .def(py::init([](T &model) {
1358 return std::make_unique<modules::AnalysisEnergyKinetic<Tvec, SPHKernel>>(model);
1359 }))
1360 .def("get_kinetic_energy", [](modules::AnalysisEnergyKinetic<Tvec, SPHKernel> &self) {
1361 return self.get_kinetic_energy();
1362 });
1363}
1364
1365template<class Tvec, template<class> class SPHKernel>
1366void add_analysisEnergyPotential_instance(py::module &m, const std::string &name_model) {
1367 using namespace shammodels::sph;
1368
1369 using Tscal = shambase::VecComponent<Tvec>;
1370 using T = Model<Tvec, SPHKernel>;
1371
1372 py::class_<modules::AnalysisEnergyPotential<Tvec, SPHKernel>>(m, name_model.c_str())
1373 .def(py::init([](T &model) {
1374 return std::make_unique<modules::AnalysisEnergyPotential<Tvec, SPHKernel>>(model);
1375 }))
1376 .def("get_potential_energy", [](modules::AnalysisEnergyPotential<Tvec, SPHKernel> &self) {
1377 return self.get_potential_energy();
1378 });
1379}
1380
1381template<class Tvec, template<class> class SPHKernel>
1382void add_analysisTotalMomentum_instance(py::module &m, const std::string &name_model) {
1383 using namespace shammodels::sph;
1384
1385 using Tscal = shambase::VecComponent<Tvec>;
1386 using T = Model<Tvec, SPHKernel>;
1387
1388 py::class_<modules::AnalysisTotalMomentum<Tvec, SPHKernel>>(m, name_model.c_str())
1389 .def(py::init([](T &model) {
1390 return std::make_unique<modules::AnalysisTotalMomentum<Tvec, SPHKernel>>(model);
1391 }))
1392 .def("get_total_momentum", [](modules::AnalysisTotalMomentum<Tvec, SPHKernel> &self) {
1393 return self.get_total_momentum();
1394 });
1395}
1396
1397template<class Tvec, template<class> class SPHKernel>
1398void add_analysisAngularMomentum_instance(py::module &m, const std::string &name_model) {
1399 using namespace shammodels::sph;
1400
1401 using Tscal = shambase::VecComponent<Tvec>;
1402 using T = Model<Tvec, SPHKernel>;
1403
1404 py::class_<modules::AnalysisAngularMomentum<Tvec, SPHKernel>>(m, name_model.c_str())
1405 .def(py::init([](T &model) {
1406 return std::make_unique<modules::AnalysisAngularMomentum<Tvec, SPHKernel>>(model);
1407 }))
1408 .def("get_angular_momentum", [](modules::AnalysisAngularMomentum<Tvec, SPHKernel> &self) {
1409 return self.get_angular_momentum();
1410 });
1411}
1412using namespace shammodels::sph;
1413
1414template<class Analysis, typename Tvec, template<class> class SPHKernel>
1415auto analysis_impl(shammodels::sph::Model<Tvec, SPHKernel> &model) -> Analysis {
1416 return Analysis(model);
1417}
1418
1419template<template<class, template<class> class> class Analysis>
1420void register_analysis_impl_for_each_kernel(py::module &msph, const char *name_class) {
1421 using namespace shammodels::sph;
1422
1423 using SPHModel_f64_3_M4 = shammodels::sph::Model<f64_3, shammath::M4>;
1424 using SPHModel_f64_3_M6 = shammodels::sph::Model<f64_3, shammath::M6>;
1425 using SPHModel_f64_3_M8 = shammodels::sph::Model<f64_3, shammath::M8>;
1426
1427 using SPHModel_f64_3_C2 = shammodels::sph::Model<f64_3, shammath::C2>;
1428 using SPHModel_f64_3_C4 = shammodels::sph::Model<f64_3, shammath::C4>;
1429 using SPHModel_f64_3_C6 = shammodels::sph::Model<f64_3, shammath::C6>;
1430
1431 msph.def(
1432 name_class,
1433 [](SPHModel_f64_3_M4 &model) {
1434 return analysis_impl<Analysis<f64_3, shammath::M4>>(model);
1435 },
1436 py::kw_only(),
1437 py::arg("model"));
1438
1439 msph.def(
1440 name_class,
1441 [](SPHModel_f64_3_M6 &model) {
1442 return analysis_impl<Analysis<f64_3, shammath::M6>>(model);
1443 },
1444 py::kw_only(),
1445 py::arg("model"));
1446
1447 msph.def(
1448 name_class,
1449 [](SPHModel_f64_3_M8 &model) {
1450 return analysis_impl<Analysis<f64_3, shammath::M8>>(model);
1451 },
1452 py::kw_only(),
1453 py::arg("model"));
1454
1455 msph.def(
1456 name_class,
1457 [](SPHModel_f64_3_C2 &model) {
1458 return analysis_impl<Analysis<f64_3, shammath::C2>>(model);
1459 },
1460 py::kw_only(),
1461 py::arg("model"));
1462
1463 msph.def(
1464 name_class,
1465 [](SPHModel_f64_3_C4 &model) {
1466 return analysis_impl<Analysis<f64_3, shammath::C4>>(model);
1467 },
1468 py::kw_only(),
1469 py::arg("model"));
1470
1471 msph.def(
1472 name_class,
1473 [](SPHModel_f64_3_C6 &model) {
1474 return analysis_impl<Analysis<f64_3, shammath::C6>>(model);
1475 },
1476 py::kw_only(),
1477 py::arg("model"));
1478}
1479
1480Register_pymod(pysphmodel) {
1481 py::module msph = m.def_submodule("model_sph", "Shamrock sph solver");
1482
1483 using namespace shammodels::sph;
1484
1485 add_instance<f64_3, shammath::M4>(msph, "SPHModel_f64_3_M4_SolverConfig", "SPHModel_f64_3_M4");
1486 add_instance<f64_3, shammath::M6>(msph, "SPHModel_f64_3_M6_SolverConfig", "SPHModel_f64_3_M6");
1487 add_instance<f64_3, shammath::M8>(msph, "SPHModel_f64_3_M8_SolverConfig", "SPHModel_f64_3_M8");
1488
1489 add_instance<f64_3, shammath::C2>(msph, "SPHModel_f64_3_C2_SolverConfig", "SPHModel_f64_3_C2");
1490 add_instance<f64_3, shammath::C4>(msph, "SPHModel_f64_3_C4_SolverConfig", "SPHModel_f64_3_C4");
1491 add_instance<f64_3, shammath::C6>(msph, "SPHModel_f64_3_C6_SolverConfig", "SPHModel_f64_3_C6");
1492
1493 using VariantSPHModelBind = std::variant<
1494 std::unique_ptr<Model<f64_3, shammath::M4>>,
1495 std::unique_ptr<Model<f64_3, shammath::M6>>,
1496 std::unique_ptr<Model<f64_3, shammath::M8>>,
1497 std::unique_ptr<Model<f64_3, shammath::C2>>,
1498 std::unique_ptr<Model<f64_3, shammath::C4>>,
1499 std::unique_ptr<Model<f64_3, shammath::C6>>>;
1500
1501 m.def(
1502 "get_Model_SPH",
1503 [](ShamrockCtx &ctx,
1504 const std::string &vector_type,
1505 const std::string &kernel) -> VariantSPHModelBind {
1506 VariantSPHModelBind ret;
1507
1508 if (vector_type == "f64_3" && kernel == "M4") {
1509 ret = std::make_unique<Model<f64_3, shammath::M4>>(ctx);
1510 } else if (vector_type == "f64_3" && kernel == "M6") {
1511 ret = std::make_unique<Model<f64_3, shammath::M6>>(ctx);
1512 } else if (vector_type == "f64_3" && kernel == "M8") {
1513 ret = std::make_unique<Model<f64_3, shammath::M8>>(ctx);
1514 } else if (vector_type == "f64_3" && kernel == "C2") {
1515 ret = std::make_unique<Model<f64_3, shammath::C2>>(ctx);
1516 } else if (vector_type == "f64_3" && kernel == "C4") {
1517 ret = std::make_unique<Model<f64_3, shammath::C4>>(ctx);
1518 } else if (vector_type == "f64_3" && kernel == "C6") {
1519 ret = std::make_unique<Model<f64_3, shammath::C6>>(ctx);
1520 } else {
1521 throw shambase::make_except_with_loc<std::invalid_argument>(
1522 "unknown combination of representation and kernel");
1523 }
1524
1525 return ret;
1526 },
1527 py::kw_only(),
1528 py::arg("context"),
1529 py::arg("vector_type"),
1530 py::arg("sph_kernel"));
1531
1532 py::class_<
1533 shammodels::sph::modules::ISPHSetupNode,
1534 std::shared_ptr<shammodels::sph::modules::ISPHSetupNode>>(msph, "ISPHSetupNode")
1535 .def("get_dot", [](std::shared_ptr<shammodels::sph::modules::ISPHSetupNode> &self) {
1536 return self->get_dot();
1537 });
1538
1539 py::class_<shammodels::sph::TimestepLog>(msph, "TimestepLog")
1540 .def(py::init<>())
1541 .def_readwrite("rank", &shammodels::sph::TimestepLog::rank)
1542 .def_readwrite("rate", &shammodels::sph::TimestepLog::rate)
1543 .def_readwrite("npart", &shammodels::sph::TimestepLog::npart)
1544 .def_readwrite("tcompute", &shammodels::sph::TimestepLog::tcompute)
1545 .def("rate_sum", &shammodels::sph::TimestepLog::rate_sum)
1546 .def("npart_sum", &shammodels::sph::TimestepLog::npart_sum);
1547
1548 add_analysisBarycenter_instance<f64_3, shammath::M4>(msph, "AnalysisBarycenter_f64_3_M4");
1549 add_analysisBarycenter_instance<f64_3, shammath::M6>(msph, "AnalysisBarycenter_f64_3_M6");
1550 add_analysisBarycenter_instance<f64_3, shammath::M8>(msph, "AnalysisBarycenter_f64_3_M8");
1551
1552 add_analysisBarycenter_instance<f64_3, shammath::C2>(msph, "AnalysisBarycenter_f64_3_C2");
1553 add_analysisBarycenter_instance<f64_3, shammath::C4>(msph, "AnalysisBarycenter_f64_3_C4");
1554 add_analysisBarycenter_instance<f64_3, shammath::C6>(msph, "AnalysisBarycenter_f64_3_C6");
1555
1556 add_analysisEnergyKinetic_instance<f64_3, shammath::M4>(msph, "AnalysisEnergyKinetic_f64_3_M4");
1557 add_analysisEnergyKinetic_instance<f64_3, shammath::M6>(msph, "AnalysisEnergyKinetic_f64_3_M6");
1558 add_analysisEnergyKinetic_instance<f64_3, shammath::M8>(msph, "AnalysisEnergyKinetic_f64_3_M8");
1559
1560 add_analysisEnergyKinetic_instance<f64_3, shammath::C2>(msph, "AnalysisEnergyKinetic_f64_3_C2");
1561 add_analysisEnergyKinetic_instance<f64_3, shammath::C4>(msph, "AnalysisEnergyKinetic_f64_3_C4");
1562 add_analysisEnergyKinetic_instance<f64_3, shammath::C6>(msph, "AnalysisEnergyKinetic_f64_3_C6");
1563
1564 add_analysisEnergyPotential_instance<f64_3, shammath::M4>(
1565 msph, "AnalysisEnergyPotential_f64_3_M4");
1566 add_analysisEnergyPotential_instance<f64_3, shammath::M6>(
1567 msph, "AnalysisEnergyPotential_f64_3_M6");
1568 add_analysisEnergyPotential_instance<f64_3, shammath::M8>(
1569 msph, "AnalysisEnergyPotential_f64_3_M8");
1570
1571 add_analysisEnergyPotential_instance<f64_3, shammath::C2>(
1572 msph, "AnalysisEnergyPotential_f64_3_C2");
1573 add_analysisEnergyPotential_instance<f64_3, shammath::C4>(
1574 msph, "AnalysisEnergyPotential_f64_3_C4");
1575 add_analysisEnergyPotential_instance<f64_3, shammath::C6>(
1576 msph, "AnalysisEnergyPotential_f64_3_C6");
1577
1578 add_analysisTotalMomentum_instance<f64_3, shammath::M4>(msph, "AnalysisTotalMomentum_f64_3_M4");
1579 add_analysisTotalMomentum_instance<f64_3, shammath::M6>(msph, "AnalysisTotalMomentum_f64_3_M6");
1580 add_analysisTotalMomentum_instance<f64_3, shammath::M8>(msph, "AnalysisTotalMomentum_f64_3_M8");
1581
1582 add_analysisTotalMomentum_instance<f64_3, shammath::C2>(msph, "AnalysisTotalMomentum_f64_3_C2");
1583 add_analysisTotalMomentum_instance<f64_3, shammath::C4>(msph, "AnalysisTotalMomentum_f64_3_C4");
1584 add_analysisTotalMomentum_instance<f64_3, shammath::C6>(msph, "AnalysisTotalMomentum_f64_3_C6");
1585
1586 add_analysisAngularMomentum_instance<f64_3, shammath::M4>(
1587 msph, "AnalysisAngularMomentum_f64_3_M4");
1588 add_analysisAngularMomentum_instance<f64_3, shammath::M6>(
1589 msph, "AnalysisAngularMomentum_f64_3_M6");
1590 add_analysisAngularMomentum_instance<f64_3, shammath::M8>(
1591 msph, "AnalysisAngularMomentum_f64_3_M8");
1592
1593 add_analysisAngularMomentum_instance<f64_3, shammath::C2>(
1594 msph, "AnalysisAngularMomentum_f64_3_C2");
1595 add_analysisAngularMomentum_instance<f64_3, shammath::C4>(
1596 msph, "AnalysisAngularMomentum_f64_3_C4");
1597 add_analysisAngularMomentum_instance<f64_3, shammath::C6>(
1598 msph, "AnalysisAngularMomentum_f64_3_C6");
1599
1600 register_analysis_impl_for_each_kernel<modules::AnalysisBarycenter>(msph, "analysisBarycenter");
1601 register_analysis_impl_for_each_kernel<modules::AnalysisEnergyKinetic>(
1602 msph, "analysisEnergyKinetic");
1603 register_analysis_impl_for_each_kernel<modules::AnalysisEnergyPotential>(
1604 msph, "analysisEnergyPotential");
1605 register_analysis_impl_for_each_kernel<modules::AnalysisTotalMomentum>(
1606 msph, "analysisTotalMomentum");
1607 register_analysis_impl_for_each_kernel<modules::AnalysisAngularMomentum>(
1608 msph, "analysisAngularMomentum");
1609}
AnalysisAngularMomentum.hpp
AnalysisAngularMomentum class.
AnalysisBarycenter.hpp
AnalysisBarycenter class with one method AnalysisBarycenter.get_barycenter()
AnalysisEnergyKinetic.hpp
AnalysisEnergyKinetic class with one method AnalysisEnergyKinetic.get_kinetic_energy()
AnalysisEnergyPotential.hpp
AnalysisEnergyPotential class with one method AnalysisEnergyPotential.get_potential_energy()
AnalysisTotalMomentum.hpp
AnalysisTotalMomentum class with one method AnalysisTotalMomentum.get_total_momentum()
PatchScheduler.hpp
MPI scheduler.
f64
double f64
Alias for double.
Definition aliases_float.hpp:20
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
u16
std::uint16_t u16
16 bit unsigned integer
Definition aliases_int.hpp:28
i32
std::int32_t i32
32 bit integer
Definition aliases_int.hpp:23
PatchScheduler
The MPI scheduler.
Definition PatchScheduler.hpp:85
shammodels::sph::Model
The shamrock SPH model.
Definition Model.hpp:55
shammodels::sph::modules::ISPHSetupNode
This class is an interface that all SPH setup nodes must implement. It describe an operation associat...
Definition ISPHSetupNode.hpp:74
exception.hpp
This header file contains utility functions related to exception handling in the code.
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::sph::modules::SetupNodePtr
std::shared_ptr< ISPHSetupNode > SetupNodePtr
Alias for a shared pointer to an ISPHSetupNode.
Definition ISPHSetupNode.hpp:141
shammodels::sph
namespace for the sph model
Definition BasicSPHGhosts.hpp:37
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...
sphkernels.hpp
sph kernels
shammath::Ray
Ray representation for intersection testing.
Definition AABB.hpp:34
shammath::RingRay
Ring ray representation for intersection testing.
Definition AABB.hpp:67
shammodels::sph::PhantomDump
Class representing a Phantom dump file.
Definition PhantomDump.hpp:347
worldInfo.hpp
Functions related to the MPI communicator.
ON_RANK_0
#define ON_RANK_0(x)
Macro to execute code only on rank 0.
Definition worldInfo.hpp:73
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8