SimBox.hpp

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
//
// -------------------------------------------------------//
 
#pragma once

 
#include "shambase/exception.hpp"
#include "shamalgs/collective/reduction.hpp"
#include "shambackends/sycl_utils.hpp"
#include "shammath/CoordRange.hpp"
#include "shamrock/patch/Patch.hpp"
#include "shamrock/patch/PatchCoordTransform.hpp"
#include "shamrock/patch/PatchDataLayerLayout.hpp"
#include <type_traits>
#include <stdexcept>
#include <tuple>
 
namespace shamrock::patch {

    class SimulationBoxInfo {
 
        static constexpr u32 dim = 3;
 
        using var_t = FieldVariant<shammath::CoordRange>;
 
        PatchDataLayerLayout &pdl;
 
        var_t bounding_box;
        PatchCoord<> patch_coord_bounding_box;
 
        public:
        inline SimulationBoxInfo(
            PatchDataLayerLayout &pdl, PatchCoord<dim> patch_coord_bounding_box)
            : pdl(pdl), patch_coord_bounding_box(std::move(patch_coord_bounding_box)),
              bounding_box(shammath::CoordRange<f32>{}) {
 
            reset_box_size();
        }
 
        void set_patch_coord_bounding_box(PatchCoord<dim> new_patch_coord_box) {
            patch_coord_bounding_box = new_patch_coord_box;
            shamlog_debug_ln(
                "SimBox",
                "changed patch coord bounds :",
                std::pair{
                    u64_3{
                        new_patch_coord_box.coord_min[0],
                        new_patch_coord_box.coord_min[1],
                        new_patch_coord_box.coord_min[2]},
                    u64_3{
                        new_patch_coord_box.coord_max[0],
                        new_patch_coord_box.coord_max[1],
                        new_patch_coord_box.coord_max[2]}});
        }

        template<class T>
        [[nodiscard]] std::tuple<T, T> get_bounding_box() const;

        template<class T>
        inline T get_bounding_box_size() const {
            auto [bmin, bmax] = get_bounding_box<T>();
            return bmax - bmin;
        }

        template<class T>
        void set_bounding_box(shammath::CoordRange<T> new_box);

        template<class T>
        inline void allreduce_set_bounding_box(shammath::CoordRange<T> new_box) {
            std::pair<T, T> reduced = shamalgs::collective::allreduce_bounds(
                std::pair<T, T>{new_box.lower, new_box.upper});
            set_bounding_box<>(shammath::CoordRange<T>(reduced));
        }

        template<class T>
        PatchCoordTransform<T> get_patch_transform() const;

        template<class T>
        std::tuple<T, T> patch_coord_to_domain(const Patch &p) const;
 
        // TODO implement box size reduction here

        void reset_box_size();

        template<class primtype>
        void clean_box(primtype tol);
 
        template<>
        inline void clean_box<f32>(f32 tol) {
 
            auto [bmin, bmax] = get_bounding_box<f32_3>();
 
            f32_3 center   = (bmin + bmax) / 2;
            f32_3 cur_delt = bmax - bmin;
            cur_delt /= 2;
 
            cur_delt *= tol;
 
            bmin = center - cur_delt;
            bmax = center + cur_delt;
 
            set_bounding_box<f32_3>({bmin, bmax});
        }
 
        template<>
        inline void clean_box<f64>(f64 tol) {
            auto [bmin, bmax] = get_bounding_box<f64_3>();
 
            f64_3 center   = (bmin + bmax) / 2;
            f64_3 cur_delt = bmax - bmin;
            cur_delt /= 2;
 
            cur_delt *= tol;
 
            bmin = center - cur_delt;
            bmax = center + cur_delt;
 
            set_bounding_box<f64_3>({bmin, bmax});
        }
 
        template<class primtype>
        inline std::tuple<sycl::vec<primtype, 3>, sycl::vec<primtype, 3>> get_box(Patch &p) {
            return patch_coord_to_domain<sycl::vec<primtype, 3>>(p);
        }
 
        template<class T>
        inline PatchCoordTransform<T> get_transform() {
            auto [bmin, bmax] = get_bounding_box<T>();
            return PatchCoordTransform<T>{
                patch_coord_bounding_box, shammath::CoordRange<T>{bmin, bmax}};
        }

        void to_json(nlohmann::json &j);

        void from_json(const nlohmann::json &j);
    };

    // out of line implementation of the simbox
 
    template<class T>
    [[nodiscard]] inline std::tuple<T, T> SimulationBoxInfo::get_bounding_box() const {
 
        if (!pdl.check_main_field_type<T>()) {
 
            throw shambase::make_except_with_loc<std::invalid_argument>(
                "the chosen type for the main field does not match the required template type\n"
                "call : "
                + std::string(__PRETTY_FUNCTION__));
        }
 
        const shammath::CoordRange<T> *pval
            = std::get_if<shammath::CoordRange<T>>(&bounding_box.value);
 
        if (!pval) {
 
            throw shambase::make_except_with_loc<std::invalid_argument>(
                "the type in SimulationBoxInfo does not match the one in the layout\n"
                "call : "
                + std::string(__PRETTY_FUNCTION__));
        }
 
        return {pval->lower, pval->upper};
    }
 
    template<class T>
    inline void SimulationBoxInfo::set_bounding_box(shammath::CoordRange<T> new_box) {
        new_box.check_throw_ranges();
        if (pdl.check_main_field_type<T>()) {
            bounding_box.value = new_box;
        } else {
            throw shambase::make_except_with_loc<std::invalid_argument>(
                "The main field is not of the required type\n"
                "call : "
                + std::string(__PRETTY_FUNCTION__));
        }
    }
 
    template<class T>
    PatchCoordTransform<T> SimulationBoxInfo::get_patch_transform() const {
 
        auto [bmin, bmax] = get_bounding_box<T>();
 
        shammath::CoordRange<T> tmp{bmin, bmax};
 
        if (tmp.is_err_mode()) {
            throw shambase::make_except_with_loc<std::runtime_error>(
                "the box size is not set, please resize the box to the domain size");
        }
 
        return PatchCoordTransform<T>{patch_coord_bounding_box.get_patch_range(), tmp};
    }
 
    template<class T>
    inline std::tuple<T, T> SimulationBoxInfo::patch_coord_to_domain(const Patch &p) const {
 
        PatchCoordTransform<T> transform = get_patch_transform<T>();
 
        auto [obj_min, obj_max] = transform.to_obj_coord(p);
 
        return {obj_min, obj_max};
    }
 
    inline void SimulationBoxInfo::reset_box_size() {
 
        if (pdl.check_main_field_type<f32_3>()) {
            bounding_box.value = shammath::CoordRange<f32_3>::max_range();
        } else if (pdl.check_main_field_type<f64_3>()) {
            bounding_box.value = shammath::CoordRange<f64_3>::max_range();
        } else if (pdl.check_main_field_type<u32_3>()) {
            bounding_box.value = shammath::CoordRange<u32_3>::max_range();
        } else if (pdl.check_main_field_type<u64_3>()) {
            bounding_box.value = shammath::CoordRange<u64_3>::max_range();
        } else if (pdl.check_main_field_type<i64_3>()) {
            bounding_box.value = shammath::CoordRange<i64_3>::max_range();
        } else {
            throw shambase::make_except_with_loc<std::runtime_error>(
                "the chosen type for the main field is not handled");
        }
    }
 
} // namespace shamrock::patch