From 2af0e650b09b2d0ee1e0c1d4b328d76eeed209a5 Mon Sep 17 00:00:00 2001 From: PhilipDeegan Date: Mon, 10 Feb 2025 13:49:38 +0100 Subject: [PATCH] 3dish --- .github/workflows/cmake_macos.yml | 2 +- pyphare/pyphare/core/box.py | 40 +++ pyphare/pyphare/pharesee/geometry.py | 129 +++++--- .../pyphare/pharesee/hierarchy/hierarchy.py | 33 ++ pyphare/pyphare/pharesee/run/utils.py | 2 +- res/amr/splitting.yml | 40 +++ res/cmake/options.cmake | 10 + .../coarsening/magnetic_field_coarsener.hpp | 238 ++++++++++---- src/amr/data/particles/refine/split.hpp | 1 + src/amr/data/particles/refine/split_1d.hpp | 24 +- src/amr/data/particles/refine/split_2d.hpp | 52 +-- src/amr/data/particles/refine/split_3d.hpp | 302 ++++++++++++++++++ src/amr/data/particles/refine/splitter.hpp | 35 +- src/core/data/grid/gridlayoutdefs.hpp | 3 +- src/core/data/ndarray/ndarray_vector.hpp | 162 ++++++++-- src/core/utilities/box/box.hpp | 2 + src/core/utilities/meta/meta_utilities.hpp | 21 +- .../data/field/refine/test_refine_field.py | 3 + .../copy/test_particles_data_copy.cpp | 90 ++++++ .../refine/input/input_3d_ratio_2.txt | 51 +++ .../amr/data/particles/refine/test_split.cpp | 2 +- .../particles/stream_pack/test_no_ghost.cpp | 297 +++++++++++++++++ tests/core/data/ndarray/test_main.cpp | 54 ++++ .../core/numerics/interpolator/test_main.cpp | 67 ++++ tests/diagnostic/CMakeLists.txt | 2 + tests/diagnostic/job_3d.py.in | 14 + tests/diagnostic/test-diagnostics_3d.cpp | 35 ++ tests/functional/harris/harris_3d.py | 188 +++++++++++ tests/simulator/__init__.py | 17 +- tests/simulator/advance/CMakeLists.txt | 5 + .../advance/test_fields_advance_3d.py | 107 +++++++ .../advance/test_particles_advance_3d.py | 51 +++ tests/simulator/initialize/CMakeLists.txt | 5 + .../initialize/test_fields_init_2d.py | 2 +- .../initialize/test_fields_init_3d.py | 49 +++ .../initialize/test_particles_init_1d.py | 3 - .../initialize/test_particles_init_2d.py | 2 +- .../initialize/test_particles_init_3d.py | 98 ++++++ tests/simulator/per_test.hpp | 8 + tests/simulator/refined_particle_nbr.py | 120 +++---- tests/simulator/test_advance.py | 56 ++-- tests/simulator/test_initialization.py | 275 +++++++--------- 42 files changed, 2240 insertions(+), 457 deletions(-) create mode 100644 src/amr/data/particles/refine/split_3d.hpp create mode 100644 tests/amr/data/particles/copy/test_particles_data_copy.cpp create mode 100644 tests/amr/data/particles/refine/input/input_3d_ratio_2.txt create mode 100644 tests/amr/data/particles/stream_pack/test_no_ghost.cpp create mode 100644 tests/diagnostic/job_3d.py.in create mode 100644 tests/diagnostic/test-diagnostics_3d.cpp create mode 100644 tests/functional/harris/harris_3d.py create mode 100644 tests/simulator/advance/test_fields_advance_3d.py create mode 100644 tests/simulator/advance/test_particles_advance_3d.py create mode 100644 tests/simulator/initialize/test_fields_init_3d.py create mode 100644 tests/simulator/initialize/test_particles_init_3d.py diff --git a/.github/workflows/cmake_macos.yml b/.github/workflows/cmake_macos.yml index 42df88708..0cbec2501 100644 --- a/.github/workflows/cmake_macos.yml +++ b/.github/workflows/cmake_macos.yml @@ -78,7 +78,7 @@ jobs: -DCMAKE_BUILD_TYPE=RelWithDebInfo \ -DENABLE_SAMRAI_TESTS=OFF -DCMAKE_C_COMPILER_LAUNCHER=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DlowResourceTests=ON \ - -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 " + -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -DPHARE_SIMULATORS=2 " - name: Build working-directory: ${{runner.workspace}}/build diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py index 9b59b3a1c..42da5084a 100644 --- a/pyphare/pyphare/core/box.py +++ b/pyphare/pyphare/core/box.py @@ -62,6 +62,9 @@ def __eq__(self, other): ) def __sub__(self, other): + if isinstance(other, (list, tuple)): + assert all([isinstance(item, Box) for item in other]) + return remove_all(self, other) assert isinstance(other, Box) return remove(self, other) @@ -179,5 +182,42 @@ def copy(arr, replace): return list(boxes.values()) +def remove_all(box, to_remove): + if len(to_remove) > 0: + remaining = box - to_remove[0] + for to_rm in to_remove[1:]: + tmp, remove = [], [] + for i, rem in enumerate(remaining): + if rem * to_rm is not None: + remove.append(i) + tmp += rem - to_rm + for rm in reversed(remove): + del remaining[rm] + remaining += tmp + return remaining + return box + + + def amr_to_local(box, ref_box): return Box(box.lower - ref_box.lower, box.upper - ref_box.lower) + + +def select(data, box): + return data[tuple([slice(l, u + 1) for l,u in zip(box.lower, box.upper)])] + +class DataSelector: + """ + can't assign return to function unless [] + usage + DataSelector(data)[box] = val + """ + def __init__(self, data): + self.data = data + def __getitem__(self, box_or_slice): + if isinstance(box_or_slice, Box): + return select(self.data, box_or_slice) + return self.data[box_or_slice] + + def __setitem__(self, box_or_slice, val): + self.__getitem__(box_or_slice)[:] = val diff --git a/pyphare/pyphare/pharesee/geometry.py b/pyphare/pyphare/pharesee/geometry.py index e22631f1c..46ce5b063 100644 --- a/pyphare/pyphare/pharesee/geometry.py +++ b/pyphare/pyphare/pharesee/geometry.py @@ -54,19 +54,24 @@ def domain_border_ghost_boxes(domain_box, patches): elif domain_box.ndim == 2: upper_x, upper_y = domain_box.upper return { - "bottom": Box( - ( - 0, - 0, - ), - (upper_x, ghost_box_width), - ), - "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)), "left": Box((0, 0), (ghost_box_width, upper_y)), "right": Box((upper_x - ghost_box_width, 0), (upper_x, upper_y)), + "bottom": Box((0, 0), (upper_x, ghost_box_width)), + "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)), } - raise ValueError("Unhandeled dimension") + else: + upper_x, upper_y, upper_z = domain_box.upper + return { + "left": Box((0, 0, 0), (ghost_box_width, upper_y, upper_z)), + "right": Box( + (upper_x - ghost_box_width, 0, 0), (upper_x, upper_y, upper_z) + ), + "bottom": Box((0, 0, 0), (upper_x, ghost_box_width, upper_z)), + "top": Box((0, upper_y - ghost_box_width, 0), (upper_x, upper_y, upper_z)), + "front": Box((0, 0, 0), (upper_x, upper_y, ghost_box_width)), + "back": Box((0, 0, upper_z - ghost_box_width), (upper_x, upper_y, upper_z)), + } def touch_domain_border(box, domain_box, border): @@ -79,21 +84,29 @@ def touch_domain_border(box, domain_box, border): def periodicity_shifts(domain_box): + shifts = {} + if domain_box.ndim == 1: shape_x = domain_box.shape - return { + shifts = { "left": shape_x, "right": -shape_x, } + shifts.update({"leftright": [shifts["left"], shifts["right"]]}) if domain_box.ndim == 2: shape_x, shape_y = domain_box.shape + if domain_box.ndim == 3: + shape_x, shape_y, shape_z = domain_box.shape + + if domain_box.ndim > 1: shifts = { "left": [(shape_x, 0)], "right": [(-shape_x, 0)], "bottom": [(0, shape_y)], "top": [(0, -shape_y)], } + shifts.update( { "bottomleft": [*shifts["left"], *shifts["bottom"], (shape_x, shape_y)], @@ -134,7 +147,7 @@ def periodicity_shifts(domain_box): shifts["topleft"][-1], shifts["topright"][-1], ], - "bottomtopleftright": [ # one patch covers domain + "leftrightbottomtop": [ # one patch covers domain *shifts["bottomleft"], *shifts["topright"], shifts["bottomright"][-1], @@ -144,7 +157,35 @@ def periodicity_shifts(domain_box): ) if domain_box.ndim == 3: - raise ValueError("Unhandeled dimension") + front = { + f"{k}front": [(v[0], v[1], shape_z) for v in l] for k, l in shifts.items() + } + back = { + f"{k}back": [([v[0], v[1], -shape_z]) for v in l] for k, l in shifts.items() + } + + shifts = {k: [([v[0], v[1], 0]) for v in l] for k, l in shifts.items()} + + shifts.update(front) + shifts.update(back) + shifts.update( + { + "back": [(0, 0, -shape_z)], + "front": [(0, 0, shape_z)], + "leftrightbottomtopfrontback": [ + *shifts["bottomleftfront"], + *shifts["bottomrightback"], + *shifts["topleftfront"], + *shifts["toprightback"], + ], + } + ) + + assert len(list(shifts.keys())) == len( + set(["".join(sorted(k)) for k in list(shifts.keys())]) + ) + + shifts = {"".join(sorted(k)): l for k, l in shifts.items()} return shifts @@ -246,7 +287,7 @@ def borders_per(patch): ] for patch_i, ref_patch in enumerate(border_patches): - in_sides = borders_per_patch[ref_patch] + in_sides = "".join(sorted(borders_per_patch[ref_patch])) assert in_sides in shifts for ref_pdname, ref_pd in ref_patch.patch_datas.items(): @@ -336,36 +377,41 @@ def get_periodic_list(patches, domain_box, n_ghosts): shift_patch(first_patch, domain_box.shape) sorted_patches.append(first_patch) + return sorted_patches + + dbu = domain_box.upper + if dim == 2: sides = { - "bottom": Box([0, 0], [domain_box.upper[0], 0]), - "top": Box( - [0, domain_box.upper[1]], [domain_box.upper[0], domain_box.upper[1]] - ), - "left": Box([0, 0], [0, domain_box.upper[1]]), - "right": Box( - [domain_box.upper[0], 0], [domain_box.upper[0], domain_box.upper[1]] - ), + "left": Box([0, 0], [0, dbu[1]]), + "right": Box([dbu[0], 0], [dbu[0], dbu[1]]), + "bottom": Box([0, 0], [dbu[0], 0]), + "top": Box([0, dbu[1]], [dbu[0], dbu[1]]), } - shifts = periodicity_shifts(domain_box) + else: + sides = { + "left": Box([0, 0, 0], [0, dbu[1], dbu[2]]), + "right": Box([dbu[0], 0, 0], [dbu[0], dbu[1], dbu[2]]), + "bottom": Box([0, 0, 0], [dbu[0], 0, dbu[2]]), + "top": Box([0, dbu[1], 0], [dbu[0], dbu[1], dbu[2]]), + "front": Box([0, 0, 0], [dbu[0], dbu[1], 0]), + "back": Box([0, 0, dbu[2]], [dbu[0], dbu[1], dbu[2]]), + } - def borders_per(box): - return "".join( - [key for key, side in sides.items() if box * side is not None] - ) + shifts = periodicity_shifts(domain_box) - for patch in patches: - in_sides = borders_per(boxm.grow(patch.box, n_ghosts)) + def borders_per(box): + return "".join([key for key, side in sides.items() if box * side is not None]) - if in_sides in shifts: # in_sides might be empty, so no borders - for shift in shifts[in_sides]: - patch_copy = copy(patch) - shift_patch(patch_copy, shift) - sorted_patches.append(patch_copy) + for patch in patches: + in_sides = "".join(sorted(borders_per(boxm.grow(patch.box, n_ghosts)))) - if dim == 3: - raise ValueError("not yet implemented") + if in_sides in shifts: # in_sides might be empty, so no borders + for shift in shifts[in_sides]: + patch_copy = copy(patch) + shift_patch(patch_copy, shift) + sorted_patches.append(patch_copy) return sorted_patches @@ -470,18 +516,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None): check_patches = patches for gabox in ghostAreaBoxes: - remaining = gabox - check_patches[0].box - - for patch in check_patches[1:]: - tmp = [] - remove = [] - for i, rem in enumerate(remaining): - if rem * patch.box is not None: - remove.append(i) - tmp += rem - patch.box - for rm in reversed(remove): - del remaining[rm] - remaining += tmp + remaining = gabox - [p.box for p in check_patches] if ilvl not in lvl_gaboxes: lvl_gaboxes[ilvl] = {} diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index a3ebc30a6..9550aec31 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -556,11 +556,44 @@ def plot2d(self, **kwargs): return fig, ax + def plot3d(self, **kwargs): + """!HAX!""" + time = kwargs.get("time", self._default_time()) + usr_lvls = kwargs.get("levels", self.levelNbrs(time)) + default_qty = None + if len(self.quantities()) == 1: + default_qty = self.quantities()[0] + qty = kwargs.get("qty", default_qty) + for lvl_nbr, lvl in self.levels(time).items(): + if lvl_nbr not in usr_lvls: + continue + for patch in self.level(lvl_nbr, time).patches: + pdat = patch.patch_datas[qty] + primals = pdat.primal_directions() + if primals[0]: + pdat._x = pdat.x[:-1] + if primals[1]: + pdat._y = pdat.y[:-1] + pdat.dataset = pdat.dataset[:, :, int(pdat.ghost_box.shape[2] / 2)] + patch.box.lower = patch.box.lower[:-1] + patch.box.upper = patch.box.upper[:-1] + patch.box.ndim = 2 + + pdat.ghost_box.lower = pdat.ghost_box.lower[:-1] + pdat.ghost_box.upper = pdat.ghost_box.upper[:-1] + pdat.ghost_box.ndim = 2 + pdat.size = np.copy(pdat.ghost_box.shape) + pdat.layout.dl = pdat.layout.dl[:-1] + + return self.plot2d(**kwargs) # ¯\_(ツ)_/¯ + def plot(self, **kwargs): if self.ndim == 1: return self.plot1d(**kwargs) elif self.ndim == 2: return self.plot2d(**kwargs) + elif self.ndim == 3: + return self.plot3d(**kwargs) def dist_plot(self, **kwargs): """ diff --git a/pyphare/pyphare/pharesee/run/utils.py b/pyphare/pyphare/pharesee/run/utils.py index e9cb5b57d..1b7a9c499 100644 --- a/pyphare/pyphare/pharesee/run/utils.py +++ b/pyphare/pyphare/pharesee/run/utils.py @@ -359,7 +359,7 @@ def _get_rank(patchdatas, patch_id, **kwargs): layout = reference_pd.layout centering = "dual" nbrGhosts = layout.nbrGhosts(layout.interp_order, centering) - shape = grow(reference_pd.box, [nbrGhosts] * 2).shape + shape = grow(reference_pd.box, [nbrGhosts] * ndim).shape if ndim == 1: pass diff --git a/res/amr/splitting.yml b/res/amr/splitting.yml index e71aa05e8..14ecdb96d 100644 --- a/res/amr/splitting.yml +++ b/res/amr/splitting.yml @@ -99,3 +99,43 @@ dimension_2: delta: 1 2 weight: .140625 .09375 .0625 .0234375 .015625 .00390625 + +dimension_3: + interp_1: + N_particles_6: + delta: .966431 + weight: .166666 + + N_particles_12: + delta: .74823 + weight: .083333 + + N_particles_27: + delta: 1 + weight: .125 .0625 + + interp_2: + N_particles_6: + delta: 1.149658 + weight: .166666 + + N_particles_12: + delta: .888184 + weight: .083333 + + N_particles_27: + delta: 1.111333 + weight: .099995 .055301 + + interp_3: + N_particles_6: + delta: 1.312622 + weight: .166666 + + N_particles_12: + delta: 1.012756 + weight: .083333 + + N_particles_27: + delta: 1.276815 + weight: .104047 .05564 diff --git a/res/cmake/options.cmake b/res/cmake/options.cmake index 20266c84a..6f132adab 100644 --- a/res/cmake/options.cmake +++ b/res/cmake/options.cmake @@ -70,10 +70,19 @@ option(withPhlop "Use phlop" OFF) # -DlowResourceTests=ON option(lowResourceTests "Disable heavy tests for CI (2d/3d/etc" OFF) +# -DhighResourceTests=ON +option(highResourceTests "Enable heavy tests for CI (3d/etc" OFF) + +if(lowResourceTests AND highResourceTests) + message(FATAL_ERROR "Both lowResourceTests and highResourceTests cannot be enabled at the same time.") +endif() + + # -DtestDuringBuild=ON enabled if devMode=ON, disabled if asan=ON (needs LD_PRELOAD) option(testDuringBuild "Runs C++ unit tests after they are built" OFF) + # -DPGO_GEN=OFF profile guided optimization generate option(PGO_GEN "profile guided optimization generate" OFF) # -DPGO_USE=OFF @@ -84,6 +93,7 @@ option(PGO_USE "profile guided optimization use" OFF) option(PHARE_COMPILER_WORKAROUNDS "Attempt silence compiler false positives" ON) # on by default + # Controlling the activation of tests if (NOT DEFINED PHARE_EXEC_LEVEL_MIN) set(PHARE_EXEC_LEVEL_MIN 1) diff --git a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp index 39d816413..a13f134bc 100644 --- a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp +++ b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp @@ -31,110 +31,206 @@ using core::dirZ; * of the enclosed fine faces * */ -template +template class MagneticFieldCoarsener { + using Point_t = core::Point; + public: - MagneticFieldCoarsener(std::array const centering, + MagneticFieldCoarsener(std::array const centering, SAMRAI::hier::Box const& sourceBox, SAMRAI::hier::Box const& destinationBox, SAMRAI::hier::IntVector const& ratio) : centering_{centering} , sourceBox_{sourceBox} , destinationBox_{destinationBox} - { } + template - void operator()(FieldT const& fineField, FieldT& coarseField, - core::Point coarseIndex) + void operator()(FieldT const& fineField, FieldT& coarseField, core::Point coarseIndex) { // For the moment we only take the case of field with the same centering TBOX_ASSERT(fineField.physicalQuantity() == coarseField.physicalQuantity()); - core::Point fineStartIndex; + coarsen(fine_start_index(coarseIndex), fineField, coarseField, + AMRToLocal(coarseIndex, destinationBox_)); + } +private: + auto fine_start_index(core::Point const coarseIndex) const + { + core::Point fineStartIndex; fineStartIndex[dirX] = coarseIndex[dirX] * this->ratio_; - - if constexpr (dimension > 1) + if constexpr (dim > 1) { fineStartIndex[dirY] = coarseIndex[dirY] * this->ratio_; - if constexpr (dimension > 2) + if constexpr (dim > 2) { fineStartIndex[dirZ] = coarseIndex[dirZ] * this->ratio_; } } + return AMRToLocal(fineStartIndex, sourceBox_); + } - fineStartIndex = AMRToLocal(fineStartIndex, sourceBox_); - coarseIndex = AMRToLocal(coarseIndex, destinationBox_); - // the following kinda assumes where B is, i.e. Yee layout centering - // as it only does faces pirmal-dual, dual-primal and dual-dual + template + typename std::enable_if::type + coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField, + Point_t const coarseIndex); - if constexpr (dimension == 1) - { - // in 1D div(B) is automatically satisfied so using this coarsening - // opertor is probably not better than the default one, but we do that - // for a kind of consistency... - // coarse flux is equal to fine flux and we're 1D so there is flux partitioned - // only for By and Bz, Bx is equal to the fine value + template + typename std::enable_if::type + coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField, + Point_t const coarseIndex); + template + typename std::enable_if::type + coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField, + Point_t const coarseIndex); - if (centering_[dirX] == core::QtyCentering::primal) // bx - { - coarseField(coarseIndex[dirX]) = fineField(fineStartIndex[dirX]); - } - else if (centering_[dirX] == core::QtyCentering::dual) // by and bz - { - coarseField(coarseIndex[dirX]) - = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX])); - } - } - - if constexpr (dimension == 2) - { - if (centering_[dirX] == core::QtyCentering::primal - and centering_[dirY] == core::QtyCentering::dual) - { - coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.5 - * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) - + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1)); - } - else if (centering_[dirX] == core::QtyCentering::dual - and centering_[dirY] == core::QtyCentering::primal) - { - coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.5 - * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) - + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY])); - } - else if (centering_[dirX] == core::QtyCentering::dual - and centering_[dirY] == core::QtyCentering::dual) - { - coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.25 - * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) - + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]) - + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1) - + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY] + 1)); - } - else - { - throw std::runtime_error("no magnetic field should end up here"); - } - } - else if constexpr (dimension == 3) - { - throw std::runtime_error("Not Implemented yet"); - } - } -private: - std::array const centering_; + std::array const centering_; SAMRAI::hier::Box const sourceBox_; SAMRAI::hier::Box const destinationBox_; static int constexpr ratio_ = 2; }; + +template +template +typename std::enable_if::type +MagneticFieldCoarsener::coarsen(Point_t const fineStartIndex, FieldT const& fineField, + FieldT& coarseField, Point_t const coarseIndex) +{ + if (centering_[dirX] == core::QtyCentering::primal) // bx + { + coarseField(coarseIndex[dirX]) = fineField(fineStartIndex[dirX]); + } + else if (centering_[dirX] == core::QtyCentering::dual) // by and bz + { + coarseField(coarseIndex[dirX]) + = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX])); + } +} + +template +template +typename std::enable_if::type +MagneticFieldCoarsener::coarsen(Point_t const fineStartIndex, FieldT const& fineField, + FieldT& coarseField, Point_t const coarseIndex) +{ + if (centering_[dirX] == core::QtyCentering::primal + and centering_[dirY] == core::QtyCentering::dual) + { + coarseField(coarseIndex[dirX], coarseIndex[dirY]) + = 0.5 + * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1)); + } + else if (centering_[dirX] == core::QtyCentering::dual + and centering_[dirY] == core::QtyCentering::primal) + { + coarseField(coarseIndex[dirX], coarseIndex[dirY]) + = 0.5 + * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY])); + } + else if (centering_[dirX] == core::QtyCentering::dual + and centering_[dirY] == core::QtyCentering::dual) + { + coarseField(coarseIndex[dirX], coarseIndex[dirY]) + = 0.25 + * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]) + + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1) + + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY] + 1)); + } + else + { + throw std::runtime_error("no magnetic field should end up here"); + } +} + +// clang-format off + +// clang-format on + +template +template +typename std::enable_if::type +MagneticFieldCoarsener::coarsen(Point_t const fineStartIndex, FieldT const& fineField, + FieldT& coarseField, Point_t const coarseIndex) +{ + auto& ci = coarseIndex; + auto& cf = coarseField; + auto& fsi = fineStartIndex; + auto& ff = fineField; + + using Fn = std::function; + using Fn_arr = std::array; + using A0_arr = std::array; + + // X = 0 = primal + // Y = 1 = dual + std::array const fns{ + A0_arr{Fn_arr{[]() { // XXX + throw std::runtime_error("no magnetic field should end up here"); + }, + [&]() { // XXY + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX], fsi[dirY], fsi[dirZ] + 1)) + * V; + }}, + Fn_arr{[&]() { // XYX + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ])) + * V; + }, + [&]() { // XYY + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ] + 1)) + * V; + }}}, + A0_arr{Fn_arr{[&]() { + // YXX + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ])) + * V; + }, + [&]() { + // YXY + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ] + 1)) + * V; + }}, + Fn_arr{[&]() { + // YYX + double constexpr static V = .5; + + cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + // + ff(fsi[dirX] + 1, fsi[dirY] + 1, fsi[dirZ])) + * V; + }, + [&]() { + // YYY + throw std::runtime_error("no magnetic field should end up here"); + }}}}; + + auto cast = [](auto const& qty) { + return static_cast>(qty); + }; + fns[cast(centering_[0])][cast(centering_[1])][cast(centering_[2])](); +} + } // namespace PHARE::amr #endif diff --git a/src/amr/data/particles/refine/split.hpp b/src/amr/data/particles/refine/split.hpp index 1d70a7689..73358bf77 100644 --- a/src/amr/data/particles/refine/split.hpp +++ b/src/amr/data/particles/refine/split.hpp @@ -3,6 +3,7 @@ #include "amr/data/particles/refine/split_1d.hpp" #include "amr/data/particles/refine/split_2d.hpp" +#include "amr/data/particles/refine/split_3d.hpp" #endif // endif PHARE_SPLIT_HPP diff --git a/src/amr/data/particles/refine/split_1d.hpp b/src/amr/data/particles/refine/split_1d.hpp index beef8ff26..c0aade255 100644 --- a/src/amr/data/particles/refine/split_1d.hpp +++ b/src/amr/data/particles/refine/split_1d.hpp @@ -17,11 +17,11 @@ using namespace PHARE::core; /**************************************************************************/ template<> -struct PinkDispatcher> : SplitPattern, RefinedParticlesConst<2>> +struct PinkPattern> : SplitPattern, RefinedParticlesConst<2>> { using Super = SplitPattern, RefinedParticlesConst<2>>; - constexpr PinkDispatcher(float const weight, float const delta) + constexpr PinkPattern(float const weight, float const delta) : Super{weight} { Super::deltas_[0] = {delta}; @@ -31,7 +31,7 @@ struct PinkDispatcher> : SplitPattern, RefinedParticlesC /**************************************************************************/ -using SplitPattern_1_1_2_Dispatcher = PatternDispatcher>>; +using SplitPattern_1_1_2_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<2>> @@ -50,7 +50,7 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<2>> /**************************************************************************/ using SplitPattern_1_1_3_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<3>> @@ -68,7 +68,7 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<3>> /**************************************************************************/ -using SplitPattern_1_2_2_Dispatcher = PatternDispatcher>>; +using SplitPattern_1_2_2_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<2>> @@ -87,7 +87,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<2>> /**************************************************************************/ using SplitPattern_1_2_3_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<3>> @@ -106,7 +106,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<3>> /**************************************************************************/ using SplitPattern_1_2_4_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<4>> @@ -124,7 +124,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<4>> /**************************************************************************/ -using SplitPattern_1_3_2_Dispatcher = PatternDispatcher>>; +using SplitPattern_1_3_2_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<2>> @@ -143,7 +143,7 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<2>> /**************************************************************************/ using SplitPattern_1_3_3_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<3>> @@ -162,7 +162,7 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<3>> /**************************************************************************/ using SplitPattern_1_3_4_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<4>> @@ -181,8 +181,8 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<4>> /**************************************************************************/ using SplitPattern_1_3_5_Dispatcher - = PatternDispatcher>, PinkDispatcher>, - PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>, + PinkPattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<5>> diff --git a/src/amr/data/particles/refine/split_2d.hpp b/src/amr/data/particles/refine/split_2d.hpp index 53a0bfa3f..f082b98ab 100644 --- a/src/amr/data/particles/refine/split_2d.hpp +++ b/src/amr/data/particles/refine/split_2d.hpp @@ -18,11 +18,11 @@ using namespace PHARE::core; /**************************************************************************/ template<> -struct PurpleDispatcher> : SplitPattern, RefinedParticlesConst<4>> +struct PurplePattern> : SplitPattern, RefinedParticlesConst<4>> { using Super = SplitPattern, RefinedParticlesConst<4>>; - constexpr PurpleDispatcher(float const weight, float const delta) + constexpr PurplePattern(float const weight, float const delta) : Super{weight} { Super::deltas_[0] = {-delta, -delta}; @@ -34,12 +34,12 @@ struct PurpleDispatcher> : SplitPattern, RefinedParticle template<> -struct BrownDispatcher> : SplitPattern, RefinedParticlesConst<8>> +struct BrownPattern> : SplitPattern, RefinedParticlesConst<8>> { using Super = SplitPattern, RefinedParticlesConst<8>>; template - constexpr BrownDispatcher(float const weight, Delta const& delta) + constexpr BrownPattern(float const weight, Delta const& delta) : Super{weight} { for (std::size_t i = 0; i < 2; i++) @@ -55,11 +55,11 @@ struct BrownDispatcher> : SplitPattern, RefinedParticles }; template<> -struct PinkDispatcher> : SplitPattern, RefinedParticlesConst<4>> +struct PinkPattern> : SplitPattern, RefinedParticlesConst<4>> { using Super = SplitPattern, RefinedParticlesConst<4>>; - constexpr PinkDispatcher(float const weight, float const delta) + constexpr PinkPattern(float const weight, float const delta) : Super{weight} { Super::deltas_[0] = {0.0f, -delta}; @@ -71,7 +71,7 @@ struct PinkDispatcher> : SplitPattern, RefinedParticlesC /**************************************************************************/ -using SplitPattern_2_1_4_Dispatcher = PatternDispatcher>>; +using SplitPattern_2_1_4_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<4>> @@ -90,7 +90,7 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<4>> /**************************************************************************/ using SplitPattern_2_1_5_Dispatcher - = PatternDispatcher>, PurpleDispatcher>>; + = PatternDispatcher>, PurplePattern>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<5>> @@ -109,7 +109,7 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<5>> /**************************************************************************/ using SplitPattern_2_1_8_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<8>> @@ -128,8 +128,8 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<8>> /**************************************************************************/ using SplitPattern_2_1_9_Dispatcher - = PatternDispatcher>, PinkDispatcher>, - PurpleDispatcher>>; + = PatternDispatcher>, PinkPattern>, + PurplePattern>>; template<> struct Splitter, InterpConst<1>, RefinedParticlesConst<9>> @@ -147,7 +147,7 @@ struct Splitter, InterpConst<1>, RefinedParticlesConst<9>> /**************************************************************************/ -using SplitPattern_2_2_4_Dispatcher = PatternDispatcher>>; +using SplitPattern_2_2_4_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<4>> @@ -166,7 +166,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<4>> /**************************************************************************/ using SplitPattern_2_2_5_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<5>> @@ -185,7 +185,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<5>> /**************************************************************************/ using SplitPattern_2_2_8_Dispatcher - = PatternDispatcher>, PurpleDispatcher>>; + = PatternDispatcher>, PurplePattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<8>> @@ -204,8 +204,8 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<8>> /**************************************************************************/ using SplitPattern_2_2_9_Dispatcher - = PatternDispatcher>, PinkDispatcher>, - PurpleDispatcher>>; + = PatternDispatcher>, PinkPattern>, + PurplePattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<9>> @@ -224,8 +224,8 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<9>> /**************************************************************************/ using SplitPattern_2_2_16_Dispatcher - = PatternDispatcher>, BrownDispatcher>, - PurpleDispatcher>>; + = PatternDispatcher>, BrownPattern>, + PurplePattern>>; template<> struct Splitter, InterpConst<2>, RefinedParticlesConst<16>> @@ -244,7 +244,7 @@ struct Splitter, InterpConst<2>, RefinedParticlesConst<16>> /**************************************************************************/ -using SplitPattern_2_3_4_Dispatcher = PatternDispatcher>>; +using SplitPattern_2_3_4_Dispatcher = PatternDispatcher>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<4>> @@ -263,7 +263,7 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<4>> /**************************************************************************/ using SplitPattern_2_3_5_Dispatcher - = PatternDispatcher>, PinkDispatcher>>; + = PatternDispatcher>, PinkPattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<5>> @@ -282,7 +282,7 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<5>> /**************************************************************************/ using SplitPattern_2_3_8_Dispatcher - = PatternDispatcher>, PurpleDispatcher>>; + = PatternDispatcher>, PurplePattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<8>> @@ -301,8 +301,8 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<8>> /**************************************************************************/ using SplitPattern_2_3_9_Dispatcher - = PatternDispatcher>, PinkDispatcher>, - PurpleDispatcher>>; + = PatternDispatcher>, PinkPattern>, + PurplePattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<9>> @@ -321,9 +321,9 @@ struct Splitter, InterpConst<3>, RefinedParticlesConst<9>> /**************************************************************************/ using SplitPattern_2_3_25_Dispatcher - = PatternDispatcher>, PinkDispatcher>, - PurpleDispatcher>, PinkDispatcher>, - BrownDispatcher>, PurpleDispatcher>>; + = PatternDispatcher>, PinkPattern>, + PurplePattern>, PinkPattern>, + BrownPattern>, PurplePattern>>; template<> struct Splitter, InterpConst<3>, RefinedParticlesConst<25>> diff --git a/src/amr/data/particles/refine/split_3d.hpp b/src/amr/data/particles/refine/split_3d.hpp new file mode 100644 index 000000000..508907abc --- /dev/null +++ b/src/amr/data/particles/refine/split_3d.hpp @@ -0,0 +1,302 @@ +/* +Splitting reference material can be found @ + https://github.com/PHAREHUB/PHARE/wiki/SplitPattern +*/ + +#ifndef PHARE_SPLIT_3D_HPP +#define PHARE_SPLIT_3D_HPP + +#include +#include +#include "core/utilities/point/point.hpp" +#include "core/utilities/types.hpp" +#include "splitter.hpp" + +namespace PHARE::amr +{ +using namespace PHARE::core; + +/**************************************************************************/ +template<> // 1 per face - centered +struct PinkPattern> : SplitPattern, RefinedParticlesConst<6>> +{ + using Super = SplitPattern, RefinedParticlesConst<6>>; + + constexpr PinkPattern(float const weight, float const delta) + : Super{weight} + { + constexpr float zero = 0; + + for (std::size_t i = 0; i < 2; i++) + { + std::size_t offset = i * 3; + float sign = i % 2 ? -1 : 1; + auto mode = delta * sign; + + Super::deltas_[0 + offset] = {mode, zero, zero}; + Super::deltas_[1 + offset] = {zero, zero, mode}; + Super::deltas_[2 + offset] = {zero, mode, zero}; + } + } +}; + +template<> // 1 per corner +struct PurplePattern> : SplitPattern, RefinedParticlesConst<8>> +{ + using Super = SplitPattern, RefinedParticlesConst<8>>; + + constexpr PurplePattern(float const weight, float const delta) + : Super{weight} + { + for (std::size_t i = 0; i < 2; i++) + { + std::size_t offset = i * 4; + float sign = i % 2 ? -1 : 1; + auto mode = delta * sign; + + Super::deltas_[0 + offset] = {mode, mode, mode}; + Super::deltas_[1 + offset] = {mode, mode, -mode}; + Super::deltas_[2 + offset] = {mode, -mode, mode}; + Super::deltas_[3 + offset] = {mode, -mode, -mode}; + } + } +}; + + +template<> // 1 per edge - centered +struct LimePattern> : SplitPattern, RefinedParticlesConst<12>> +{ + using Super = SplitPattern, RefinedParticlesConst<12>>; + + constexpr LimePattern(float const weight, float const delta) + : Super{weight} + { + constexpr float zero = 0; + + auto addSquare = [&delta, this](size_t offset, float y) { + Super::deltas_[0 + offset] = {zero, y, delta}; + Super::deltas_[1 + offset] = {zero, y, -delta}; + Super::deltas_[2 + offset] = {delta, y, zero}; + Super::deltas_[3 + offset] = {-delta, y, zero}; + }; + + addSquare(0, delta); // top + addSquare(4, -delta); // bottom + addSquare(8, 0); // middle + } +}; + + +template<> // 2 per edge - equidistant from center +class WhitePattern> : SplitPattern, RefinedParticlesConst<24>> +{ + using Super = SplitPattern, RefinedParticlesConst<24>>; + + template + constexpr WhitePattern(float const weight, Delta const& delta) + : Super{weight} + { + auto addSquare = [&delta, this](size_t offset, float x, float y, float z) { + Super::deltas_[0 + offset] = {x, y, z}; + Super::deltas_[1 + offset] = {x, -y, z}; + Super::deltas_[2 + offset] = {-x, y, z}; + Super::deltas_[3 + offset] = {-x, -y, z}; + }; + + auto addSquares = [addSquare, &delta, this](size_t offset, float x, float y, float z) { + addSquare(offset, x, y, z); + addSquare(offset + 4, x, y, -z); + }; + + addSquares(0, delta[0], delta[0], delta[1]); + addSquares(8, delta[0], delta[1], delta[0]); + addSquares(16, delta[1], delta[0], delta[0]); + } +}; + + +/**************************************************************************/ +/****************************** INTERP == 1 *******************************/ +/**************************************************************************/ +using SplitPattern_3_1_6_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<1>, RefinedParticlesConst<6>> + : public ASplitter, InterpConst<1>, RefinedParticlesConst<6>>, + SplitPattern_3_1_6_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_1_6_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {.966431}; + static constexpr std::array weight = {.166666}; +}; + + +/**************************************************************************/ +using SplitPattern_3_1_12_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<1>, RefinedParticlesConst<12>> + : public ASplitter, InterpConst<1>, RefinedParticlesConst<12>>, + SplitPattern_3_1_12_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_1_12_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {.74823}; + static constexpr std::array weight = {.083333}; +}; + + +/**************************************************************************/ +using SplitPattern_3_1_27_Dispatcher + = PatternDispatcher>, PinkPattern>, + PurplePattern>, LimePattern>>; + +template<> +struct Splitter, InterpConst<1>, RefinedParticlesConst<27>> + : public ASplitter, InterpConst<1>, RefinedParticlesConst<27>>, + SplitPattern_3_1_27_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_1_27_Dispatcher{ + {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}} + { + } + + static constexpr std::array delta = {1}; + static constexpr std::array weight = {0.125, 0.0625}; +}; + + + +/**************************************************************************/ +/****************************** INTERP == 2 *******************************/ +/**************************************************************************/ +using SplitPattern_3_2_6_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<2>, RefinedParticlesConst<6>> + : public ASplitter, InterpConst<2>, RefinedParticlesConst<6>>, + SplitPattern_3_2_6_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_2_6_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {1.149658}; + static constexpr std::array weight = {.166666}; +}; + + +/**************************************************************************/ +using SplitPattern_3_2_12_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<2>, RefinedParticlesConst<12>> + : public ASplitter, InterpConst<2>, RefinedParticlesConst<12>>, + SplitPattern_3_2_12_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_2_12_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {.888184}; + static constexpr std::array weight = {.083333}; +}; + + +/**************************************************************************/ +using SplitPattern_3_2_27_Dispatcher + = PatternDispatcher>, PinkPattern>, + PurplePattern>, LimePattern>>; + +template<> +struct Splitter, InterpConst<2>, RefinedParticlesConst<27>> + : public ASplitter, InterpConst<2>, RefinedParticlesConst<27>>, + SplitPattern_3_2_27_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_2_27_Dispatcher{ + {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}} + { + } + + static constexpr std::array delta = {1.111333}; + static constexpr std::array weight = {.099995, .055301}; +}; + + +/**************************************************************************/ +/****************************** INTERP == 3 *******************************/ +/**************************************************************************/ +using SplitPattern_3_3_6_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<3>, RefinedParticlesConst<6>> + : public ASplitter, InterpConst<3>, RefinedParticlesConst<6>>, + SplitPattern_3_3_6_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_3_6_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {1.312622}; + static constexpr std::array weight = {.166666}; +}; + + +/**************************************************************************/ +using SplitPattern_3_3_12_Dispatcher = PatternDispatcher>>; + +template<> +struct Splitter, InterpConst<3>, RefinedParticlesConst<12>> + : public ASplitter, InterpConst<3>, RefinedParticlesConst<12>>, + SplitPattern_3_3_12_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_3_12_Dispatcher{{weight[0], delta[0]}} + { + } + + static constexpr std::array delta = {1.012756}; + static constexpr std::array weight = {.083333}; +}; + + +/**************************************************************************/ +using SplitPattern_3_3_27_Dispatcher + = PatternDispatcher>, PinkPattern>, + PurplePattern>, LimePattern>>; + +template<> +struct Splitter, InterpConst<3>, RefinedParticlesConst<27>> + : public ASplitter, InterpConst<3>, RefinedParticlesConst<27>>, + SplitPattern_3_3_27_Dispatcher +{ + constexpr Splitter() + : SplitPattern_3_3_27_Dispatcher{ + {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}} + { + } + + static constexpr std::array delta = {1.276815}; + static constexpr std::array weight = {.104047, .05564}; +}; + + +/**************************************************************************/ + + +} // namespace PHARE::amr + + +#endif /* PHARE_SPLIT_3D_HPP */ diff --git a/src/amr/data/particles/refine/splitter.hpp b/src/amr/data/particles/refine/splitter.hpp index accd4d387..cf252f40e 100644 --- a/src/amr/data/particles/refine/splitter.hpp +++ b/src/amr/data/particles/refine/splitter.hpp @@ -6,7 +6,6 @@ #include #include #include -#include #include #include #include @@ -67,17 +66,17 @@ class PatternDispatcher using FineParticle = decltype(particles[0]); // may be a reference core::apply(patterns, [&](auto const& pattern) { + auto const weight = static_cast(pattern.weight_); for (size_t rpIndex = 0; rpIndex < pattern.deltas_.size(); rpIndex++) { FineParticle fineParticle = particles[idx++]; - fineParticle.weight = particle.weight * static_cast(pattern.weight_) - * power[dimension - 1]; - fineParticle.charge = particle.charge; - fineParticle.iCell = particle.iCell; - fineParticle.delta = particle.delta; - fineParticle.v = particle.v; - - for (size_t iDim = 0; iDim < dimension; iDim++) + fineParticle.weight = particle.weight * weight * power[dimension - 1]; + fineParticle.charge = particle.charge; + fineParticle.iCell = particle.iCell; + fineParticle.delta = particle.delta; + fineParticle.v = particle.v; + + for (size_t iDim = 0; iDim < dimension; ++iDim) { fineParticle.delta[iDim] += static_cast(pattern.deltas_[rpIndex][iDim]); @@ -114,25 +113,33 @@ class Splitter : public ASplitter }; template -struct BlackDispatcher : SplitPattern> +struct BlackPattern : SplitPattern> { using Super = SplitPattern>; - constexpr BlackDispatcher(float const weight) + constexpr BlackPattern(float const weight) : Super{weight} { } }; template -struct PurpleDispatcher +struct PinkPattern +{ +}; +template +struct PurplePattern +{ +}; +template +struct BrownPattern { }; template -struct BrownDispatcher +struct LimePattern { }; template -struct PinkDispatcher +struct WhitePattern { }; diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp index 23920396e..f5cef3bba 100644 --- a/src/core/data/grid/gridlayoutdefs.hpp +++ b/src/core/data/grid/gridlayoutdefs.hpp @@ -2,6 +2,7 @@ #define PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP #include +#include #include "core/hybrid/hybrid_quantities.hpp" #include "core/utilities/types.hpp" @@ -14,7 +15,7 @@ namespace core enum class Direction { X, Y, Z }; - enum class QtyCentering { primal = 0, dual = 1 }; + enum class QtyCentering : std::uint16_t { primal = 0, dual = 1 }; template diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index 7d52d1d01..98898b796 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -3,6 +3,7 @@ #include "core/def.hpp" #include +#include #include #include #include @@ -123,6 +124,10 @@ class MaskedView NO_DISCARD auto yend() const { return shape_[1] - 1 - mask_.max(); } + auto zstart() const { return mask_.min(); } + auto zend() const { return shape_[2] - 1 - mask_.max(); } + + private: Array& array_; std::array shape_; @@ -344,10 +349,10 @@ class NdArrayMask { auto shape = array.shape(); - for (std::size_t i = min_; i <= max_; ++i) + for (auto i = min_; i <= max_; ++i) array(i) = val; - for (std::size_t i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i) + for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i) array(i) = val; } @@ -357,24 +362,24 @@ class NdArrayMask auto shape = array.shape(); // left border - for (std::size_t i = min_; i <= max_; ++i) - for (std::size_t j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto i = min_; i <= max_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) array(i, j) = val; // right border - for (std::size_t i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i) - for (std::size_t j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) array(i, j) = val; - for (std::size_t i = min_; i <= shape[0] - 1 - min_; ++i) + for (auto i = min_; i <= shape[0] - 1 - min_; ++i) { // bottom border - for (std::size_t j = min_; j <= max_; ++j) + for (auto j = min_; j <= max_; ++j) array(i, j) = val; // top border - for (std::size_t j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j) + for (auto j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j) array(i, j) = val; } } @@ -382,7 +387,44 @@ class NdArrayMask template void fill3D(Array& array, typename Array::type val) const { - throw std::runtime_error("3d not implemented"); + auto shape = array.shape(); + + // left border + for (auto i = min_; i <= shape[0] - 1 - max_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto k = min_; k <= max_; ++k) + array(i, j, k) = val; + + // // right border + for (auto i = min_; i <= shape[0] - 1 - max_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto k = shape[2] - 1 - max_; k <= shape[2] - 1 - min_; ++k) + array(i, j, k) = val; + + for (auto i = min_; i <= shape[0] - 1 - min_; ++i) + { + // bottom border + for (auto j = min_; j <= max_; ++j) + for (auto k = min_; k <= shape[2] - 1 - min_; ++k) + array(i, j, k) = val; + + // top border + for (auto j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j) + for (auto k = min_; k <= shape[2] - 1 - min_; ++k) + array(i, j, k) = val; + } + + // front + for (auto i = min_; i <= max_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto k = min_; k <= shape[2] - 1 - min_; ++k) + array(i, j, k) = val; + + // back + for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i) + for (auto j = min_; j <= shape[1] - 1 - max_; ++j) + for (auto k = min_; k <= shape[2] - 1 - min_; ++k) + array(i, j, k) = val; } template @@ -392,17 +434,23 @@ class NdArrayMask std::size_t cells = 0; - if constexpr (Array::dimension == 1) - for (std::size_t i = min_; i <= max_; ++i) + for (auto i = min_; i <= max_; ++i) + { + if constexpr (Array::dimension == 1) cells += 2; - if constexpr (Array::dimension == 2) - for (std::size_t i = min_; i <= max_; ++i) + if constexpr (Array::dimension == 2) cells += (shape[0] - (i * 2) - 2) * 2 + (shape[1] - (i * 2) - 2) * 2 + 4; - if constexpr (Array::dimension == 3) - throw std::runtime_error("Not implemented dimension"); - + if constexpr (Array::dimension == 3) + { + auto [x, y, z] = shape; + x -= i * 2; + y -= i * 2; + z -= i * 2; + cells += (x * y * 2) + (y * (z - 2) * 2) + ((z - 2) * (x - 2) * 2); + } + } return cells; } @@ -422,20 +470,21 @@ void operator>>(MaskedView&& inner, MaskedView&& outer { using MaskedView_t = MaskedView; + assert(inner.xstart() > outer.xstart() and inner.xend() < outer.xend()); + if constexpr (MaskedView_t::dimension == 1) { - assert(inner.xstart() > outer.xstart()); - assert(inner.xend() < outer.xend()); outer(outer.xstart()) = inner(inner.xstart()); outer(outer.xend()) = inner(inner.xend()); } + if constexpr (MaskedView_t::dimension > 1) + { + assert(inner.ystart() > outer.ystart() and inner.yend() < outer.yend()); + } if constexpr (MaskedView_t::dimension == 2) { - assert(inner.xstart() > outer.xstart() and inner.xend() < outer.xend() - and inner.ystart() > outer.ystart() and inner.yend() < outer.yend()); - for (auto ix = inner.xstart(); ix <= inner.xend(); ++ix) { outer(ix, outer.ystart()) = inner(ix, inner.ystart()); // bottom @@ -452,7 +501,7 @@ void operator>>(MaskedView&& inner, MaskedView&& outer for (auto ix = outer.xstart(); ix < inner.xstart(); ++ix) outer(ix, outer.ystart()) = inner(inner.xstart(), inner.ystart()); - for (std::size_t iy = outer.ystart(); iy < inner.ystart(); ++iy) + for (auto iy = outer.ystart(); iy < inner.ystart(); ++iy) outer(outer.xstart(), iy) = inner(inner.xstart(), inner.ystart()); @@ -481,7 +530,72 @@ void operator>>(MaskedView&& inner, MaskedView&& outer if constexpr (MaskedView_t::dimension == 3) { - throw std::runtime_error("3d not implemented"); + assert(inner.zstart() > outer.zstart() and inner.zend() < outer.zend()); + + for (auto i = inner.xstart(); i <= inner.xend(); ++i) + { + for (auto k = inner.zstart(); k <= inner.zend(); ++k) + { + outer(i, outer.ystart(), k) = inner(i, inner.ystart(), k); + outer(i, outer.yend(), k) = inner(i, inner.yend(), k); + } + for (auto j = inner.ystart(); j <= inner.yend(); ++j) + { + outer(i, j, outer.zstart()) = inner(i, j, inner.zstart()); + outer(i, j, outer.zend()) = inner(i, j, inner.zend()); + } + } + + for (auto j = inner.ystart(); j <= inner.yend(); ++j) + { + for (auto k = inner.zstart(); k <= inner.zend(); ++k) + { + outer(outer.xstart(), j, k) = inner(inner.xstart(), j, k); + outer(outer.xend(), j, k) = inner(inner.xend(), j, k); + } + } + + for (auto k = inner.zstart(); k <= inner.zend(); ++k) + { + outer(outer.xstart(), outer.ystart(), k) = inner(outer.xstart(), inner.ystart(), k); + outer(outer.xstart(), outer.yend(), k) = inner(outer.xstart(), inner.yend(), k); + + outer(outer.xend(), outer.ystart(), k) = inner(outer.xend(), inner.ystart(), k); + outer(outer.xend(), outer.yend(), k) = inner(outer.xend(), inner.yend(), k); + } + + for (auto j = inner.ystart(); j <= inner.yend(); ++j) + { + outer(outer.xstart(), j, outer.zstart()) = inner(outer.xstart(), j, inner.zstart()); + outer(outer.xstart(), j, outer.zend()) = inner(outer.xstart(), j, inner.zend()); + + outer(outer.xend(), j, outer.zstart()) = inner(outer.xend(), j, inner.zstart()); + outer(outer.xend(), j, outer.zend()) = inner(outer.xend(), j, inner.zend()); + } + + for (auto i = inner.xstart(); i <= inner.xend(); ++i) + { + outer(i, outer.ystart(), outer.zstart()) = inner(i, outer.ystart(), inner.zstart()); + outer(i, outer.ystart(), outer.zend()) = inner(i, outer.ystart(), inner.zend()); + + outer(i, outer.yend(), outer.zstart()) = inner(i, outer.yend(), inner.zstart()); + outer(i, outer.yend(), outer.zend()) = inner(i, outer.yend(), inner.zend()); + } + + auto corner = [&](auto xouter, auto xinner) { + outer(xouter, outer.ystart(), outer.zstart()) + = inner(xinner, outer.ystart(), inner.zstart()); + + outer(xouter, outer.ystart(), outer.zend()) + = inner(xinner, outer.ystart(), inner.zend()); + + outer(xouter, outer.yend(), outer.zstart()) + = inner(xinner, outer.yend(), inner.zstart()); + outer(xouter, outer.yend(), outer.zend()) = inner(xinner, outer.yend(), inner.zend()); + }; + + corner(outer.xstart(), inner.xstart()); + corner(outer.xend(), inner.xend()); } } diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index 0ff5d1a8a..ccb225fa0 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -148,6 +148,8 @@ struct Box } }; + + template class box_iterator { diff --git a/src/core/utilities/meta/meta_utilities.hpp b/src/core/utilities/meta/meta_utilities.hpp index 3bea76ba3..00c13dbc3 100644 --- a/src/core/utilities/meta/meta_utilities.hpp +++ b/src/core/utilities/meta/meta_utilities.hpp @@ -6,6 +6,10 @@ #include "core/utilities/types.hpp" +#if !defined(PHARE_SIMULATORS) +#define PHARE_SIMULATORS 3 +#endif + namespace PHARE { namespace core @@ -74,11 +78,24 @@ namespace core // inner tuple = dim, interp, list[possible nbrParticles for dim/interp] return std::tuple, InterpConst<1>, 2, 3>, SimulatorOption, InterpConst<2>, 2, 3, 4>, - SimulatorOption, InterpConst<3>, 2, 3, 4, 5>, + SimulatorOption, InterpConst<3>, 2, 3, 4, 5> +#if PHARE_SIMULATORS > 1 + , SimulatorOption, InterpConst<1>, 4, 5, 8, 9>, SimulatorOption, InterpConst<2>, 4, 5, 8, 9, 16>, - SimulatorOption, InterpConst<3>, 4, 5, 8, 9, 25>>{}; + SimulatorOption, InterpConst<3>, 4, 5, 8, 9, 25> +#endif + + // TODO add in the rest of 3d nbrParticles permutations + // possibly consider compile time activation for uncommon cases +#if PHARE_SIMULATORS > 2 + , + SimulatorOption, InterpConst<1>, 6, 12 /*, 27*/>, + SimulatorOption, InterpConst<2>, 6, 12>, + SimulatorOption, InterpConst<3>, 6, 12> +#endif +>{}; } diff --git a/tests/amr/data/field/refine/test_refine_field.py b/tests/amr/data/field/refine/test_refine_field.py index 188f80191..a3e66f756 100644 --- a/tests/amr/data/field/refine/test_refine_field.py +++ b/tests/amr/data/field/refine/test_refine_field.py @@ -102,6 +102,9 @@ def refine_electric(field, **kwargs): fine_data[::2, 1::2] = 0.5 * (coarse_data[:, 1:] + coarse_data[:, :-1]) fine_data[1::2, 1::2] = 0.5 * (coarse_data[:, 1:] + coarse_data[:, :-1]) + elif field.box.ndim == 3: + assert False # fix + return cropToFieldData(fine_data, field) diff --git a/tests/amr/data/particles/copy/test_particles_data_copy.cpp b/tests/amr/data/particles/copy/test_particles_data_copy.cpp new file mode 100644 index 000000000..83f5be411 --- /dev/null +++ b/tests/amr/data/particles/copy/test_particles_data_copy.cpp @@ -0,0 +1,90 @@ + +#include "core/def/phare_mpi.hpp" + + +#include "amr/data/particles/particles_data.hpp" +#include +#include + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + + +using testing::DoubleEq; +using testing::Eq; +using testing::Pointwise; + +using namespace PHARE::core; +using namespace PHARE::amr; + + +template +struct AParticlesDataND : public testing::Test +{ + static constexpr auto dim = dimType{}(); + + SAMRAI::tbox::Dimension dimension{dim}; + SAMRAI::hier::BlockId blockId{0}; + + SAMRAI::hier::Box sourceDomain{SAMRAI::hier::Index{dimension, 6}, + SAMRAI::hier::Index{dimension, 11}, blockId}; + + SAMRAI::hier::Box destDomain{SAMRAI::hier::Index{dimension, 1}, + SAMRAI::hier::Index{dimension, 5}, blockId}; + + SAMRAI::hier::IntVector ghost{SAMRAI::hier::IntVector::getOne(dimension)}; + + ParticlesData> destData{destDomain, ghost, "name"}; + ParticlesData> sourceData{sourceDomain, ghost, "name"}; + typename ParticleArray::Particle_t particle; + + AParticlesDataND() + { + particle.weight = 1.0; + particle.charge = 1.0; + particle.v = {1.0, 1.0, 1.0}; + } +}; + + + +using WithAllDim = testing::Types, DimConst<2>, DimConst<3>>; + +TYPED_TEST_SUITE(AParticlesDataND, WithAllDim); + + +TYPED_TEST(AParticlesDataND, copy_test) +{ + // particle is in the ghost of the source patchdata + static constexpr auto dim = TypeParam{}(); + + this->particle.iCell = ConstArray(5); + + this->sourceData.patchGhostParticles.push_back(this->particle); + this->destData.copy(this->sourceData); + + ASSERT_THAT(this->destData.domainParticles.size(), Eq(1)); + ASSERT_THAT(this->destData.patchGhostParticles.size(), Eq(0)); +} + + + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + + SAMRAI::tbox::SAMRAI_MPI::init(&argc, &argv); + SAMRAI::tbox::SAMRAIManager::initialize(); + SAMRAI::tbox::SAMRAIManager::startup(); + + int testResult = RUN_ALL_TESTS(); + + // Finalize + + SAMRAI::tbox::SAMRAIManager::shutdown(); + SAMRAI::tbox::SAMRAIManager::finalize(); + SAMRAI::tbox::SAMRAI_MPI::finalize(); + + return testResult; +} diff --git a/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt new file mode 100644 index 000000000..064cf8fb5 --- /dev/null +++ b/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt @@ -0,0 +1,51 @@ +CartesianGridGeometry +{ + domain_boxes = [ (0, 0, 0) , (8, 8, 8) ] + x_lo = 0.0, 0.0, 0.0 + x_up = 1.0, 1.0, 1.0 + periodic_dimension = 1, 1, 1 +} +Main +{ + dim = 1 +} +PatchHierarchy +{ + max_levels = 2 + proper_nesting_buffer = 2 + // vector of coarse ratio with dim dimension + ratio_to_coarser + { + level_1 = 2, 2, 2 + } + smallest_patch_size + { + level_0 = 5, 5, 5 + // All finer level will use same values in as level_0 + } +} +ChopAndPackLoadBalancer +{ + bin_pack_method = "SPATIAL" +} +StandardTagAndInitialize +{ + at_0 + { + cycle = 0 + tag_0 + { + tagging_method = "REFINE_BOXES" + level_0 + { + boxes = [ (1, 1, 1) , (3, 3, 3)] + } + } + } +} +TileClustering +{ +} +GriddingAlgorithm +{ +} diff --git a/tests/amr/data/particles/refine/test_split.cpp b/tests/amr/data/particles/refine/test_split.cpp index 8027c35a2..77651dd75 100644 --- a/tests/amr/data/particles/refine/test_split.cpp +++ b/tests/amr/data/particles/refine/test_split.cpp @@ -22,7 +22,7 @@ struct SplitterTest : public ::testing::Test SplitterTest() { Splitter splitter; } }; -using Splitters = testing::Types, Splitter<2, 1, 8> /*, Splitter<3, 1, 27>*/>; +using Splitters = testing::Types, Splitter<2, 1, 8>, Splitter<3, 1, 27>>; TYPED_TEST_SUITE(SplitterTest, Splitters); diff --git a/tests/amr/data/particles/stream_pack/test_no_ghost.cpp b/tests/amr/data/particles/stream_pack/test_no_ghost.cpp new file mode 100644 index 000000000..77718194e --- /dev/null +++ b/tests/amr/data/particles/stream_pack/test_no_ghost.cpp @@ -0,0 +1,297 @@ + +#include "core/def/phare_mpi.hpp" + + +#include +#include + +#include "amr/data/particles/particles_data.hpp" +#include "amr/data/particles/particles_data_factory.hpp" +#include +#include +#include +#include +#include +#include +#include + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + + +#include "core/utilities/types.hpp" + +using testing::Eq; + +using namespace PHARE::core; +using namespace PHARE::amr; + + +template +struct AParticlesData +{ + static constexpr auto dimension = dim; + + SAMRAI::tbox::Dimension amr_dimension{dim}; + SAMRAI::hier::BlockId blockId{0}; + + SAMRAI::hier::Box destDomain{SAMRAI::hier::Index{amr_dimension, 0}, + SAMRAI::hier::Index{amr_dimension, 5}, blockId}; + + SAMRAI::hier::Box sourceDomain{SAMRAI::hier::Index{amr_dimension, 10}, + SAMRAI::hier::Index{amr_dimension, 15}, blockId}; + + SAMRAI::hier::IntVector ghost{SAMRAI::hier::IntVector::getOne(amr_dimension)}; + + std::shared_ptr patchDescriptor{ + std::make_shared()}; + + SAMRAI::hier::Patch destPatch{destDomain, patchDescriptor}; + SAMRAI::hier::Patch sourcePatch{sourceDomain, patchDescriptor}; + + ParticlesData> destData{destDomain, ghost, "name"}; + ParticlesData> sourceData{sourceDomain, ghost, "name"}; + + std::shared_ptr destGeom{ + std::make_shared(destPatch.getBox(), ghost)}; + + std::shared_ptr sourceGeom{ + std::make_shared(sourcePatch.getBox(), ghost)}; + + + SAMRAI::hier::Box srcMask{sourceData.getGhostBox()}; + SAMRAI::hier::Box fillBox{destData.getGhostBox()}; + + bool overwriteInterior{true}; + + SAMRAI::hier::Index oneIndex{SAMRAI::hier::IntVector::getOne(amr_dimension)}; + + SAMRAI::hier::Transformation transformation{destDomain.lower() - sourceDomain.upper() + - oneIndex}; + + + std::shared_ptr cellOverlap{ + std::dynamic_pointer_cast(destGeom->calculateOverlap( + *sourceGeom, srcMask, fillBox, overwriteInterior, transformation))}; + + + typename ParticleArray::Particle_t particle; + + + AParticlesData() + { + particle.weight = 1.0; + particle.charge = 1.0; + particle.v = {1.0, 1.0, 1.0}; + } +}; + + +template +struct StreamPackTest : public ::testing::Test +{ +}; + +using ParticlesDatas = testing::Types, AParticlesData<2>, AParticlesData<3>>; +TYPED_TEST_SUITE(StreamPackTest, ParticlesDatas); + +TYPED_TEST(StreamPackTest, PreserveVelocityWhenPackStreamWithPeriodics) +{ + using ParticlesData = TypeParam; + constexpr auto dim = ParticlesData::dimension; + + ParticlesData param; + auto& particle = param.particle; + auto& sourceData = param.sourceData; + auto& cellOverlap = param.cellOverlap; + auto& destData = param.destData; + + particle.iCell = ConstArray(15); + sourceData.domainParticles.push_back(particle); + + SAMRAI::tbox::MessageStream particlesWriteStream; + + sourceData.packStream(particlesWriteStream, *cellOverlap); + + SAMRAI::tbox::MessageStream particlesReadStream{particlesWriteStream.getCurrentSize(), + SAMRAI::tbox::MessageStream::Read, + particlesWriteStream.getBufferStart()}; + + destData.unpackStream(particlesReadStream, *cellOverlap); + + ASSERT_THAT(destData.patchGhostParticles.size(), Eq(1)); + ASSERT_THAT(destData.patchGhostParticles[0].v, Eq(particle.v)); +} + + + + +TYPED_TEST(StreamPackTest, ShiftTheiCellWhenPackStreamWithPeriodics) +{ + using ParticlesData = TypeParam; + constexpr auto dim = ParticlesData::dimension; + + ParticlesData param; + auto& particle = param.particle; + auto& sourceData = param.sourceData; + auto& cellOverlap = param.cellOverlap; + auto& destData = param.destData; + + particle.iCell = ConstArray(15); + + sourceData.domainParticles.push_back(particle); + + SAMRAI::tbox::MessageStream particlesWriteStream; + + sourceData.packStream(particlesWriteStream, *cellOverlap); + + SAMRAI::tbox::MessageStream particlesReadStream{particlesWriteStream.getCurrentSize(), + SAMRAI::tbox::MessageStream::Read, + particlesWriteStream.getBufferStart()}; + + destData.unpackStream(particlesReadStream, *cellOverlap); + + // patch0 start at 0 , patch1 start at 10 + // with periodics condition, we have 0 equivalent to 15 + auto expectediCell = ConstArray(-1); + + + ASSERT_THAT(destData.patchGhostParticles.size(), Eq(1)); + ASSERT_THAT(destData.patchGhostParticles[0].iCell, Eq(expectediCell)); +} + + + +TYPED_TEST(StreamPackTest, PackInTheCorrectBufferWithPeriodics) +{ + using ParticlesData = TypeParam; + constexpr auto dim = ParticlesData::dimension; + + ParticlesData param; + auto& particle = param.particle; + auto& sourceData = param.sourceData; + auto& cellOverlap = param.cellOverlap; + auto& destData = param.destData; + + particle.iCell = ConstArray(15); + + sourceData.domainParticles.push_back(particle); + + SAMRAI::tbox::MessageStream particlesWriteStream; + + sourceData.packStream(particlesWriteStream, *cellOverlap); + + SAMRAI::tbox::MessageStream particlesReadStream{particlesWriteStream.getCurrentSize(), + SAMRAI::tbox::MessageStream::Read, + particlesWriteStream.getBufferStart()}; + + destData.unpackStream(particlesReadStream, *cellOverlap); + + auto expectediCell = ConstArray(-1); + + ASSERT_THAT(destData.patchGhostParticles.size(), Eq(1)); + ASSERT_THAT(destData.patchGhostParticles[0].iCell, Eq(expectediCell)); +} + + + +TYPED_TEST(StreamPackTest, + PreserveParticleAttributesWhenPackingWithPeriodicsFromGhostSrcToDomainDest) +{ + using ParticlesData = TypeParam; + constexpr auto dim = ParticlesData::dimension; + + ParticlesData param; + auto& particle = param.particle; + auto& sourceData = param.sourceData; + auto& cellOverlap = param.cellOverlap; + auto& destData = param.destData; + + particle.iCell = ConstArray(16); + + sourceData.domainParticles.push_back(particle); + + SAMRAI::tbox::MessageStream particlesWriteStream; + + sourceData.packStream(particlesWriteStream, *cellOverlap); + + SAMRAI::tbox::MessageStream particlesReadStream{particlesWriteStream.getCurrentSize(), + SAMRAI::tbox::MessageStream::Read, + particlesWriteStream.getBufferStart()}; + + destData.unpackStream(particlesReadStream, *cellOverlap); + + auto expectediCell = ConstArray(0); + + EXPECT_THAT(destData.domainParticles[0].v, Eq(particle.v)); + EXPECT_THAT(destData.domainParticles[0].iCell, Eq(expectediCell)); + EXPECT_THAT(destData.domainParticles[0].delta, Eq(particle.delta)); + EXPECT_THAT(destData.domainParticles[0].weight, Eq(particle.weight)); + EXPECT_THAT(destData.domainParticles[0].charge, Eq(particle.charge)); +} + + + +TYPED_TEST(StreamPackTest, + PreserveParticleAttributesWhenPackingWithPeriodicsFromDomainSrcToGhostDest) +{ + using ParticlesData = TypeParam; + constexpr auto dim = ParticlesData::dimension; + + ParticlesData param; + auto& particle = param.particle; + auto& sourceData = param.sourceData; + auto& cellOverlap = param.cellOverlap; + auto& destData = param.destData; + + particle.iCell = ConstArray(15); + + sourceData.domainParticles.push_back(particle); + + SAMRAI::tbox::MessageStream particlesWriteStream; + + sourceData.packStream(particlesWriteStream, *cellOverlap); + + SAMRAI::tbox::MessageStream particlesReadStream{particlesWriteStream.getCurrentSize(), + SAMRAI::tbox::MessageStream::Read, + particlesWriteStream.getBufferStart()}; + + destData.unpackStream(particlesReadStream, *cellOverlap); + + auto expectediCell = ConstArray(-1); + + EXPECT_THAT(destData.patchGhostParticles[0].v, Eq(particle.v)); + EXPECT_THAT(destData.patchGhostParticles[0].iCell, Eq(expectediCell)); + EXPECT_THAT(destData.patchGhostParticles[0].delta, Eq(particle.delta)); + EXPECT_THAT(destData.patchGhostParticles[0].weight, Eq(particle.weight)); + EXPECT_THAT(destData.patchGhostParticles[0].charge, Eq(particle.charge)); +} + + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + + + + SAMRAI::tbox::SAMRAI_MPI::init(&argc, &argv); + + SAMRAI::tbox::SAMRAIManager::initialize(); + + SAMRAI::tbox::SAMRAIManager::startup(); + + + int testResult = RUN_ALL_TESTS(); + + + // Finalize + + SAMRAI::tbox::SAMRAIManager::shutdown(); + + SAMRAI::tbox::SAMRAIManager::finalize(); + + SAMRAI::tbox::SAMRAI_MPI::finalize(); + + return testResult; +} diff --git a/tests/core/data/ndarray/test_main.cpp b/tests/core/data/ndarray/test_main.cpp index be0bebd9f..b034d1b08 100644 --- a/tests/core/data/ndarray/test_main.cpp +++ b/tests/core/data/ndarray/test_main.cpp @@ -392,6 +392,60 @@ TEST(MaskedView2d, maskOps2) + Mask{0u}.nCells(array)); } +TEST(MaskedView3d, maskOps3) +{ + constexpr std::size_t dim = 3; + constexpr std::uint32_t size0 = 10; + constexpr std::uint32_t sizeCu = size0 * size0 * size0; + using Mask = PHARE::core::NdArrayMask; + + auto sum = [](auto const& array) { return std::accumulate(array.begin(), array.end(), 0); }; + + { + NdArrayVector array{size0, size0, size0}; + EXPECT_EQ(sum(array), 0); + std::fill(array.begin(), array.end(), 1); + EXPECT_EQ(sum(array), sizeCu); + } + + { + NdArrayVector array{size0, size0, size0}; + EXPECT_EQ(std::accumulate(array.begin(), array.end(), 0), 0); + array[Mask{0}] = 1; + EXPECT_EQ(std::accumulate(array.begin(), array.end(), 0), 488); + + // outter cells of a 10**3 cube = + // (10 * 10 * 2) + (10 * 8 * 2) + (8 * 8 * 2); + // or + // (8 * 8 * 6) + (10 * 4) + (8 * 8); + // = 488 + } + + std::uint32_t ten = 10; + PHARE::core::NdArrayVector<3> array(ten, ten, ten); + + array[Mask{0}] = 1; + EXPECT_EQ(sum(array), 488); + array[Mask{1}] >> array[Mask{0}]; + EXPECT_EQ(sum(array), 0); + + array[Mask{2}] = 1; + EXPECT_EQ(sum(array), 152); + array[Mask{1}] = 1; + EXPECT_EQ(sum(array), 448); + array[Mask{1}] = 0; + EXPECT_EQ(sum(array), 152); + + array[Mask{2}] >> array[Mask{1}]; + EXPECT_EQ(sum(array), 448); + array[Mask{2}] = 0; + EXPECT_EQ(sum(array), 296); + + EXPECT_EQ(Mask{1}.nCells(array), 296); + EXPECT_EQ(Mask{2}.nCells(array), 152); +} + + int main(int argc, char** argv) { ::testing::InitGoogleTest(&argc, argv); diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp index 28e2fff14..b69025e7a 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -721,6 +721,73 @@ using My2dTypes = ::testing::Types, Interpolator<2, 2>, Inter INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_2d, My2dTypes); + +/*********************************************************************************************/ +template +struct ACollectionOfParticles_3d : public ::testing::Test +{ + static constexpr auto interp_order = Interpolator::interp_order; + static constexpr std::size_t dim = 3; + static constexpr std::uint32_t nx = 15, ny = 15, nz = 15; + static constexpr int start = 0, end = 5; + constexpr static auto safeLayer = static_cast(1 + ghostWidthForParticles()); + + using PHARE_TYPES = PHARE::core::PHARE_Types; + using ParticleArray_t = typename PHARE_TYPES::ParticleArray_t; + using GridLayout_t = typename PHARE_TYPES::GridLayout_t; + using Grid_t = typename PHARE_TYPES::Grid_t; + using UsableVecFieldND = UsableVecField; + + GridLayout_t layout{ConstArray(.1), {nx, ny}, ConstArray(0)}; + ParticleArray_t particles; + Grid_t rho; + UsableVecFieldND v; + Interpolator interpolator; + + ACollectionOfParticles_3d() + : particles{grow(layout.AMRBox(), safeLayer)} + , rho{"field", HybridQuantity::Scalar::rho, nx, ny, nz} + , v{"v", layout, HybridQuantity::Vector::V} + { + double weight = [](auto const& meshSize) { + return std::accumulate(meshSize.begin(), meshSize.end(), 1.0, + std::multiplies()); + }(layout.meshSize()); + + for (int i = start; i < end; i++) + for (int j = start; j < end; j++) + for (int k = start; k < end; k++) + { + auto& part = particles.emplace_back(); + part.iCell = {i, j, k}; + part.delta = ConstArray(.5); + part.weight = weight; + part.v[0] = +2.; + part.v[1] = -1.; + part.v[2] = +1.; + } + + interpolator(particles, rho, v, layout); + } +}; +TYPED_TEST_SUITE_P(ACollectionOfParticles_3d); + + +TYPED_TEST_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d) +{ + // auto const& [vx, vy, vz] = this->v(); + // EXPECT_DOUBLE_EQ(this->rho(7, 7, 7), 1.0); + // EXPECT_DOUBLE_EQ(vx(7, 7, 7), 2.0); + // EXPECT_DOUBLE_EQ(vy(7, 7, 7), -1.0); + // EXPECT_DOUBLE_EQ(vz(7, 7, 7), 1.0); +} +REGISTER_TYPED_TEST_SUITE_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d); + + +using My3dTypes = ::testing::Types, Interpolator<3, 2>, Interpolator<3, 3>>; +INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_3d, My3dTypes); +/*********************************************************************************************/ + int main(int argc, char** argv) { ::testing::InitGoogleTest(&argc, argv); diff --git a/tests/diagnostic/CMakeLists.txt b/tests/diagnostic/CMakeLists.txt index 1554c1020..8ce79e4be 100644 --- a/tests/diagnostic/CMakeLists.txt +++ b/tests/diagnostic/CMakeLists.txt @@ -34,9 +34,11 @@ if(HighFive) _add_diagnostics_test(test-diagnostics_1d) _add_diagnostics_test(test-diagnostics_2d) + _add_diagnostics_test(test-diagnostics_3d) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_1d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_1d.py @ONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_2d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_2d.py @ONLY) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_3d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_3d.py @ONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/__init__.py ${CMAKE_CURRENT_BINARY_DIR}/__init__.py @ONLY) message(STATUS "diagnostic working directory " ${PHARE_PROJECT_DIR}) diff --git a/tests/diagnostic/job_3d.py.in b/tests/diagnostic/job_3d.py.in new file mode 100644 index 000000000..27e45c041 --- /dev/null +++ b/tests/diagnostic/job_3d.py.in @@ -0,0 +1,14 @@ +#!/usr/bin/env python3 + +import pyphare.pharein as ph +from pyphare.pharein import ElectronModel +from tests.simulator import basicSimulatorArgs, makeBasicModel +from tests.diagnostic import dump_all_diags + +out = "phare_outputs/diags_3d/" +simInput = {"diag_options": {"format": "phareh5", "options": {"dir": out, "mode" : "overwrite"}}} + +ph.Simulation(**basicSimulatorArgs(dim = 3, interp = 1, **simInput)) +model = makeBasicModel() +ElectronModel(closure="isothermal",Te = 0.12) +dump_all_diags(model.populations) diff --git a/tests/diagnostic/test-diagnostics_3d.cpp b/tests/diagnostic/test-diagnostics_3d.cpp new file mode 100644 index 000000000..54e5741eb --- /dev/null +++ b/tests/diagnostic/test-diagnostics_3d.cpp @@ -0,0 +1,35 @@ + +#include "core/def/phare_mpi.hpp" + +#include "test_diagnostics.ipp" + +static std::string const job_file = "job_3d"; +static std::string const out_dir = "phare_outputs/diags_3d/"; + +TYPED_TEST(Simulator3dTest, fluid) +{ + fluid_test(TypeParam{job_file}, out_dir); +} + +TYPED_TEST(Simulator3dTest, particles) +{ + particles_test(TypeParam{job_file}, out_dir); +} + +TYPED_TEST(Simulator3dTest, electromag) +{ + electromag_test(TypeParam{job_file}, out_dir); +} + +TYPED_TEST(Simulator3dTest, allFromPython) +{ + allFromPython_test(TypeParam{job_file}, out_dir); +} + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + PHARE::SamraiLifeCycle samsam(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/functional/harris/harris_3d.py b/tests/functional/harris/harris_3d.py new file mode 100644 index 000000000..ed6cecf9b --- /dev/null +++ b/tests/functional/harris/harris_3d.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python3 +import os +import numpy as np +from pathlib import Path + +import pyphare.pharein as ph +from pyphare.cpp import cpp_lib +from pyphare.pharesee.run import Run +from pyphare.simulator.simulator import Simulator +from pyphare.simulator.simulator import startMPI + +os.environ["PHARE_SCOPE_TIMING"] = "0" # turn on scope timing + + +ph.NO_GUI() +cpp = cpp_lib() +startMPI() + +cells = (40, 40, 40) +dl = (0.2, 0.2, 0.2) + +diag_outputs = "phare_outputs/test/harris/3d" +time_step_nbr = 1000 +time_step = 0.001 +final_time = time_step * time_step_nbr + +timestamps = [0, final_time] +if time_step_nbr > 100: + dt = 100 * time_step + nt = final_time / dt + 1 + timestamps = dt * np.arange(nt) + +plot_dir = Path(f"{diag_outputs}_plots") +plot_dir.mkdir(parents=True, exist_ok=True) + + +def config(): + sim = ph.Simulation( + time_step=time_step, + time_step_nbr=time_step_nbr, + dl=dl, + cells=cells, + refinement="tagging", + max_nbr_levels=1, + hyper_resistivity=0.001, + resistivity=0.001, + diag_options={ + "format": "phareh5", + "options": {"dir": diag_outputs, "mode": "overwrite"}, + }, + # strict=True, + ) + + def density(x, y, z): + L = sim.simulation_domain()[1] + return ( + 0.2 + + 1.0 / np.cosh((y - L * 0.3) / 0.5) ** 2 + + 1.0 / np.cosh((y - L * 0.7) / 0.5) ** 2 + ) + + def S(y, y0, l): + return 0.5 * (1.0 + np.tanh((y - y0) / l)) + + def by(x, y, z): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + w1 = 0.2 + w2 = 1.0 + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) + w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) + w5 = 2.0 * w1 / w2 + return (w5 * x0 * w3) + (-w5 * x0 * w4) + + def bx(x, y, z): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + w1 = 0.2 + w2 = 1.0 + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) + w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) + w5 = 2.0 * w1 / w2 + v1 = -1 + v2 = 1.0 + return ( + v1 + + (v2 - v1) * (S(y, Ly * 0.3, 0.5) - S(y, Ly * 0.7, 0.5)) + + (-w5 * y1 * w3) + + (+w5 * y2 * w4) + ) + + def bz(x, y, z): + return bx(x, y, z) * -1 + + def b2(x, y, z): + return bx(x, y, z) ** 2 + by(x, y, z) ** 2 + bz(x, y, z) ** 2 + + def T(x, y, z): + K = 1 + temp = 1.0 / density(x, y, z) * (K - b2(x, y, z) * 0.5) + assert np.all(temp > 0) + return temp + + def vxyz(x, y, z): + return 0.0 + + def vthxyz(x, y, z): + return np.sqrt(T(x, y, z)) + + C = "xyz" + vvv = { + **{f"vbulk{c}": vxyz for c in C}, + **{f"vth{c}": vthxyz for c in C}, + "nbr_part_per_cell": 100, + } + protons = {"charge": 1, "density": density, **vvv, "init": {"seed": 12334}} + ph.MaxwellianFluidModel(bx=bx, by=by, bz=bz, protons=protons) + ph.ElectronModel(closure="isothermal", Te=0.0) + + for quantity in ["density", "bulkVelocity"]: + ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) + ph.FluidDiagnostics( + quantity="density", write_timestamps=timestamps, population_name="protons" + ) + for quantity in ["E", "B"]: + ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) + ph.InfoDiagnostics(quantity="particle_count") # defaults all coarse time steps + + return sim + + +def plot_file_for_qty(qty, time): + return f"{plot_dir}/harris_{qty}_t{time}.png" + + +def plot(diag_dir): + run = Run(diag_dir) + for time in timestamps: + # run.GetDivB(time).plot( + # filename=plot_file_for_qty("divb", time), + # plot_patches=True, + # vmin=1e-11, + # vmax=2e-10, + # ) + # run.GetRanks(time).plot( + # filename=plot_file_for_qty("Ranks", time), + # plot_patches=True, + # ) + run.GetN(time, pop_name="protons").plot( + filename=plot_file_for_qty("N", time), + plot_patches=True, + ) + for c in ["x", "y", "z"]: + run.GetB(time).plot( + filename=plot_file_for_qty(f"b{c}", time), + qty=f"B{c}", + plot_patches=True, + ) + # run.GetJ(time).plot( + # filename=plot_file_for_qty("jz", time), + # qty="Jz", + # plot_patches=True, + # vmin=-2, + # vmax=2, + # ) + + +def main(): + Simulator(config()).run() + if cpp.mpi_rank() == 0: + plot(diag_outputs) + # try: + # from tools.python3 import plotting as m_plotting + + # m_plotting.plot_run_timer_data(diag_outputs, cpp.mpi_rank()) + # except ImportError: + # print("Phlop not found - install with: `pip install phlop`") + cpp.mpi_barrier() + + +if __name__ == "__main__": + main() diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index 5c100634c..74a65f1d4 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -94,12 +94,17 @@ def density_2d_periodic(sim, x, y): ) -# def density_3d_periodic(sim, x, y, z): -# xmax, ymax, zmax = sim.simulation_domain() -# background_particles = 0.3 # avoids 0 density -# xx, yy, zz = meshify(x, y, z) -# r = np.exp(-(xx-0.5*xmax)**2)*np.exp(-(yy-ymax/2.)**2)*np.exp(-(zz-zmax/2.)**2) + background_particles -# return r +def density_3d_periodic(sim, x, y, z): + xmax, ymax, zmax = sim.simulation_domain() + background_particles = 0.3 # avoids 0 density + xx, yy, zz = meshify(x, y, z) + r = ( + np.exp(-((xx - 0.5 * xmax) ** 2)) + * np.exp(-((yy - ymax / 2.0) ** 2)) + * np.exp(-((zz - zmax / 2.0) ** 2)) + + background_particles + ) + return r def defaultPopulationSettings(sim, density_fn, vbulk_fn): diff --git a/tests/simulator/advance/CMakeLists.txt b/tests/simulator/advance/CMakeLists.txt index 0b4ad77a5..4304d8aa0 100644 --- a/tests/simulator/advance/CMakeLists.txt +++ b/tests/simulator/advance/CMakeLists.txt @@ -18,6 +18,11 @@ if(HighFive) phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-2d-fields test_fields_advance_2d.py ${CMAKE_CURRENT_BINARY_DIR}) phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-2d-particles test_particles_advance_2d.py ${CMAKE_CURRENT_BINARY_DIR}) + if(highResourceTests) # off by default as it's quite intensive + phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-3d-fields test_fields_advance_3d.py ${CMAKE_CURRENT_BINARY_DIR}) + phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-3d-particles test_particles_advance_3d.py ${CMAKE_CURRENT_BINARY_DIR}) + endif(highResourceTests) + endif() endif() diff --git a/tests/simulator/advance/test_fields_advance_3d.py b/tests/simulator/advance/test_fields_advance_3d.py new file mode 100644 index 000000000..272ff40a8 --- /dev/null +++ b/tests/simulator/advance/test_fields_advance_3d.py @@ -0,0 +1,107 @@ +""" + This file exists independently from test_advance.py to isolate dimension + test cases and allow each to be overridden in some way if required. +""" + +import unittest + +import matplotlib +from ddt import data, ddt, unpack +from pyphare.core.box import Box3D + +from tests.simulator.test_advance import AdvanceTestBase + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc, cells = 10, 30 + + +def per_interp(dic): + return [(interp, dic) for interp in interp_orders] + + +@ddt +class AdvanceTest(AdvanceTestBase): + @data( + *per_interp({}), + *per_interp({"L0": [Box3D(10, 19)]}), + *per_interp({"L0": [Box3D(8, 20)]}), + ) + @unpack + def test_overlaped_fields_are_equal(self, interp_order, refinement_boxes): + print(f"{self._testMethodName}_{ndim}d") + time_step_nbr = 3 + time_step = 0.001 + datahier = self.getHierarchy( + ndim, + interp_order, + refinement_boxes, + "eb", + cells=cells, + time_step=time_step, + time_step_nbr=time_step_nbr, + nbr_part_per_cell=ppc, + ) + self._test_overlaped_fields_are_equal(datahier, time_step_nbr, time_step) + + @data( + *per_interp({}), + *per_interp({"L0": [Box3D(5, 14)]}), + ) + @unpack + def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts( + self, interp_order, refinement_boxes + ): + print(f"{self._testMethodName}_{ndim}d") + time_step_nbr = 3 + time_step = 0.001 + from pyphare.pharein.simulation import check_patch_size + + largest_patch_size, smallest_patch_size = check_patch_size( + ndim, interp_order=interp_order, cells=[cells] * ndim + ) + datahier = self.getHierarchy( + ndim, + interp_order, + refinement_boxes, + "eb", + cells=cells, + smallest_patch_size=smallest_patch_size, + largest_patch_size=smallest_patch_size, + time_step=time_step, + time_step_nbr=time_step_nbr, + nbr_part_per_cell=ppc, + ) + self._test_overlaped_fields_are_equal(datahier, time_step_nbr, time_step) + + # @data( + # *per_interp(({"L0": {"B0": Box3D(10, 14)}})), + # *per_interp(({"L0": {"B0": Box3D(10, 14), "B1": Box3D(15, 19)}})), + # *per_interp(({"L0": {"B0": Box3D(6, 23)}})), + # *per_interp(({"L0": {"B0": Box3D( 2, 12), "B1": Box3D(13, 25)}})), + # *per_interp(({"L0": {"B0": Box3D( 5, 20)}, "L1": {"B0": Box3D(15, 19)}})), + # *per_interp(({"L0": {"B0": Box3D( 5, 20)}, "L1": {"B0": Box3D(12, 38)}, "L2": {"B0": Box3D(30, 52)} })), + # ) + # @unpack + # def test_field_coarsening_via_subcycles(self, interp_order, refinement_boxes): + # print(f"{self._testMethodName}_{ndim}d") + # self._test_field_coarsening_via_subcycles(ndim, interp_order, refinement_boxes, dl=.3, cells=cells) + + # @unittest.skip("should change to work on moments") + # @data( # only supports a hierarchy with 2 levels + # *per_interp(({"L0": [Box3D(0, 4)]})), + # *per_interp(({"L0": [Box3D(10, 14)]})), + # *per_interp(({"L0": [Box3D(0, 4), Box3D(10, 14)]})), + # *per_interp(({"L0": [Box3D(0, 4), Box3D(5, 9), Box3D(10, 14)]})), + # *per_interp(({"L0": [Box3D(20, 24)]})), + # ) + # @unpack + # def test_field_level_ghosts_via_subcycles_and_coarser_interpolation(self, interp_order, refinement_boxes): + # print(f"{self._testMethodName}_{ndim}d") + # self._test_field_level_ghosts_via_subcycles_and_coarser_interpolation(ndim, interp_order, refinement_boxes) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulator/advance/test_particles_advance_3d.py b/tests/simulator/advance/test_particles_advance_3d.py new file mode 100644 index 000000000..a095986d9 --- /dev/null +++ b/tests/simulator/advance/test_particles_advance_3d.py @@ -0,0 +1,51 @@ +""" + This file exists independently from test_advance.py to isolate dimension + test cases and allow each to be overridden in some way if required. +""" + +import unittest + +import matplotlib +from ddt import data, ddt, unpack +from pyphare.core.box import Box3D + +from tests.simulator.test_advance import AdvanceTestBase + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc = 5 + + +def per_interp(dic): + return [(interp, dic) for interp in interp_orders] + + +@ddt +class AdvanceTest(AdvanceTestBase): + @data( + *per_interp({}), + *per_interp({"L0": [Box3D(10, 19)]}), + *per_interp({"L0": [Box3D(5, 9), Box3D(10, 14)]}), + ) + @unpack + def test_overlapped_particledatas_have_identical_particles( + self, interp_order, refinement_boxes + ): + self._test_overlapped_particledatas_have_identical_particles( + ndim, + interp_order, + refinement_boxes, + ppc=ppc, + cells=40, + largest_patch_size=20, + ) + + @data(*interp_orders) + def test_L0_particle_number_conservation(self, interp): + self._test_L0_particle_number_conservation(ndim, interp, ppc=ppc, cells=30) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulator/initialize/CMakeLists.txt b/tests/simulator/initialize/CMakeLists.txt index cbaf44d02..1f1c93827 100644 --- a/tests/simulator/initialize/CMakeLists.txt +++ b/tests/simulator/initialize/CMakeLists.txt @@ -18,6 +18,11 @@ if(HighFive) phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-2d-fields test_fields_init_2d.py ${CMAKE_CURRENT_BINARY_DIR}) phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-2d-particles test_particles_init_2d.py ${CMAKE_CURRENT_BINARY_DIR}) + if(highResourceTests) + phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-3d-fields test_fields_init_3d.py ${CMAKE_CURRENT_BINARY_DIR}) + phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-3d-particles test_particles_init_3d.py ${CMAKE_CURRENT_BINARY_DIR}) + endif() + endif() endif() diff --git a/tests/simulator/initialize/test_fields_init_2d.py b/tests/simulator/initialize/test_fields_init_2d.py index af47a5d99..22e4587d1 100644 --- a/tests/simulator/initialize/test_fields_init_2d.py +++ b/tests/simulator/initialize/test_fields_init_2d.py @@ -23,7 +23,7 @@ class Initialization2DTest(InitializationTest): @data(*interp_orders) def test_B_is_as_provided_by_user(self, interp_order): print(f"\n{self._testMethodName}_{ndim}d") - self._test_B_is_as_provided_by_user(ndim, interp_order, nbr_part_per_cell=ppc) + self._test_B_is_as_provided_by_user(ndim, interp_order, ppc=ppc) @data(*interp_orders) def test_bulkvel_is_as_provided_by_user(self, interp_order): diff --git a/tests/simulator/initialize/test_fields_init_3d.py b/tests/simulator/initialize/test_fields_init_3d.py new file mode 100644 index 000000000..05d8f5a53 --- /dev/null +++ b/tests/simulator/initialize/test_fields_init_3d.py @@ -0,0 +1,49 @@ +""" + This file exists independently from test_initialization.py to isolate dimension + test cases and allow each to be overridden in some way if required. +""" + +import unittest + +import numpy as np +import matplotlib +from ddt import data, ddt + +from tests.simulator.test_initialization import InitializationTest + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc, cells = 10, 20 + + +@ddt +class Initialization3DTest(InitializationTest): + @data(*interp_orders) + def test_B_is_as_provided_by_user(self, interp_order): + print(f"\n{self._testMethodName}_{ndim}d") + self._test_B_is_as_provided_by_user(ndim, interp_order, ppc=ppc, cells=cells) + + @data(*interp_orders) + def test_bulkvel_is_as_provided_by_user(self, interp_order): + print(f"\n{self._testMethodName}_{ndim}d") + self._test_bulkvel_is_as_provided_by_user( + ndim, interp_order, ppc=ppc, cells=cells + ) + + @data(*interp_orders) + def test_density_is_as_provided_by_user(self, interp_order): + print(f"\n{self._testMethodName}_{ndim}d") + self._test_density_is_as_provided_by_user(ndim, interp_order, cells=cells) + + @data(*interp_orders) # uses too much RAM - to isolate somehow + def test_density_decreases_as_1overSqrtN(self, interp_order): + print(f"\n{self._testMethodName}_{ndim}d") + self._test_density_decreases_as_1overSqrtN( + ndim, interp_order, np.asarray([10, 25, 50, 75]), cells=cells + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulator/initialize/test_particles_init_1d.py b/tests/simulator/initialize/test_particles_init_1d.py index 60c44d692..708995b33 100644 --- a/tests/simulator/initialize/test_particles_init_1d.py +++ b/tests/simulator/initialize/test_particles_init_1d.py @@ -8,12 +8,9 @@ import matplotlib from ddt import data, ddt, unpack from pyphare.core.box import Box1D -from pyphare.cpp import cpp_lib from tests.simulator.test_initialization import InitializationTest -cpp = cpp_lib() - matplotlib.use("Agg") # for systems without GUI ndim = 1 diff --git a/tests/simulator/initialize/test_particles_init_2d.py b/tests/simulator/initialize/test_particles_init_2d.py index cc56392f8..758d3b0a2 100644 --- a/tests/simulator/initialize/test_particles_init_2d.py +++ b/tests/simulator/initialize/test_particles_init_2d.py @@ -23,7 +23,7 @@ def per_interp(dic): @ddt -class Initialization1DTest(InitializationTest): +class Initialization2DTest(InitializationTest): @data(*interp_orders) def test_nbr_particles_per_cell_is_as_provided(self, interp_order): print(f"{self._testMethodName}_{ndim}d") diff --git a/tests/simulator/initialize/test_particles_init_3d.py b/tests/simulator/initialize/test_particles_init_3d.py new file mode 100644 index 000000000..a87c5328b --- /dev/null +++ b/tests/simulator/initialize/test_particles_init_3d.py @@ -0,0 +1,98 @@ +""" + This file exists independently from test_initialization.py to isolate dimension + test cases and allow each to be overridden in some way if required. +""" + +import unittest + +import matplotlib +from ddt import data, ddt, unpack +from pyphare.core.box import Box3D + +from tests.simulator.test_initialization import InitializationTest + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc, cells = 10, 30 + + +def per_interp(dic): + return [(interp, dic) for interp in interp_orders] + + +@ddt +class Initialization3DTest(InitializationTest): + @data(*interp_orders) + def test_nbr_particles_per_cell_is_as_provided(self, interp_order): + print(f"{self._testMethodName}_{ndim}d") + self._test_nbr_particles_per_cell_is_as_provided( + ndim, interp_order, ppc, cells=cells + ) + + @data( + *per_interp(({"L0": {"B0": Box3D(10, 14)}})), + *per_interp(({"L0": {"B0": Box3D(10, 14)}, "L1": {"B0": Box3D(22, 26)}})), + *per_interp(({"L0": {"B0": Box3D(2, 6), "B1": Box3D(7, 11)}})), + ) + @unpack + def test_levelghostparticles_have_correct_split_from_coarser_particle( + self, interp_order, refinement_boxes + ): + print(f"\n{self._testMethodName}_{ndim}d") + now = self.datetime_now() + self._test_levelghostparticles_have_correct_split_from_coarser_particle( + self.getHierarchy( + ndim, + interp_order, + refinement_boxes, + "particles", + cells=cells, + nbr_part_per_cell=ppc, + ) + ) + print( + f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds" + ) + + @data( + *per_interp(({"L0": {"B0": Box3D(10, 14)}})), + *per_interp(({"L0": {"B0": Box3D(5, 14)}, "L1": {"B0": Box3D(15, 19)}})), + *per_interp(({"L0": {"B0": Box3D(2, 12), "B1": Box3D(13, 25)}})), + ) + @unpack + def test_domainparticles_have_correct_split_from_coarser_particle( + self, interp_order, refinement_boxes + ): + print(f"\n{self._testMethodName}_{ndim}d") + now = self.datetime_now() + self._test_domainparticles_have_correct_split_from_coarser_particle( + ndim, interp_order, refinement_boxes, nbr_part_per_cell=ppc + ) + print( + f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds" + ) + + # @data({"cells": 40, "smallest_patch_size": 20, "largest_patch_size": 20, "nbr_part_per_cell" : ppc}) + # def test_no_patch_ghost_on_refined_level_case(self, simInput): + # print(f"\n{self._testMethodName}_{ndim}d") + # now = self.datetime_now() + # self._test_patch_ghost_on_refined_level_case(ndim, False, **simInput) + # print(f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds") + + # @data({"cells": 40, "interp_order": 1, "nbr_part_per_cell" : ppc}) + # def test_has_patch_ghost_on_refined_level_case(self, simInput): + # print(f"\n{self._testMethodName}_{ndim}d") + # from pyphare.pharein.simulation import check_patch_size + # diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts_{ndim}_{self.ddt_test_id()}" + # _, smallest_patch_size = check_patch_size(ndim, **simInput) + # simInput["smallest_patch_size"] = smallest_patch_size + # simInput["largest_patch_size"] = smallest_patch_size + # now = self.datetime_now() + # self._test_patch_ghost_on_refined_level_case(ndim, True, **simInput) + # print(f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds") + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulator/per_test.hpp b/tests/simulator/per_test.hpp index d4a8be4e0..901c7c428 100644 --- a/tests/simulator/per_test.hpp +++ b/tests/simulator/per_test.hpp @@ -108,6 +108,14 @@ using Simulators2d = testing::Types< TYPED_TEST_SUITE(Simulator2dTest, Simulators2d); +template +struct Simulator3dTest : public ::testing::Test +{ +}; +using Simulator3d = testing::Types, SimulatorTestParam<3, 2, 6>, + SimulatorTestParam<3, 3, 6>>; +TYPED_TEST_SUITE(Simulator3dTest, Simulator3d); + #endif /* PHARE_TEST_SIMULATOR_PER_TEST_H */ diff --git a/tests/simulator/refined_particle_nbr.py b/tests/simulator/refined_particle_nbr.py index 95e21dda9..c111738f0 100644 --- a/tests/simulator/refined_particle_nbr.py +++ b/tests/simulator/refined_particle_nbr.py @@ -30,27 +30,21 @@ def __init__(self, *args, **kwargs): print(exc) sys.exit(1) + # needed in case exception is raised in test and Simulator not reset properly + def tearDown(self): + if self.simulator is not None: + self.simulator.reset() + # while splitting particles may leave the domain area # so we remove the particles from the border cells of each patch def _less_per_dim(self, dim, refined_particle_nbr, patch): if dim == 1: return refined_particle_nbr * 2 cellNbr = patch.upper - patch.lower + 1 + dim2 = refined_particle_nbr * ((cellNbr[0] * 2 + (cellNbr[1] * 2))) if dim == 2: - return refined_particle_nbr * ((cellNbr[0] * 2 + (cellNbr[1] * 2))) - raise ValueError("Unhandled dimension for function") - - def _check_deltas_and_weights(self, dim, interp, refined_particle_nbr): - yaml_dim = self.yaml_root["dimension_" + str(dim)] - yaml_interp = yaml_dim["interp_" + str(interp)] - yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)] - - yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")] - yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")] - - splitter_t = splitter_type(dim, interp, refined_particle_nbr) - np.testing.assert_allclose(yaml_delta, splitter_t.delta) - np.testing.assert_allclose(yaml_weight, splitter_t.weight) + return dim2 + return dim2 * (cellNbr[2] * 2) def _do_dim(self, dim, min_diff, max_diff): from pyphare.pharein.simulation import valid_refined_particle_nbr @@ -58,11 +52,7 @@ def _do_dim(self, dim, min_diff, max_diff): for interp in range(1, 4): prev_split_particle_max = 0 for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]: - self._check_deltas_and_weights(dim, interp, refined_particle_nbr) - - simInput = NoOverwriteDict( - {"refined_particle_nbr": refined_particle_nbr} - ) + simInput = {"refined_particle_nbr": refined_particle_nbr} self.simulator = Simulator(populate_simulation(dim, interp, **simInput)) self.simulator.initialize() dw = self.simulator.data_wrangler() @@ -78,9 +68,6 @@ def _do_dim(self, dim, min_diff, max_diff): per_pop += patch.data.size() max_per_pop = max(max_per_pop, per_pop) prev_min_diff = prev_split_particle_max * min_diff - - # while splitting particles may leave the domain area - # so we remove the particles from the border cells of each patch self.assertTrue(max_per_pop > prev_min_diff - (leaving_particles)) if prev_split_particle_max > 0: prev_max_diff = prev_min_diff * dim * max_diff @@ -88,45 +75,70 @@ def _do_dim(self, dim, min_diff, max_diff): prev_split_particle_max = max_per_pop self.simulator = None - """ 1d - refine 10 cells in 1d, ppc 100 - 10 * 2 * ppc = 200 - 10 * 3 * ppc = 300 300 / 200 = 1.5 - 10 * 4 * ppc = 400 500 / 400 = 1.33 - 10 * 5 * ppc = 500 500 / 400 = 1.25 - taking the minimul diff across permutations - current to previous should be at least this times bigger - """ - PREVIOUS_ITERATION_MIN_DIFF_1d = 1.25 - PREVIOUS_ITERATION_MAX_DIFF_1d = 1.50 + # 1d + # refine 10 cells in 1d, ppc 100 + # 10 * 2 * ppc = 200 + # 10 * 3 * ppc = 300 300 / 200 = 1.5 + # 10 * 4 * ppc = 400 500 / 400 = 1.33 + # 10 * 5 * ppc = 500 500 / 400 = 1.25 + # taking the minimul diff across permutations + # current to previous should be at least this times bigger + PRIOR_MIN_DIFF_1d = 1.25 + PRIOR_MAX_DIFF_1d = 1.50 def test_1d(self): This = type(self) - self._do_dim( - 1, This.PREVIOUS_ITERATION_MIN_DIFF_1d, This.PREVIOUS_ITERATION_MAX_DIFF_1d - ) - - """ 2d - refine 10x10 cells in 2d, ppc 100 - 10 * 10 * 4 * ppc = 400 - 10 * 10 * 8 * ppc = 800 800 / 400 = 1.5 - 10 * 10 * 9 * ppc = 900 900 / 800 = 1.125 - """ - PREVIOUS_ITERATION_MIN_DIFF_2d = 1.125 - PREVIOUS_ITERATION_MAX_DIFF_2d = 1.50 + self._do_dim(1, This.PRIOR_MIN_DIFF_1d, This.PRIOR_MAX_DIFF_1d) + + # 2d + # refine 10x10 cells in 2d, ppc 100 + # 10 * 10 * 4 * ppc = 400 + # 10 * 10 * 8 * ppc = 800 800 / 400 = 1.5 + # 10 * 10 * 9 * ppc = 900 900 / 800 = 1.125 + PRIOR_MIN_DIFF_2d = 1.125 + PRIOR_MAX_DIFF_2d = 1.50 def test_2d(self): This = type(self) - self._do_dim( - 2, This.PREVIOUS_ITERATION_MIN_DIFF_2d, This.PREVIOUS_ITERATION_MAX_DIFF_2d - ) + self._do_dim(2, This.PRIOR_MIN_DIFF_2d, This.PRIOR_MAX_DIFF_2d) - def tearDown(self): - # needed in case exception is raised in test and Simulator - # not reset properly - if self.simulator is not None: - self.simulator.reset() + # 3d + # refine 10x10x10 cells in 3d, ppc 100 + # 10 * 10 * 10 * 6 * ppc = 6000 + # 10 * 10 * 10 * 12 * ppc = 12000 - 12000 / 6000 = 2 + # 10 * 10 * 10 * 27 * ppc = 27000 - 27000 / 12000 = 2.25 + PRIOR_MIN_DIFF_3d = 2 + PRIOR_MAX_DIFF_3d = 2.25 + + def test_3d(self): + This = type(self) + self._do_dim(3, This.PRIOR_MIN_DIFF_3d, This.PRIOR_MAX_DIFF_3d) + + def _check_deltas_and_weights(self, dim, interp, refined_particle_nbr): + yaml_dim = self.yaml_root["dimension_" + str(dim)] + yaml_interp = yaml_dim["interp_" + str(interp)] + yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)] + + yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")] + yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")] + + splitter_t = splitter_type(dim, interp, refined_particle_nbr) + np.testing.assert_allclose(yaml_delta, splitter_t.delta) + np.testing.assert_allclose(yaml_weight, splitter_t.weight) + + def test_values(self): + from pyphare.pharein.simulation import valid_refined_particle_nbr + + for dim in range(1, 4): + for interp in range(1, 4): + for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]: + self._check_deltas_and_weights(dim, interp, refined_particle_nbr) if __name__ == "__main__": - unittest.main() + # the following ensures the order of execution so delta/weights are verified first + suite = unittest.TestSuite() + suite.addTest(SimulatorRefinedParticleNbr("test_values")) + for dim in range(1, 4): + suite.addTest(SimulatorRefinedParticleNbr("test_" + str(dim) + "d")) + unittest.TextTestRunner(failfast=True).run(suite) diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 569f2c873..ecb8db623 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -86,6 +86,9 @@ def getHierarchy( strict=True, ) + def S(x, x0, l): + return 0.5 * (1 + np.tanh((x - x0) / l)) + def bx(*xyz): return 1.0 @@ -259,29 +262,15 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time): # this is because the overlap box has been calculated from # the intersection of possibly shifted patch data ghost boxes - loc_b1 = boxm.amr_to_local( - box, boxm.shift(pd1.ghost_box, offsets[0]) + slice1 = boxm.select( + pd1.dataset, + boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0])), ) - loc_b2 = boxm.amr_to_local( - box, boxm.shift(pd2.ghost_box, offsets[1]) + slice2 = boxm.select( + pd2.dataset, + boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1])), ) - - data1 = pd1.dataset - data2 = pd2.dataset - - if box.ndim == 1: - slice1 = data1[loc_b1.lower[0] : loc_b1.upper[0] + 1] - slice2 = data2[loc_b2.lower[0] : loc_b2.upper[0] + 1] - - if box.ndim == 2: - slice1 = data1[ - loc_b1.lower[0] : loc_b1.upper[0] + 1, - loc_b1.lower[1] : loc_b1.upper[1] + 1, - ] - slice2 = data2[ - loc_b2.lower[0] : loc_b2.upper[0] + 1, - loc_b2.lower[1] : loc_b2.upper[1] + 1, - ] + assert slice1.dtype == np.float64 try: # empirical max absolute observed 5.2e-15 @@ -298,12 +287,9 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time): print(pd1.y.mean()) print(pd2.x.mean()) print(pd2.y.mean()) - print(loc_b1) - print(loc_b2) print(coarsest_time) print(slice1) print(slice2) - print(data1[:]) if self.rethrow_: raise e return diff_boxes(slice1, slice2, box) @@ -381,8 +367,9 @@ def _test_overlapped_particledatas_have_identical_particles( part2.iCells = part2.iCells + offsets[1] self.assertEqual(part1, part2) - def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100): - cells = 120 + def _test_L0_particle_number_conservation( + self, ndim, interp_order, ppc=100, cells=120 + ): time_step_nbr = 10 time_step = 0.001 @@ -409,7 +396,7 @@ def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100): self.assertEqual(n_particles, n_particles_at_t) def _test_field_coarsening_via_subcycles( - self, dim, interp_order, refinement_boxes, **kwargs + self, dim, interp_order, refinement_boxes, cells=60, **kwargs ): print( "test_field_coarsening_via_subcycles for dim/interp : {}/{}".format( @@ -423,12 +410,14 @@ def _test_field_coarsening_via_subcycles( time_step_nbr = 3 + diag_outputs = f"subcycle_coarsening/{dim}/{interp_order}/{self.ddt_test_id()}" datahier = self.getHierarchy( dim, interp_order, refinement_boxes, "fields", - cells=60, + cells=cells, + diag_outputs=diag_outputs, time_step=0.001, extra_diag_options={"fine_dump_lvl_max": 10}, time_step_nbr=time_step_nbr, @@ -501,16 +490,7 @@ def _test_field_coarsening_via_subcycles( afterCoarse = np.copy(coarse_pdDataset) # change values that should be updated to make failure obvious - assert dim < 3 # update - if dim == 1: - afterCoarse[ - dataBox.lower[0] : dataBox.upper[0] + 1 - ] = -144123 - if dim == 2: - afterCoarse[ - dataBox.lower[0] : dataBox.upper[0] + 1, - dataBox.lower[1] : dataBox.upper[1] + 1, - ] = -144123 + boxm.DataSelector(afterCoarse)[dataBox] = -144123 coarsen( qty, diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 2ea135909..70937fbf6 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -54,7 +54,6 @@ def getHierarchy( diag_outputs = self.unique_diag_dir_for_test_case( "phare_outputs/init", ndim, interp_order, diag_outputs ) - from pyphare.pharein import global_vars global_vars.sim = None @@ -252,19 +251,26 @@ def vthz(*xyz): ) return mom_hier - def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs): + def _test_B_is_as_provided_by_user(self, dim, interp_order, ppc=100, **kwargs): print( "test_B_is_as_provided_by_user : dim {} interp_order : {}".format( dim, interp_order ) ) + now = self.datetime_now() hier = self.getHierarchy( dim, interp_order, refinement_boxes=None, qty="b", + nbr_part_per_cell=ppc, + diag_outputs=f"test_b/{dim}/{interp_order}/{self.ddt_test_id()}", **kwargs, ) + print( + f"\n{self._testMethodName}_{dim}d init took {self.datetime_diff(now)} seconds" + ) + now = self.datetime_now() from pyphare.pharein import global_vars @@ -320,16 +326,44 @@ def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs): ) if dim == 3: - raise ValueError("Unsupported dimension") + zbx = bx_pd.z[:] + zby = by_pd.z[:] + zbz = bz_pd.z[:] - def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order): + xbx, ybx, zbx = [ + a.flatten() for a in np.meshgrid(xbx, ybx, zbx, indexing="ij") + ] + xby, yby, zby = [ + a.flatten() for a in np.meshgrid(xby, yby, zby, indexing="ij") + ] + xbz, ybz, zbz = [ + a.flatten() for a in np.meshgrid(xbz, ybz, zbz, indexing="ij") + ] + + np.testing.assert_allclose( + bx, bx_fn(xbx, ybx, zbx), atol=1e-16, rtol=0 + ) + np.testing.assert_allclose( + by, by_fn(xby, yby, zby).reshape(by.shape), atol=1e-16, rtol=0 + ) + np.testing.assert_allclose( + bz, bz_fn(xbz, ybz, zbz).reshape(bz.shape), atol=1e-16, rtol=0 + ) + + print(f"\n{self._testMethodName}_{dim}d took {self.datetime_diff(now)} seconds") + + def _test_bulkvel_is_as_provided_by_user( + self, dim, interp_order, ppc=100, **kwargs + ): hier = self.getHierarchy( dim, interp_order, {"L0": {"B0": nDBox(dim, 10, 19)}}, "moments", - nbr_part_per_cell=100, + nbr_part_per_cell=ppc, beam=True, + diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}", + **kwargs, ) from pyphare.pharein import global_vars @@ -347,117 +381,58 @@ def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order): for ip, patch in enumerate(level.patches): print("patch {}".format(ip)) - layout = patch.patch_datas["protons_Fx"].layout - centering = layout.centering["X"][ - patch.patch_datas["protons_Fx"].field_name - ] - nbrGhosts = layout.nbrGhosts(interp_order, centering) - - if dim == 1: - x = patch.patch_datas["protons_Fx"].x[nbrGhosts:-nbrGhosts] - - fpx = patch.patch_datas["protons_Fx"].dataset[nbrGhosts:-nbrGhosts] - fpy = patch.patch_datas["protons_Fy"].dataset[nbrGhosts:-nbrGhosts] - fpz = patch.patch_datas["protons_Fz"].dataset[nbrGhosts:-nbrGhosts] - - fbx = patch.patch_datas["beam_Fx"].dataset[nbrGhosts:-nbrGhosts] - fby = patch.patch_datas["beam_Fy"].dataset[nbrGhosts:-nbrGhosts] - fbz = patch.patch_datas["beam_Fz"].dataset[nbrGhosts:-nbrGhosts] - - ni = patch.patch_datas["rho"].dataset[nbrGhosts:-nbrGhosts] - - vxact = (fpx + fbx) / ni - vyact = (fpy + fby) / ni - vzact = (fpz + fbz) / ni - - vxexp = (nprot(x) * vx_fn(x) + nbeam(x) * vx_fn(x)) / ( - nprot(x) + nbeam(x) - ) - vyexp = (nprot(x) * vy_fn(x) + nbeam(x) * vy_fn(x)) / ( - nprot(x) + nbeam(x) - ) - vzexp = (nprot(x) * vz_fn(x) + nbeam(x) * vz_fn(x)) / ( - nprot(x) + nbeam(x) - ) - - for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)): - std = np.std(vexp - vact) - print("sigma(user v - actual v) = {}".format(std)) - self.assertTrue( - std < 1e-2 - ) # empirical value obtained from print just above - - def reshape(patch_data, nGhosts): + pdatas = patch.patch_datas + layout = pdatas["protons_Fx"].layout + centering = layout.centering["X"][pdatas["protons_Fx"].field_name] + nbrGhosts = layout.nbrGhosts( + interp_order, centering + ) # primal in all directions + select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(dim)]) + + def domain(patch_data): + if dim == 1: + return patch_data.dataset[select] return patch_data.dataset[:].reshape( - patch.box.shape + (nGhosts * 2) + 1 - ) - - if dim == 2: - xx, yy = np.meshgrid( - patch.patch_datas["protons_Fx"].x, - patch.patch_datas["protons_Fx"].y, - indexing="ij", - ) - - density = reshape(patch.patch_datas["rho"], nbrGhosts) - - protons_Fx = reshape(patch.patch_datas["protons_Fx"], nbrGhosts) - protons_Fy = reshape(patch.patch_datas["protons_Fy"], nbrGhosts) - protons_Fz = reshape(patch.patch_datas["protons_Fz"], nbrGhosts) - - beam_Fx = reshape(patch.patch_datas["beam_Fx"], nbrGhosts) - beam_Fy = reshape(patch.patch_datas["beam_Fy"], nbrGhosts) - beam_Fz = reshape(patch.patch_datas["beam_Fz"], nbrGhosts) - - x = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - y = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - - fpx = protons_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fpy = protons_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fpz = protons_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - - fbx = beam_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fby = beam_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fbz = beam_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - - ni = density[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - - vxact = (fpx + fbx) / ni - vyact = (fpy + fby) / ni - vzact = (fpz + fbz) / ni - - vxexp = (nprot(x, y) * vx_fn(x, y) + nbeam(x, y) * vx_fn(x, y)) / ( - nprot(x, y) + nbeam(x, y) - ) - vyexp = (nprot(x, y) * vy_fn(x, y) + nbeam(x, y) * vy_fn(x, y)) / ( - nprot(x, y) + nbeam(x, y) - ) - vzexp = (nprot(x, y) * vz_fn(x, y) + nbeam(x, y) * vz_fn(x, y)) / ( - nprot(x, y) + nbeam(x, y) - ) - - for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)): - self.assertTrue(np.std(vexp - vact) < 1e-2) + patch.box.shape + (nbrGhosts * 2) + 1 + )[select] + + ni = domain(pdatas["rho"]) + vxact = (domain(pdatas["protons_Fx"]) + domain(pdatas["beam_Fx"])) / ni + vyact = (domain(pdatas["protons_Fy"]) + domain(pdatas["beam_Fy"])) / ni + vzact = (domain(pdatas["protons_Fz"]) + domain(pdatas["beam_Fz"])) / ni + + select = pdatas["protons_Fx"].meshgrid(select) + vxexp = ( + nprot(*select) * vx_fn(*select) + nbeam(*select) * vx_fn(*select) + ) / (nprot(*select) + nbeam(*select)) + vyexp = ( + nprot(*select) * vy_fn(*select) + nbeam(*select) * vy_fn(*select) + ) / (nprot(*select) + nbeam(*select)) + vzexp = ( + nprot(*select) * vz_fn(*select) + nbeam(*select) * vz_fn(*select) + ) / (nprot(*select) + nbeam(*select)) + for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)): + self.assertTrue(np.std(vexp - vact) < 1e-2) + + def _test_density_is_as_provided_by_user(self, ndim, interp_order, **kwargs): + print( + f"test_density_is_as_provided_by_user : dim {ndim} interp_order {interp_order}" + ) - def _test_density_is_as_provided_by_user(self, dim, interp_order): empirical_dim_devs = { 1: 6e-3, 2: 3e-2, + 3: 2e-1, } - - nbParts = {1: 10000, 2: 1000} - print( - "test_density_is_as_provided_by_user : interp_order : {}".format( - interp_order - ) - ) + nbParts = {1: 10000, 2: 1000, 3: 20} hier = self.getHierarchy( - dim, + ndim, interp_order, - {"L0": {"B0": nDBox(dim, 10, 20)}}, + {"L0": {"B0": nDBox(ndim, 5, 14)}}, qty="moments", - nbr_part_per_cell=nbParts[dim], + nbr_part_per_cell=nbParts[ndim], beam=True, + **kwargs, ) from pyphare.pharein import global_vars @@ -479,34 +454,16 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order): layout = patch.patch_datas["rho"].layout centering = layout.centering["X"][patch.patch_datas["rho"].field_name] nbrGhosts = layout.nbrGhosts(interp_order, centering) + select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(ndim)]) - if dim == 1: - protons_expected = proton_density_fn(x[nbrGhosts:-nbrGhosts]) - beam_expected = beam_density_fn(x[nbrGhosts:-nbrGhosts]) - ion_expected = protons_expected + beam_expected + mesh = patch.patch_datas["rho"].meshgrid(select) + protons_expected = proton_density_fn(*mesh) + beam_expected = beam_density_fn(*mesh) + ion_expected = protons_expected + beam_expected - ion_actual = ion_density[nbrGhosts:-nbrGhosts] - beam_actual = beam_density[nbrGhosts:-nbrGhosts] - protons_actual = proton_density[nbrGhosts:-nbrGhosts] - - if dim == 2: - y = patch.patch_datas["rho"].y - xx, yy = np.meshgrid(x, y, indexing="ij") - - x0 = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - y0 = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - - protons_expected = proton_density_fn(x0, y0) - beam_expected = beam_density_fn(x0, y0) - ion_expected = protons_expected + beam_expected - - ion_actual = ion_density[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - beam_actual = beam_density[ - nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts - ] - protons_actual = proton_density[ - nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts - ] + protons_actual = proton_density[select] + beam_actual = beam_density[select] + ion_actual = ion_density[select] names = ("ions", "protons", "beam") expected = (ion_expected, protons_expected, beam_expected) @@ -519,11 +476,11 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order): for name, dev in devs.items(): print(f"sigma(user density - {name} density) = {dev}") self.assertLess( - dev, empirical_dim_devs[dim], f"{name} has dev = {dev}" + dev, empirical_dim_devs[ndim], f"{name} has dev = {dev}" ) def _test_density_decreases_as_1overSqrtN( - self, dim, interp_order, nbr_particles=None, cells=960 + self, ndim, interp_order, nbr_particles=None, cells=960 ): import matplotlib.pyplot as plt @@ -536,12 +493,12 @@ def _test_density_decreases_as_1overSqrtN( for inbr, nbrpart in enumerate(nbr_particles): hier = self.getHierarchy( - dim, + ndim, interp_order, None, "moments", nbr_part_per_cell=nbrpart, - diag_outputs=f"{nbrpart}", + diag_outputs=f"1overSqrtN/{ndim}/{interp_order}/{nbrpart}", density=lambda *xyz: np.zeros(tuple(_.shape[0] for _ in xyz)) + 1.0, smallest_patch_size=int(cells / 2), largest_patch_size=int(cells / 2), @@ -559,7 +516,7 @@ def _test_density_decreases_as_1overSqrtN( centering = layout.centering["X"][patch.patch_datas["rho"].field_name] nbrGhosts = layout.nbrGhosts(interp_order, centering) - select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(dim)]) + select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(ndim)]) ion_density = patch.patch_datas["rho"].dataset[:] mesh = patch.patch_datas["rho"].meshgrid(select) @@ -568,14 +525,14 @@ def _test_density_decreases_as_1overSqrtN( noise[inbr] = np.std(expected - actual) print(f"noise is {noise[inbr]} for {nbrpart} particles per cell") - if dim == 1: + if ndim == 1: x = patch.patch_datas["rho"].x plt.figure() - plt.plot(x[nbrGhosts:-nbrGhosts], actual, label="actual") - plt.plot(x[nbrGhosts:-nbrGhosts], expected, label="expected") + plt.plot(x[select], actual, label="actual") + plt.plot(x[select], expected, label="expected") plt.legend() plt.title(r"$\sigma =$ {}".format(noise[inbr])) - plt.savefig(f"noise_{nbrpart}_interp_{dim}_{interp_order}.png") + plt.savefig(f"noise_{nbrpart}_interp_{ndim}_{interp_order}.png") plt.close("all") plt.figure() @@ -587,7 +544,7 @@ def _test_density_decreases_as_1overSqrtN( ) plt.xlabel("nbr_particles") plt.legend() - plt.savefig(f"noise_nppc_interp_{dim}_{interp_order}.png") + plt.savefig(f"noise_nppc_interp_{ndim}_{interp_order}.png") plt.close("all") noiseMinusTheory = noise / noise[0] - 1 / np.sqrt( @@ -601,25 +558,32 @@ def _test_density_decreases_as_1overSqrtN( ) plt.xlabel("nbr_particles") plt.legend() - plt.savefig(f"noise_nppc_minus_theory_interp_{dim}_{interp_order}.png") + plt.savefig(f"noise_nppc_minus_theory_interp_{ndim}_{interp_order}.png") plt.close("all") self.assertGreater(3e-2, noiseMinusTheory[1:].mean()) def _test_nbr_particles_per_cell_is_as_provided( - self, dim, interp_order, default_ppc=100 + self, ndim, interp_order, ppc=100, **kwargs ): + ddt_test_id = self.ddt_test_id() datahier = self.getHierarchy( - dim, + ndim, interp_order, - {"L0": {"B0": nDBox(dim, 10, 20)}}, + {}, "particles", + diag_outputs=f"ppc/{ndim}/{interp_order}/{ddt_test_id}", + nbr_part_per_cell=ppc, + **kwargs, ) - for patch in datahier.level(0).patches: + if cpp.mpi_rank() > 0: + return + + for pi, patch in enumerate(datahier.level(0).patches): pd = patch.patch_datas["protons_particles"] icells = pd.dataset[patch.box].iCells - H, _ = np.histogramdd(icells) - self.assertTrue((H == default_ppc).all()) + H, edges = np.histogramdd(icells, bins=patch.box.shape) + self.assertTrue((H == ppc).all()) def _domainParticles_for(self, datahier, ilvl): patch0 = datahier.levels()[ilvl].patches[0] @@ -654,6 +618,9 @@ def _test_domainparticles_have_correct_split_from_coarser_particle( **kwargs, ) + if cpp.mpi_rank() > 0: + return + from pyphare.pharein.global_vars import sim assert sim is not None and len(sim.cells) == ndim @@ -691,15 +658,15 @@ def _test_domainparticles_have_correct_split_from_coarser_particle( part2 = coarse_split_particles[pop_name].select(patch.box) self.assertEqual(part1, part2) - def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs): + def _test_patch_ghost_on_refined_level_case(self, ndim, has_patch_ghost, **kwargs): import pyphare.pharein as ph out = "phare_outputs" - refinement_boxes = {"L0": [nDBox(dim, 10, 19)]} + refinement_boxes = {"L0": [nDBox(ndim, 10, 19)]} kwargs["interp_order"] = kwargs.get("interp_order", 1) kwargs["diag_outputs"] = f"{has_patch_ghost}" datahier = self.getHierarchy( - ndim=dim, + ndim, refinement_boxes=refinement_boxes, qty="particles_patch_ghost", **kwargs, @@ -722,12 +689,12 @@ def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs def _test_levelghostparticles_have_correct_split_from_coarser_particle( self, datahier ): - dim = datahier.level(0).patches[0].box.ndim + ndim = datahier.level(0).patches[0].box.ndim from pyphare.pharein.global_vars import sim assert sim is not None - assert len(sim.cells) == dim + assert len(sim.cells) == ndim particle_level_ghost_boxes_per_level = level_ghost_boxes(datahier, "particles")