diff --git a/.gitignore b/.gitignore index 085d61c97..2dc1404d2 100644 --- a/.gitignore +++ b/.gitignore @@ -8,7 +8,8 @@ Testing .idea *.png *.jpg -data* +data_out* .sublime* compile_commands.json .clangd +phare_outputs diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 18b4db1aa..588a1a9a6 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -317,9 +317,11 @@ def get_max_ghosts(): return max(grid.nbrGhosts(kwargs["interp_order"], x) for x in ['primal','dual']) max_ghosts = get_max_ghosts() - small_invalid_patch_size = phare_utilities.np_array_ify(max_ghosts - 1, ndim) + small_invalid_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim) largest_patch_size = kwargs.get("largest_patch_size", None) - smallest_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim) + + # to prevent primal ghost overlaps of non adjacent patches, we need smallest_patch_size+=1 + smallest_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim) + 1 if "smallest_patch_size" in kwargs and kwargs["smallest_patch_size"] is not None: smallest_patch_size = phare_utilities.np_array_ify(kwargs["smallest_patch_size"], ndim) diff --git a/src/amr/CMakeLists.txt b/src/amr/CMakeLists.txt index d9c6782b8..62f1bcf23 100644 --- a/src/amr/CMakeLists.txt +++ b/src/amr/CMakeLists.txt @@ -64,6 +64,7 @@ set( SOURCES_CPP resources_manager/amr_utils.cpp data/field/coarsening/field_coarsen.cpp messengers/messenger_factory.cpp + data/field/field_variable_fill_pattern.cpp ) add_library(${PROJECT_NAME} ${SOURCES_INC} ${SOURCES_CPP}) diff --git a/src/amr/data/field/field_geometry.h b/src/amr/data/field/field_geometry.h index a1f5747f8..1543d25da 100644 --- a/src/amr/data/field/field_geometry.h +++ b/src/amr/data/field/field_geometry.h @@ -9,6 +9,24 @@ #include "core/data/grid/gridlayout_impl.h" #include "field_overlap.h" + + +namespace PHARE::amr +{ +class AFieldGeometry : public SAMRAI::hier::BoxGeometry +{ +public: + virtual ~AFieldGeometry() {} + AFieldGeometry(SAMRAI::hier::Box const& box) + : patchBox{box} + { + } + + SAMRAI::hier::Box const patchBox; +}; + +} // namespace PHARE::amr + namespace PHARE { namespace amr @@ -17,7 +35,7 @@ namespace amr /** * @brief The FieldGeometry class */ - class FieldGeometry : public SAMRAI::hier::BoxGeometry + class FieldGeometry : public AFieldGeometry { public: static constexpr std::size_t dimension = GridLayoutT::dimension; @@ -27,7 +45,8 @@ namespace amr * with a temporary gridlayout */ FieldGeometry(SAMRAI::hier::Box const& box, GridLayoutT layout, PhysicalQuantity qty) - : ghostBox_{toFieldBox(box, qty, layout)} + : AFieldGeometry{box} + , ghostBox_{toFieldBox(box, qty, layout)} , interiorBox_{toFieldBox(box, qty, layout, false)} , layout_{std::move(layout)} , quantity_{qty} @@ -222,7 +241,6 @@ namespace amr PhysicalQuantity quantity_; - /*** \brief Compute destination box representing the intersection of two geometry * * \param destinationBoxes BoxContainer that will be filled of box @@ -286,10 +304,10 @@ namespace amr SAMRAI::hier::Box const together(destinationBox * sourceBox * fillField); - // if the interesection is not empty we either push it into the container // if we don't want to fill the interior we remove it from the intersection // which may add multiple boxes to the container. + if (!together.empty()) { if (overwriteInterior) diff --git a/src/amr/data/field/field_variable_fill_pattern.cpp b/src/amr/data/field/field_variable_fill_pattern.cpp new file mode 100644 index 000000000..6aa3d800d --- /dev/null +++ b/src/amr/data/field/field_variable_fill_pattern.cpp @@ -0,0 +1,10 @@ + + +#include "field_variable_fill_pattern.h" + +namespace PHARE::amr +{ +// the value of this string cannot change, SAMRAI breaks :| +const std::string FieldFillPattern::s_name_id = "BOX_GEOMETRY_FILL_PATTERN"; + +} // namespace PHARE::amr diff --git a/src/amr/data/field/field_variable_fill_pattern.h b/src/amr/data/field/field_variable_fill_pattern.h new file mode 100644 index 000000000..732bc9a0b --- /dev/null +++ b/src/amr/data/field/field_variable_fill_pattern.h @@ -0,0 +1,97 @@ +#ifndef PHARE_SRC_AMR_FIELD_FIELD_VARIABLE_FILL_PATTERN_H +#define PHARE_SRC_AMR_FIELD_FIELD_VARIABLE_FILL_PATTERN_H + +#include +#include "SAMRAI/xfer/VariableFillPattern.h" +#include "amr/data/field/field_geometry.h" + +namespace PHARE::amr +{ +class FieldFillPattern : public SAMRAI::xfer::VariableFillPattern +{ +public: + FieldFillPattern() {} + virtual ~FieldFillPattern() {} + + std::shared_ptr + calculateOverlap(const SAMRAI::hier::BoxGeometry& dst_geometry, + const SAMRAI::hier::BoxGeometry& src_geometry, + const SAMRAI::hier::Box& dst_patch_box, const SAMRAI::hier::Box& src_mask, + const SAMRAI::hier::Box& fill_box, const bool overwrite_interior_, + const SAMRAI::hier::Transformation& transformation) const + { +#ifndef DEBUG_CHECK_DIM_ASSERTIONS + NULL_USE(dst_patch_box); +#endif + TBOX_ASSERT_OBJDIM_EQUALITY2(dst_patch_box, src_mask); + + bool overwrite_interior = true; // replace func param + assert(overwrite_interior_ == overwrite_interior); + + auto& dst_cast = dynamic_cast(dst_geometry); + auto& src_cast = dynamic_cast(src_geometry); + + // for shared border node value sync + if (src_cast.patchBox.getGlobalId().getOwnerRank() + != dst_cast.patchBox.getGlobalId().getOwnerRank()) + overwrite_interior = src_cast.patchBox.getGlobalId() > dst_cast.patchBox.getGlobalId(); + // overwrite_interior = src_cast.patchBox.getLocalId() > dst_cast.patchBox.getLocalId(); + + return dst_geometry.calculateOverlap(src_geometry, src_mask, fill_box, overwrite_interior, + transformation); + } + + const std::string& getPatternName() const { return s_name_id; } + +private: + FieldFillPattern(const FieldFillPattern&); // not implemented + FieldFillPattern& operator=(const FieldFillPattern&); // not implemented + + static const std::string s_name_id; // = "GHOST_ONLY_FILL_PATTERN"; + + const SAMRAI::hier::IntVector& getStencilWidth() + { + TBOX_ERROR("getStencilWidth() should not be\n" + << "called. This pattern creates overlaps based on\n" + << "the BoxGeometry objects and is not restricted to a\n" + << "specific stencil.\n"); + + /* + * Dummy return value that will never get reached. + */ + return SAMRAI::hier::IntVector::getZero(SAMRAI::tbox::Dimension(1)); + } + + /* + ************************************************************************* + * + * Compute BoxOverlap that specifies data to be filled by refinement + * operator. + * + ************************************************************************* + */ + std::shared_ptr + computeFillBoxesOverlap(const SAMRAI::hier::BoxContainer& fill_boxes, + const SAMRAI::hier::BoxContainer& node_fill_boxes, + const SAMRAI::hier::Box& patch_box, const SAMRAI::hier::Box& data_box, + const SAMRAI::hier::PatchDataFactory& pdf) const + { + NULL_USE(node_fill_boxes); + + /* + * For this (default) case, the overlap is simply the intersection of + * fill_boxes and data_box. + */ + SAMRAI::hier::Transformation transformation( + SAMRAI::hier::IntVector::getZero(patch_box.getDim())); + + SAMRAI::hier::BoxContainer overlap_boxes(fill_boxes); + overlap_boxes.intersectBoxes(data_box); + + return pdf.getBoxGeometry(patch_box)->setUpOverlap(overlap_boxes, transformation); + } +}; + +} // namespace PHARE::amr + +#endif /* PHARE_SRC_AMR_FIELD_FIELD_VARIABLE_FILL_PATTERN_H */ diff --git a/src/amr/messengers/quantity_communicator.h b/src/amr/messengers/quantity_communicator.h index 71b8b5fea..e7852edf1 100644 --- a/src/amr/messengers/quantity_communicator.h +++ b/src/amr/messengers/quantity_communicator.h @@ -4,6 +4,7 @@ #include "amr/messengers/hybrid_messenger_info.h" +#include "amr/data/field/field_variable_fill_pattern.h" #include "amr/data/field/coarsening/field_coarsen_operator.h" #include "amr/data/field/refine/field_refine_operator.h" #include "amr/data/field/time_interpolate/field_linear_time_interpolate.h" @@ -24,6 +25,7 @@ + namespace PHARE { namespace amr @@ -48,6 +50,18 @@ namespace amr { }; + class XFieldFillPattern : public FieldFillPattern + { + }; + + class YFieldFillPattern : public FieldFillPattern + { + }; + + class ZFieldFillPattern : public FieldFillPattern + { + }; + struct Refiner { using schedule_type = SAMRAI::xfer::RefineSchedule; @@ -125,6 +139,9 @@ namespace amr */ void add(std::shared_ptr schedule, int levelNumber) { + // for shared border node value sync + schedule->setDeterministicUnpackOrderingFlag(true); + schedules_[levelNumber] = std::move(schedule); } }; @@ -160,11 +177,12 @@ namespace amr std::shared_ptr timeOp) { std::shared_ptr xVariableFillPattern - = std::make_shared(); + = std::make_shared(); std::shared_ptr yVariableFillPattern - = std::make_shared(); + = std::make_shared(); std::shared_ptr zVariableFillPattern - = std::make_shared(); + = std::make_shared(); + Communicator com; auto registerRefine diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index 9672698b6..a024c4e24 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -18,15 +18,18 @@ def __setitem__(self, k, v): def basicSimulatorArgs(dim: int, interp: int, **kwargs): from pyphare.pharein.simulation import valid_refined_particle_nbr + from pyphare.pharein.simulation import check_patch_size cells = kwargs.get("cells", [20 for i in range(dim)]) if not isinstance(cells, (list, tuple)): cells = [cells] * dim + + _, smallest_patch_size = check_patch_size(dim, interp_order=interp, cells=cells) dl = [1.0 / v for v in cells] b0 = [[3] * dim, [8] * dim] args = { "interp_order": interp, - "smallest_patch_size": [5] * dim, + "smallest_patch_size": smallest_patch_size, "largest_patch_size": [20] * dim, "time_step_nbr": 1000, "final_time": 1.0, diff --git a/tests/simulator/refined_particle_nbr.py b/tests/simulator/refined_particle_nbr.py index c41b60aa2..ff1078a09 100644 --- a/tests/simulator/refined_particle_nbr.py +++ b/tests/simulator/refined_particle_nbr.py @@ -20,6 +20,7 @@ class SimulatorRefinedParticleNbr(unittest.TestCase): def __init__(self, *args, **kwargs): super(SimulatorRefinedParticleNbr, self).__init__(*args, **kwargs) + self.simulator = None with open(os.path.join(project_root, "res/amr/splitting.yml"), 'r') as stream: try: self.yaml_root = yaml.safe_load(stream) diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 31b78cdcc..a3c0aa39c 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -2,6 +2,7 @@ cpp = cpp_lib() from pyphare.simulator.simulator import Simulator, startMPI +from pyphare.core.phare_utilities import np_array_ify from pyphare.pharesee.hierarchy import hierarchy_from, merge_particles from pyphare.pharein import MaxwellianFluidModel from pyphare.pharein.diagnostics import ParticleDiagnostics, FluidDiagnostics, ElectromagDiagnostics @@ -26,12 +27,18 @@ def ddt_test_id(self): def getHierarchy(self, interp_order, refinement_boxes, qty, nbr_part_per_cell=100, diag_outputs="phare_outputs", - smallest_patch_size=5, largest_patch_size=20, + smallest_patch_size=None, largest_patch_size=20, cells=120, time_step=0.001, model_init={}, - dl=0.1, extra_diag_options={}, time_step_nbr=1, timestamps=None): + dl=0.1, extra_diag_options={}, time_step_nbr=1, timestamps=None, ndim=1): + cells = np_array_ify(cells, ndim) from pyphare.pharein import global_vars global_vars.sim = None + + if smallest_patch_size is None: + from pyphare.pharein.simulation import check_patch_size + _, smallest_patch_size = check_patch_size(ndim, interp_order=interp_order, cells=cells) + startMPI() extra_diag_options["mode"] = "overwrite" extra_diag_options["dir"] = diag_outputs @@ -41,8 +48,8 @@ def getHierarchy(self, interp_order, refinement_boxes, qty, nbr_part_per_cell=10 time_step_nbr=time_step_nbr, time_step=time_step, boundary_types="periodic", - cells=cells, - dl=dl, + cells=np_array_ify(cells, ndim), + dl=np_array_ify(dl, ndim), interp_order=interp_order, refinement_boxes=refinement_boxes, diag_options={"format": "phareh5", @@ -210,9 +217,9 @@ def _test_overlaped_fields_are_equal(self, time_step, time_step_nbr, datahier): slice2 = data2[loc_b2.lower[0]:loc_b2.upper[0] + 1] try: - np.testing.assert_allclose(slice1, slice2, atol=1e-6) + np.testing.assert_equal(slice1, slice2) except AssertionError as e: - print("error", coarsest_time, overlap) + print("error", pd1.name, coarsest_time, overlap) raise e self.assertGreater(check, time_step_nbr) @@ -235,16 +242,15 @@ def test_overlaped_fields_are_equal(self, refinement_boxes): @data( + {}, {"L0": [Box1D(10, 19)]}, - # {"L0": [Box2D(10, 19)]}, - # {"L0": [Box3D(10, 19)]}, ) def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self, refinement_boxes): from pyphare.pharein.simulation import check_patch_size - dim = refinement_boxes["L0"][0].ndim + dim = 1 - cells = [30] * dim + cells = [60] * dim time_step_nbr=3 time_step=0.001 diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts{self.ddt_test_id()}" @@ -252,13 +258,12 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self, largest_patch_size, smallest_patch_size = check_patch_size(dim, interp_order=interp_order, cells=cells) datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size, - time_step=time_step, time_step_nbr=time_step_nbr) + time_step=time_step, time_step_nbr=time_step_nbr, cells=cells) self._test_overlaped_fields_are_equal(time_step, time_step_nbr, datahier) - def _test_patch_ghost_particle_are_clone_of_overlaped_patch_domain_particles(self, dim, interp_order, refinement_boxes): print("test_patch_ghost_particle_are_clone_of_overlaped_patch_domain_particles") print("interporder : {}".format(interp_order)) @@ -452,7 +457,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=30, diag_outputs=diag_outputs, time_step=0.001, extra_diag_options={"fine_dump_lvl_max": 10}, - time_step_nbr=time_step_nbr, smallest_patch_size=5, + time_step_nbr=time_step_nbr, largest_patch_size=30) lvl_steps = global_vars.sim.level_time_steps @@ -556,7 +561,7 @@ def _test_field_level_ghosts_via_subcycles_and_coarser_interpolation(self, ndim, def _getHier(diag_dir, boxes=[]): return self.getHierarchy(interp_order, boxes, "eb", cells=30, - time_step_nbr=1, smallest_patch_size=5, largest_patch_size=30, + time_step_nbr=1, largest_patch_size=30, diag_outputs=diag_dir, extra_diag_options={"fine_dump_lvl_max": 10}, time_step=0.001, model_init={"seed": rando} ) @@ -657,7 +662,7 @@ def test_hierarchy_timestamp_cadence(self, refinement_boxes): diag_outputs=f"phare_outputs_hierarchy_timestamp_cadence_{self.ddt_test_id()}_{i}" hier = self.getHierarchy(interp_order=1, refinement_boxes=refinement_boxes, qty="eb", cells=30, diag_outputs=diag_outputs, time_step=time_step, - time_step_nbr=time_step_nbr, smallest_patch_size=5, + time_step_nbr=time_step_nbr, largest_patch_size=30, timestamps=timestamps) time_hier_keys = list(hier.time_hier.keys()) diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 75243a74d..e48a36eeb 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -528,8 +528,8 @@ def _test_patch_ghost_on_refined_level_case(self, has_patch_ghost, **kwargs): local_out = f"{out}/dim{dim}_interp{interp}_mpi_n_{cpp.mpi_size()}_id{test_id}/{str(has_patch_ghost)}" kwargs["diag_outputs"] = local_out - - datahier = self.getHierarchy(interp, refinement_boxes, "particles_patch_ghost", **kwargs) + kwargs["interp_order"] = kwargs.get("interp_order", 1) + datahier = self.getHierarchy(refinement_boxes=refinement_boxes, qty="particles_patch_ghost", **kwargs) self.assertTrue(any([diagInfo.quantity.endswith("patchGhost") for diagInfo in ph.global_vars.sim.diagnostics])) @@ -551,11 +551,15 @@ def test_no_patch_ghost_on_refined_level_case(self, simInput): _has_patch_ghost_on_refined_level_case = ( { "cells": 40, - "smallest_patch_size": 5, - "largest_patch_size": 5}, + "interp_order": 1 + }, ) @data(*_has_patch_ghost_on_refined_level_case) def test_has_patch_ghost_on_refined_level_case(self, simInput): + from pyphare.pharein.simulation import check_patch_size + _, smallest_patch_size = check_patch_size(ndim=1, **simInput) + simInput["smallest_patch_size"] = smallest_patch_size + simInput["largest_patch_size"] = smallest_patch_size self._test_patch_ghost_on_refined_level_case(True, **simInput)