diff --git a/doc/sphinx/examples/sph/run_phantom_run_emulation.py b/doc/sphinx/examples/sph/run_phantom_run_emulation.py index 89a1bdfcce..96cb2a82fd 100644 --- a/doc/sphinx/examples/sph/run_phantom_run_emulation.py +++ b/doc/sphinx/examples/sph/run_phantom_run_emulation.py @@ -38,4 +38,49 @@ shamrock.utils.url.download_file(input_file_url, input_file_path) shamrock.utils.url.download_file(dump_file_url, dump_file_path) -shamrock.utils.phantom.run_phantom_simulation(dump_folder, "disc") + +def plot_that_rippa(ctx, model, dump_number): + pixel_x = 1920 + pixel_y = 1080 + radius = 5 + center = (0.0, 0.0, 0.0) + + aspect = pixel_x / pixel_y + pic_range = [-radius * aspect, radius * aspect, -radius, radius] + delta_x = (radius * 2 * aspect, 0.0, 0.0) + delta_y = (0.0, radius * 2, 0.0) + + arr_rho = model.render_cartesian_column_integ( + "rho", + "f64", + center=(0.0, 0.0, 0.0), + delta_x=delta_x, + delta_y=delta_y, + nx=pixel_x, + ny=pixel_y, + ) + + import copy + + import matplotlib + import matplotlib.pyplot as plt + + my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat")) # copy the default cmap + my_cmap.set_bad(color="black") + + plt.figure(figsize=(16 / 2, 9 / 2)) + res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9) + + cbar = plt.colorbar(res, extend="both") + cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]") + # or r"$\rho$ [code unit]" for slices + + plt.title("t = {:0.3f} [code unit]".format(model.get_time())) + plt.xlabel("x") + plt.ylabel("z") + plt.show() + + +ctx, model, in_params = shamrock.utils.phantom.run_phantom_simulation( + dump_folder, "disc", callback=plot_that_rippa +) diff --git a/src/pylib/shamrock/utils/phantom/__init__.py b/src/pylib/shamrock/utils/phantom/__init__.py index 79c7e23738..ff50bb74d4 100644 --- a/src/pylib/shamrock/utils/phantom/__init__.py +++ b/src/pylib/shamrock/utils/phantom/__init__.py @@ -138,7 +138,7 @@ def load_simulation(simulation_path, dump_file_name=None, in_file_name=None, do_ return ctx, model, in_params -def run_phantom_simulation(simulation_folder, sim_name): +def run_phantom_simulation(simulation_folder, sim_name, callback=None): """ Run a Phantom simulation in Shamrock. """ @@ -191,6 +191,12 @@ def do_dump(dump_number): ): os.remove(os.path.join(simulation_folder, f"{sim_name}_00000.tmp")) + if callback is not None: + callback(ctx, model, dump_number) + + model.set_next_dt(0) + model.timestep() + do_dump(dump_number) dump_number += 1 @@ -204,7 +210,7 @@ def do_dump(dump_number): next_time = tmax last_step = True - last_step = model.evolve_until(next_time) + model.evolve_until(next_time) do_dump(dump_number) dump_number += 1 diff --git a/src/shambackends/include/shambackends/DeviceBufferRef.hpp b/src/shambackends/include/shambackends/DeviceBufferRef.hpp new file mode 100644 index 0000000000..5fd1d3e65a --- /dev/null +++ b/src/shambackends/include/shambackends/DeviceBufferRef.hpp @@ -0,0 +1,26 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +#pragma once + +/** + * @file DeviceBufferRef.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/DeviceBuffer.hpp" + +namespace sham { + + template + using DeviceBufferRef = std::reference_wrapper>; + +} diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp index 118dadc459..3d493487af 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp @@ -70,6 +70,9 @@ namespace shammodels::gsph { using RTree = shamtree::CompressedLeafBVH; + std::shared_ptr exchange_gz_node; + std::shared_ptr exchange_gz_positions; + /// Particle counts per patch std::shared_ptr> part_counts; std::shared_ptr> part_counts_with_ghost; diff --git a/src/shammodels/gsph/src/Solver.cpp b/src/shammodels/gsph/src/Solver.cpp index 6ff34dbc3f..62629005a5 100644 --- a/src/shammodels/gsph/src/Solver.cpp +++ b/src/shammodels/gsph/src/Solver.cpp @@ -93,6 +93,10 @@ void shammodels::gsph::Solver::init_solver_graph() { storage.soundspeed = std::make_shared>(1, "soundspeed", "c_s"); + storage.exchange_gz_node + = std::make_shared(storage.ghost_layout); + storage.exchange_gz_positions + = std::make_shared(storage.xyzh_ghost_layout); // Initialize gradient fields for MUSCL reconstruction // These are only used when reconstruct_config.is_muscl() == true storage.grad_density @@ -190,6 +194,9 @@ template class Kern> void shammodels::gsph::Solver::merge_position_ghost() { StackEntry stack_loc{}; + std::shared_ptr exchange_gz_node + = std::make_shared(storage.xyzh_ghost_layout); + storage.merged_xyzh.set( storage.ghost_handler.get().build_comm_merge_positions(storage.ghost_patch_cache.get())); @@ -618,7 +625,7 @@ void shammodels::gsph::Solver::init_ghost_layout() { storage.ghost_layout = std::make_shared(); shamrock::patch::PatchDataLayerLayout &ghost_layout - = shambase::get_check_ref(storage.ghost_layout.get()); + = shambase::get_check_ref(storage.ghost_layout); solver_config.set_ghost_layout(ghost_layout); } @@ -746,6 +753,9 @@ void shammodels::gsph::Solver::communicate_merge_ghosts_fields() { } }); + std::shared_ptr exchange_gz_node + = std::make_shared(storage.ghost_layout); + // Communicate ghost data across MPI ranks shambase::DistributedDataShared interf_pdat = ghost_handle.communicate_pdat(ghost_layout_ptr, std::move(pdat_interf)); diff --git a/src/shammodels/sph/src/io/Phantom2Shamrock.cpp b/src/shammodels/sph/src/io/Phantom2Shamrock.cpp index 29e73ab6e5..1d3f7369ff 100644 --- a/src/shammodels/sph/src/io/Phantom2Shamrock.cpp +++ b/src/shammodels/sph/src/io/Phantom2Shamrock.cpp @@ -151,6 +151,20 @@ namespace shammodels::sph { f64 utime = phdump.read_header_float("utime"); f64 umagfd = phdump.read_header_float("umagfd"); + if (udist == 1.0 || umass == 1.0 || utime == 1.0 || umagfd == 3.54491) { + logger::warn_ln("SPH", "phantom dump units are not set, defaulting to SI"); + logger::warn_ln("SPH", "udist =", udist); + logger::warn_ln("SPH", "umass =", umass); + logger::warn_ln("SPH", "utime =", utime); + logger::warn_ln("SPH", "umagfd =", umagfd); + + return shamunits::UnitSystem(); + } + + // convert from CGS to SI + udist /= 100.0; + umass /= 1000.0; + return shamunits::UnitSystem( utime, udist, umass // unit_current = 1 , @@ -169,9 +183,10 @@ namespace shammodels::sph { dump.table_header_f64.add("umass", units->kg_inv); dump.table_header_f64.add("utime", units->s_inv); - f64 umass = units->template to(); + // Back to freakin CGS (worst units system ever, well no ... there is imperial) + f64 umass = units->template to() / 1000.0; f64 utime = units->template to(); - f64 udist = units->template to(); + f64 udist = units->template to() / 100.0; shamunits::Constants ctes{*units}; f64 ccst = ctes.c(); diff --git a/src/shammodels/sph/src/modules/SPHSetup.cpp b/src/shammodels/sph/src/modules/SPHSetup.cpp index 087fee9fc0..d6c2549eb8 100644 --- a/src/shammodels/sph/src/modules/SPHSetup.cpp +++ b/src/shammodels/sph/src/modules/SPHSetup.cpp @@ -444,6 +444,9 @@ void shammodels::sph::modules::SPHSetup::apply_setup_new( has_been_limited = false; using namespace shamrock::patch; + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + sham::DeviceBuffer mask_get_ids_where(0, dev_sched); + // inject in local domains first PatchCoordTransform ptransf = sched.get_sim_box().get_patch_transform(); sched.for_each_local_patchdata([&](const Patch p, PatchDataLayer &pdat) {