From 79c323657b1a45e862be19cba2a13b2b55679ae9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?No=C3=A9=20Brucy?= Date: Fri, 20 Feb 2026 17:18:49 +0100 Subject: [PATCH 1/5] First attempt to add analytical gravity --- src/shammodels/ramses/CMakeLists.txt | 1 + .../shammodels/ramses/SolverConfig.hpp | 24 +++++- .../ramses/config/enum_GravityMode.hpp | 10 +++ .../modules/ComputeAnalyticalGravity.hpp | 65 ++++++++++++++++ .../ramses/modules/SolverStorage.hpp | 1 + src/shammodels/ramses/src/Solver.cpp | 31 ++++++-- .../src/modules/ComputeAnalyticalGravity.cpp | 77 +++++++++++++++++++ .../ramses/src/modules/TimeIntegrator.cpp | 24 ++++++ 8 files changed, 225 insertions(+), 8 deletions(-) create mode 100644 src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp create mode 100644 src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp diff --git a/src/shammodels/ramses/CMakeLists.txt b/src/shammodels/ramses/CMakeLists.txt index d9e291cec5..1fc604d97d 100644 --- a/src/shammodels/ramses/CMakeLists.txt +++ b/src/shammodels/ramses/CMakeLists.txt @@ -45,6 +45,7 @@ set(Sources src/modules/FuseGhostLayer.cpp src/modules/render/GridRender.cpp src/modules/ComputeCoordinates.cpp + src/modules/ComputeAnalyticalGravity.cpp src/SolverConfig.cpp ) diff --git a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp index 9e2f1ab285..4ae9877294 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp @@ -80,10 +80,24 @@ namespace shammodels::basegodunov { struct GravityConfig { using Tscal = shambase::VecComponent; GravityMode gravity_mode = NoGravity; - bool analytical_gravity = false; // whether to use an external analytical gravity Tscal tol = 1e-6; inline Tscal get_tolerance() { return tol; } inline bool is_gravity_on() { return gravity_mode != NoGravity; } + Tscal point_mass_GM = 1; + }; + + template + struct AnalyticalGravityConfig { + using Tscal = shambase::VecComponent; + + AnalyticalGravityMode analytical_gravity_mode = POINTMASS; + // parameters for analytical gravity can be added here + Tscal point_mass_GM = 1; + Tvec point_mass_pos = Tvec(0.5, 0.5, 0.5); + Tscal epsilon_softening = 1e-10; + inline bool is_analytical_gravity_on() { + return analytical_gravity_mode != NoAnalyticalGravity; + } }; template @@ -208,9 +222,15 @@ struct shammodels::basegodunov::SolverConfig { auto scal_G = get_constant_G(); return 4 * M_PI * scal_G; } + AnalyticalGravityConfig analytical_gravity_config{}; inline Tscal get_grav_tol() { return gravity_config.get_tolerance(); } inline bool is_gravity_on() { return gravity_config.is_gravity_on(); } - inline bool is_coordinate_field_required() { return gravity_config.analytical_gravity; } + inline bool is_coordinate_field_required() { + return analytical_gravity_config.is_analytical_gravity_on(); + } + inline bool is_analytical_gravity_on() { + return analytical_gravity_config.is_analytical_gravity_on(); + } ////////////////////////////////////////////////////////////////////////////////////////////// // Gravity config (END) diff --git a/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp b/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp index 463714c93b..535a644e39 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/config/enum_GravityMode.hpp @@ -38,4 +38,14 @@ namespace shammodels::basegodunov { {GravityMode::BICGSTAB, "bicgstab"}, {GravityMode::MULTIGRID, "multigrid"}}); + enum AnalyticalGravityMode { + NoAnalyticalGravity = 0, + POINTMASS = 1, // point mass analytical gravity + }; + + SHAMROCK_JSON_SERIALIZE_ENUM( + AnalyticalGravityMode, + {{AnalyticalGravityMode::NoAnalyticalGravity, "no_analytical_gravity"}, + {AnalyticalGravityMode::POINTMASS, "point_mass"}}); + } // namespace shammodels::basegodunov diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp new file mode 100644 index 0000000000..ec5b1ad868 --- /dev/null +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeAnalyticalGravity.hpp @@ -0,0 +1,65 @@ +// -------------------------------------------------------// +// +// 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 ComputeAnalyticalGravity.hpp + * @author Noé Brucy (noe.brucy@ens-lyon.fr) + * @brief Compute a gravitational acceleration at the center of each cell, + using an analytical formula. + */ + +#include "shambackends/vec.hpp" +#include "shammodels/ramses/SolverConfig.hpp" +#include "shamrock/solvergraph/Field.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarsEdge.hpp" + +namespace shammodels::basegodunov::modules { + + template + class NodeComputeAnalyticalGravity : public shamrock::solvergraph::INode { + + using AGConfig = typename shammodels::basegodunov::AnalyticalGravityConfig; + using Tscal = shambase::VecComponent; + + u32 block_size; + AGConfig config; + + public: + NodeComputeAnalyticalGravity(u32 block_size, AGConfig config) + : block_size(block_size), config(config) {} + +#define NODE_COMPUTE_ANALYTICAL_GRAVITY(X_RO, X_RW) \ + /* inputs */ \ + X_RO( \ + shamrock::solvergraph::Indexes, \ + sizes) /* number of blocks per patch for all patches on the current MPI process*/ \ + X_RO( \ + shamrock::solvergraph::IFieldSpan, \ + spans_coordinates) /* center coordinates of each cell */ \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_gravitational_force) + + EXPAND_NODE_EDGES(NODE_COMPUTE_ANALYTICAL_GRAVITY) + +#undef NODE_COMPUTE_ANALYTICAL_GRAVITY + + void _impl_evaluate_internal(); + + void _impl_reset_internal() {}; + + inline virtual std::string _impl_get_label() const { return "ComputeAnalyticalGravity"; }; + + virtual std::string _impl_get_tex() const { return "ComputeAnalyticalGravity"; }; + }; + +} // namespace shammodels::basegodunov::modules diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp index fd7ab07fcf..0ccabc4b47 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp @@ -94,6 +94,7 @@ namespace shammodels::basegodunov { std::shared_ptr> block_cell_sizes; std::shared_ptr> cell0block_aabb_lower; std::shared_ptr> coordinates; + std::shared_ptr> gravitational_force; std::shared_ptr> grad_rho; std::shared_ptr> dx_v; diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index 899604ce67..a46c1645da 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -25,6 +25,7 @@ #include "shammodels/ramses/SolverConfig.hpp" #include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp" #include "shammodels/ramses/modules/BlockNeighToCellNeigh.hpp" +#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp" #include "shammodels/ramses/modules/ComputeCFL.hpp" #include "shammodels/ramses/modules/ComputeCellAABB.hpp" #include "shammodels/ramses/modules/ComputeCoordinates.hpp" @@ -423,6 +424,12 @@ void shammodels::basegodunov::Solver::init_solver_graph() { "level0_amr", "level0_amr"); } + if (solver_config.is_analytical_gravity_on()) { + // will be filled by NodeComputeAnalyticalGravity + storage.gravitational_force = std::make_shared>( + AMRBlock::block_size, "gravitational_force", "\\mathbf{f}_a"); + } + storage.grad_rho = std::make_shared>( AMRBlock::block_size, "grad_rho", "\\nabla \\rho"); storage.dx_v = std::make_shared>( @@ -707,7 +714,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { cfg_bc_to_geom(solver_config.bc_config.get_y()), cfg_bc_to_geom(solver_config.bc_config.get_z())}; - // if outflow we want zero gradient so we skip the vector transformation in TransformGhostLayer + // if outflow we want zero gradient so we skip the vector transformation in + // TransformGhostLayer bool transform_vec_x = solver_config.bc_config.get_x() != BCConfig::GhostType::Outflow; bool transform_vec_y = solver_config.bc_config.get_y() != BCConfig::GhostType::Outflow; bool transform_vec_z = solver_config.bc_config.get_z() != BCConfig::GhostType::Outflow; @@ -792,7 +800,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { // enable this to debug GZ //{ - // auto filename_edge = std::make_shared>( + // auto filename_edge = + // std::make_shared>( // "debug_fuse.vtk", "debug_fuse.vtk"); // filename_edge->data = "debug_fuse.vtk"; // @@ -954,8 +963,8 @@ void shammodels::basegodunov::Solver::init_solver_graph() { std::make_shared(std::move(node_coordinates))); } - if (solver_config.amr_mode.need_level_zero_compute()) { // compute level0 sizes in patch (to be - // enabled later when needed) + if (solver_config.amr_mode.need_level_zero_compute()) { // compute level0 sizes in patch (to + // be enabled later when needed) modules::ComputeLevel0CellSize node_level0_sizes{}; node_level0_sizes.set_edges( graph.get_edge_ptr>>("global_patch_boxes"), @@ -963,6 +972,15 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.level0_size); solver_sequence.push_back( std::make_shared(std::move(node_level0_sizes))); + if (solver_config.is_analytical_gravity_on()) { + modules::NodeComputeAnalyticalGravity node_analytical_gravity{ + AMRBlock::block_size, solver_config.analytical_gravity_config}; + node_analytical_gravity.set_edges( + storage.block_counts, storage.coordinates, storage.gravitational_force); + solver_sequence.push_back( + std::make_shared( + std::move(node_analytical_gravity))); + } } if (solver_config.should_compute_rho_mean()) { @@ -1707,8 +1725,9 @@ void shammodels::basegodunov::Solver::do_debug_vtk_dump(std::str */ /* - std::unique_ptr> dtrho = storage.dtrho.get().rankgather_computefield(sched); - writer.write_field("dtrho", dtrho, num_obj * block_size); + std::unique_ptr> dtrho = + storage.dtrho.get().rankgather_computefield(sched); writer.write_field("dtrho", dtrho, + num_obj * block_size); std::unique_ptr> dtrhov = storage.dtrhov.get().rankgather_computefield(sched); diff --git a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp new file mode 100644 index 0000000000..c1f544b58c --- /dev/null +++ b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp @@ -0,0 +1,77 @@ +// -------------------------------------------------------// +// +// 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 +// +// -------------------------------------------------------// + +/** + * @file ComputeAnalyticalGravity.cpp + * @author Noé Brucy (noe.brucy@ens-lyon.fr) + * @brief Compute a gravitational acceleration at the center of each cell, + using an analytical formula. + */ + +#include "shammodels/ramses/modules/ComputeAnalyticalGravity.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace shammodels::basegodunov::modules { + + template + void NodeComputeAnalyticalGravity::_impl_evaluate_internal() { + using Tscal = shambase::VecComponent; + + auto edges = get_edges(); + + edges.spans_coordinates.check_sizes(edges.sizes.indexes); + edges.spans_gravitational_force.ensure_sizes(edges.sizes.indexes); + + auto &coordinates_spans = edges.spans_coordinates.get_spans(); + auto &gravitational_force_spans = edges.spans_gravitational_force.get_spans(); + + shambase::DistributedData cell_counts + = edges.sizes.indexes.template map([&](u64 id, u32 block_count) { + u32 cell_count = block_count * block_size; + return cell_count; + }); + + static constexpr u32 dim = shambase::VectorProperties::dimension; + + if (config.analytical_gravity_mode == POINTMASS) { + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{coordinates_spans}, + sham::DDMultiRef{gravitational_force_spans}, + cell_counts, + [this]( + u32 i, const Tvec *__restrict cell_coor, Tvec *__restrict gravitational_force) { + u32 block_id = i / block_size; + Tvec cell_pos = cell_coor[i]; + + // compute vector from the point mass to the cell + Tvec r_vec = cell_pos - config.point_mass_pos; + + // compute distance to the point mass + Tscal r = 0; + for (u32 d = 0; d < dim; d++) { + r += r_vec[d] * r_vec[d]; + } + r = sqrt(r); + + // compute gravitational acceleration using softened point mass potential + // add small term to avoid division by zero + Tvec g_vec + = -config.point_mass_GM * r_vec / (r * r * r + config.epsilon_softening); + + gravitational_force[i] = g_vec; + }); + } + } + +} // namespace shammodels::basegodunov::modules + +template class shammodels::basegodunov::modules::NodeComputeAnalyticalGravity; diff --git a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp index 94ac21b428..9493cdbfd8 100644 --- a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp +++ b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp @@ -27,6 +27,8 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e shamrock::solvergraph::Field &cfield_dtrho = shambase::get_check_ref(storage.dtrho); shamrock::solvergraph::Field &cfield_dtrhov = shambase::get_check_ref(storage.dtrhov); shamrock::solvergraph::Field &cfield_dtrhoe = shambase::get_check_ref(storage.dtrhoe); + shamrock::solvergraph::Field &cfield_gravitational_force + = shambase::get_check_ref(storage.gravitational_force); // load layout info PatchDataLayerLayout &pdl = scheduler().pdl_old(); @@ -47,6 +49,8 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e sham::DeviceBuffer &dt_rho_patch = cfield_dtrho.get_buf(id); sham::DeviceBuffer &dt_rhov_patch = cfield_dtrhov.get_buf(id); sham::DeviceBuffer &dt_rhoe_patch = cfield_dtrhoe.get_buf(id); + sham::DeviceBuffer &gravitational_force_patch + = cfield_gravitational_force.get_buf(id); u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; @@ -58,11 +62,30 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e auto acc_dt_rho_patch = dt_rho_patch.get_read_access(depends_list); auto acc_dt_rhov_patch = dt_rhov_patch.get_read_access(depends_list); auto acc_dt_rhoe_patch = dt_rhoe_patch.get_read_access(depends_list); + auto acc_gravitational_force_patch + = gravitational_force_patch.get_read_access(depends_list); auto rho = buf_rho.get_write_access(depends_list); auto rhov = buf_rhov.get_write_access(depends_list); auto rhoe = buf_rhoe.get_write_access(depends_list); + if (solver_config.is_analytical_gravity_on()) { + auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { + shambase::parallel_for( + cgh, cell_count, "add analytical gravity", [=](u32 id_a) { + rhov[id_a] += dt * acc_gravitational_force_patch[id_a] * rho[id_a]; + Tscal force_squared = 0; + for (u32 d = 0; d < shambase::VectorProperties::dimension; d++) { + force_squared += acc_gravitational_force_patch[id_a][d] + * acc_gravitational_force_patch[id_a][d]; + } + + rhoe[id_a] += 0.5 * dt * dt * force_squared * rho[id_a]; + ; + }); + }); + } + auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { shambase::parallel_for(cgh, cell_count, "accumulate fluxes", [=](u32 id_a) { const u32 cell_global_id = (u32) id_a; @@ -76,6 +99,7 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e dt_rho_patch.complete_event_state(e); dt_rhov_patch.complete_event_state(e); dt_rhoe_patch.complete_event_state(e); + gravitational_force_patch.complete_event_state(e); buf_rho.complete_event_state(e); buf_rhov.complete_event_state(e); From d7b49077242460f53ec8608455b9ad095fb5dffc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?No=C3=A9=20Brucy?= Date: Fri, 20 Feb 2026 18:08:56 +0100 Subject: [PATCH 2/5] solve formula issues --- .../shammodels/ramses/SolverConfig.hpp | 3 +- .../src/modules/ComputeAnalyticalGravity.cpp | 15 ++++------ .../ramses/src/modules/TimeIntegrator.cpp | 28 ++++++++++++++----- 3 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp index 4ae9877294..8a7dc4aca6 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp @@ -83,7 +83,6 @@ namespace shammodels::basegodunov { Tscal tol = 1e-6; inline Tscal get_tolerance() { return tol; } inline bool is_gravity_on() { return gravity_mode != NoGravity; } - Tscal point_mass_GM = 1; }; template @@ -94,7 +93,7 @@ namespace shammodels::basegodunov { // parameters for analytical gravity can be added here Tscal point_mass_GM = 1; Tvec point_mass_pos = Tvec(0.5, 0.5, 0.5); - Tscal epsilon_softening = 1e-10; + Tscal epsilon_softening = 0.1; inline bool is_analytical_gravity_on() { return analytical_gravity_mode != NoAnalyticalGravity; } diff --git a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp index c1f544b58c..1bd9504a14 100644 --- a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp +++ b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp @@ -47,7 +47,7 @@ namespace shammodels::basegodunov::modules { sham::DDMultiRef{coordinates_spans}, sham::DDMultiRef{gravitational_force_spans}, cell_counts, - [this]( + [config = this->config, block_size = this->block_size]( u32 i, const Tvec *__restrict cell_coor, Tvec *__restrict gravitational_force) { u32 block_id = i / block_size; Tvec cell_pos = cell_coor[i]; @@ -56,18 +56,15 @@ namespace shammodels::basegodunov::modules { Tvec r_vec = cell_pos - config.point_mass_pos; // compute distance to the point mass - Tscal r = 0; + Tscal r_squared = 0; for (u32 d = 0; d < dim; d++) { - r += r_vec[d] * r_vec[d]; + r_squared += r_vec[d] * r_vec[d]; } - r = sqrt(r); - // compute gravitational acceleration using softened point mass potential - // add small term to avoid division by zero - Tvec g_vec - = -config.point_mass_GM * r_vec / (r * r * r + config.epsilon_softening); + Tscal one_over_r_cube_soft = pow(r_squared + config.epsilon_softening, -1.5f); - gravitational_force[i] = g_vec; + // compute gravitational acceleration using softened point mass potential + gravitational_force[i] = -config.point_mass_GM * r_vec * one_over_r_cube_soft; }); } } diff --git a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp index 9493cdbfd8..f2c0a31d04 100644 --- a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp +++ b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp @@ -70,20 +70,34 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e auto rhoe = buf_rhoe.get_write_access(depends_list); if (solver_config.is_analytical_gravity_on()) { - auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { + auto grav_event = q.submit(depends_list, [&, dt](sycl::handler &cgh) { shambase::parallel_for( cgh, cell_count, "add analytical gravity", [=](u32 id_a) { - rhov[id_a] += dt * acc_gravitational_force_patch[id_a] * rho[id_a]; - Tscal force_squared = 0; + Tvec g_vec = acc_gravitational_force_patch[id_a]; + Tvec rhov_old = rhov[id_a]; + + Tscal ekin_old = 0; for (u32 d = 0; d < shambase::VectorProperties::dimension; d++) { - force_squared += acc_gravitational_force_patch[id_a][d] - * acc_gravitational_force_patch[id_a][d]; + ekin_old += rhov_old[d] * rhov_old[d]; } + ekin_old *= 0.5 / rho[id_a]; + + Tscal e_other = rhoe[id_a] - ekin_old; - rhoe[id_a] += 0.5 * dt * dt * force_squared * rho[id_a]; - ; + // Update momentum + rhov[id_a] += dt * g_vec * rho[id_a]; + + Tscal ekin_new = 0; + for (u32 d = 0; d < shambase::VectorProperties::dimension; d++) { + ekin_new += rhov[id_a][d] * rhov[id_a][d]; + } + ekin_new *= 0.5 / rho[id_a]; + + rhoe[id_a] = e_other + ekin_new; }); }); + + depends_list.add_event(grav_event); } auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { From 1435cf7c2ee19f72c7ec4999bc7d0f41994dab1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?No=C3=A9=20Brucy?= Date: Mon, 23 Feb 2026 11:38:47 +0100 Subject: [PATCH 3/5] add disc example --- exemples/ramses/godunov_disc.py | 209 ++++++++++++++++++++++++++++++++ 1 file changed, 209 insertions(+) create mode 100644 exemples/ramses/godunov_disc.py diff --git a/exemples/ramses/godunov_disc.py b/exemples/ramses/godunov_disc.py new file mode 100644 index 0000000000..a617472910 --- /dev/null +++ b/exemples/ramses/godunov_disc.py @@ -0,0 +1,209 @@ +import os + +import matplotlib +import matplotlib.animation as animation +import matplotlib.pyplot as plt +import numpy as np + +import shamrock + +### Parameters + +## Grid setup +resolution = 128 + +## Initial conditions +h_over_r = 0.1 +contrast_factor = 10 +ref_density = 1 +ref_sound_speed = 1 +cs_expo = 0.5 +point_mass_Gmass = 1 +gravity_softening = 0.1 +radius = 1 +gamma0 = 1.6 + +center = 0.5 + + +## Initial conditions + + +def get_rc_soft(rmin, rmax): + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + rc_soft = np.sqrt(x**2 + y**2 + gravity_softening**2) + return rc_soft + + +def get_cs(rc_soft): + cs = ref_sound_speed * (rc_soft / radius) ** -cs_expo + return cs + + +def get_omega(rs_soft, rc_soft, cs_expo): + + cs = get_cs(rc_soft) + if cs_expo == 0.5: + omega = point_mass_Gmass / (rc_soft * rc_soft * rs_soft) + omega -= (4.0 - cs_expo) * (cs * cs / rc_soft * rc_soft) + omega = max(omega, 0) + omega = np.sqrt(omega) + else: + omega = np.sqrt( + max( + point_mass_Gmass / pow(rs_soft, 3) + - (3.0 - cs_expo) * cs * cs / (rc_soft * rc_soft), + 0.0, + ) + ) + + return omega + + +def rho_map(rmin, rmax): + """ + Initial density + ---- + rmin : tuple + Coordinates of the lower corner of the cell + rmax : tuple + Coordinates of the upper corner of the cell + """ + rc_soft = get_rc_soft(rmin, rmax) + dens = ref_density * (radius / rc_soft) ** ((5 - 2 * cs_expo) / 2.0) + + return dens + + +def rhovel_map(rmin, rmax): + rc_soft = get_rc_soft(rmin, rmax) + rho = rho_map(rmin, rmax) + + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + rs = np.sqrt(x**2 + y**2 + z**2) + rs_soft = np.sqrt(x**2 + y**2 + z**2 + gravity_softening**2) + x_soft = x * (rs_soft / rs) + y_soft = y * (rs_soft / rs) + omega = get_omega(rs_soft, rc_soft, cs_expo) + + vx = -rho * omega * y_soft + vy = rho * omega * x_soft + vz = 0.0 + + return (vx * rho, vy * rho, vz * rho) + + +def rhoe_map(rmin, rmax): + """ + Initial total energy + """ + rho = rho_map(rmin, rmax) + rho_vx, rho_vy, rho_vz = rhovel_map(rmin, rmax) + + rc_soft = get_rc_soft(rmin, rmax) + cs = get_cs(rc_soft) + + eint = rho * cs * cs / (gamma0 - 1) + ekin = 0.5 * (rho_vx**2 + rho_vy**2 + rho_vz**2) / rho + + return eint + ekin + + +def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y): + # Create the cylindrical coordinate grid + x_vals = np.linspace(min_x, max_x, nx) + y_vals = np.linspace(min_y, max_y, ny) + + # Create meshgrid + x_grid, y_grid = np.meshgrid(x_vals, y_vals) + + # Convert to Cartesian coordinates (z = 0 for a disc in the xy-plane) + z_grid = z_val * np.ones_like(x_grid) + + # Flatten and stack to create list of positions + positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()]) + + return [tuple(pos) for pos in positions] + + +nx, ny = resolution, resolution + + +positions = make_cartesian_coords(nx, ny, 0.5, 0, 1 - 1e-6, 0, 1 - 1e-6) + + +def plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name, dpi=200): + ext = metadata["extent"] + + my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap + my_cmap.set_bad(color="black") + + arr_rho_pos = np.array(arr_rho_pos).reshape(nx, ny) + + plt.figure(dpi=dpi) + res = plt.imshow( + arr_rho_pos, + cmap=my_cmap, + origin="lower", + extent=ext, + aspect="auto", + ) + plt.xlabel("x") + plt.ylabel("y") + plt.title(f"t = {metadata['time']:0.3f} [seconds]") + cbar = plt.colorbar(res, extend="both") + cbar.set_label(r"$\rho$ [code unit]") + plt.savefig(os.path.join(".", f"rho_{case_name}_{iplot:04d}.png")) + plt.close() + + +def plot(t, iplot): + metadata = {"extent": [0, 1, 0, 1], "time": t} + arr_rho_pos = model.render_slice("rho", "f64", positions) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk") + + +### Initialisation + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + +cell_size = 2 +base = resolution + +cfg = model.gen_default_config() + +scale_fact = 1 / (cell_size * base) +cfg.set_scale_factor(scale_fact) +cfg.set_eos_gamma(1.4) # Can we have something which is NOT an EOS? + +cfg.set_riemann_solver_hll() +cfg.set_slope_lim_minmod() +cfg.set_face_time_interpolation(True) + +model.set_solver_config(cfg) + +model.init_scheduler(int(1e7), 1) # (crit_split - patches, crit_merge - patches) +model.make_base_grid( + (0, 0, 0), (cell_size, cell_size, cell_size), (base, base, base) +) # What is this doing? + + +model.set_field_value_lambda_f64("rho", rho_map) +model.set_field_value_lambda_f64("rhoetot", rhoe_map) +model.set_field_value_lambda_f64_3("rhovel", rhovel_map) + +dt = 0.01 +for i in range(3): + model.dump_vtk("disk" + str(i) + ".vtk") + plot(dt * i, i) + model.evolve_until(dt) + +plot(dt * (i + 1), i) From 5fb4135df2fab302ec11d8846c2dcbb92e9d4553 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?No=C3=A9=20Brucy?= Date: Mon, 23 Feb 2026 14:15:01 +0100 Subject: [PATCH 4/5] fix merging bug --- src/shammodels/ramses/src/Solver.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index a46c1645da..e138d6b42c 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -972,15 +972,16 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.level0_size); solver_sequence.push_back( std::make_shared(std::move(node_level0_sizes))); - if (solver_config.is_analytical_gravity_on()) { - modules::NodeComputeAnalyticalGravity node_analytical_gravity{ - AMRBlock::block_size, solver_config.analytical_gravity_config}; - node_analytical_gravity.set_edges( - storage.block_counts, storage.coordinates, storage.gravitational_force); - solver_sequence.push_back( - std::make_shared( - std::move(node_analytical_gravity))); - } + } + + if (solver_config.is_analytical_gravity_on()) { + modules::NodeComputeAnalyticalGravity node_analytical_gravity{ + AMRBlock::block_size, solver_config.analytical_gravity_config}; + node_analytical_gravity.set_edges( + storage.block_counts, storage.coordinates, storage.gravitational_force); + solver_sequence.push_back( + std::make_shared( + std::move(node_analytical_gravity))); } if (solver_config.should_compute_rho_mean()) { From dee93901ffe9e3e042decbb8950517bd7eaffb39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?No=C3=A9=20Brucy?= Date: Mon, 23 Feb 2026 22:37:32 +0100 Subject: [PATCH 5/5] Separated the kick, fix IC --- exemples/ramses/godunov_disc.py | 187 +++++++++++------- .../shammodels/ramses/SolverConfig.hpp | 7 +- .../ramses/modules/TimeIntegrator.hpp | 2 +- src/shammodels/ramses/src/Solver.cpp | 6 + src/shammodels/ramses/src/SolverConfig.cpp | 1 + .../src/modules/ComputeAnalyticalGravity.cpp | 3 +- .../ramses/src/modules/TimeIntegrator.cpp | 165 +++++++++------- 7 files changed, 228 insertions(+), 143 deletions(-) diff --git a/exemples/ramses/godunov_disc.py b/exemples/ramses/godunov_disc.py index a617472910..195117b750 100644 --- a/exemples/ramses/godunov_disc.py +++ b/exemples/ramses/godunov_disc.py @@ -4,6 +4,7 @@ import matplotlib.animation as animation import matplotlib.pyplot as plt import numpy as np +from numba import njit import shamrock @@ -13,102 +14,121 @@ resolution = 128 ## Initial conditions -h_over_r = 0.1 -contrast_factor = 10 +contrast_factor = 1000 ref_density = 1 -ref_sound_speed = 1 +ref_sound_speed = 0.2 cs_expo = 0.5 point_mass_Gmass = 1 -gravity_softening = 0.1 -radius = 1 +gravity_softening = 0.05 +outer_radius = 0.25 +inner_radius_factor = 0.1 +inner_iso_z_flaring = 2 gamma0 = 1.6 - +smallr = 1e-4 center = 0.5 +h_over_r = 1 ## Initial conditions -def get_rc_soft(rmin, rmax): - rmin = np.array(rmin) - center - rmax = np.array(rmax) - center - x, y, z = (rmin + rmax) / 2 - rc_soft = np.sqrt(x**2 + y**2 + gravity_softening**2) - return rc_soft +@njit +def get_cs(rc_soft, z): + inner_radius = np.sqrt( + outer_radius + * inner_radius_factor + * (outer_radius * inner_radius_factor + inner_iso_z_flaring * np.abs(z)) + ) + if rc_soft > inner_radius: + cs = ref_sound_speed * (rc_soft / outer_radius) ** (-cs_expo) + else: + cs = ref_sound_speed * (inner_radius / outer_radius) ** (-cs_expo) -def get_cs(rc_soft): - cs = ref_sound_speed * (rc_soft / radius) ** -cs_expo - return cs + return cs, inner_radius -def get_omega(rs_soft, rc_soft, cs_expo): +@njit +def get_omega(rs_soft, rc_soft, cs, inner_radius): - cs = get_cs(rc_soft) - if cs_expo == 0.5: - omega = point_mass_Gmass / (rc_soft * rc_soft * rs_soft) - omega -= (4.0 - cs_expo) * (cs * cs / rc_soft * rc_soft) - omega = max(omega, 0) - omega = np.sqrt(omega) - else: - omega = np.sqrt( - max( - point_mass_Gmass / pow(rs_soft, 3) - - (3.0 - cs_expo) * cs * cs / (rc_soft * rc_soft), - 0.0, - ) + if cs_expo == 0.5 or rc_soft > inner_radius: + omega2 = point_mass_Gmass / (rc_soft * rc_soft * rs_soft) - (4.0 - cs_expo) * (cs * cs) / ( + rc_soft * rc_soft ) + omega = np.sqrt(np.maximum(omega2, 0.0)) + else: + omega2 = point_mass_Gmass / (rs_soft**3) - (3.0 - cs_expo) * (cs * cs) / (rc_soft * rc_soft) + omega = np.sqrt(np.maximum(omega2, 0.0)) return omega +@njit def rho_map(rmin, rmax): - """ - Initial density - ---- - rmin : tuple - Coordinates of the lower corner of the cell - rmax : tuple - Coordinates of the upper corner of the cell - """ - rc_soft = get_rc_soft(rmin, rmax) - dens = ref_density * (radius / rc_soft) ** ((5 - 2 * cs_expo) / 2.0) - return dens + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + rs_soft = np.sqrt(x * x + y * y + z * z + gravity_softening**2) + + cs, _ = get_cs(rc_soft, z) + + # Radial power law + dens = ref_density * (outer_radius / rc_soft) ** ((5.0 - 2.0 * cs_expo) / 2.0) + + # Vertical hydrostatic equilibrium + dens *= np.exp((point_mass_Gmass / (cs * cs)) * (1.0 / rs_soft - 1.0 / rc_soft)) + # Outer taper + if rc_soft > outer_radius or np.abs(z) > 0.5 * outer_radius: + dens /= contrast_factor + return np.maximum(dens, smallr) + + +@njit def rhovel_map(rmin, rmax): - rc_soft = get_rc_soft(rmin, rmax) + rho = rho_map(rmin, rmax) rmin = np.array(rmin) - center rmax = np.array(rmax) - center x, y, z = (rmin + rmax) / 2 - rs = np.sqrt(x**2 + y**2 + z**2) - rs_soft = np.sqrt(x**2 + y**2 + z**2 + gravity_softening**2) + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + rs = np.sqrt(x * x + y * y + z * z) + rs_soft = np.sqrt(x * x + y * y + z * z + gravity_softening**2) + + cs, inner_radius = get_cs(rc_soft, z) + omega = get_omega(rs_soft, rc_soft, cs, inner_radius) + x_soft = x * (rs_soft / rs) y_soft = y * (rs_soft / rs) - omega = get_omega(rs_soft, rc_soft, cs_expo) - vx = -rho * omega * y_soft - vy = rho * omega * x_soft + vx = -omega * y_soft + vy = omega * x_soft vz = 0.0 - return (vx * rho, vy * rho, vz * rho) + return (rho * vx, rho * vy, rho * vz) +@njit def rhoe_map(rmin, rmax): - """ - Initial total energy - """ + rho = rho_map(rmin, rmax) rho_vx, rho_vy, rho_vz = rhovel_map(rmin, rmax) - rc_soft = get_rc_soft(rmin, rmax) - cs = get_cs(rc_soft) + rmin = np.array(rmin) - center + rmax = np.array(rmax) - center + x, y, z = (rmin + rmax) / 2 + + rc_soft = np.sqrt(x * x + y * y + gravity_softening**2) + cs, _ = get_cs(rc_soft, z) - eint = rho * cs * cs / (gamma0 - 1) - ekin = 0.5 * (rho_vx**2 + rho_vy**2 + rho_vz**2) / rho + eint = rho * cs * cs / (gamma0 - 1.0) + ekin = 0.5 * (rho_vx * rho_vx + rho_vy * rho_vy + rho_vz * rho_vz) / np.maximum(rho, smallr) return eint + ekin @@ -130,7 +150,7 @@ def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y): return [tuple(pos) for pos in positions] -nx, ny = resolution, resolution +nx, ny = 128, 128 positions = make_cartesian_coords(nx, ny, 0.5, 0, 1 - 1e-6, 0, 1 - 1e-6) @@ -140,7 +160,7 @@ def plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name, dpi=200): ext = metadata["extent"] my_cmap = matplotlib.colormaps["gist_heat"].copy() # copy the default cmap - my_cmap.set_bad(color="black") + my_cmap.set_bad(color="yellow") arr_rho_pos = np.array(arr_rho_pos).reshape(nx, ny) @@ -151,11 +171,13 @@ def plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name, dpi=200): origin="lower", extent=ext, aspect="auto", + norm=matplotlib.colors.LogNorm(vmin=1e-3, vmax=6), ) plt.xlabel("x") plt.ylabel("y") - plt.title(f"t = {metadata['time']:0.3f} [seconds]") + plt.title(f"t = {metadata['time']:0.5f} [seconds]") cbar = plt.colorbar(res, extend="both") + plt.gca().set_aspect("equal") cbar.set_label(r"$\rho$ [code unit]") plt.savefig(os.path.join(".", f"rho_{case_name}_{iplot:04d}.png")) plt.close() @@ -165,10 +187,26 @@ def plot(t, iplot): metadata = {"extent": [0, 1, 0, 1], "time": t} arr_rho_pos = model.render_slice("rho", "f64", positions) plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk") + pos_XZ = [(x, z, y) for (x, y, z) in positions] + arr_rho_pos = model.render_slice("rho", "f64", pos_XZ) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk_xz") + pos_YZ = [(z, y, x) for (x, y, z) in positions] + arr_rho_pos = model.render_slice("rho", "f64", pos_YZ) + plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, "disk_yz") + + +def plot_force(t, iplot): + metadata = {"extent": [0, 1, 0, 1], "time": t} + force = model.render_slice("gravitational_force", "f64_3", positions) + force = np.array(force).reshape(nx, ny, 3) + force = np.linalg.norm(force, axis=2) + plot_rho_slice_cartesian(metadata, force, iplot, "force") ### Initialisation +print("initialising model...") + ctx = shamrock.Context() ctx.pdata_layout_new() @@ -182,7 +220,7 @@ def plot(t, iplot): scale_fact = 1 / (cell_size * base) cfg.set_scale_factor(scale_fact) -cfg.set_eos_gamma(1.4) # Can we have something which is NOT an EOS? +cfg.set_eos_gamma(1.00001) cfg.set_riemann_solver_hll() cfg.set_slope_lim_minmod() @@ -190,20 +228,33 @@ def plot(t, iplot): model.set_solver_config(cfg) -model.init_scheduler(int(1e7), 1) # (crit_split - patches, crit_merge - patches) -model.make_base_grid( - (0, 0, 0), (cell_size, cell_size, cell_size), (base, base, base) -) # What is this doing? - +model.init_scheduler(int(1e9), 1) # (crit_split - patches, crit_merge - patches) +model.make_base_grid((0, 0, 0), (cell_size, cell_size, cell_size), (base, base, base)) +print("initialising fields...", flush=True) +print("initial density...", flush=True) model.set_field_value_lambda_f64("rho", rho_map) -model.set_field_value_lambda_f64("rhoetot", rhoe_map) +print("initial momentum...", flush=True) model.set_field_value_lambda_f64_3("rhovel", rhovel_map) +print("initial total energy...", flush=True) +model.set_field_value_lambda_f64("rhoetot", rhoe_map) +print("initialisation done, plotting initial conditions...", flush=True) dt = 0.01 -for i in range(3): - model.dump_vtk("disk" + str(i) + ".vtk") - plot(dt * i, i) - model.evolve_until(dt) +for i in range(0, 30): + t = dt * i + model.evolve_until(dt * i) + d = ctx.collect_data() + print( + f"Step {i}, t = {t:0.5f} seconds, min rho = {d['rho'].min():0.5e}, max rho = {d['rho'].max():0.5e}" + ) + print( + f"Step {i}, t = {t:0.5f} seconds, min rhoe = {d['rhoetot'].min():0.5e}, max rhov = {d['rhoetot'].max():0.5e}" + ) + print( + f"Step {i}, t = {t:0.5f} seconds, min rhov = {d['rhovel'].min():0.5e}, max rhov = {d['rhovel'].max():0.5e}" + ) + plot(t, i) + plot(dt * (i + 1), i) diff --git a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp index 8a7dc4aca6..c2d84a17dd 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp @@ -91,9 +91,10 @@ namespace shammodels::basegodunov { AnalyticalGravityMode analytical_gravity_mode = POINTMASS; // parameters for analytical gravity can be added here - Tscal point_mass_GM = 1; - Tvec point_mass_pos = Tvec(0.5, 0.5, 0.5); - Tscal epsilon_softening = 0.1; + Tscal point_mass_GM = 1; + Tvec point_mass_pos = Tvec(0.5, 0.5, 0.5); + Tscal epsilon_softening = 0.05; + Tscal min_rho_for_analytical_gravity = 1e-5; // minimum density to apply analytical gravity inline bool is_analytical_gravity_on() { return analytical_gravity_mode != NoAnalyticalGravity; } diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp index 538e960993..c8649b1962 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/TimeIntegrator.hpp @@ -46,7 +46,7 @@ namespace shammodels::basegodunov::modules { TimeIntegrator(ShamrockCtx &context, Config &solver_config, Storage &storage) : context(context), solver_config(solver_config), storage(storage) {} - void forward_euler(Tscal dt); + void forward_euler(Tscal dt, bool precall = false); private: inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index e138d6b42c..f3417ce1df 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -1547,6 +1547,12 @@ void shammodels::basegodunov::Solver::evolve_once() { shambase::get_check_ref(storage.solver_sequence).evaluate(); } + // apply gravity source term + { + modules::TimeIntegrator dt_integ(context, solver_config, storage); + dt_integ.forward_euler(dt_input, true); + } + // compute dt fields modules::ComputeTimeDerivative dt_compute(context, solver_config, storage); dt_compute.compute_dt_fields(); diff --git a/src/shammodels/ramses/src/SolverConfig.cpp b/src/shammodels/ramses/src/SolverConfig.cpp index 3d9f0fc068..2418d758d9 100644 --- a/src/shammodels/ramses/src/SolverConfig.cpp +++ b/src/shammodels/ramses/src/SolverConfig.cpp @@ -35,6 +35,7 @@ namespace shammodels::basegodunov { pdl.add_field("rho", AMRBlock::block_size); pdl.add_field("rhovel", AMRBlock::block_size); pdl.add_field("rhoetot", AMRBlock::block_size); + pdl.add_field("gravitational_force", AMRBlock::block_size); if (is_dust_on()) { u32 ndust = dust_config.ndust; diff --git a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp index 1bd9504a14..3b6b26d787 100644 --- a/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp +++ b/src/shammodels/ramses/src/modules/ComputeAnalyticalGravity.cpp @@ -61,7 +61,8 @@ namespace shammodels::basegodunov::modules { r_squared += r_vec[d] * r_vec[d]; } - Tscal one_over_r_cube_soft = pow(r_squared + config.epsilon_softening, -1.5f); + Tscal one_over_r_cube_soft = pow( + r_squared + config.epsilon_softening * config.epsilon_softening, -1.5f); // compute gravitational acceleration using softened point mass potential gravitational_force[i] = -config.point_mass_GM * r_vec * one_over_r_cube_soft; diff --git a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp index f2c0a31d04..b0e872b3da 100644 --- a/src/shammodels/ramses/src/modules/TimeIntegrator.cpp +++ b/src/shammodels/ramses/src/modules/TimeIntegrator.cpp @@ -16,7 +16,8 @@ #include "shammodels/ramses/modules/TimeIntegrator.hpp" template -void shammodels::basegodunov::modules::TimeIntegrator::forward_euler(Tscal dt) { +void shammodels::basegodunov::modules::TimeIntegrator::forward_euler( + Tscal dt, bool precall) { StackEntry stack_loc{}; @@ -33,47 +34,51 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e // load layout info PatchDataLayerLayout &pdl = scheduler().pdl_old(); - const u32 icell_min = pdl.get_field_idx("cell_min"); - const u32 icell_max = pdl.get_field_idx("cell_max"); - const u32 irho = pdl.get_field_idx("rho"); - const u32 irhoetot = pdl.get_field_idx("rhoetot"); - const u32 irhovel = pdl.get_field_idx("rhovel"); + const u32 icell_min = pdl.get_field_idx("cell_min"); + const u32 icell_max = pdl.get_field_idx("cell_max"); + const u32 irho = pdl.get_field_idx("rho"); + const u32 irhoetot = pdl.get_field_idx("rhoetot"); + const u32 irhovel = pdl.get_field_idx("rhovel"); + const u32 igravforce = pdl.get_field_idx("gravitational_force"); - scheduler().for_each_patchdata_nonempty( - [&, dt](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) { - shamlog_debug_ln("[AMR Flux]", "forward euler integration patch", p.id_patch); + Tscal min_rho_for_analytical_gravity + = solver_config.analytical_gravity_config.min_rho_for_analytical_gravity; - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - u32 id = p.id_patch; + scheduler().for_each_patchdata_nonempty([&, dt, min_rho_for_analytical_gravity, precall]( + const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + shamlog_debug_ln("[AMR Flux]", "forward euler integration patch", p.id_patch); - sham::DeviceBuffer &dt_rho_patch = cfield_dtrho.get_buf(id); - sham::DeviceBuffer &dt_rhov_patch = cfield_dtrhov.get_buf(id); - sham::DeviceBuffer &dt_rhoe_patch = cfield_dtrhoe.get_buf(id); - sham::DeviceBuffer &gravitational_force_patch - = cfield_gravitational_force.get_buf(id); + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u32 id = p.id_patch; + + u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; - u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; + sham::DeviceBuffer &buf_rho = pdat.get_field_buf_ref(irho); + sham::DeviceBuffer &buf_rhov = pdat.get_field_buf_ref(irhovel); + sham::DeviceBuffer &buf_rhoe = pdat.get_field_buf_ref(irhoetot); - sham::DeviceBuffer &buf_rho = pdat.get_field_buf_ref(irho); - sham::DeviceBuffer &buf_rhov = pdat.get_field_buf_ref(irhovel); - sham::DeviceBuffer &buf_rhoe = pdat.get_field_buf_ref(irhoetot); + sham::EventList depends_list; - sham::EventList depends_list; - auto acc_dt_rho_patch = dt_rho_patch.get_read_access(depends_list); - auto acc_dt_rhov_patch = dt_rhov_patch.get_read_access(depends_list); - auto acc_dt_rhoe_patch = dt_rhoe_patch.get_read_access(depends_list); + auto rho = buf_rho.get_write_access(depends_list); + auto rhov = buf_rhov.get_write_access(depends_list); + auto rhoe = buf_rhoe.get_write_access(depends_list); + + if (precall && solver_config.is_analytical_gravity_on()) { + sham::DeviceBuffer &gravitational_force_patch + = cfield_gravitational_force.get_buf(id); + sham::DeviceBuffer &buf_gravforce = pdat.get_field_buf_ref(igravforce); auto acc_gravitational_force_patch = gravitational_force_patch.get_read_access(depends_list); + auto gravforce = buf_gravforce.get_write_access(depends_list); - auto rho = buf_rho.get_write_access(depends_list); - auto rhov = buf_rhov.get_write_access(depends_list); - auto rhoe = buf_rhoe.get_write_access(depends_list); - - if (solver_config.is_analytical_gravity_on()) { - auto grav_event = q.submit(depends_list, [&, dt](sycl::handler &cgh) { + auto grav_event = q.submit( + depends_list, [&, dt, min_rho_for_analytical_gravity](sycl::handler &cgh) { shambase::parallel_for( cgh, cell_count, "add analytical gravity", [=](u32 id_a) { - Tvec g_vec = acc_gravitational_force_patch[id_a]; + Tvec g_vec = acc_gravitational_force_patch[id_a]; + + rho[id_a] = std::max(min_rho_for_analytical_gravity, rho[id_a]); Tvec rhov_old = rhov[id_a]; Tscal ekin_old = 0; @@ -94,11 +99,27 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e ekin_new *= 0.5 / rho[id_a]; rhoe[id_a] = e_other + ekin_new; + + gravforce[id_a] = g_vec; }); }); - depends_list.add_event(grav_event); - } + gravitational_force_patch.complete_event_state(grav_event); + + buf_rho.complete_event_state(grav_event); + buf_rhov.complete_event_state(grav_event); + buf_rhoe.complete_event_state(grav_event); + buf_gravforce.complete_event_state(grav_event); + } + + if (!precall) { + + sham::DeviceBuffer &dt_rho_patch = cfield_dtrho.get_buf(id); + sham::DeviceBuffer &dt_rhov_patch = cfield_dtrhov.get_buf(id); + sham::DeviceBuffer &dt_rhoe_patch = cfield_dtrhoe.get_buf(id); + auto acc_dt_rho_patch = dt_rho_patch.get_read_access(depends_list); + auto acc_dt_rhov_patch = dt_rhov_patch.get_read_access(depends_list); + auto acc_dt_rhoe_patch = dt_rhoe_patch.get_read_access(depends_list); auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { shambase::parallel_for(cgh, cell_count, "accumulate fluxes", [=](u32 id_a) { @@ -113,60 +134,64 @@ void shammodels::basegodunov::modules::TimeIntegrator::forward_e dt_rho_patch.complete_event_state(e); dt_rhov_patch.complete_event_state(e); dt_rhoe_patch.complete_event_state(e); - gravitational_force_patch.complete_event_state(e); buf_rho.complete_event_state(e); buf_rhov.complete_event_state(e); buf_rhoe.complete_event_state(e); - }); + } + }); - if (solver_config.is_dust_on()) { + if (!precall) { + if (solver_config.is_dust_on()) { - shamrock::solvergraph::Field &cfield_dtrho_dust - = shambase::get_check_ref(storage.dtrho_dust); - shamrock::solvergraph::Field &cfield_dtrhov_dust - = shambase::get_check_ref(storage.dtrhov_dust); + shamrock::solvergraph::Field &cfield_dtrho_dust + = shambase::get_check_ref(storage.dtrho_dust); + shamrock::solvergraph::Field &cfield_dtrhov_dust + = shambase::get_check_ref(storage.dtrhov_dust); - const u32 irho_dust = pdl.get_field_idx("rho_dust"); - const u32 irhovel_dust = pdl.get_field_idx("rhovel_dust"); + const u32 irho_dust = pdl.get_field_idx("rho_dust"); + const u32 irhovel_dust = pdl.get_field_idx("rhovel_dust"); - scheduler().for_each_patchdata_nonempty([&, dt]( - const shamrock::patch::Patch p, - shamrock::patch::PatchDataLayer &pdat) { - shamlog_debug_ln( - "[AMR Flux]", "forward euler integration patch for dust fields", p.id_patch); + scheduler().for_each_patchdata_nonempty([&, dt]( + const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + shamlog_debug_ln( + "[AMR Flux]", "forward euler integration patch for dust fields", p.id_patch); - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - u32 id = p.id_patch; + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u32 id = p.id_patch; - sham::DeviceBuffer &dt_rho_dust_patch = cfield_dtrho_dust.get_buf(id); - sham::DeviceBuffer &dt_rhov_dust_patch = cfield_dtrhov_dust.get_buf(id); + sham::DeviceBuffer &dt_rho_dust_patch = cfield_dtrho_dust.get_buf(id); + sham::DeviceBuffer &dt_rhov_dust_patch = cfield_dtrhov_dust.get_buf(id); - u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; - u32 ndust = solver_config.dust_config.ndust; + u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size; + u32 ndust = solver_config.dust_config.ndust; - sham::DeviceBuffer &buf_rho_dust = pdat.get_field_buf_ref(irho_dust); - sham::DeviceBuffer &buf_rhov_dust = pdat.get_field_buf_ref(irhovel_dust); + sham::DeviceBuffer &buf_rho_dust = pdat.get_field_buf_ref(irho_dust); + sham::DeviceBuffer &buf_rhov_dust + = pdat.get_field_buf_ref(irhovel_dust); - sham::EventList depends_list; - auto acc_dt_rho_dust_patch = dt_rho_dust_patch.get_read_access(depends_list); - auto acc_dt_rhov_dust_patch = dt_rhov_dust_patch.get_read_access(depends_list); + sham::EventList depends_list; + auto acc_dt_rho_dust_patch = dt_rho_dust_patch.get_read_access(depends_list); + auto acc_dt_rhov_dust_patch = dt_rhov_dust_patch.get_read_access(depends_list); - auto rho_dust = buf_rho_dust.get_write_access(depends_list); - auto rhov_dust = buf_rhov_dust.get_write_access(depends_list); + auto rho_dust = buf_rho_dust.get_write_access(depends_list); + auto rhov_dust = buf_rhov_dust.get_write_access(depends_list); - auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { - shambase::parallel_for(cgh, ndust * cell_count, "accumulate fluxes", [=](u32 id_a) { - rho_dust[id_a] += dt * acc_dt_rho_dust_patch[id_a]; - rhov_dust[id_a] += dt * acc_dt_rhov_dust_patch[id_a]; + auto e = q.submit(depends_list, [&, dt](sycl::handler &cgh) { + shambase::parallel_for( + cgh, ndust * cell_count, "accumulate fluxes", [=](u32 id_a) { + rho_dust[id_a] += dt * acc_dt_rho_dust_patch[id_a]; + rhov_dust[id_a] += dt * acc_dt_rhov_dust_patch[id_a]; + }); }); - }); - dt_rho_dust_patch.complete_event_state(e); - dt_rhov_dust_patch.complete_event_state(e); - buf_rho_dust.complete_event_state(e); - buf_rhov_dust.complete_event_state(e); - }); + dt_rho_dust_patch.complete_event_state(e); + dt_rhov_dust_patch.complete_event_state(e); + buf_rho_dust.complete_event_state(e); + buf_rhov_dust.complete_event_state(e); + }); + } } }