From cd3225272e4c2a68f6b3e1e53a247ff6d729b739 Mon Sep 17 00:00:00 2001 From: Leodasce Sewanou Date: Mon, 16 Mar 2026 14:29:08 +0100 Subject: [PATCH 1/5] [Ramses][AMR][Global PR] 2:1 refinement check This commit is the first in a series leading to a PR that updates the AMR algorithms. Specifically, it implements the 2:1 refinement consistency check. --- .../modules/AMRGridRefinementHandler.hpp | 15 +- .../src/modules/AMRGridRefinementHandler.cpp | 176 +++++++++++------- 2 files changed, 119 insertions(+), 72 deletions(-) diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp index bab7476bb7..a0a6938fb9 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp @@ -23,6 +23,9 @@ #include "shamrock/amr/AMRCell.hpp" namespace shammodels::basegodunov::modules { + using namespace shamrock::patch; + using Direction_ = shammodels::basegodunov::modules::Direction; + using AMRGraphLinkiterator = shammodels::basegodunov::modules::AMRGraph::ro_access; template class AMRGridRefinementHandler { @@ -43,6 +46,8 @@ namespace shammodels::basegodunov::modules { using BlockCoord = shamrock::amr::AMRBlockCoord; using OrientedAMRGraph = OrientedAMRGraph; + using TgridUint = typename std::make_unsigned>::type; + ShamrockCtx &context; Config &solver_config; Storage &storage; @@ -69,10 +74,16 @@ namespace shammodels::basegodunov::modules { */ template void gen_refine_block_changes( - shambase::DistributedData> &refine_list, - shambase::DistributedData> &derefine_list, + shambase::DistributedData> &dd_refine_flags, + shambase::DistributedData> &dd_derefine_flags, T &&...args); + /** + * @brief + */ + void enforce_two_to_one_refinement( + shambase::DistributedData> &&dd_refine_flags); + template bool internal_refine_grid(shambase::DistributedData> &&refine_list); diff --git a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp index d43f70b7c9..3bdd7275ed 100644 --- a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp +++ b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp @@ -24,15 +24,10 @@ template template void shammodels::basegodunov::modules::AMRGridRefinementHandler:: gen_refine_block_changes( - shambase::DistributedData> &refine_list, - shambase::DistributedData> &derefine_list, + shambase::DistributedData> &dd_refine_flags, + shambase::DistributedData> &dd_derefine_flags, T &&...args) { - using namespace shamrock::patch; - - u64 tot_refine = 0; - u64 tot_derefine = 0; - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); @@ -79,81 +74,122 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: uacc.finalize(resulting_events, id_patch, cur_p, pdat, args...); } - sham::DeviceBuffer &buf_cell_min = pdat.get_field_buf_ref(0); - sham::DeviceBuffer &buf_cell_max = pdat.get_field_buf_ref(1); - - sham::EventList depends_list; - auto acc_min = buf_cell_min.get_read_access(depends_list); - auto acc_max = buf_cell_max.get_read_access(depends_list); - auto acc_merge_flag = derefine_flags.get_write_access(depends_list); - - // keep only derefine flags on only if the eight cells want to merge and if they can - auto e = q.submit(depends_list, [&](sycl::handler &cgh) { - cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) { - u32 id = gid.get_linear_id(); - - std::array blocks; - bool do_merge = true; - - // This avoid the case where we are in the last block of the buffer to avoid the - // out-of-bound read - if (id + split_count <= obj_cnt) { - bool all_want_to_merge = true; - - for (u32 lid = 0; lid < split_count; lid++) { - blocks[lid] = BlockCoord{acc_min[gid + lid], acc_max[gid + lid]}; - all_want_to_merge = all_want_to_merge && acc_merge_flag[gid + lid]; - } - - do_merge = all_want_to_merge && BlockCoord::are_mergeable(blocks); + dd_refine_flags.add_obj(id_patch, std::move(refine_flags)); + dd_derefine_flags.add_obj(id_patch, std::move(derefine_flags)); + }); +} - } else { - do_merge = false; - } +/** + * @brief check and enforce 2:1 rule for refinement + * @tparam Tvec + * @tparam TgridVec + * @param dd_refine_flags + */ +template +void shammodels::basegodunov::modules::AMRGridRefinementHandler:: + enforce_two_to_one_refinement( + shambase::DistributedData> &&dd_refine_flags) { - acc_merge_flag[gid] = do_merge; + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u64 id_patch = cur_p.id_patch; + sham::DeviceBuffer &patch_refine_flags = dd_refine_flags.get(id_patch); + u32 obj_cnt = pdat.get_obj_cnt(); + + // blocks graph in each direction for the current patch + AMRGraph &block_graph_neighs_xp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::xp) + .get(id_patch); + AMRGraph &block_graph_neighs_xm = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::xm) + .get(id_patch); + AMRGraph &block_graph_neighs_yp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::yp) + .get(id_patch); + AMRGraph &block_graph_neighs_ym = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::ym) + .get(id_patch); + AMRGraph &block_graph_neighs_zp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::zp) + .get(id_patch); + AMRGraph &block_graph_neighs_zm = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::zm) + .get(id_patch); + // get levels in the current patch + sham::DeviceBuffer &buf_amr_block_levels + = shambase::get_check_ref(storage.amr_block_levels).get_buf(id_patch); + + for (u32 pass = 0; pass < 100; pass++) { + sham::EventList depend_list; + AMRGraphLinkiterator block_graph_xp + = block_graph_neighs_xp.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_xm + = block_graph_neighs_xm.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_yp + = block_graph_neighs_yp.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_ym + = block_graph_neighs_ym.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_zp + = block_graph_neighs_zp.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_zm + = block_graph_neighs_zm.get_read_access(depend_list); + auto acc_amr_levels = buf_amr_block_levels.get_read_access(depend_list); + auto acc_ref_flags = patch_refine_flags.get_write_access(depend_list); + + auto e = q.submit(depend_list, [&](sycl::handler &cgh) { + cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) { + u32 block_id = gid.get_linear_id(); + // get refinement flag and amr level of the current block + u32 cur_ref_flag = acc_ref_flags[block_id]; + auto cur_block_level = acc_amr_levels[block_id]; + + auto check_2To1_ref = [&](u32 nid) { + if (nid < obj_cnt) { + // get refinement flag and amr level of the neighborh block + u32 neigh_ref_flag = acc_ref_flags[nid]; + auto neigh_block_level = acc_amr_levels[nid]; + if (cur_ref_flag + && sycl::abs(cur_block_level + 1 - neigh_block_level) > 1) { + sycl::atomic_ref< + u32, + sycl::memory_order::relaxed, + sycl::memory_scope::system> + atomic_flag(acc_ref_flags[nid]); + atomic_flag.store(1); + } + } + }; + block_graph_xp.for_each_object_link(block_id, check_2To1_ref); + block_graph_xm.for_each_object_link(block_id, check_2To1_ref); + block_graph_yp.for_each_object_link(block_id, check_2To1_ref); + block_graph_ym.for_each_object_link(block_id, check_2To1_ref); + block_graph_zp.for_each_object_link(block_id, check_2To1_ref); + block_graph_zm.for_each_object_link(block_id, check_2To1_ref); + }); }); - }); - - buf_cell_min.complete_event_state(e); - buf_cell_max.complete_event_state(e); - derefine_flags.complete_event_state(e); + block_graph_neighs_xp.complete_event_state(e); + block_graph_neighs_xm.complete_event_state(e); + block_graph_neighs_yp.complete_event_state(e); + block_graph_neighs_ym.complete_event_state(e); + block_graph_neighs_zp.complete_event_state(e); + block_graph_neighs_zm.complete_event_state(e); + buf_amr_block_levels.complete_event_state(e); + patch_refine_flags.complete_event_state(e); + } //////////////////////////////////////////////////////////////////////////////// // refinement //////////////////////////////////////////////////////////////////////////////// // perform stream compactions on the refinement flags - auto buf_refine = shamalgs::numeric::stream_compact(dev_sched, refine_flags, obj_cnt); - - shamlog_debug_ln( - "AMRGrid", "patch ", id_patch, "refine block count = ", buf_refine.get_size()); - - tot_refine += buf_refine.get_size(); - - // add the results to the map - refine_list.add_obj(id_patch, std::move(buf_refine)); - - //////////////////////////////////////////////////////////////////////////////// - // derefinement - //////////////////////////////////////////////////////////////////////////////// - - // perform stream compactions on the derefinement flags - auto buf_derefine = shamalgs::numeric::stream_compact(dev_sched, derefine_flags, obj_cnt); - + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + auto dev_buf_ref + = shamalgs::numeric::stream_compact(dev_sched, patch_refine_flags, obj_cnt); shamlog_debug_ln( - "AMRGrid", "patch ", id_patch, "merge block count = ", buf_derefine.get_size()); - - tot_derefine += buf_derefine.get_size(); - - // add the results to the map - derefine_list.add_obj(id_patch, std::move(buf_derefine)); + "AMRGrid", "patch ", id_patch, dev_buf_ref.get_size(), "marked for refinement + 2:1"); }); - - logger::info_ln("AMRGrid", "on this process", tot_refine, "blocks were refined"); - logger::info_ln( - "AMRGrid", "on this process", tot_derefine * split_count, "blocks were derefined"); } + template template bool shammodels::basegodunov::modules::AMRGridRefinementHandler:: From 0d4af23a05cc8b35b939dbf382ed0223d9261669 Mon Sep 17 00:00:00 2001 From: Leodasce Sewanou Date: Mon, 16 Mar 2026 17:36:01 +0100 Subject: [PATCH 2/5] updates --- .../modules/AMRGridRefinementHandler.hpp | 12 +- .../ramses/modules/SolverStorage.hpp | 1 + .../src/modules/AMRGridRefinementHandler.cpp | 600 ++++++++++++++---- 3 files changed, 493 insertions(+), 120 deletions(-) diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp index a0a6938fb9..80acfba8a8 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/AMRGridRefinementHandler.hpp @@ -84,12 +84,20 @@ namespace shammodels::basegodunov::modules { void enforce_two_to_one_refinement( shambase::DistributedData> &&dd_refine_flags); + /** + * @brief + */ + void enforce_two_to_one_derefinement( + shambase::DistributedData> &&dd_derefine_flags, + shambase::DistributedData> &&dd_refine_flags); + template - bool internal_refine_grid(shambase::DistributedData> &&refine_list); + bool internal_refine_grid( + shambase::DistributedData> &&dd_refine_flags); template bool internal_derefine_grid( - shambase::DistributedData> &&derefine_list); + shambase::DistributedData> &&dd_derefine_flags); template void internal_update_refinement(); diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp index 4e78f8b371..3a68e700fb 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/SolverStorage.hpp @@ -91,6 +91,7 @@ namespace shammodels::basegodunov { std::shared_ptr> vel; std::shared_ptr> press; std::shared_ptr> vel_dust; + std::shared_ptr> rho_primitive; std::shared_ptr> block_cell_sizes; std::shared_ptr> cell0block_aabb_lower; diff --git a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp index 3bdd7275ed..beb12297b8 100644 --- a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp +++ b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp @@ -14,9 +14,10 @@ * */ -#include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp" +#include "shambase/memory.hpp" #include "shamalgs/details/algorithm/algorithm.hpp" #include "shamcomm/logs.hpp" +#include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp" #include "shammodels/ramses/modules/AMRSortBlocks.hpp" #include @@ -190,37 +191,275 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: }); } +/** + * @brief check and enforce 2:1 rule for derefinement + * @tparam Tvec + * @tparam TgridVec + * @param dd_derefine_flags + * @param dd_refine_flags + */ + +template +void shammodels::basegodunov::modules::AMRGridRefinementHandler:: + enforce_two_to_one_derefinement( + shambase::DistributedData> &&dd_derefine_flags, + shambase::DistributedData> &&dd_refine_flags) { + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + u64 id_patch = cur_p.id_patch; + + sham::DeviceBuffer &patch_derefine_flag = dd_derefine_flags.get(id_patch); + sham::DeviceBuffer &patch_refine_flag = dd_refine_flags.get(id_patch); + + u32 obj_cnt = pdat.get_obj_cnt(); + + // blocks graph in each direction for the current patch + AMRGraph &block_graph_neighs_xp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::xp) + .get(id_patch); + AMRGraph &block_graph_neighs_xm = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::xm) + .get(id_patch); + AMRGraph &block_graph_neighs_yp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::yp) + .get(id_patch); + AMRGraph &block_graph_neighs_ym = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::ym) + .get(id_patch); + AMRGraph &block_graph_neighs_zp = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::zp) + .get(id_patch); + AMRGraph &block_graph_neighs_zm = shambase::get_check_ref(storage.block_graph_edge) + .get_refs_dir(Direction_::zm) + .get(id_patch); + // get the current buffer of block levels in the current patch + sham::DeviceBuffer &buf_amr_block_levels + = shambase::get_check_ref(storage.amr_block_levels).get_buf(id_patch); + + //////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////// + + auto dev_buf_deref_0 + = shamalgs::numeric::stream_compact(dev_sched, patch_derefine_flag, obj_cnt); + + logger::raw_ln( + " Count block's flag for derefinement [No geometry validity check and no 2:1 check] \t " + ": ", + dev_buf_deref_0.get_size(), + "\n"); + //////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////////// + // keep derefine flags on only if the eight cells want to merge and if they can + sham::DeviceBuffer &buf_cell_min = pdat.get_field_buf_ref(0); + sham::DeviceBuffer &buf_cell_max = pdat.get_field_buf_ref(1); + + sham::EventList depends_list; + auto acc_min = buf_cell_min.get_read_access(depends_list); + auto acc_max = buf_cell_max.get_read_access(depends_list); + auto acc_amr_levels = buf_amr_block_levels.get_read_access(depends_list); + + auto acc_merge_flag = patch_derefine_flag.get_write_access(depends_list); + auto acc_refine_flag = patch_refine_flag.get_read_access(depends_list); + + auto e = q.submit(depends_list, [&](sycl::handler &cgh) { + cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) { + u32 id = gid.get_linear_id(); + + std::array blocks; + bool do_merge = true; + bool all_same_level = true; + + // This avoid the case where we are in the last block of the buffer to + // avoid the out-of-bound read + if (id + split_count <= obj_cnt) { + bool all_want_to_merge = true; + + auto get_coord = [](u32 i) -> std::array { + constexpr u32 NsideBlockPow = 1; + constexpr u32 Nside = 1U << NsideBlockPow; + constexpr u32 side_size = Nside; + constexpr u32 block_size = shambase::pow_constexpr(Nside); + + if constexpr (dim == 3) { + const u32 tmp = i >> NsideBlockPow; + // This line is why derefinement never happens + // return {i % Nside, (tmp) % Nside, (tmp ) >> NsideBlockPow}; + return {(tmp) >> NsideBlockPow, (tmp) % Nside, i % Nside}; + } + }; + + auto get_split + = [=](BlockCoord target_block) -> std::array { + std::array ret; + auto bmin = target_block.bmin; + auto bmax = target_block.bmax; + auto split = bmin + (bmax - bmin) / 2; + std::array szs = {bmin, split, bmax}; + for (u32 i = 0; i < split_count; i++) { + auto [lx, ly, lz] = get_coord(i); + + ret[i].bmin = TgridVec{szs[lx].x(), szs[ly].y(), szs[lz].z()}; + ret[i].bmax + = TgridVec{szs[lx + 1].x(), szs[ly + 1].y(), szs[lz + 1].z()}; + } + + return ret; + }; + + for (u32 b_lid = 0; b_lid < split_count; b_lid++) { + blocks[b_lid] = BlockCoord{acc_min[id + b_lid], acc_max[id + b_lid]}; + all_want_to_merge = all_want_to_merge && acc_merge_flag[id + b_lid]; + all_same_level + = all_same_level && (acc_amr_levels[id] == acc_amr_levels[id + b_lid]); + } + + BlockCoord merged = BlockCoord::get_merge(blocks); + std::array splitted = get_split(merged); + for (u32 lid = 0; lid < split_count; lid++) { + do_merge = do_merge && sham::equals(blocks[lid].bmin, splitted[lid].bmin) + && sham::equals(blocks[lid].bmax, splitted[lid].bmax); + } + + do_merge = do_merge && all_want_to_merge && all_same_level; + if (acc_refine_flag[id] && do_merge) { + do_merge = false; + } + + } else { + do_merge = false; + } + acc_merge_flag[id] = do_merge; + }); + }); + buf_cell_min.complete_event_state(e); + buf_cell_max.complete_event_state(e); + buf_amr_block_levels.complete_event_state(e); + patch_derefine_flag.complete_event_state(e); + patch_refine_flag.complete_event_state(e); + + /////////////////////////////////////////////////// + ////////////////////////////////////////////////// + auto buf_derefine_1 + = shamalgs::numeric::stream_compact(dev_sched, patch_derefine_flag, obj_cnt); + logger::raw_ln( + " Count block's flag for derefinement [After geometry validity check and before 2:1 " + "check] " + "\t : ", + buf_derefine_1.get_size(), + "\n"); + ///////////////////////////////////////////////// + //////////////////////////////////////////////// + + // //////////////////////////////////////////////////////////////////////////////////// + // // // enforce 2:1 at parent level + // /////////////////////////////////////////////////////////////////////////////////// + + for (int it = 0; it < 3; it++) { + + sham::EventList depend_list; + AMRGraphLinkiterator block_graph_xp + = block_graph_neighs_xp.get_read_access(depend_list); + + AMRGraphLinkiterator block_graph_xm + = block_graph_neighs_xm.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_yp + = block_graph_neighs_yp.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_ym + = block_graph_neighs_ym.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_zp + = block_graph_neighs_zp.get_read_access(depend_list); + AMRGraphLinkiterator block_graph_zm + = block_graph_neighs_zm.get_read_access(depend_list); + + auto acc_amr_levels = buf_amr_block_levels.get_read_access(depend_list); + auto acc_deref_flag = patch_derefine_flag.get_write_access(depend_list); + auto acc_ref_flag = patch_refine_flag.get_read_access(depend_list); + + auto e_2to1 = q.submit(depend_list, [&](sycl::handler &cgh) { + cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) { + auto lid = gid.get_linear_id(); + + auto check_2To1_der = [&](u32 nid) { + if (nid < obj_cnt && (acc_amr_levels[lid] < acc_amr_levels[nid])) { + acc_deref_flag[lid] = 0; + } + if (nid < obj_cnt && (acc_amr_levels[lid] == acc_amr_levels[nid]) + && (acc_ref_flag[nid] == 1)) { + acc_deref_flag[lid] = 0; + } + }; + + if (acc_deref_flag[lid]) { + + for (u32 i = 0; i < AMRBlock::block_size; i++) { + block_graph_xp.for_each_object_link((lid + i), check_2To1_der); + block_graph_xm.for_each_object_link((lid + i), check_2To1_der); + block_graph_yp.for_each_object_link((lid + i), check_2To1_der); + block_graph_ym.for_each_object_link((lid + i), check_2To1_der); + block_graph_zp.for_each_object_link((lid + i), check_2To1_der); + block_graph_zm.for_each_object_link((lid + i), check_2To1_der); + } + } + }); + }); + block_graph_neighs_xp.complete_event_state(e_2to1); + block_graph_neighs_xm.complete_event_state(e_2to1); + block_graph_neighs_yp.complete_event_state(e_2to1); + block_graph_neighs_ym.complete_event_state(e_2to1); + block_graph_neighs_zp.complete_event_state(e_2to1); + block_graph_neighs_zm.complete_event_state(e_2to1); + buf_amr_block_levels.complete_event_state(e_2to1); + patch_derefine_flag.complete_event_state(e_2to1); + patch_refine_flag.complete_event_state(e_2to1); + } + // //////////////////////////////////////////////////////////////////////////////// + // // derefinement + // //////////////////////////////////////////////////////////////////////////////// + // perform stream compactions on the derefinement flags + auto buf_derefine + = shamalgs::numeric::stream_compact(dev_sched, patch_derefine_flag, obj_cnt); + + logger::raw_ln( + " Count block's flag for derefinement [After geometry validity check and after 2:1 " + "check] \t : ", + buf_derefine.get_size(), + "\n"); + + shamlog_debug_ln( + "AMRGrid", "patch ", id_patch, buf_derefine.get_size(), "marked for derefinement "); + }); +} + template template bool shammodels::basegodunov::modules::AMRGridRefinementHandler:: - internal_refine_grid(shambase::DistributedData> &&refine_list) { - - using namespace shamrock::patch; + internal_refine_grid(shambase::DistributedData> &&dd_refine_flags) { u64 sum_block_count = 0; bool new_cell_were_added = false; + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); scheduler().for_each_patch_data([&](u64 id_patch, Patch cur_p, PatchDataLayer &pdat) { - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - - u32 old_obj_cnt = pdat.get_obj_cnt(); - - sham::DeviceBuffer &refine_flags = refine_list.get(id_patch); - - if (refine_flags.get_size() > 0) { + u32 old_obj_cnt = pdat.get_obj_cnt(); + auto stream_compaction_result = shamalgs::numeric::stream_compact( + dev_sched, dd_refine_flags.get(id_patch), old_obj_cnt); + if (stream_compaction_result.get_size() > 0) { // alloc memory for the new blocks to be created - pdat.expand(refine_flags.get_size() * (split_count - 1)); - + pdat.expand(static_cast(stream_compaction_result.get_size()) * (split_count - 1)); sham::DeviceBuffer &buf_cell_min = pdat.get_field_buf_ref(0); sham::DeviceBuffer &buf_cell_max = pdat.get_field_buf_ref(1); - sham::EventList depends_list; + auto block_bound_low = buf_cell_min.get_write_access(depends_list); auto block_bound_high = buf_cell_max.get_write_access(depends_list); UserAcc uacc(depends_list, pdat); - auto index_to_ref = refine_flags.get_read_access(depends_list); + auto index_to_ref = stream_compaction_result.get_read_access(depends_list); // Refine the block (set the positions) and fill the corresponding fields auto e = q.submit(depends_list, [&](sycl::handler &cgh) { @@ -228,54 +467,53 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: constexpr u32 new_splits = split_count - 1; - cgh.parallel_for(sycl::range<1>(refine_flags.get_size()), [=](sycl::item<1> gid) { - u32 tid = gid.get_linear_id(); + cgh.parallel_for( + sycl::range<1>(stream_compaction_result.get_size()), [=](sycl::item<1> gid) { + u32 tid = gid.get_linear_id(); - u32 idx_to_refine = index_to_ref[tid]; + u32 idx_to_refine = index_to_ref[gid]; - // gen splits coordinates - BlockCoord cur_block{ - block_bound_low[idx_to_refine], block_bound_high[idx_to_refine]}; + // gen splits coordinates + BlockCoord cur_block{ + block_bound_low[idx_to_refine], block_bound_high[idx_to_refine]}; - std::array block_coords - = BlockCoord::get_split(cur_block.bmin, cur_block.bmax); + std::array block_coords + = BlockCoord::get_split(cur_block.bmin, cur_block.bmax); - // generate index for the refined blocks - std::array blocks_ids; - blocks_ids[0] = idx_to_refine; + // generate index for the refined blocks + std::array blocks_ids; + blocks_ids[0] = idx_to_refine; // generate index for the new blocks (the current index is reused for the first // new block, the others are pushed at the end of the patchdata) #pragma unroll - for (u32 pid = 0; pid < new_splits; pid++) { - blocks_ids[pid + 1] = start_index_push + tid * new_splits + pid; - } + for (u32 pid = 0; pid < new_splits; pid++) { + blocks_ids[pid + 1] = start_index_push + tid * new_splits + pid; + } // write coordinates #pragma unroll - for (u32 pid = 0; pid < split_count; pid++) { - block_bound_low[blocks_ids[pid]] = block_coords[pid].bmin; - block_bound_high[blocks_ids[pid]] = block_coords[pid].bmax; - } + for (u32 pid = 0; pid < split_count; pid++) { + block_bound_low[blocks_ids[pid]] = block_coords[pid].bmin; + block_bound_high[blocks_ids[pid]] = block_coords[pid].bmax; + } - // user lambda to fill the fields - uacc.apply_refine(idx_to_refine, cur_block, blocks_ids, block_coords, uacc); - }); + // user lambda to fill the fields + uacc.apply_refine(idx_to_refine, cur_block, blocks_ids, block_coords, uacc); + }); }); - sham::EventList resulting_events{e}; + sham::EventList resulting_events; buf_cell_min.complete_event_state(resulting_events); buf_cell_max.complete_event_state(resulting_events); - + stream_compaction_result.complete_event_state(e); uacc.finalize(resulting_events, pdat); - - refine_flags.complete_event_state(resulting_events); } - + shamlog_debug_ln("AMRGrid", "patch ", id_patch, "new block count = ", pdat.get_obj_cnt()); sum_block_count += pdat.get_obj_cnt(); - new_cell_were_added = new_cell_were_added || refine_flags.get_size() > 0; + new_cell_were_added = new_cell_were_added || (stream_compaction_result.get_size() > 0); }); logger::info_ln("AMRGrid", "process block count =", sum_block_count); @@ -286,7 +524,7 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: template template bool shammodels::basegodunov::modules::AMRGridRefinementHandler:: - internal_derefine_grid(shambase::DistributedData> &&derefine_list) { + internal_derefine_grid(shambase::DistributedData> &&dd_derefine_flags) { using namespace shamrock::patch; @@ -296,77 +534,76 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); scheduler().for_each_patch_data([&](u64 id_patch, Patch cur_p, PatchDataLayer &pdat) { - u32 old_obj_cnt = pdat.get_obj_cnt(); - - sham::DeviceBuffer &derefine_flags = derefine_list.get(id_patch); - - if (derefine_flags.get_size() > 0) { - + u32 old_obj_cnt = pdat.get_obj_cnt(); + u32 old_obj_cnt_before_refinement = dd_derefine_flags.get(id_patch).get_size(); + auto stream_compact_results = shamalgs::numeric::stream_compact( + dev_sched, dd_derefine_flags.get(id_patch), old_obj_cnt_before_refinement); + if (stream_compact_results.get_size() > 0) { // init flag table + sham::DeviceBuffer keep_block_flag(old_obj_cnt, dev_sched); keep_block_flag.fill(1); sham::DeviceBuffer &buf_cell_min = pdat.get_field_buf_ref(0); sham::DeviceBuffer &buf_cell_max = pdat.get_field_buf_ref(1); - sham::EventList depends_list; auto block_bound_low = buf_cell_min.get_write_access(depends_list); auto block_bound_high = buf_cell_max.get_write_access(depends_list); UserAcc uacc(depends_list, pdat); - auto index_to_deref = derefine_flags.get_read_access(depends_list); + auto index_to_deref = stream_compact_results.get_read_access(depends_list); auto flag_keep = keep_block_flag.get_write_access(depends_list); // edit block content + make flag of blocks to keep auto e = q.submit(depends_list, [&](sycl::handler &cgh) { - cgh.parallel_for(sycl::range<1>(derefine_flags.get_size()), [=](sycl::item<1> gid) { - u32 tid = gid.get_linear_id(); + cgh.parallel_for( + sycl::range<1>(stream_compact_results.get_size()), [=](sycl::item<1> gid) { + u32 tid = gid.get_linear_id(); - u32 idx_to_derefine = index_to_deref[gid]; + u32 idx_to_derefine = index_to_deref[gid]; - // compute old block indexes - std::array old_indexes; + // compute old block indexes + std::array old_indexes; #pragma unroll - for (u32 pid = 0; pid < split_count; pid++) { - old_indexes[pid] = idx_to_derefine + pid; - } + for (u32 pid = 0; pid < split_count; pid++) { + old_indexes[pid] = idx_to_derefine + pid; + } - // load block coords - std::array block_coords; + // load block coords + std::array block_coords; #pragma unroll - for (u32 pid = 0; pid < split_count; pid++) { - block_coords[pid] = BlockCoord{ - block_bound_low[old_indexes[pid]], block_bound_high[old_indexes[pid]]}; - } + for (u32 pid = 0; pid < split_count; pid++) { + block_coords[pid] = BlockCoord{ + block_bound_low[old_indexes[pid]], + block_bound_high[old_indexes[pid]]}; + } - // make new block coord - BlockCoord merged_block_coord = BlockCoord::get_merge(block_coords); + // make new block coord + BlockCoord merged_block_coord = BlockCoord::get_merge(block_coords); - // write new coord - block_bound_low[idx_to_derefine] = merged_block_coord.bmin; - block_bound_high[idx_to_derefine] = merged_block_coord.bmax; + // write new coord + block_bound_low[idx_to_derefine] = merged_block_coord.bmin; + block_bound_high[idx_to_derefine] = merged_block_coord.bmax; // flag the old blocks for removal #pragma unroll - for (u32 pid = 1; pid < split_count; pid++) { - flag_keep[idx_to_derefine + pid] = 0; - } + for (u32 pid = 1; pid < split_count; pid++) { + flag_keep[idx_to_derefine + pid] = 0; + } - // user lambda to fill the fields - uacc.apply_derefine( - old_indexes, block_coords, idx_to_derefine, merged_block_coord, uacc); - }); + // user lambda to fill the fields + + uacc.apply_derefine( + old_indexes, block_coords, idx_to_derefine, merged_block_coord, uacc); + }); }); - sham::EventList resulting_events{e}; + sham::EventList resulting_events; buf_cell_min.complete_event_state(resulting_events); buf_cell_max.complete_event_state(resulting_events); - uacc.finalize(resulting_events, pdat); - + stream_compact_results.complete_event_state(resulting_events); keep_block_flag.complete_event_state(resulting_events); - derefine_flags.complete_event_state(resulting_events); - // stream compact the flags auto buf_keep = shamalgs::numeric::stream_compact(dev_sched, keep_block_flag, old_obj_cnt); @@ -376,18 +613,17 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: "patch", id_patch, "derefine block count ", - old_obj_cnt, - "->", + old_obj_cnt - buf_keep.get_size(), + "new block count = ", buf_keep.get_size()); if (buf_keep.get_size() == 0) { throw std::runtime_error("buf keep must contain something at this point"); } - // remap pdat according to stream compact pdat.index_remap_resize(buf_keep, buf_keep.get_size()); - cell_were_removed = cell_were_removed || derefine_flags.get_size() > 0; + cell_were_removed = cell_were_removed || stream_compact_results.get_size() > 0; } }); @@ -460,11 +696,8 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: Tscal dxfact, Tscal wanted_mass) { - sham::DeviceBuffer &buf_cell_low_bound = pdat.get_field(0).get_buf(); - sham::DeviceBuffer &buf_cell_high_bound = pdat.get_field(1).get_buf(); - - buf_cell_low_bound.complete_event_state(resulting_events); - buf_cell_high_bound.complete_event_state(resulting_events); + pdat.get_field(0).get_buf().complete_event_state(resulting_events); + pdat.get_field(1).get_buf().complete_event_state(resulting_events); pdat.get_field(pdat.pdl().get_field_idx("rho")) .get_buf() .complete_event_state(resulting_events); @@ -550,13 +783,19 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: }; auto get_gid_write = [&](std::array &glid) -> u32 { + // First, get the block id (it's the block to be refine) in wich the new cell glid + // is located. std::array bid = {glid[0] >> AMRBlock::NsideBlockPow, glid[1] >> AMRBlock::NsideBlockPow, glid[2] >> AMRBlock::NsideBlockPow}; - // logger::raw_ln(glid,bid); - return new_blocks[get_index_block(bid)] * AMRBlock::block_size + // get the new global block id + auto new_glob_id = new_blocks[get_index_block(bid)] * AMRBlock::block_size; + + // then added to new_glob_id the local index (between 0 and 7) of the generated + // cells to get. This give the global ids of the new generated cells. + return new_glob_id + AMRBlock::get_index( {glid[0] % AMRBlock::Nside, glid[1] % AMRBlock::Nside, @@ -592,20 +831,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: std::array glid = {lx * 2 + sx, ly * 2 + sy, lz * 2 + sz}; u32 new_cell_idx = get_gid_write(glid); - /* - if (1627 == cur_idx) { - logger::raw_ln( - cur_idx, - "set cell ", - new_cell_idx, - " from cell", - old_cell_idx, - "old", - rho_block, - rho_vel_block, - rhoE_block); - } - */ + acc.rho[new_cell_idx] = rho_block; acc.rho_vel[new_cell_idx] = rho_vel_block; acc.rhoE[new_cell_idx] = rhoE_block; @@ -631,20 +857,20 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: rhoE_block[cell_id] = {}; } + // for each siblings block, perform restriction from its 8 children cells for (u32 pid = 0; pid < 8; pid++) { + auto rho_pid = rho_block[pid]; + auto rho_vel_pid = rho_vel_block[pid]; + auto rhoe_pid = rhoE_block[pid]; + for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) { - rho_block[cell_id] += acc.rho[old_blocks[pid] * AMRBlock::block_size + cell_id]; - rho_vel_block[cell_id] - += acc.rho_vel[old_blocks[pid] * AMRBlock::block_size + cell_id]; - rhoE_block[cell_id] - += acc.rhoE[old_blocks[pid] * AMRBlock::block_size + cell_id]; + rho_pid += acc.rho[old_blocks[pid] * AMRBlock::block_size + cell_id]; + rho_vel_pid += acc.rho_vel[old_blocks[pid] * AMRBlock::block_size + cell_id]; + rhoe_pid += acc.rhoE[old_blocks[pid] * AMRBlock::block_size + cell_id]; } - } - - for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) { - rho_block[cell_id] /= 8; - rho_vel_block[cell_id] /= 8; - rhoE_block[cell_id] /= 8; + rho_block[pid] = rho_pid * (1. / 8.); + rho_vel_block[pid] = rho_vel_pid * (1. / 8.); + rhoE_block[pid] = rhoe_pid * (1. / 8.); } for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) { @@ -656,6 +882,144 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: } }; + class RefineCritPseudoGradientAccessor { + public: + Tscal one_over_Nside = 1. / AMRBlock::Nside; + const TgridVec *block_low_bound; + const TgridVec *block_high_bound; + const Tscal *block_rho; + const f64 *block_pressure; + const f64_3 *block_velocity; + + Tscal error_min; + Tscal error_max; + + AMRGraphLinkiterator cell_graph_xp; + AMRGraphLinkiterator cell_graph_xm; + AMRGraphLinkiterator cell_graph_yp; + AMRGraphLinkiterator cell_graph_ym; + AMRGraphLinkiterator cell_graph_zp; + AMRGraphLinkiterator cell_graph_zm; + + RefineCritPseudoGradientAccessor( + sham::EventList &depends_list, + Storage &storage, + u64 id_patch, + shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat, + Tscal err_min, + Tscal err_max) + : error_min(err_min), error_max(err_max), + cell_graph_xp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_xm( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xm) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_yp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::yp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_ym( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::ym) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_zp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_zm( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zm) + .get(id_patch) + .get() + .get_read_access(depends_list)) + + { + block_low_bound = pdat.get_field(0).get_buf().get_read_access(depends_list); + block_high_bound = pdat.get_field(1).get_buf().get_read_access(depends_list); + block_rho = shambase::get_check_ref(storage.rho_primitive) + .get_buf(id_patch) + .get_read_access(depends_list); + block_pressure = shambase::get_check_ref(storage.press) + .get_buf(id_patch) + .get_read_access(depends_list); + block_velocity = shambase::get_check_ref(storage.vel) + .get_buf(id_patch) + .get_read_access(depends_list); + } + + void finalize( + sham::EventList &resulting_events, + Storage &storage, + u64 id_patch, + shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat, + Tscal err_min, + Tscal err_max) { + + pdat.get_field(0).get_buf().complete_event_state(resulting_events); + pdat.get_field(1).get_buf().complete_event_state(resulting_events); + + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xm) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::yp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::ym) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zm) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + + shambase::get_check_ref(storage.rho_primitive) + .get_buf(id_patch) + .complete_event_state(resulting_events); + + shambase::get_check_ref(storage.press) + .get_buf(id_patch) + .complete_event_state(resulting_events); + + shambase::get_check_ref(storage.vel) + .get_buf(id_patch) + .complete_event_state(resulting_events); + } + }; + using AMRmode_None = typename AMRMode::None; using AMRmode_DensityBased = typename AMRMode::DensityBased; From 30aae1e545c2c3161ff8fa66401fff171bed5857 Mon Sep 17 00:00:00 2001 From: Leodasce Sewanou Date: Thu, 19 Mar 2026 19:33:20 +0100 Subject: [PATCH 3/5] [Ramses][AMR] RefinementHandler modules + first tests --- .../TO_MIGRATE/ramses/sod_amr_reflective.py | 237 ++++++++++++++++ .../shammodels/ramses/SolverConfig.hpp | 14 +- .../ramses/modules/ComputeLevel0CellSize.hpp | 2 +- .../ramses/modules/ConsToPrimGas.hpp | 13 +- .../modules/SlopeLimitedGradientUtilities.hpp | 157 +++++++++++ src/shammodels/ramses/src/Solver.cpp | 34 ++- .../src/modules/AMRGridRefinementHandler.cpp | 254 ++++++++++++++++-- .../ramses/src/modules/ComputeAMRLevel.cpp | 2 + .../ramses/src/modules/ConsToPrimGas.cpp | 23 +- src/shammodels/ramses/src/pyRamsesModel.cpp | 8 + 10 files changed, 704 insertions(+), 40 deletions(-) create mode 100644 examples/TO_MIGRATE/ramses/sod_amr_reflective.py diff --git a/examples/TO_MIGRATE/ramses/sod_amr_reflective.py b/examples/TO_MIGRATE/ramses/sod_amr_reflective.py new file mode 100644 index 0000000000..de8a1b065b --- /dev/null +++ b/examples/TO_MIGRATE/ramses/sod_amr_reflective.py @@ -0,0 +1,237 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np + +import shamrock + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + +multx = 1 +multy = 1 +multz = 1 +max_amr_lev = 2 +cell_size = 2 << max_amr_lev # refinement is limited to cell_size = 2 +base = 16 + +cfg = model.gen_default_config() +scale_fact = 1 / (cell_size * base * multx) +cfg.set_scale_factor(scale_fact) + +gamma = 1.4 +cfg.set_eos_gamma(gamma) +cfg.set_Csafe(0.3) +cfg.set_boundary_condition("x", "reflective") +cfg.set_boundary_condition("y", "reflective") +cfg.set_boundary_condition("z", "reflective") +cfg.set_riemann_solver_hllc() + + +cfg.set_slope_lim_minmod() +cfg.set_face_time_interpolation(True) + +err_min = 0.25 +err_max = 0.15 + +cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max) + +mass_crit = 1e-6 * 5 * 2 * 2 +# cfg.set_amr_mode_density_based(crit_mass=mass_crit) + + +crit_refin = 0.1 +crit_coars = 0.2 +# cfg.set_amr_mode_second_order_derivative_based(crit_min=crit_refin, crit_max=crit_coars) +model.set_solver_config(cfg) + + +model.init_scheduler(int(1e7), 1) +model.make_base_grid( + (0, 0, 0), (cell_size, cell_size, cell_size), (base * multx, base * multy, base * multz) +) + + +def rho_map(rmin, rmax): + + x, y, z = rmin + if y < 0.5: + return 1 + else: + return 0.125 + + +etot_L = 1.0 / (gamma - 1) +etot_R = 0.1 / (gamma - 1) + + +def rhoetot_map(rmin, rmax): + + rho = rho_map(rmin, rmax) + + x, y, z = rmin + if y < 0.5: + return etot_L + else: + return etot_R + + +def rhovel_map(rmin, rmax): + rho = rho_map(rmin, rmax) + + return (0, 0, 0) + + +model.set_field_value_lambda_f64("rho", rho_map) +model.set_field_value_lambda_f64("rhoetot", rhoetot_map) +model.set_field_value_lambda_f64_3("rhovel", rhovel_map) + + +def convert_to_cell_coords(dic): + + cmin = dic["cell_min"] + cmax = dic["cell_max"] + + xmin = [] + ymin = [] + zmin = [] + xmax = [] + ymax = [] + zmax = [] + + for i in range(len(cmin)): + m, M = cmin[i], cmax[i] + + mx, my, mz = m + Mx, My, Mz = M + + for j in range(8): + a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j) + + x, y, z = a + xmin.append(x) + ymin.append(y) + zmin.append(z) + + x, y, z = b + xmax.append(x) + ymax.append(y) + zmax.append(z) + + dic["xmin"] = np.array(xmin) + dic["ymin"] = np.array(ymin) + dic["zmin"] = np.array(zmin) + dic["xmax"] = np.array(xmax) + dic["ymax"] = np.array(ymax) + dic["zmax"] = np.array(zmax) + + return dic + + +t_target = 0.245 + +dt = 0 +t = 0 +freq = 1 +dX0 = [] +for i in range(5): + next_dt = model.evolve_once_override_time(t, dt) + if i == 0: + dic0 = convert_to_cell_coords(ctx.collect_data()) + dX0.append(dic0["ymax"][i] - dic0["ymin"][i]) + + t += dt + dt = next_dt + + if i % freq == 0: + model.dump_vtk(f"test{i:04d}.vtk") + + if t_target < t + next_dt: + dt = t_target - t + if t == t_target: + model.dump_vtk(f"test{i:04d}.vtk") + break + +xref = 0.5 +xrange = 0.5 +sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1) +sodanalysis = model.make_analysis_sodtube(sod, (0, 1, 0), t_target, xref, 0.0, 1.0) + + +################# +### Plot +################# +# do plot or not +if True: + dic = convert_to_cell_coords(ctx.collect_data()) + + X = [] + dX = [] + rho = [] + rhovelx = [] + rhoetot = [] + + for i in range(len(dic["ymin"])): + X.append(dic["ymin"][i]) + dX.append(dic["ymax"][i] - dic["ymin"][i]) + rho.append(dic["rho"][i]) + rhovelx.append(dic["rhovel"][i][1]) + rhoetot.append(dic["rhoetot"][i]) + + X = np.array(X) + dX = np.array(dX) + dX0 = np.array(dX0) + rho = np.array(rho) + rhovelx = np.array(rhovelx) + rhoetot = np.array(rhoetot) + + vx = rhovelx / rho + + fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(9, 6), dpi=125) + + ax1 = plt.gca() + ax2 = ax1.twinx() + + l = -np.log2(dX / max(dX0.max(), dX.max())) + + ax1.scatter(X, rho, rasterized=True, s=12 * np.ones(X.shape), label="rho") + ax1.scatter(X, vx, rasterized=True, s=12 * np.ones(X.shape), label="v") + ax1.scatter( + X, + (rhoetot - 0.5 * rho * (vx**2)) * (gamma - 1), + rasterized=True, + s=12 * np.ones(X.shape), + label="P", + ) + idx = np.argsort(X) + ax2.plot(X[idx], l[idx], color="purple", marker="D", linewidth=2.0, ls="-.", label="AMR level") + # plt.scatter(X,rhoetot, rasterized=True,label="rhoetot") + ax1.legend(loc=0) + ax2.legend(loc=0) + ax1.grid() + + #### add analytical soluce + arr_x = np.linspace(xref - xrange, xref + xrange, 1000) + + arr_rho = [] + arr_P = [] + arr_vx = [] + + for i in range(len(arr_x)): + x_ = arr_x[i] - xref + + _rho, _vx, _P = sod.get_value(t_target, x_) + arr_rho.append(_rho) + arr_vx.append(_vx) + arr_P.append(_P) + + ax1.plot(arr_x, arr_rho, ls="--", lw=2.0, color="black", label="analytic") + ax1.plot(arr_x, arr_vx, ls="--", lw=2.0, color="black") + ax1.plot(arr_x, arr_P, ls="--", lw=2.0, color="black") + ax2.set_ylabel("AMR level") + plt.title(f"Threshold = {err_max}, derefinement factor = {err_min}") + plt.savefig("sod_tube.png") + ####### diff --git a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp index 31831e8356..874de25b05 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp @@ -100,14 +100,22 @@ namespace shammodels::basegodunov { Tscal crit_mass; }; - using mode = std::variant; + struct PseudoGradientBased { + Tscal error_min; + Tscal error_max; + }; + + using mode = std::variant; mode config = None{}; void set_refine_none() { config = None{}; } void set_refine_density_based(Tscal crit_mass) { config = DensityBased{crit_mass}; } + void set_refine_pseudo_gradient_based(Tscal error_min, Tscal error_max) { + config = PseudoGradientBased{error_min, error_max}; + } - bool need_level_zero_compute() { return false; } - bool need_amr_level_compute() { return false; } + bool need_level_zero_compute() { return true; } + bool need_amr_level_compute() { return true; } }; struct BCConfig { diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeLevel0CellSize.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeLevel0CellSize.hpp index 008b744194..af6afa8298 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeLevel0CellSize.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeLevel0CellSize.hpp @@ -12,7 +12,7 @@ /** * @file ComputeLevel0CellSize.hpp * @author Léodasce Sewanou (leodasce.sewanou@ens-lyon.fr) - * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Timothée David--Cléris (tim.shamrock@proton.me) --no git blame-- * @brief * */ diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/ConsToPrimGas.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/ConsToPrimGas.hpp index f129e772bc..15a9f4c58c 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/ConsToPrimGas.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/ConsToPrimGas.hpp @@ -36,6 +36,9 @@ namespace shammodels::basegodunov::modules { const shamrock::solvergraph::IFieldSpan &spans_rho; const shamrock::solvergraph::IFieldSpan &spans_rhov; const shamrock::solvergraph::IFieldSpan &spans_rhoe; + /**/ + shamrock::solvergraph::IFieldSpan &spans_rho_fields; + /**/ shamrock::solvergraph::IFieldSpan &spans_vel; shamrock::solvergraph::IFieldSpan &spans_P; }; @@ -45,10 +48,13 @@ namespace shammodels::basegodunov::modules { std::shared_ptr> spans_rho, std::shared_ptr> spans_rhov, std::shared_ptr> spans_rhoe, + /**/ + std::shared_ptr> spans_rho_fields, + /**/ std::shared_ptr> spans_vel, std::shared_ptr> spans_P) { __internal_set_ro_edges({sizes, spans_rho, spans_rhov, spans_rhoe}); - __internal_set_rw_edges({spans_vel, spans_P}); + __internal_set_rw_edges({spans_rho_fields, spans_vel, spans_P}); } inline Edges get_edges() { @@ -57,8 +63,9 @@ namespace shammodels::basegodunov::modules { get_ro_edge>(1), get_ro_edge>(2), get_ro_edge>(3), - get_rw_edge>(0), - get_rw_edge>(1), + get_rw_edge>(0), + get_rw_edge>(1), + get_rw_edge>(2), }; } diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp index 1dd4f519ab..5014be9986 100644 --- a/src/shammodels/ramses/include/shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp +++ b/src/shammodels/ramses/include/shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp @@ -280,4 +280,161 @@ namespace { return {lim_slope_W_x, lim_slope_W_y, lim_slope_W_z}; } + + template + inline T get_pseudo_grad( + const u32 cell_global_id, + const AMRGraphLinkiterator &graph_iter_xp, + const AMRGraphLinkiterator &graph_iter_xm, + const AMRGraphLinkiterator &graph_iter_yp, + const AMRGraphLinkiterator &graph_iter_ym, + const AMRGraphLinkiterator &graph_iter_zp, + const AMRGraphLinkiterator &graph_iter_zm, + ACCField &&field_access) + + { + + using namespace sham; + using namespace sham::details; + + auto get_avg_neigh = [&](auto &graph_links, u32 dir) -> T { + T acc = shambase::VectorProperties::get_zero(); + u32 cnt = graph_links.for_each_object_link_cnt(cell_global_id, [&](u32 id_b) { + acc += field_access(id_b); + }); + + return (cnt > 0) ? acc / cnt : shambase::VectorProperties::get_zero(); + }; + + auto epsilon = shambase::get_epsilon(); + T u_cur = field_access(cell_global_id); + T u_xp = get_avg_neigh(graph_iter_xp, 0); + T u_xm = get_avg_neigh(graph_iter_xm, 1); + T u_yp = get_avg_neigh(graph_iter_yp, 2); + T u_ym = get_avg_neigh(graph_iter_ym, 3); + T u_zp = get_avg_neigh(graph_iter_zp, 4); + T u_zm = get_avg_neigh(graph_iter_zm, 5); + + // RAMSES LIKE + + T x_scal = 2 + * g_sycl_max( + g_sycl_abs((u_cur - u_xm) / (epsilon + u_cur + u_xm)), + g_sycl_abs((u_cur - u_xp) / (epsilon + u_cur + u_xp))); + + T y_scal = 2 + * g_sycl_max( + g_sycl_abs((u_cur - u_ym) / (epsilon + u_cur + u_ym)), + g_sycl_abs((u_cur - u_yp) / (epsilon + u_cur + u_yp))); + T z_scal = 2 + * g_sycl_max( + g_sycl_abs((u_cur - u_zm) / (epsilon + u_cur + u_zm)), + g_sycl_abs((u_cur - u_zp) / (epsilon + u_cur + u_zp))); + + T res = g_sycl_max(x_scal, g_sycl_max(y_scal, z_scal)); + return res; + } + + /** + */ + template + inline T baryonic_normalized_slope_criterion( + const u32 cell_global_id, + const AMRGraphLinkiterator &graph_iter_xp, + const AMRGraphLinkiterator &graph_iter_xm, + const AMRGraphLinkiterator &graph_iter_yp, + const AMRGraphLinkiterator &graph_iter_ym, + const AMRGraphLinkiterator &graph_iter_zp, + const AMRGraphLinkiterator &graph_iter_zm, + ACCField &&field_access) + + { + + using namespace sham; + using namespace sham::details; + + auto get_avg_neigh = [&](auto &graph_links, u32 dir) -> T { + T acc = shambase::VectorProperties::get_zero(); + u32 cnt = graph_links.for_each_object_link_cnt(cell_global_id, [&](u32 id_b) { + acc += field_access(id_b); + }); + + return (cnt > 0) ? acc / cnt : shambase::VectorProperties::get_zero(); + }; + + auto epsilon = shambase::get_epsilon(); + T u_cur = field_access(cell_global_id); + T u_xp = get_avg_neigh(graph_iter_xp, 0); + T u_xm = get_avg_neigh(graph_iter_xm, 1); + T u_yp = get_avg_neigh(graph_iter_yp, 2); + T u_ym = get_avg_neigh(graph_iter_ym, 3); + T u_zp = get_avg_neigh(graph_iter_zp, 4); + T u_zm = get_avg_neigh(graph_iter_zm, 5); + + T norm_slope_x = g_sycl_abs((u_xm - u_xp) / (2 * u_cur + epsilon)); + T norm_slope_y = g_sycl_abs((u_ym - u_yp) / (2 * u_cur + epsilon)); + T norm_slope_z = g_sycl_abs((u_zm - u_zp) / (2 * u_cur + epsilon)); + + T res = g_sycl_max(norm_slope_x, g_sycl_max(norm_slope_y, norm_slope_z)); + + return res; + } + + /** + */ + template + inline T modif_second_derivative( + const u32 cell_global_id, + const AMRGraphLinkiterator &graph_iter_xp, + const AMRGraphLinkiterator &graph_iter_xm, + const AMRGraphLinkiterator &graph_iter_yp, + const AMRGraphLinkiterator &graph_iter_ym, + const AMRGraphLinkiterator &graph_iter_zp, + const AMRGraphLinkiterator &graph_iter_zm, + ACCField &&field_access) { + using namespace sham; + using namespace sham::details; + + auto get_avg_neigh = [&](auto &graph_links) -> T { + T acc = shambase::VectorProperties::get_zero(); + u32 cnt = graph_links.for_each_object_link_cnt(cell_global_id, [&](u32 id_b) { + acc += field_access(id_b); + }); + return (cnt > 0) ? acc / cnt : shambase::VectorProperties::get_zero(); + }; + + auto eps_ref = 0.01; + auto epsilon = shambase::get_epsilon(); + T u_cur = field_access(cell_global_id); + T u_xp = get_avg_neigh(graph_iter_xp); + T u_xm = get_avg_neigh(graph_iter_xm); + T u_yp = get_avg_neigh(graph_iter_yp); + T u_ym = get_avg_neigh(graph_iter_ym); + T u_zp = get_avg_neigh(graph_iter_zp); + T u_zm = get_avg_neigh(graph_iter_zm); + + T delta_u_xp = u_xp - u_cur; + T delta_u_xm = u_xm - u_cur; + T delta_u_yp = u_yp - u_cur; + T delta_u_ym = u_ym - u_cur; + T delta_u_zp = u_zp - u_cur; + T delta_u_zm = u_zm - u_cur; + + T scalar_x = g_sycl_abs(u_xp) + g_sycl_abs(u_xm) + 2 * g_sycl_abs(u_cur); + T scalar_y = g_sycl_abs(u_yp) + g_sycl_abs(u_ym) + 2 * g_sycl_abs(u_cur); + T scalar_z = g_sycl_abs(u_zp) + g_sycl_abs(u_zm) + 2 * g_sycl_abs(u_cur); + + T res_x + = g_sycl_abs(delta_u_xm + delta_u_xp) + / (g_sycl_abs(delta_u_xm) + g_sycl_abs(delta_u_xp) + eps_ref * scalar_x + epsilon); + T res_y + = g_sycl_abs(delta_u_ym + delta_u_yp) + / (g_sycl_abs(delta_u_ym) + g_sycl_abs(delta_u_yp) + eps_ref * scalar_y + epsilon); + T res_z + = g_sycl_abs(delta_u_zm + delta_u_zp) + / (g_sycl_abs(delta_u_zm) + g_sycl_abs(delta_u_zp) + eps_ref * scalar_z + epsilon); + + return (res_x + res_y + res_z); + } + } // namespace diff --git a/src/shammodels/ramses/src/Solver.cpp b/src/shammodels/ramses/src/Solver.cpp index 923ac97f10..037c1eab62 100644 --- a/src/shammodels/ramses/src/Solver.cpp +++ b/src/shammodels/ramses/src/Solver.cpp @@ -388,6 +388,9 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.press = std::make_shared>(AMRBlock::block_size, "P", "P"); + storage.rho_primitive = std::make_shared>( + AMRBlock::block_size, "rho-prim", "rho-prim"); + if (solver_config.is_dust_on()) { u32 ndust = solver_config.dust_config.ndust; @@ -425,6 +428,12 @@ void shammodels::basegodunov::Solver::init_solver_graph() { "level0_amr", "level0_amr"); } + if (solver_config.amr_mode.need_amr_level_compute()) { + using TgridUint = typename std::make_unsigned>::type; + storage.amr_block_levels + = std::make_shared>(1, "", ""); + } + storage.grad_rho = std::make_shared>( AMRBlock::block_size, "grad_rho", "\\nabla \\rho"); storage.dx_v = std::make_shared>( @@ -1001,6 +1010,7 @@ void shammodels::basegodunov::Solver::init_solver_graph() { storage.refs_rho, storage.refs_rhov, storage.refs_rhoe, + storage.rho_primitive, storage.vel, storage.press); @@ -1614,8 +1624,28 @@ void shammodels::basegodunov::Solver::evolve_once() { shambase::throw_unimplemented(); } - modules::AMRGridRefinementHandler refinement(context, solver_config, storage); - refinement.update_refinement(); + { + + modules::NodeConsToPrimGas node_ctp_after_updated{ + AMRBlock::block_size, solver_config.eos_gamma}; + node_ctp_after_updated.set_edges( + storage.block_counts_with_ghost, + storage.refs_rho, + storage.refs_rhov, + storage.refs_rhoe, + // /**/ + storage.rho_primitive, + // /**/ + storage.vel, + storage.press); + + node_ctp_after_updated.evaluate(); + } + + if (dt_input > 0) { + modules::AMRGridRefinementHandler refinement(context, solver_config, storage); + refinement.update_refinement(); + } modules::ComputeCFL cfl_compute(context, solver_config, storage); f64 new_dt = cfl_compute.compute_cfl(); diff --git a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp index beb12297b8..b483bfe1ee 100644 --- a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp +++ b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp @@ -17,8 +17,10 @@ #include "shambase/memory.hpp" #include "shamalgs/details/algorithm/algorithm.hpp" #include "shamcomm/logs.hpp" +#include "shammodels/common/amr/NeighGraph.hpp" #include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp" #include "shammodels/ramses/modules/AMRSortBlocks.hpp" +#include "shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp" #include template @@ -44,7 +46,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: { sham::EventList depends_list; - UserAcc uacc(depends_list, id_patch, cur_p, pdat, args...); + UserAcc uacc(depends_list, storage, id_patch, cur_p, pdat, args...); auto refine_acc = refine_flags.get_write_access(depends_list); auto derefine_acc = derefine_flags.get_write_access(depends_list); @@ -72,7 +74,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: refine_flags.complete_event_state(resulting_events); derefine_flags.complete_event_state(resulting_events); - uacc.finalize(resulting_events, id_patch, cur_p, pdat, args...); + uacc.finalize(resulting_events, storage, id_patch, cur_p, pdat, args...); } dd_refine_flags.add_obj(id_patch, std::move(refine_flags)); @@ -458,7 +460,7 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: auto block_bound_low = buf_cell_min.get_write_access(depends_list); auto block_bound_high = buf_cell_max.get_write_access(depends_list); - UserAcc uacc(depends_list, pdat); + UserAcc uacc(depends_list, storage, id_patch, pdat); auto index_to_ref = stream_compaction_result.get_read_access(depends_list); // Refine the block (set the positions) and fill the corresponding fields @@ -509,7 +511,7 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: buf_cell_min.complete_event_state(resulting_events); buf_cell_max.complete_event_state(resulting_events); stream_compaction_result.complete_event_state(e); - uacc.finalize(resulting_events, pdat); + uacc.finalize(resulting_events, storage, id_patch, pdat); } shamlog_debug_ln("AMRGrid", "patch ", id_patch, "new block count = ", pdat.get_obj_cnt()); sum_block_count += pdat.get_obj_cnt(); @@ -540,21 +542,26 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: dev_sched, dd_derefine_flags.get(id_patch), old_obj_cnt_before_refinement); if (stream_compact_results.get_size() > 0) { // init flag table + sycl::buffer keep_block_flag + = shamalgs::algorithm::gen_buffer_device(q.q, old_obj_cnt, [](u32 i) -> u32 { + return 1; + }); - sham::DeviceBuffer keep_block_flag(old_obj_cnt, dev_sched); - keep_block_flag.fill(1); + // sham::DeviceBuffer keep_block_flag(old_obj_cnt, dev_sched); + // keep_block_flag.fill(1); sham::DeviceBuffer &buf_cell_min = pdat.get_field_buf_ref(0); sham::DeviceBuffer &buf_cell_max = pdat.get_field_buf_ref(1); sham::EventList depends_list; auto block_bound_low = buf_cell_min.get_write_access(depends_list); auto block_bound_high = buf_cell_max.get_write_access(depends_list); - UserAcc uacc(depends_list, pdat); + UserAcc uacc(depends_list, storage, id_patch, pdat); auto index_to_deref = stream_compact_results.get_read_access(depends_list); - auto flag_keep = keep_block_flag.get_write_access(depends_list); + // auto flag_keep = keep_block_flag.get_write_access(depends_list); // edit block content + make flag of blocks to keep auto e = q.submit(depends_list, [&](sycl::handler &cgh) { + sycl::accessor flag_keep{keep_block_flag, cgh, sycl::read_write}; cgh.parallel_for( sycl::range<1>(stream_compact_results.get_size()), [=](sycl::item<1> gid) { u32 tid = gid.get_linear_id(); @@ -601,29 +608,54 @@ bool shammodels::basegodunov::modules::AMRGridRefinementHandler: buf_cell_min.complete_event_state(resulting_events); buf_cell_max.complete_event_state(resulting_events); - uacc.finalize(resulting_events, pdat); + uacc.finalize(resulting_events, storage, id_patch, pdat); + stream_compact_results.complete_event_state(resulting_events); - keep_block_flag.complete_event_state(resulting_events); - // stream compact the flags - auto buf_keep - = shamalgs::numeric::stream_compact(dev_sched, keep_block_flag, old_obj_cnt); - shamlog_debug_ln( + // stream compact the flags (get new block ids map after merged) + auto [opt_buf, len] + = shamalgs::numeric::stream_compact(q.q, keep_block_flag, old_obj_cnt); + + logger::info_ln( "AMR Grid", "patch", id_patch, - "derefine block count ", - old_obj_cnt - buf_keep.get_size(), + "derefine block count = ", + old_obj_cnt - len, "new block count = ", - buf_keep.get_size()); + len); - if (buf_keep.get_size() == 0) { - throw std::runtime_error("buf keep must contain something at this point"); + if (!opt_buf) { + throw std::runtime_error("opt buf must contain something at this point"); } - // remap pdat according to stream compact - pdat.index_remap_resize(buf_keep, buf_keep.get_size()); + + // remap pdat according to stream compact (for each field in patchdataleyer resize + // according to new block ids map) + pdat.index_remap_resize(*opt_buf, len); cell_were_removed = cell_were_removed || stream_compact_results.get_size() > 0; + + // // keep_block_flag.complete_event_state(resulting_events); + // // stream compact the flags + // auto buf_keep + // = shamalgs::numeric::stream_compact(dev_sched, keep_block_flag, old_obj_cnt); + + // shamlog_debug_ln( + // "AMR Grid", + // "patch", + // id_patch, + // "derefine block count ", + // old_obj_cnt - buf_keep.get_size(), + // "new block count = ", + // buf_keep.get_size()); + + // if (buf_keep.get_size() == 0) { + // throw std::runtime_error("buf keep must contain something at this point"); + // } + // // remap pdat according to stream compact + // pdat.index_remap_resize(buf_keep, buf_keep.get_size()); + + // cell_were_removed = cell_were_removed || stream_compact_results.get_size() > 0; } }); @@ -674,6 +706,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: RefineCritBlock( sham::EventList &depends_list, + Storage &storage, u64 id_patch, shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat, @@ -690,6 +723,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: void finalize( sham::EventList &resulting_events, + Storage &storage, u64 id_patch, shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat, @@ -742,18 +776,104 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: f64 *rho; f64_3 *rho_vel; f64 *rhoE; + u64 p_id; - RefineCellAccessor(sham::EventList &depends_list, shamrock::patch::PatchDataLayer &pdat) { + // this will be needed for interpolation during refinement + AMRGraphLinkiterator cell_graph_xp; + AMRGraphLinkiterator cell_graph_xm; + AMRGraphLinkiterator cell_graph_yp; + AMRGraphLinkiterator cell_graph_ym; + AMRGraphLinkiterator cell_graph_zp; + AMRGraphLinkiterator cell_graph_zm; + RefineCellAccessor( + sham::EventList &depends_list, + Storage &storage, + u64 &id_patch, + shamrock::patch::PatchDataLayer &pdat) + : cell_graph_xp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::xp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_xm( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::xm) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_yp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::yp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_ym( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::ym) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_zp( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::zp) + .get(id_patch) + .get() + .get_read_access(depends_list)), + cell_graph_zm( + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction::zm) + .get(id_patch) + .get() + .get_read_access(depends_list)) + + { + p_id = id_patch; rho = pdat.get_field(2).get_buf().get_write_access(depends_list); rho_vel = pdat.get_field(3).get_buf().get_write_access(depends_list); rhoE = pdat.get_field(4).get_buf().get_write_access(depends_list); } - void finalize(sham::EventList &resulting_events, shamrock::patch::PatchDataLayer &pdat) { + void finalize( + sham::EventList &resulting_events, + Storage &storage, + u64 &id_patch, + shamrock::patch::PatchDataLayer &pdat) { pdat.get_field(2).get_buf().complete_event_state(resulting_events); pdat.get_field(3).get_buf().complete_event_state(resulting_events); pdat.get_field(4).get_buf().complete_event_state(resulting_events); + + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::xm) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::yp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::ym) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zp) + .get(id_patch) + .get() + .complete_event_state(resulting_events); + shambase::get_check_ref(storage.cell_graph_edge) + .get_refs_dir(Direction_::zm) + .get(id_patch) + .get() + .complete_event_state(resulting_events); } void apply_refine( @@ -1018,10 +1138,69 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: .get_buf(id_patch) .complete_event_state(resulting_events); } + + void refine_criterion( + u32 block_id, + RefineCritPseudoGradientAccessor acc, + bool &should_refine, + bool &should_derefine) const { + TgridVec low_bound = acc.block_low_bound[block_id]; + TgridVec high_bound = acc.block_high_bound[block_id]; + + Tscal block_rho_slope = shambase::VectorProperties::get_zero(); + for (u32 i = 0; i < AMRBlock::block_size; i++) { + block_rho_slope = sham::details::g_sycl_max( + block_rho_slope, + baryonic_normalized_slope_criterion( + i + block_id * AMRBlock::block_size, + cell_graph_xp, + cell_graph_xm, + cell_graph_yp, + cell_graph_ym, + cell_graph_zp, + cell_graph_zm, + [=](u32 id) { + return block_rho[id]; + })); + } + + Tscal block_press_grad = shambase::VectorProperties::get_zero(); + for (u32 i = 0; i < AMRBlock::block_size; i++) { + block_press_grad = sham::details::g_sycl_max( + block_press_grad, + get_pseudo_grad( + i + block_id * AMRBlock::block_size, + cell_graph_xp, + cell_graph_xm, + cell_graph_yp, + cell_graph_ym, + cell_graph_zp, + cell_graph_zm, + [=](u32 id) { + return block_pressure[id]; + })); + } + + Tscal error = sham::details::g_sycl_max( + block_press_grad, sham::details::g_sycl_max(block_rho_slope, 0.0)); + + should_refine = false; + should_derefine = false; + if (error > error_max) { + should_refine = true; + } else if (error < (error_min * error_max)) { + should_derefine = true; + } + + should_refine = should_refine && (high_bound.x() - low_bound.x() > AMRBlock::Nside); + should_refine = should_refine && (high_bound.y() - low_bound.y() > AMRBlock::Nside); + should_refine = should_refine && (high_bound.z() - low_bound.z() > AMRBlock::Nside); + } }; - using AMRmode_None = typename AMRMode::None; - using AMRmode_DensityBased = typename AMRMode::DensityBased; + using AMRmode_None = typename AMRMode::None; + using AMRmode_DensityBased = typename AMRMode::DensityBased; + using AMRmode_PseudoGradientBased = typename AMRMode::PseudoGradientBased; bool has_cell_order_changed = false; @@ -1053,6 +1232,31 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler: has_cell_order_changed = has_cell_order_changed || (change_refine || change_derefine); } + else if ( + AMRmode_PseudoGradientBased *cfg + = std::get_if(&solver_config.amr_mode.config)) { + + shambase::DistributedData> refine_list; + shambase::DistributedData> derefine_list; + + gen_refine_block_changes( + refine_list, derefine_list, cfg->error_min, cfg->error_max); + + ///// enforce 2:1 for refinement /////// + enforce_two_to_one_refinement(std::move(refine_list)); + + /////// enforce 2:1 for derefinement ////// + enforce_two_to_one_derefinement(std::move(derefine_list), std::move(refine_list)); + + //////// apply refine //////// + // Note that this only add new blocks at the end of the patchdata + bool change_refine = internal_refine_grid(std::move(refine_list)); + + bool change_derefine = internal_derefine_grid(std::move(derefine_list)); + + has_cell_order_changed = has_cell_order_changed || (change_refine) || (change_derefine); + } + if (has_cell_order_changed) { // Ensure that the blocks are sorted before refinement AMRSortBlocks block_sorter(context, solver_config, storage); diff --git a/src/shammodels/ramses/src/modules/ComputeAMRLevel.cpp b/src/shammodels/ramses/src/modules/ComputeAMRLevel.cpp index 553746974a..71520568f5 100644 --- a/src/shammodels/ramses/src/modules/ComputeAMRLevel.cpp +++ b/src/shammodels/ramses/src/modules/ComputeAMRLevel.cpp @@ -20,6 +20,7 @@ #include "shambackends/kernel_call.hpp" #include "shambackends/kernel_call_distrib.hpp" #include "shambackends/math.hpp" +#include "shamcomm/logs.hpp" #include "shammodels/ramses/modules/ComputeAMRLevel.hpp" #include "shamrock/patch/PatchDataField.hpp" #include "shamsys/NodeInstance.hpp" @@ -28,6 +29,7 @@ namespace shammodels::basegodunov::modules { template void ComputeAMRLevel::_impl_evaluate_internal() { + auto edges = get_edges(); edges.block_min.check_sizes(edges.sizes.indexes); diff --git a/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp b/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp index 30469afa8e..ee89fe36db 100644 --- a/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp +++ b/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp @@ -31,7 +31,9 @@ namespace { const shambase::DistributedData> &spans_rho, const shambase::DistributedData> &spans_rhov, const shambase::DistributedData> &spans_rhoe, - + // /**/ + shambase::DistributedData> &spans_rho_fields, + // /**/ shambase::DistributedData> &spans_vel, shambase::DistributedData> &spans_P, const shambase::DistributedData &sizes, @@ -47,13 +49,17 @@ namespace { sham::distributed_data_kernel_call( shamsys::instance::get_compute_scheduler_ptr(), sham::DDMultiRef{spans_rho, spans_rhov, spans_rhoe}, - sham::DDMultiRef{spans_vel, spans_P}, + sham::DDMultiRef{spans_rho_fields, spans_vel, spans_P}, cell_counts, [gamma]( u32 i, const Tscal *__restrict rho, const Tvec *__restrict rhov, const Tscal *__restrict rhoe, + + // /**/ + Tscal *__restrict rho_field, + // /**/ Tvec *__restrict vel, Tscal *__restrict P) { auto conststate = shammath::ConsState{rho[i], rhoe[i], rhov[i]}; @@ -62,8 +68,9 @@ namespace { SHAM_ASSERT(prim_state.press >= 0.0); - vel[i] = prim_state.vel; - P[i] = prim_state.press; + vel[i] = prim_state.vel; + P[i] = prim_state.press; + rho_field[i] = rho[i]; }); } }; @@ -80,6 +87,7 @@ namespace shammodels::basegodunov::modules { edges.spans_rhov.check_sizes(edges.sizes.indexes); edges.spans_rhoe.check_sizes(edges.sizes.indexes); + edges.spans_rho_fields.ensure_sizes(edges.sizes.indexes); edges.spans_vel.ensure_sizes(edges.sizes.indexes); edges.spans_P.ensure_sizes(edges.sizes.indexes); @@ -87,6 +95,9 @@ namespace shammodels::basegodunov::modules { edges.spans_rho.get_spans(), edges.spans_rhov.get_spans(), edges.spans_rhoe.get_spans(), + // /**/ + edges.spans_rho_fields.get_spans(), + // /**/ edges.spans_vel.get_spans(), edges.spans_P.get_spans(), edges.sizes.indexes, @@ -101,8 +112,8 @@ namespace shammodels::basegodunov::modules { auto rho = get_ro_edge_base(1).get_tex_symbol(); auto rhov = get_ro_edge_base(2).get_tex_symbol(); auto rhoe = get_ro_edge_base(3).get_tex_symbol(); - auto vel = get_rw_edge_base(0).get_tex_symbol(); - auto P = get_rw_edge_base(1).get_tex_symbol(); + auto vel = get_rw_edge_base(1).get_tex_symbol(); + auto P = get_rw_edge_base(2).get_tex_symbol(); std::string tex = R"tex( Conservative to primitive variable (gas) diff --git a/src/shammodels/ramses/src/pyRamsesModel.cpp b/src/shammodels/ramses/src/pyRamsesModel.cpp index de4da61a4d..5a86171605 100644 --- a/src/shammodels/ramses/src/pyRamsesModel.cpp +++ b/src/shammodels/ramses/src/pyRamsesModel.cpp @@ -199,6 +199,14 @@ namespace shammodels::basegodunov { }, py::kw_only(), py::arg("crit_mass")) + .def( + "set_amr_mode_pseudo_gradient_based", + [](TConfig &self, Tscal error_min, Tscal error_max) { + self.amr_mode.set_refine_pseudo_gradient_based(error_min, error_max); + }, + py::kw_only(), + py::arg("error_min"), + py::arg("error_max")) .def( "set_gravity_mode_no_gravity", [](TConfig &self) { From 9c01cfcb458a906201d49a7ad67c158407e5afac Mon Sep 17 00:00:00 2001 From: Leodasce Sewanou Date: Fri, 20 Mar 2026 15:27:46 +0100 Subject: [PATCH 4/5] sod and sedov_blast test --- examples/TO_MIGRATE/ramses/sedo_amr.py | 110 ++++++++++++++++++ .../TO_MIGRATE/ramses/sod_amr_reflective.py | 12 +- 2 files changed, 117 insertions(+), 5 deletions(-) create mode 100644 examples/TO_MIGRATE/ramses/sedo_amr.py diff --git a/examples/TO_MIGRATE/ramses/sedo_amr.py b/examples/TO_MIGRATE/ramses/sedo_amr.py new file mode 100644 index 0000000000..5f89198a69 --- /dev/null +++ b/examples/TO_MIGRATE/ramses/sedo_amr.py @@ -0,0 +1,110 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np + +import shamrock + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + +multx = 1 +multy = 1 +multz = 1 +max_amr_lev = 3 +cell_size = 2 << max_amr_lev # refinement is limited to cell_size = 2 +base = 16 + +cfg = model.gen_default_config() +scale_fact = 1 / (cell_size * base * multx) +cfg.set_scale_factor(scale_fact) + +center = (0.5 * base * scale_fact, 0.5 * base * scale_fact, 0.5 * base * scale_fact) +xc, yc, zc = center +Rstart = 1.0 / (2.0 * base) + 1e-4 +# Rstart = 0.01 +gamma = 5.0 / 3.0 + +cfg.set_eos_gamma(gamma) +cfg.set_Csafe(0.3) +cfg.set_boundary_condition("x", "periodic") +cfg.set_boundary_condition("y", "periodic") +cfg.set_boundary_condition("z", "periodic") +cfg.set_riemann_solver_hllc() +cfg.set_slope_lim_minmod() +cfg.set_face_time_interpolation(True) + +err_min = 0.30 +err_max = 0.10 + +cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max) + +model.set_solver_config(cfg) + + +model.init_scheduler(int(1e7), 1) +model.make_base_grid( + (0, 0, 0), (cell_size, cell_size, cell_size), (base * multx, base * multy, base * multz) +) + + +def rho_map(rmin, rmax): + x_min, y_min, z_min = rmin + x_max, y_max, z_max = rmax + x = (x_min + x_max) * 0.5 - 0.5 + y = (y_min + y_max) * 0.5 - 0.5 + z = (z_min + z_max) * 0.5 - 0.5 + r = np.sqrt(x * x + y * y + z * z) + # if r < Rstart: + # return 1. + # else: + # return 1.2 + + return 1.0 + + +def rhoe_map(rmin, rmax): + x_min, y_min, z_min = rmin + x_max, y_max, z_max = rmax + x = (x_min + x_max) * 0.5 - 0.5 + y = (y_min + y_max) * 0.5 - 0.5 + z = (z_min + z_max) * 0.5 - 0.5 + r = np.sqrt(x * x + y * y + z * z) + if r < Rstart: + return 10.0 / (gamma - 1.0) + else: + return 1e-5 / (gamma - 1.0) + + +def rhovel_map(rmin, rmax): + return (0.0, 0.0, 0.0) + + +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) + +tmax = 0.2 + + +dt = 0 +t = 0 +freq = 1 +dX0 = [] +for i in range(10000): + next_dt = model.evolve_once_override_time(t, dt) + + t += dt + dt = next_dt + + if i % freq == 0: + model.dump_vtk(f"test{i:04d}_ref_new.vtk") + + if tmax < t + next_dt: + dt = tmax - t + if t == tmax: + model.dump_vtk(f"test{i:04d}_ref_new.vtk") + break diff --git a/examples/TO_MIGRATE/ramses/sod_amr_reflective.py b/examples/TO_MIGRATE/ramses/sod_amr_reflective.py index de8a1b065b..1ea343681a 100644 --- a/examples/TO_MIGRATE/ramses/sod_amr_reflective.py +++ b/examples/TO_MIGRATE/ramses/sod_amr_reflective.py @@ -14,7 +14,7 @@ multx = 1 multy = 1 multz = 1 -max_amr_lev = 2 +max_amr_lev = 3 cell_size = 2 << max_amr_lev # refinement is limited to cell_size = 2 base = 16 @@ -135,9 +135,9 @@ def convert_to_cell_coords(dic): dt = 0 t = 0 -freq = 1 +freq = 100 dX0 = [] -for i in range(5): +for i in range(10000): next_dt = model.evolve_once_override_time(t, dt) if i == 0: dic0 = convert_to_cell_coords(ctx.collect_data()) @@ -195,7 +195,9 @@ def convert_to_cell_coords(dic): ax1 = plt.gca() ax2 = ax1.twinx() - l = -np.log2(dX / max(dX0.max(), dX.max())) + l_0 = np.log2(base * 2) + + l = -np.log2(dX / max(dX0.max(), dX.max())) + l_0 ax1.scatter(X, rho, rasterized=True, s=12 * np.ones(X.shape), label="rho") ax1.scatter(X, vx, rasterized=True, s=12 * np.ones(X.shape), label="v") @@ -233,5 +235,5 @@ def convert_to_cell_coords(dic): ax1.plot(arr_x, arr_P, ls="--", lw=2.0, color="black") ax2.set_ylabel("AMR level") plt.title(f"Threshold = {err_max}, derefinement factor = {err_min}") - plt.savefig("sod_tube.png") + plt.savefig("sod_tube_3_1_baryonic_density.png") ####### From 0f590649e16f0fc5a2feb9812902cf312853965e Mon Sep 17 00:00:00 2001 From: Leodasce Sewanou Date: Fri, 20 Mar 2026 15:46:59 +0100 Subject: [PATCH 5/5] updates --- examples/TO_MIGRATE/ramses/sedo_amr.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/examples/TO_MIGRATE/ramses/sedo_amr.py b/examples/TO_MIGRATE/ramses/sedo_amr.py index 5f89198a69..2fe8c7eeee 100644 --- a/examples/TO_MIGRATE/ramses/sedo_amr.py +++ b/examples/TO_MIGRATE/ramses/sedo_amr.py @@ -23,9 +23,7 @@ cfg.set_scale_factor(scale_fact) center = (0.5 * base * scale_fact, 0.5 * base * scale_fact, 0.5 * base * scale_fact) -xc, yc, zc = center Rstart = 1.0 / (2.0 * base) + 1e-4 -# Rstart = 0.01 gamma = 5.0 / 3.0 cfg.set_eos_gamma(gamma) @@ -52,17 +50,6 @@ def rho_map(rmin, rmax): - x_min, y_min, z_min = rmin - x_max, y_max, z_max = rmax - x = (x_min + x_max) * 0.5 - 0.5 - y = (y_min + y_max) * 0.5 - 0.5 - z = (z_min + z_max) * 0.5 - 0.5 - r = np.sqrt(x * x + y * y + z * z) - # if r < Rstart: - # return 1. - # else: - # return 1.2 - return 1.0