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