AMRGridRefinementHandler.cpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//

 
#include "shambase/memory.hpp"
#include "shamalgs/details/algorithm/algorithm.hpp"
#include "shamcomm/logs.hpp"
#include "shammodels/common/amr/NeighGraph.hpp"
#include "shammodels/ramses/config/enum_SlopeMode.hpp"
#include "shammodels/ramses/modules/AMRGridRefinementHandler.hpp"
#include "shammodels/ramses/modules/AMRSortBlocks.hpp"
#include "shammodels/ramses/modules/SlopeLimitedGradientUtilities.hpp"
#include <stdexcept>
#include <variant>
 
template<class Tvec, class TgridVec>
template<class UserAcc, class... T>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    gen_refine_block_changes_old(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &dd_refine_list,
        shambase::DistributedData<sham::DeviceBuffer<u32>> &dd_derefine_list,
        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();
 
    scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
        u64 id_patch = cur_p.id_patch;
 
        // create the refine and derefine flags buffers
        u32 obj_cnt = pdat.get_obj_cnt();
 
        sham::DeviceBuffer<u32> refine_flags(obj_cnt, dev_sched);
        sham::DeviceBuffer<u32> derefine_flags(obj_cnt, dev_sched);
 
        {
            sham::EventList depends_list;
 
            UserAcc uacc(depends_list, 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);
 
            // fill in the flags
            auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
                cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) {
                    bool flag_refine   = false;
                    bool flag_derefine = false;
                    uacc.refine_criterion(gid.get_linear_id(), uacc, flag_refine, flag_derefine);
 
                    // This is just a safe guard to avoid this nonsensicall case
                    if (flag_refine && flag_derefine) {
                        flag_derefine = false;
                    }
 
                    refine_acc[gid]   = (flag_refine) ? 1 : 0;
                    derefine_acc[gid] = (flag_derefine) ? 1 : 0;
                });
            });
 
            sham::EventList resulting_events;
            resulting_events.add_event(e);
 
            refine_flags.complete_event_state(resulting_events);
            derefine_flags.complete_event_state(resulting_events);
 
            uacc.finalize(resulting_events, id_patch, cur_p, pdat, args...);
        }
 
        sham::DeviceBuffer<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
        sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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<BlockCoord, split_count> 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);
 
                } else {
                    do_merge = false;
                }
 
                acc_merge_flag[gid] = do_merge;
            });
        });
 
        buf_cell_min.complete_event_state(e);
        buf_cell_max.complete_event_state(e);
        derefine_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
        dd_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);
 
        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
        dd_derefine_list.add_obj(id_patch, std::move(buf_derefine));
    });
 
    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<class Tvec, class TgridVec>
template<class UserAcc>
bool shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_refine_grid_old(shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_refine_list) {
 
    using namespace shamrock::patch;
 
    u64 sum_block_count = 0;
 
    bool new_cell_were_added = false;
 
    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<u32> &refine_flags = dd_refine_list.get(id_patch);
 
        if (refine_flags.get_size() > 0) {
 
            // alloc memory for the new blocks to be created
            pdat.expand(refine_flags.get_size() * (split_count - 1));
 
            sham::DeviceBuffer<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
            sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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);
 
            // Refine the block (set the positions) and fill the corresponding fields
            auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
                u32 start_index_push = old_obj_cnt;
 
                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();
 
                    u32 idx_to_refine = index_to_ref[tid];
 
                    // gen splits coordinates
                    BlockCoord cur_block{
                        block_bound_low[idx_to_refine], block_bound_high[idx_to_refine]};
 
                    std::array<BlockCoord, split_count> block_coords
                        = BlockCoord::get_split(cur_block.bmin, cur_block.bmax);
 
                    // generate index for the refined blocks
                    std::array<u32, split_count> 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;
                    }
 
                    // 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;
                    }
 
                    // user lambda to fill the fields
                    uacc.apply_refine(idx_to_refine, cur_block, blocks_ids, block_coords, uacc);
                });
            });
 
            sham::EventList resulting_events{e};
 
            buf_cell_min.complete_event_state(resulting_events);
            buf_cell_max.complete_event_state(resulting_events);
 
            uacc.finalize(resulting_events, pdat);
 
            refine_flags.complete_event_state(resulting_events);
        }
 
        sum_block_count += pdat.get_obj_cnt();
        new_cell_were_added = new_cell_were_added || refine_flags.get_size() > 0;
    });
 
    logger::info_ln("AMRGrid", "process block count =", sum_block_count);
 
    return new_cell_were_added;
}
 
template<class Tvec, class TgridVec>
template<class UserAcc>
bool shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_derefine_grid_old(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_derefine_list) {
 
    using namespace shamrock::patch;
 
    bool cell_were_removed = 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) {
        u32 old_obj_cnt = pdat.get_obj_cnt();
 
        sham::DeviceBuffer<u32> &derefine_flags = dd_derefine_list.get(id_patch);
 
        if (derefine_flags.get_size() > 0) {
 
            // init flag table
            sham::DeviceBuffer<u32> keep_block_flag(old_obj_cnt, dev_sched);
            keep_block_flag.fill(1);
 
            sham::DeviceBuffer<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
            sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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 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();
 
                    u32 idx_to_derefine = index_to_deref[gid];
 
                    // compute old block indexes
                    std::array<u32, split_count> old_indexes;
#pragma unroll
                    for (u32 pid = 0; pid < split_count; pid++) {
                        old_indexes[pid] = idx_to_derefine + pid;
                    }
 
                    // load block coords
                    std::array<BlockCoord, split_count> 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]]};
                    }
 
                    // 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;
 
// flag the old blocks for removal
#pragma unroll
                    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);
                });
            });
 
            sham::EventList resulting_events{e};
 
            buf_cell_min.complete_event_state(resulting_events);
            buf_cell_max.complete_event_state(resulting_events);
 
            uacc.finalize(resulting_events, pdat);
 
            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);
 
            shamlog_debug_ln(
                "AMR Grid",
                "patch",
                id_patch,
                "derefine block count ",
                old_obj_cnt,
                "->",
                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;
        }
    });
 
    return cell_were_removed;
}
 
template<class Tvec, class TgridVec>
template<class UserAccCrit, class UserAccSplit, class UserAccMerge>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_update_refinement_old() {
 
    // Ensure that the blocks are sorted before refinement
    AMRSortBlocks block_sorter(context, solver_config, storage);
    block_sorter.reorder_amr_blocks();
 
    // get refine and derefine list
    shambase::DistributedData<sham::DeviceBuffer<u32>> dd_refine_list;
    shambase::DistributedData<sham::DeviceBuffer<u32>> dd_derefine_list;
 
    gen_refine_block_changes<UserAccCrit>(dd_refine_list, dd_derefine_list);

    // Note that this only add new blocks at the end of the patchdata
    internal_refine_grid<UserAccSplit>(std::move(dd_refine_list));

    // Note that this will perform the merge then remove the old blocks
    // This is ok to call straight after the refine without edditing the index list in derefine_list
    // since no permutations were applied in internal_refine_grid and no cells can be both refined
    // and derefined in the same pass
    internal_derefine_grid<UserAccMerge>(std::move(dd_derefine_list));
}
 
template<class Tvec, class TgridVec>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    update_refinement_old() {
 
    class RefineCritBlock {
        public:
        const TgridVec *block_low_bound;
        const TgridVec *block_high_bound;
        const Tscal *block_density_field;
 
        Tscal one_over_Nside = 1. / AMRBlock::Nside;
 
        Tscal dxfact;
        Tscal wanted_mass;
 
        RefineCritBlock(
            sham::EventList &depends_list,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal dxfact,
            Tscal wanted_mass)
            : dxfact(dxfact), wanted_mass(wanted_mass) {
 
            block_low_bound  = pdat.get_field<TgridVec>(0).get_buf().get_read_access(depends_list);
            block_high_bound = pdat.get_field<TgridVec>(1).get_buf().get_read_access(depends_list);
            block_density_field = pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                                      .get_buf()
                                      .get_read_access(depends_list);
        }
 
        void finalize(
            sham::EventList &resulting_events,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal dxfact,
            Tscal wanted_mass) {
 
            sham::DeviceBuffer<i64_3> &buf_cell_low_bound  = pdat.get_field<i64_3>(0).get_buf();
            sham::DeviceBuffer<i64_3> &buf_cell_high_bound = pdat.get_field<i64_3>(1).get_buf();
 
            buf_cell_low_bound.complete_event_state(resulting_events);
            buf_cell_high_bound.complete_event_state(resulting_events);
            pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                .get_buf()
                .complete_event_state(resulting_events);
        }
 
        void refine_criterion(
            u32 block_id, RefineCritBlock 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];
 
            Tvec lower_flt = low_bound.template convert<Tscal>() * dxfact;
            Tvec upper_flt = high_bound.template convert<Tscal>() * dxfact;
 
            Tvec block_cell_size = (upper_flt - lower_flt) * one_over_Nside;
 
            Tscal sum_mass = 0;
            for (u32 i = 0; i < AMRBlock::block_size; i++) {
                sum_mass += acc.block_density_field[i + block_id * AMRBlock::block_size];
            }
            sum_mass *= block_cell_size.x() * block_cell_size.y() * block_cell_size.z();
 
            if (sum_mass > wanted_mass * 8) {
                should_refine   = true;
                should_derefine = false;
            } else if (sum_mass < wanted_mass) {
                should_refine   = false;
                should_derefine = true;
            } else {
                should_refine   = false;
                should_derefine = false;
            }
 
            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);
        }
    };
 
    class RefineCellAccessor {
        public:
        f64 *rho;
        f64_3 *rho_vel;
        f64 *rhoE;
 
        RefineCellAccessor(sham::EventList &depends_list, shamrock::patch::PatchDataLayer &pdat) {
 
            rho     = pdat.get_field<f64>(2).get_buf().get_write_access(depends_list);
            rho_vel = pdat.get_field<f64_3>(3).get_buf().get_write_access(depends_list);
            rhoE    = pdat.get_field<f64>(4).get_buf().get_write_access(depends_list);
        }
 
        void finalize(sham::EventList &resulting_events, shamrock::patch::PatchDataLayer &pdat) {
            pdat.get_field<f64>(2).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64_3>(3).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64>(4).get_buf().complete_event_state(resulting_events);
        }
 
        void apply_refine(
            u32 cur_idx,
            BlockCoord cur_coords,
            std::array<u32, 8> new_blocks,
            std::array<BlockCoord, 8> new_block_coords,
            RefineCellAccessor acc) const {
 
            auto get_coord_ref = [](u32 i) -> std::array<u32, dim> {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    const u32 tmp = i >> NsideBlockPow;
                    return {i % Nside, (tmp) % Nside, (tmp) >> NsideBlockPow};
                }
            };
 
            auto get_index_block = [](std::array<u32, dim> coord) -> u32 {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    return coord[0] + Nside * coord[1] + Nside * Nside * coord[2];
                }
            };
 
            auto get_gid_write = [&](std::array<u32, dim> &glid) -> u32 {
                std::array<u32, dim> 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
                       + AMRBlock::get_index(
                           {glid[0] % AMRBlock::Nside,
                            glid[1] % AMRBlock::Nside,
                            glid[2] % AMRBlock::Nside});
            };
 
            std::array<f64, AMRBlock::block_size> old_rho_block;
            std::array<f64_3, AMRBlock::block_size> old_rho_vel_block;
            std::array<f64, AMRBlock::block_size> old_rhoE_block;
 
            // save old block
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz]         = get_coord_ref(loc_id);
                u32 old_cell_idx          = cur_idx * AMRBlock::block_size + loc_id;
                old_rho_block[loc_id]     = acc.rho[old_cell_idx];
                old_rho_vel_block[loc_id] = acc.rho_vel[old_cell_idx];
                old_rhoE_block[loc_id]    = acc.rhoE[old_cell_idx];
            }
 
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz] = get_coord_ref(loc_id);
                u32 old_cell_idx  = cur_idx * AMRBlock::block_size + loc_id;
 
                Tscal rho_block    = old_rho_block[loc_id];
                Tvec rho_vel_block = old_rho_vel_block[loc_id];
                Tscal rhoE_block   = old_rhoE_block[loc_id];
                for (u32 subdiv_lid = 0; subdiv_lid < 8; subdiv_lid++) {
 
                    auto [sx, sy, sz] = get_coord_ref(subdiv_lid);
 
                    std::array<u32, 3> 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;
                }
            }
        }
 
        void apply_derefine(
            std::array<u32, 8> old_blocks,
            std::array<BlockCoord, 8> old_coords,
            u32 new_cell,
            BlockCoord new_coord,
 
            RefineCellAccessor acc) const {
 
            std::array<f64, AMRBlock::block_size> rho_block;
            std::array<f64_3, AMRBlock::block_size> rho_vel_block;
            std::array<f64, AMRBlock::block_size> rhoE_block;
 
            for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) {
                rho_block[cell_id]     = {};
                rho_vel_block[cell_id] = {};
                rhoE_block[cell_id]    = {};
            }
 
            for (u32 pid = 0; pid < 8; 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];
                }
            }
 
            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;
            }
 
            for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) {
                u32 newcell_idx          = new_cell * AMRBlock::block_size + cell_id;
                acc.rho[newcell_idx]     = rho_block[cell_id];
                acc.rho_vel[newcell_idx] = rho_vel_block[cell_id];
                acc.rhoE[newcell_idx]    = rhoE_block[cell_id];
            }
        }
    };
 
    using AMRmode_None         = typename AMRMode<Tvec, TgridVec>::None;
    using AMRmode_DensityBased = typename AMRMode<Tvec, TgridVec>::DensityBased;
 
    bool has_cell_order_changed = false;
 
    if (AMRmode_None *cfg = std::get_if<AMRmode_None>(&solver_config.amr_mode.config)) {
        // no refinment here turn around there is nothing to see
    } else if (
        AMRmode_DensityBased *cfg
        = std::get_if<AMRmode_DensityBased>(&solver_config.amr_mode.config)) {
        Tscal dxfact(solver_config.grid_coord_to_pos_fact);
 
        // get refine and derefine list
        shambase::DistributedData<sham::DeviceBuffer<u32>> refine_list;
        shambase::DistributedData<sham::DeviceBuffer<u32>> derefine_list;
 
        gen_refine_block_changes_old<RefineCritBlock>(
            refine_list, derefine_list, dxfact, cfg->crit_mass);

        // Note that this only add new blocks at the end of the patchdata
        bool change_refine = internal_refine_grid_old<RefineCellAccessor>(std::move(refine_list));

        // Note that this will perform the merge then remove the old blocks
        // This is ok to call straight after the refine without edditing the index list in
        // derefine_list since no permutations were applied in internal_refine_grid and no cells can
        // be both refined and derefined in the same pass
        bool change_derefine
            = internal_derefine_grid_old<RefineCellAccessor>(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);
        block_sorter.reorder_amr_blocks();
    }
}
 
template<class Tvec, class TgridVec>
template<class UserAcc, class... T>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    gen_refine_block_changes_new(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &dd_refine_flags,
        shambase::DistributedData<sham::DeviceBuffer<u32>> &dd_derefine_flags,
        T &&...args) {
 
    sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
    auto dev_sched       = shamsys::instance::get_compute_scheduler_ptr();
 
    scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
        u64 id_patch = cur_p.id_patch;
 
        // create the refine and derefine flags buffers
        u32 obj_cnt = pdat.get_obj_cnt();
 
        sham::DeviceBuffer<u32> refine_flags(obj_cnt, dev_sched);
        sham::DeviceBuffer<u32> derefine_flags(obj_cnt, dev_sched);
 
        {
            sham::EventList depends_list;
 
            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);
 
            // fill in the flags
            auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
                cgh.parallel_for(sycl::range<1>(obj_cnt), [=](sycl::item<1> gid) {
                    bool flag_refine   = false;
                    bool flag_derefine = false;
                    uacc.refine_criterion_new(
                        gid.get_linear_id(), uacc, flag_refine, flag_derefine);
 
                    // This is just a safe guard to avoid this nonsensicall case
                    if (flag_refine && flag_derefine) {
                        flag_derefine = false;
                    }
 
                    refine_acc[gid]   = (flag_refine) ? 1 : 0;
                    derefine_acc[gid] = (flag_derefine) ? 1 : 0;
                });
            });
 
            sham::EventList resulting_events;
            resulting_events.add_event(e);
 
            refine_flags.complete_event_state(resulting_events);
            derefine_flags.complete_event_state(resulting_events);
 
            uacc.finalize_new(resulting_events, storage, id_patch, cur_p, pdat, args...);
        }
 
        dd_refine_flags.add_obj(id_patch, std::move(refine_flags));
        dd_derefine_flags.add_obj(id_patch, std::move(derefine_flags));
    });
}

template<class Tvec, class TgridVec>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    enforce_two_to_one_refinement_new(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_refine_flags) {
 
    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<u32> &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<TgridUint> &buf_amr_block_levels
            = shambase::get_check_ref(storage.amr_block_levels).get_buf(id_patch);
 
        std::shared_ptr<sham::DeviceScheduler> dev_sched
            = shamsys::instance::get_compute_scheduler_ptr();
 
        sham::DeviceBuffer<u32> changed_buf(1, dev_sched);
 
        for (u32 pass = 0; pass < 100; pass++) {
 
            changed_buf.set_val_at_idx(0, 0);
 
            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_changed    = changed_buf.get_write_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();
 
                    u32 cur_ref_flag = acc_ref_flags[block_id];
 
                    if (!cur_ref_flag)
                        return;
 
                    auto cur_block_level = acc_amr_levels[block_id];
 
                    auto check_2To1_ref = [&](u32 nid) {
                        if (nid >= obj_cnt)
                            return;
 
                        // 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];
 
                        auto cur_future = cur_block_level + (cur_ref_flag ? 1 : 0);
 
                        auto neigh_future = neigh_block_level + (neigh_ref_flag ? 1 : 0);
 
                        if (cur_ref_flag && (cur_future > neigh_future + 1)) {
 
                            if (!neigh_ref_flag) {
                                sycl::atomic_ref<
                                    u32,
                                    sycl::memory_order::relaxed,
                                    sycl::memory_scope::system>
                                    atomic_neigh_flag(acc_ref_flags[nid]);
                                atomic_neigh_flag.exchange(1);
 
                                sycl::atomic_ref<
                                    u32,
                                    sycl::memory_order::relaxed,
                                    sycl::memory_scope::system>
                                    atomic_changed(acc_changed[0]);
                                atomic_changed.exchange(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);
                });
            });
            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);
 
            changed_buf.complete_event_state(e);
 
            patch_refine_flags.complete_event_state(e);
 
            e.wait();
 
            if (changed_buf.get_val_at_idx(0) == 0) {
                logger::raw_ln("Refinement 2:1 balance converged in ", pass + 1, " sweeps");
                break;
            }
        }

        // refinement
 
        // perform stream compactions on the refinement flags
        auto dev_buf_ref
            = shamalgs::numeric::stream_compact(dev_sched, patch_refine_flags, obj_cnt);
 
        shamlog_debug_ln(
            "AMRGrid", "patch ", id_patch, dev_buf_ref.get_size(), "marked for refinement + 2:1");
    });
}

 
template<class Tvec, class TgridVec>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    enforce_two_to_one_derefinement_new(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_derefine_flags,
        shambase::DistributedData<sham::DeviceBuffer<u32>> &&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<u32> &patch_derefine_flag = dd_derefine_flags.get(id_patch);
        sham::DeviceBuffer<u32> &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<TgridUint> &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<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
        sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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<BlockCoord, split_count> 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<u32, dim> {
                        constexpr u32 NsideBlockPow = 1;
                        constexpr u32 Nside         = 1U << NsideBlockPow;
                        constexpr u32 side_size     = Nside;
                        constexpr u32 block_size    = shambase::pow_constexpr<dim>(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<BlockCoord, split_count> {
                        std::array<BlockCoord, split_count> ret;
                        auto bmin                   = target_block.bmin;
                        auto bmax                   = target_block.bmax;
                        auto split                  = bmin + (bmax - bmin) / 2;
                        std::array<TgridVec, 3> 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<BlockCoord, split_count> 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
        //     ///////////////////////////////////////////////////////////////////////////////////
 
        std::shared_ptr<sham::DeviceScheduler> dev_sched
            = shamsys::instance::get_compute_scheduler_ptr();
 
        sham::DeviceBuffer<u32> changed_buf(1, dev_sched);
 
        // copy old deref buffer to avoid race condition
        sham::DeviceBuffer<u32> patch_derefine_flag_old(obj_cnt, dev_sched);
        patch_derefine_flag.copy_range(0, obj_cnt, patch_derefine_flag_old);
 
        //
        sham::DeviceBuffer<u32> patch_derefine_flag_new(obj_cnt, dev_sched);
 
        for (int it = 0; it < 100; it++) {
            changed_buf.set_val_at_idx(0, 0);
 
            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_flag   = patch_refine_flag.get_read_access(depend_list);
            auto acc_changed    = changed_buf.get_write_access(depend_list);
 
            auto acc_deref_old = patch_derefine_flag_old.get_read_access(depend_list);
            auto acc_deref_new = patch_derefine_flag_new.get_write_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 old_flag = acc_deref_old[lid];
                    auto new_flag = old_flag;
 
                    auto check_2To1_der = [&](u32 nid) {
                        if (nid < obj_cnt)
 
                        {
 
                            auto neigh_future = acc_amr_levels[nid] + (acc_ref_flag[nid] ? 1 : 0)
                                                - (acc_deref_old[nid] ? 1 : 0);
 
                            auto my_future = acc_amr_levels[lid] - 1;
 
                            if (neigh_future > my_future + 1) {
                                new_flag = 0;
                            }
                        }
                    };
 
                    if (old_flag) {
 
                        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);
                        }
 
                        if (old_flag != new_flag) {
                            sycl::atomic_ref<
                                u32,
                                sycl::memory_order::relaxed,
                                sycl::memory_scope::system>
                                atomic_changed(acc_changed[0]);
                            atomic_changed.exchange(1);
                        }
                    }
 
                    acc_deref_new[lid] = new_flag;
                });
            });
            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_refine_flag.complete_event_state(e_2to1);
            changed_buf.complete_event_state(e_2to1);
 
            patch_derefine_flag_old.complete_event_state(e_2to1);
            patch_derefine_flag_new.complete_event_state(e_2to1);
            e_2to1.wait();
 
            std::swap(patch_derefine_flag_old, patch_derefine_flag_new);
 
            if (changed_buf.get_val_at_idx(0) == 0) {
            logger:
                shamcomm::logs::raw_ln(
                    "Derefinement 2:1 balance converge in \t ", it + 1, "\t sweeps \n\n");
                break;
            }
        }
 
        // copy back to ..
        patch_derefine_flag_old.copy_range(0, obj_cnt, patch_derefine_flag);
 
        // ////////////////////////////////////////////////////////////////////////////////
        // // 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<class Tvec, class TgridVec>
template<class UserAcc>
bool shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_refine_grid_new(shambase::DistributedData<sham::DeviceBuffer<u32>> &&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) {
        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
 
            logger::raw_ln(
                "Will refine \t", stream_compaction_result.get_size(), " \t blocks \n\n");
 
            pdat.expand(static_cast<u32>(stream_compaction_result.get_size()) * (split_count - 1));
            sham::DeviceBuffer<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
            sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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, 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
            auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
                u32 start_index_push = old_obj_cnt;
 
                constexpr u32 new_splits = split_count - 1;
 
                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[gid];
 
                        // gen splits coordinates
                        BlockCoord cur_block{
                            block_bound_low[idx_to_refine], block_bound_high[idx_to_refine]};
 
                        std::array<BlockCoord, split_count> block_coords
                            = BlockCoord::get_split(cur_block.bmin, cur_block.bmax);
 
                        // generate index for the refined blocks
                        std::array<u32, split_count> 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;
                        }
 
                    // 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;
                        }
 
                        // user lambda to fill the fields
                        uacc.apply_refine_new(
                            idx_to_refine, cur_block, blocks_ids, block_coords, uacc);
                    });
            });
 
            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_new(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();
        new_cell_were_added = new_cell_were_added || (stream_compaction_result.get_size() > 0);
    });
 
    logger::info_ln("AMRGrid", "process block count =", sum_block_count);
 
    return new_cell_were_added;
}
 
template<class Tvec, class TgridVec>
template<class UserAcc>
bool shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_derefine_grid_new(
        shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_derefine_flags) {
 
    using namespace shamrock::patch;
 
    bool cell_were_removed = 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) {
        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
            sycl::buffer<u32> keep_block_flag
                = shamalgs::algorithm::gen_buffer_device(q.q, old_obj_cnt, [](u32 i) -> u32 {
                      return 1;
                  });
 
            sham::DeviceBuffer<TgridVec> &buf_cell_min = pdat.get_field_buf_ref<TgridVec>(0);
            sham::DeviceBuffer<TgridVec> &buf_cell_max = pdat.get_field_buf_ref<TgridVec>(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, storage, id_patch, pdat);
            auto index_to_deref = stream_compact_results.get_read_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();
 
                        u32 idx_to_derefine = index_to_deref[gid];
 
                        // compute old block indexes
                        std::array<u32, split_count> old_indexes;
#pragma unroll
                        for (u32 pid = 0; pid < split_count; pid++) {
                            old_indexes[pid] = idx_to_derefine + pid;
                        }
 
                        // load block coords
                        std::array<BlockCoord, split_count> 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]]};
                        }
 
                        // 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;
 
// flag the old blocks for removal
#pragma unroll
                        for (u32 pid = 1; pid < split_count; pid++) {
                            flag_keep[idx_to_derefine + pid] = 0;
                        }
 
                        // user lambda to fill the fields
 
                        uacc.apply_derefine_new(
                            old_indexes, block_coords, idx_to_derefine, merged_block_coord, uacc);
                    });
            });
 
            sham::EventList resulting_events;
 
            buf_cell_min.complete_event_state(resulting_events);
            buf_cell_max.complete_event_state(resulting_events);
            uacc.finalize_new(resulting_events, storage, id_patch, pdat);
 
            stream_compact_results.complete_event_state(resulting_events);
 
            // 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 - len,
                "new block count = ",
                len);
 
            if (!opt_buf) {
                throw std::runtime_error("opt buf must contain something at this point");
            }
 
            // 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;
        }
    });
 
    return cell_were_removed;
}
 
template<class Tvec, class TgridVec>
template<class UserAccCrit, class UserAccSplit, class UserAccMerge>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    internal_update_refinement_new() {
 
    // Ensure that the blocks are sorted before refinement
    AMRSortBlocks block_sorter(context, solver_config, storage);
    block_sorter.reorder_amr_blocks();
 
    // get refine and derefine list
    shambase::DistributedData<sham::DeviceBuffer<u32>> dd_refine_list;
    shambase::DistributedData<sham::DeviceBuffer<u32>> dd_derefine_list;
 
    gen_refine_block_changes_new<UserAccCrit>(dd_refine_list, dd_derefine_list);

    // Note that this only add new blocks at the end of the patchdata
    internal_refine_grid_new<UserAccSplit>(std::move(dd_refine_list));

    // Note that this will perform the merge then remove the old blocks
    // This is ok to call straight after the refine without edditing the index list in derefine_list
    // since no permutations were applied in internal_refine_grid_new and no cells can be both
    // refined and derefined in the same pass
    internal_derefine_grid_new<UserAccMerge>(std::move(dd_derefine_list));
}
 
template<class Tvec, class TgridVec>
void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>::
    update_refinement_new() {
 
    class RefineCritBlock {
        public:
        const TgridVec *block_low_bound;
        const TgridVec *block_high_bound;
        const Tscal *block_density_field;
 
        Tscal one_over_Nside = 1. / AMRBlock::Nside;
 
        Tscal dxfact;
        Tscal wanted_mass;
 
        RefineCritBlock(
            sham::EventList &depends_list,
            Storage &storage,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal dxfact,
            Tscal wanted_mass)
            : dxfact(dxfact), wanted_mass(wanted_mass) {
 
            block_low_bound  = pdat.get_field<TgridVec>(0).get_buf().get_read_access(depends_list);
            block_high_bound = pdat.get_field<TgridVec>(1).get_buf().get_read_access(depends_list);
            block_density_field = pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                                      .get_buf()
                                      .get_read_access(depends_list);
        }
 
        void finalize_new(
            sham::EventList &resulting_events,
            Storage &storage,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal dxfact,
            Tscal wanted_mass) {
 
            pdat.get_field<TgridVec>(0).get_buf().complete_event_state(resulting_events);
            pdat.get_field<TgridVec>(1).get_buf().complete_event_state(resulting_events);
            pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                .get_buf()
                .complete_event_state(resulting_events);
        }
 
        void refine_criterion_new(
            u32 block_id, RefineCritBlock 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];
 
            Tvec lower_flt = low_bound.template convert<Tscal>() * dxfact;
            Tvec upper_flt = high_bound.template convert<Tscal>() * dxfact;
 
            Tvec block_cell_size = (upper_flt - lower_flt) * one_over_Nside;
 
            Tscal sum_mass = 0;
            for (u32 i = 0; i < AMRBlock::block_size; i++) {
                sum_mass += acc.block_density_field[i + block_id * AMRBlock::block_size];
            }
            sum_mass *= block_cell_size.x() * block_cell_size.y() * block_cell_size.z();
 
            if (sum_mass > wanted_mass * 8) {
                should_refine   = true;
                should_derefine = false;
            } else if (sum_mass < wanted_mass) {
                should_refine   = false;
                should_derefine = true;
            } else {
                should_refine   = false;
                should_derefine = false;
            }
 
            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);
        }
    };
 
    class RefineCellAccessor {
        public:
        f64 *rho;
        f64_3 *rho_vel;
        f64 *rhoE;
        u64 p_id;
        // f64* cell_sizes;
 
        // 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<f64>(2).get_buf().get_write_access(depends_list);
            rho_vel = pdat.get_field<f64_3>(3).get_buf().get_write_access(depends_list);
            rhoE    = pdat.get_field<f64>(4).get_buf().get_write_access(depends_list);
            // cell_sizes = shambase::get_check_ref(storage.block_cell_sizes)
            //                 .get_buf(id_patch)
            //                 .get_write_access(depends_list);
        }
 
        void finalize_new(
            sham::EventList &resulting_events,
            Storage &storage,
            u64 &id_patch,
            shamrock::patch::PatchDataLayer &pdat) {
            pdat.get_field<f64>(2).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64_3>(3).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64>(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);
 
            // shambase::get_check_ref(storage.block_cell_sizes)
            //     .get_buf(id_patch)
            //     .complete_event_state(resulting_events);
        }
 
        void apply_refine_new(
            u32 cur_idx,
            BlockCoord cur_coords,
            std::array<u32, 8> new_blocks,
            std::array<BlockCoord, 8> new_block_coords,
            RefineCellAccessor acc) const {
 
            auto get_coord_ref = [](u32 i) -> std::array<u32, dim> {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    const u32 tmp = i >> NsideBlockPow;
                    return {i % Nside, (tmp) % Nside, (tmp) >> NsideBlockPow};
                }
            };
 
            auto get_index_block = [](std::array<u32, dim> coord) -> u32 {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    return coord[0] + Nside * coord[1] + Nside * Nside * coord[2];
                }
            };
 
            auto get_gid_write = [&](std::array<u32, dim> &glid) -> u32 {
                // First, get the block id (it's the block to be refine) in wich the new cell glid
                // is located.
                std::array<u32, dim> bid
                    = {glid[0] >> AMRBlock::NsideBlockPow,
                       glid[1] >> AMRBlock::NsideBlockPow,
                       glid[2] >> AMRBlock::NsideBlockPow};
 
                // 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,
                            glid[2] % AMRBlock::Nside});
            };
 
            std::array<f64, AMRBlock::block_size> old_rho_block;
            std::array<f64_3, AMRBlock::block_size> old_rho_vel_block;
            std::array<f64, AMRBlock::block_size> old_rhoE_block;
 
            // save old block
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz]         = get_coord_ref(loc_id);
                u32 old_cell_idx          = cur_idx * AMRBlock::block_size + loc_id;
                old_rho_block[loc_id]     = acc.rho[old_cell_idx];
                old_rho_vel_block[loc_id] = acc.rho_vel[old_cell_idx];
                old_rhoE_block[loc_id]    = acc.rhoE[old_cell_idx];
            }
 
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz] = get_coord_ref(loc_id);
                u32 old_cell_idx  = cur_idx * AMRBlock::block_size + loc_id;
 
                // // // cell size in the refined block
                // // Tscal delta_cell = cell_sizes[cur_idx];
                // // Tscal c_offset   = delta_cell * 0.25;
                // // std::array<f64_3, AMRBlock::block_size> child_center_offsets;
                // //  child_center_offsets[0] = {-c_offset, -c_offset, -c_offset}; /*(0,0,0) */
                // // child_center_offsets[1] = {c_offset, -c_offset, -c_offset};   /*(1,0,0)*/
                // // child_center_offsets[2] = {-c_offset, c_offset, -c_offset};  /* (0,1,0)*/
                // // child_center_offsets[3] = {c_offset, c_offset, -c_offset};   /*(1,1,0)*/
                // // child_center_offsets[4] = {-c_offset, -c_offset, c_offset}; /*(0,0,1)*/
                // // child_center_offsets[5] = {c_offset, -c_offset, c_offset}; /*(1,0,1)*/
                // // child_center_offsets[6] = {-c_offset, c_offset, c_offset};  /*(0,1,1)*/
                // // child_center_offsets[7] = {c_offset, c_offset, c_offset};   /*(1,1,1)*/
 
                // auto cons_var_slopes = get_3d_grad_cons<Tvec, Minmod>(
                //     old_cell_idx,
                //     delta_cell,
                //     cell_graph_xp,
                //     cell_graph_xm,
                //     cell_graph_yp,
                //     cell_graph_ym,
                //     cell_graph_zp,
                //     cell_graph_zm,
                //     [=](u32 id){
                //         return acc.rho[id];
                //     },
                //     [=](u32 id){
                //         return acc.rho_vel[id];
                //     },
                //     [=](u32 id){
                //         return acc.rhoE[id];
                //     });
 
                Tscal rho_block    = old_rho_block[loc_id];
                Tvec rho_vel_block = old_rho_vel_block[loc_id];
                Tscal rhoE_block   = old_rhoE_block[loc_id];
                for (u32 subdiv_lid = 0; subdiv_lid < 8; subdiv_lid++) {
 
                    auto [sx, sy, sz] = get_coord_ref(subdiv_lid);
 
                    std::array<u32, 3> glid = {lx * 2 + sx, ly * 2 + sy, lz * 2 + sz};
 
                    u32 new_cell_idx = get_gid_write(glid);
 
                    // shammath::ConsState<Tvec> cons_var_interp =
                    // child_center_offsets[subdiv_lid][0] * cons_var_slopes[0] +
                    // child_center_offsets[subdiv_lid][1] * cons_var_slopes[1] +
                    // child_center_offsets[subdiv_lid][2] * cons_var_slopes[2];
 
                    // // acc.rho[new_cell_idx]     = rho_block + cons_var_interp.rho ;
                    // // acc.rho_vel[new_cell_idx] = rho_vel_block + cons_var_interp.rhovel;
                    // // acc.rhoE[new_cell_idx]    = rhoE_block + cons_var_interp.rhoe;
 
                    acc.rho[new_cell_idx]     = rho_block;
                    acc.rho_vel[new_cell_idx] = rho_vel_block;
                    acc.rhoE[new_cell_idx]    = rhoE_block;
                }
            }
        }
 
        void apply_derefine_new(
            std::array<u32, 8> old_blocks,
            std::array<BlockCoord, 8> old_coords,
            u32 new_cell,
            BlockCoord new_coord,
 
            RefineCellAccessor acc) const {
 
            std::array<f64, AMRBlock::block_size> rho_block;
            std::array<f64_3, AMRBlock::block_size> rho_vel_block;
            std::array<f64, AMRBlock::block_size> rhoE_block;
 
            for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) {
                rho_block[cell_id]     = {};
                rho_vel_block[cell_id] = {};
                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_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];
                }
                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++) {
                u32 newcell_idx          = new_cell * AMRBlock::block_size + cell_id;
                acc.rho[newcell_idx]     = rho_block[cell_id];
                acc.rho_vel[newcell_idx] = rho_vel_block[cell_id];
                acc.rhoE[newcell_idx]    = rhoE_block[cell_id];
            }
        }
    };

    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;
 
        const Tscal *rho_cons;
        Tscal error_min;
        Tscal error_max;
        u32 nblock_per_patch;
 
        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<TgridVec>(0).get_buf().get_read_access(depends_list);
            block_high_bound = pdat.get_field<TgridVec>(1).get_buf().get_read_access(depends_list);
            rho_cons         = pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                                   .get_buf()
                                   .get_read_access(depends_list);
 
            nblock_per_patch = pdat.get_obj_cnt();
 
            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_new(
            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<i64_3>(0).get_buf().complete_event_state(resulting_events);
            pdat.get_field<i64_3>(1).get_buf().complete_event_state(resulting_events);
            pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("rho"))
                .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);
        }
 
        void refine_criterion_new(
            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<Tscal>::get_zero();
            // for (u32 i = 0; i < AMRBlock::block_size; i++) {
            //     block_rho_slope = sham::details::g_sycl_max(
            //         block_rho_slope,
            //         modif_second_derivative<Tscal, Tvec>(
            //             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 acc.block_rho[id];
            //             }));
            // }
 
            Tscal block_rho_slope = shambase::VectorProperties<Tscal>::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<Tscal>
            //         // get_pseudo_grad<Tscal, Tvec>
            //         (
            //             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 rho_cons[id];
            //                 return  acc.block_rho[id];
            //             }));
            // }
 
            Tscal block_press_grad = shambase::VectorProperties<Tscal>::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<Tscal, Tvec>(
                        i + block_id * AMRBlock::block_size,
                        cell_graph_xp,
                        cell_graph_xm,
                        cell_graph_yp,
                        cell_graph_ym,
                        cell_graph_zp,
                        cell_graph_zm,
                        [=, this](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);
        }
    };

    class RefineCritShearAccessor {
        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 threshold;
        Tscal gamma;
        Tscal dxfact;
 
        AMRGraphLinkiterator cell_graph_xp;
        AMRGraphLinkiterator cell_graph_xm;
        AMRGraphLinkiterator cell_graph_yp;
        AMRGraphLinkiterator cell_graph_ym;
        AMRGraphLinkiterator cell_graph_zp;
        AMRGraphLinkiterator cell_graph_zm;
 
        RefineCritShearAccessor(
            sham::EventList &depends_list,
            Storage &storage,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal threshold,
            Tscal gamma,
            Tscal dxfact)
            : threshold(threshold), gamma(gamma), dxfact(dxfact),
              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<TgridVec>(0).get_buf().get_read_access(depends_list);
            block_high_bound = pdat.get_field<TgridVec>(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_new(
            sham::EventList &resulting_events,
            Storage &storage,
            u64 id_patch,
            shamrock::patch::Patch p,
            shamrock::patch::PatchDataLayer &pdat,
            Tscal threshold,
            Tscal gamma,
            Tscal dxfact) {
 
            pdat.get_field<i64_3>(0).get_buf().complete_event_state(resulting_events);
            pdat.get_field<i64_3>(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);
        }
 
        void refine_criterion_new(
            u32 block_id,
            RefineCritShearAccessor 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];
 
            Tvec lower_flt = low_bound.template convert<Tscal>() * dxfact;
            Tvec upper_flt = high_bound.template convert<Tscal>() * dxfact;
 
            Tvec block_cell_size = (upper_flt - lower_flt) * one_over_Nside;
 
            Tscal block_normalized_shear = shambase::VectorProperties<Tscal>::get_zero();
            for (u32 i = 0; i < AMRBlock::block_size; i++) {
                auto cell_id = i + block_id * AMRBlock::block_size;
                auto cs = sycl::sqrt(gamma * acc.block_pressure[cell_id] / acc.block_rho[cell_id]);
                block_normalized_shear = sham::details::g_sycl_max(
                    block_normalized_shear,
                    normalized_shear<Tvec>(
                        cell_id,
                        cs,
                        block_cell_size,
 
                        cell_graph_xp,
                        cell_graph_xm,
                        cell_graph_yp,
                        cell_graph_ym,
                        cell_graph_zp,
                        cell_graph_zm,
                        [=](u32 id) {
                            return acc.block_velocity[id];
                        }));
            }
            should_refine   = false;
            should_derefine = false;
            if (block_normalized_shear > threshold * threshold) {
                should_refine = true;
            } else if (block_normalized_shear < 0.25 * threshold * threshold) {
                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);
        }
    };

    class RefineCellAccessorAutogravity {
        public:
        f64 *rho;
        f64_3 *rho_vel;
        f64 *rhoE;
        f64 *phi_old;
        f64 *phi_new;
 
        u64 p_id;
        // f64* cell_sizes;
 
        // 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;
 
        RefineCellAccessorAutogravity(
            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<f64>(2).get_buf().get_write_access(depends_list);
            rho_vel = pdat.get_field<f64_3>(3).get_buf().get_write_access(depends_list);
            rhoE    = pdat.get_field<f64>(4).get_buf().get_write_access(depends_list);
            phi_old = pdat.get_field<f64>(pdat.pdl().get_field_idx<Tscal>("phi_old"))
                          .get_buf()
                          .get_write_access(depends_list);
            phi_new = pdat.get_field<f64>(pdat.pdl().get_field_idx<Tscal>("phi"))
                          .get_buf()
                          .get_write_access(depends_list);
        }
 
        void finalize_new(
            sham::EventList &resulting_events,
            Storage &storage,
            u64 &id_patch,
            shamrock::patch::PatchDataLayer &pdat) {
            pdat.get_field<f64>(2).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64_3>(3).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64>(4).get_buf().complete_event_state(resulting_events);
            pdat.get_field<f64>(pdat.pdl().get_field_idx<Tscal>("phi_old"))
                .get_buf()
                .complete_event_state(resulting_events);
            pdat.get_field<f64>(pdat.pdl().get_field_idx<Tscal>("phi"))
                .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_new(
            u32 cur_idx,
            BlockCoord cur_coords,
            std::array<u32, 8> new_blocks,
            std::array<BlockCoord, 8> new_block_coords,
            RefineCellAccessorAutogravity acc) const {
 
            auto get_coord_ref = [](u32 i) -> std::array<u32, dim> {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    const u32 tmp = i >> NsideBlockPow;
                    return {i % Nside, (tmp) % Nside, (tmp) >> NsideBlockPow};
                }
            };
 
            auto get_index_block = [](std::array<u32, dim> coord) -> u32 {
                constexpr u32 NsideBlockPow = 1;
                constexpr u32 Nside         = 1U << NsideBlockPow;
 
                if constexpr (dim == 3) {
                    return coord[0] + Nside * coord[1] + Nside * Nside * coord[2];
                }
            };
 
            auto get_gid_write = [&](std::array<u32, dim> &glid) -> u32 {
                // First, get the block id (it's the block to be refine) in wich the new cell glid
                // is located.
                std::array<u32, dim> bid
                    = {glid[0] >> AMRBlock::NsideBlockPow,
                       glid[1] >> AMRBlock::NsideBlockPow,
                       glid[2] >> AMRBlock::NsideBlockPow};
 
                // 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,
                            glid[2] % AMRBlock::Nside});
            };
 
            std::array<f64, AMRBlock::block_size> old_rho_block;
            std::array<f64_3, AMRBlock::block_size> old_rho_vel_block;
            std::array<f64, AMRBlock::block_size> old_rhoE_block;
            std::array<f64, AMRBlock::block_size> old_phi_old_block;
            std::array<f64, AMRBlock::block_size> old_phi_new_block;
 
            // save old block
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz]         = get_coord_ref(loc_id);
                u32 old_cell_idx          = cur_idx * AMRBlock::block_size + loc_id;
                old_rho_block[loc_id]     = acc.rho[old_cell_idx];
                old_rho_vel_block[loc_id] = acc.rho_vel[old_cell_idx];
                old_rhoE_block[loc_id]    = acc.rhoE[old_cell_idx];
                old_phi_old_block[loc_id] = acc.phi_old[old_cell_idx];
                old_phi_new_block[loc_id] = acc.phi_new[old_cell_idx];
            }
 
            for (u32 loc_id = 0; loc_id < AMRBlock::block_size; loc_id++) {
 
                auto [lx, ly, lz] = get_coord_ref(loc_id);
                u32 old_cell_idx  = cur_idx * AMRBlock::block_size + loc_id;
 
                // // // cell size in the refined block
                // // Tscal delta_cell = cell_sizes[cur_idx];
                // // Tscal c_offset   = delta_cell * 0.25;
                // // std::array<f64_3, AMRBlock::block_size> child_center_offsets;
                // //  child_center_offsets[0] = {-c_offset, -c_offset, -c_offset}; /*(0,0,0) */
                // // child_center_offsets[1] = {c_offset, -c_offset, -c_offset};   /*(1,0,0)*/
                // // child_center_offsets[2] = {-c_offset, c_offset, -c_offset};  /* (0,1,0)*/
                // // child_center_offsets[3] = {c_offset, c_offset, -c_offset};   /*(1,1,0)*/
                // // child_center_offsets[4] = {-c_offset, -c_offset, c_offset}; /*(0,0,1)*/
                // // child_center_offsets[5] = {c_offset, -c_offset, c_offset}; /*(1,0,1)*/
                // // child_center_offsets[6] = {-c_offset, c_offset, c_offset};  /*(0,1,1)*/
                // // child_center_offsets[7] = {c_offset, c_offset, c_offset};   /*(1,1,1)*/
 
                // auto cons_var_slopes = get_3d_grad_cons<Tvec, Minmod>(
                //     old_cell_idx,
                //     delta_cell,
                //     cell_graph_xp,
                //     cell_graph_xm,
                //     cell_graph_yp,
                //     cell_graph_ym,
                //     cell_graph_zp,
                //     cell_graph_zm,
                //     [=](u32 id){
                //         return acc.rho[id];
                //     },
                //     [=](u32 id){
                //         return acc.rho_vel[id];
                //     },
                //     [=](u32 id){
                //         return acc.rhoE[id];
                //     });
 
                Tscal rho_block     = old_rho_block[loc_id];
                Tvec rho_vel_block  = old_rho_vel_block[loc_id];
                Tscal rhoE_block    = old_rhoE_block[loc_id];
                Tscal phi_old_block = old_phi_old_block[loc_id];
                Tscal phi_new_block = old_phi_new_block[loc_id];
 
                for (u32 subdiv_lid = 0; subdiv_lid < 8; subdiv_lid++) {
 
                    auto [sx, sy, sz] = get_coord_ref(subdiv_lid);
 
                    std::array<u32, 3> glid = {lx * 2 + sx, ly * 2 + sy, lz * 2 + sz};
 
                    u32 new_cell_idx = get_gid_write(glid);
 
                    // shammath::ConsState<Tvec> cons_var_interp =
                    // child_center_offsets[subdiv_lid][0] * cons_var_slopes[0] +
                    // child_center_offsets[subdiv_lid][1] * cons_var_slopes[1] +
                    // child_center_offsets[subdiv_lid][2] * cons_var_slopes[2];
 
                    // // acc.rho[new_cell_idx]     = rho_block + cons_var_interp.rho ;
                    // // acc.rho_vel[new_cell_idx] = rho_vel_block + cons_var_interp.rhovel;
                    // // acc.rhoE[new_cell_idx]    = rhoE_block + cons_var_interp.rhoe;
 
                    acc.rho[new_cell_idx]     = rho_block;
                    acc.rho_vel[new_cell_idx] = rho_vel_block;
                    acc.rhoE[new_cell_idx]    = rhoE_block;
                    acc.phi_old[new_cell_idx] = phi_old_block;
                    acc.phi_new[new_cell_idx] = phi_new_block;
                }
            }
        }
 
        void apply_derefine_new(
            std::array<u32, 8> old_blocks,
            std::array<BlockCoord, 8> old_coords,
            u32 new_cell,
            BlockCoord new_coord,
 
            RefineCellAccessorAutogravity acc) const {
 
            std::array<f64, AMRBlock::block_size> rho_block;
            std::array<f64_3, AMRBlock::block_size> rho_vel_block;
            std::array<f64, AMRBlock::block_size> rhoE_block;
            std::array<f64, AMRBlock::block_size> phi_old_block;
            std::array<f64, AMRBlock::block_size> phi_new_block;
 
            for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) {
                rho_block[cell_id]     = {};
                rho_vel_block[cell_id] = {};
                rhoE_block[cell_id]    = {};
                phi_old_block[cell_id] = {};
                phi_new_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];
                auto phi_old_pid = phi_old_block[pid];
                auto phi_new_pid = phi_new_block[pid];
 
                for (u32 cell_id = 0; cell_id < 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];
                    phi_old_pid += acc.phi_old[old_blocks[pid] * AMRBlock::block_size + cell_id];
                    phi_new_pid += acc.phi_new[old_blocks[pid] * AMRBlock::block_size + cell_id];
                }
                rho_block[pid]     = rho_pid * (1. / 8.);
                rho_vel_block[pid] = rho_vel_pid * (1. / 8.);
                rhoE_block[pid]    = rhoe_pid * (1. / 8.);
                // phi_old_block[pid]    = phi_old_pid * (1. / 8.);
                // phi_new_block[pid]    = phi_new_pid * (1. / 8.);
            }
 
            for (u32 cell_id = 0; cell_id < AMRBlock::block_size; cell_id++) {
                u32 newcell_idx          = new_cell * AMRBlock::block_size + cell_id;
                acc.rho[newcell_idx]     = rho_block[cell_id];
                acc.rho_vel[newcell_idx] = rho_vel_block[cell_id];
                acc.rhoE[newcell_idx]    = rhoE_block[cell_id];
 
                // acc.phi_old[newcell_idx] = phi_old_block[cell_id];
                // acc.phi_new[newcell_idx] = phi_new_block[cell_id];
            }
        }
    };
 
    using AMRmode_None                = typename AMRMode<Tvec, TgridVec>::None;
    using AMRmode_DensityBased        = typename AMRMode<Tvec, TgridVec>::DensityBased;
    using AMRmode_PseudoGradientBased = typename AMRMode<Tvec, TgridVec>::PseudoGradientBased;
    using AMRmode_JeansLengthBased    = typename AMRMode<Tvec, TgridVec>::JeansLengthBased;
    using AMRmode_ShearBased          = typename AMRMode<Tvec, TgridVec>::ShearBased;
 
    bool has_cell_order_changed = false;
 
    // get refine and derefine list
    shambase::DistributedData<sham::DeviceBuffer<u32>> refine_list;
    shambase::DistributedData<sham::DeviceBuffer<u32>> derefine_list;
 
    if (AMRmode_None *cfg = std::get_if<AMRmode_None>(&solver_config.amr_mode.config)) {
        // no refinment here turn around there is nothing to see
    } else {
        if (AMRmode_DensityBased *cfg
            = std::get_if<AMRmode_DensityBased>(&solver_config.amr_mode.config)) {
 
            Tscal dxfact(solver_config.grid_coord_to_pos_fact);
            gen_refine_block_changes_new<RefineCritBlock>(
                refine_list, derefine_list, dxfact, cfg->crit_mass);
        }
 
        else if (
            AMRmode_PseudoGradientBased *cfg
            = std::get_if<AMRmode_PseudoGradientBased>(&solver_config.amr_mode.config)) {
 
            gen_refine_block_changes_new<RefineCritPseudoGradientAccessor>(
                refine_list, derefine_list, cfg->error_min, cfg->error_max);
        }
 
        else if (
            AMRmode_ShearBased *cfg
            = std::get_if<AMRmode_ShearBased>(&solver_config.amr_mode.config)) {
            Tscal dxfact(solver_config.grid_coord_to_pos_fact);
            Tscal gamma(solver_config.eos_gamma);
 
            gen_refine_block_changes_new<RefineCritShearAccessor>(
                refine_list, derefine_list, cfg->threshold, gamma, dxfact);
        }

        enforce_two_to_one_refinement_new(std::move(refine_list));
        enforce_two_to_one_derefinement_new(std::move(derefine_list), std::move(refine_list));
        // Note that this only add new blocks at the end of the patchdata
        bool change_refine = internal_refine_grid_new<RefineCellAccessor>(std::move(refine_list));

        // Note that this will perform the merge then remove the old blocks
        // This is ok to call straight after the refine without edditing the index list in
        // derefine_list since no permutations were applied in internal_refine_grid_new and no cells
        // can be both refined and derefined in the same pass
        bool change_derefine
            = internal_derefine_grid_new<RefineCellAccessor>(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);
            block_sorter.reorder_amr_blocks();
        }
    }
}
 
template class shammodels::basegodunov::modules::AMRGridRefinementHandler<f64_3, i64_3>;