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