From 2c907d89c4172db1328620a9c9a59587f03066d6 Mon Sep 17 00:00:00 2001 From: dekken Date: Wed, 6 Dec 2023 17:37:31 +0100 Subject: [PATCH] 3d(mostly) --- .github/workflows/cmake_macos.yml | 2 +- .github/workflows/cmake_ubuntu.yml | 2 +- pyphare/pyphare/core/box.py | 40 +++ pyphare/pyphare/core/gridlayout.py | 14 + pyphare/pyphare/pharein/simulation.py | 10 +- pyphare/pyphare/pharesee/geometry.py | 113 ++++--- pyphare/pyphare/pharesee/hierarchy.py | 23 +- res/amr/splitting.yml | 40 +++ res/cmake/def.cmake | 16 + res/cmake/options.cmake | 11 + 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 | 19 +- .../hybrid_hybrid_messenger_strategy.hpp | 4 + src/core/data/grid/gridlayoutimplyee.hpp | 4 +- src/core/data/ndarray/ndarray_vector.hpp | 162 ++++++++-- src/core/utilities/box/box.hpp | 23 +- src/core/utilities/meta/meta_utilities.hpp | 21 +- src/core/utilities/point/point.hpp | 14 +- src/phare/phare.cpp | 13 +- src/phare/phare_init.py | 2 +- src/python3/CMakeLists.txt | 18 +- src/python3/cpp_simulator.cpp | 6 + src/simulator/simulator.hpp | 1 + .../data/field/refine/test_refine_field.py | 3 + .../refine/input/input_1d_ratio_2.txt | 2 +- .../refine/input/input_2d_ratio_2.txt | 2 +- .../refine/input/input_3d_ratio_2.txt | 51 +++ .../amr/data/particles/refine/test_split.cpp | 2 +- tests/core/data/ndarray/test_main.cpp | 54 ++++ .../core/numerics/interpolator/test_main.cpp | 74 +++++ tests/diagnostic/CMakeLists.txt | 2 + tests/diagnostic/job_3d.py.in | 14 + tests/diagnostic/test-diagnostics_3d.cpp | 35 ++ tests/diagnostic/test_diagnostics.hpp | 1 + tests/simulator/__init__.py | 47 ++- tests/simulator/advance/CMakeLists.txt | 5 + .../advance/test_fields_advance_2d.py | 9 +- .../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/refinement/test_2d_10_core.py | 6 +- tests/simulator/test_advance.py | 58 ++-- tests/simulator/test_initialization.py | 275 +++++++--------- tests/simulator/test_python_concurrent.py | 77 ----- tools/python3/__init__.py | 37 ++- 55 files changed, 1665 insertions(+), 471 deletions(-) create mode 100644 src/amr/data/particles/refine/split_3d.hpp create mode 100644 tests/amr/data/particles/refine/input/input_3d_ratio_2.txt create mode 100644 tests/diagnostic/job_3d.py.in create mode 100644 tests/diagnostic/test-diagnostics_3d.cpp 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 delete mode 100644 tests/simulator/test_python_concurrent.py diff --git a/.github/workflows/cmake_macos.yml b/.github/workflows/cmake_macos.yml index 9b7741231..74883fee9 100644 --- a/.github/workflows/cmake_macos.yml +++ b/.github/workflows/cmake_macos.yml @@ -77,7 +77,7 @@ jobs: cmake $GITHUB_WORKSPACE -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \ -DENABLE_SAMRAI_TESTS=OFF -DCMAKE_C_COMPILER_LAUNCHER=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DlowResourceTests=ON \ - -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -O2" + -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -O2 -DPHARE_SIMULATORS=2" - name: Build working-directory: ${{runner.workspace}}/build diff --git a/.github/workflows/cmake_ubuntu.yml b/.github/workflows/cmake_ubuntu.yml index fed78d6dc..e1fc7005c 100644 --- a/.github/workflows/cmake_ubuntu.yml +++ b/.github/workflows/cmake_ubuntu.yml @@ -80,7 +80,7 @@ jobs: cmake $GITHUB_WORKSPACE -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON --fresh \ -DCMAKE_C_COMPILER_LAUNCHER=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \ - -DlowResourceTests=ON -DdevMode=ON \ + -DlowResourceTests=ON -DdevMode=ON -Dphare_configurator=ON \ -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -O2" - name: Build diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py index 15d5e571d..9c6201292 100644 --- a/pyphare/pyphare/core/box.py +++ b/pyphare/pyphare/core/box.py @@ -67,6 +67,9 @@ def __eq__(self, other): return isinstance(other, Box) and (self.lower == other.lower).all() and (self.upper == other.upper).all() 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) @@ -183,5 +186,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/core/gridlayout.py b/pyphare/pyphare/core/gridlayout.py index 1f033d828..dc50d1d68 100644 --- a/pyphare/pyphare/core/gridlayout.py +++ b/pyphare/pyphare/core/gridlayout.py @@ -331,6 +331,20 @@ def yeeCoords(self, knode, iStart, centering, direction, ds, origin, derivOrder) return x + def meshCoords(self, qty): + ndim = self.ndim + assert ndim > 0 and ndim < 4 + x = self.yeeCoordsFor(qty, "x") + if ndim == 1: + return x + y = self.yeeCoordsFor(qty, "y") + if ndim == 2: + X,Y = np.meshgrid(x,y,indexing="ij") + return np.array([X.flatten(),Y.flatten()]).T.reshape((len(x), len(y), ndim)) + z = self.yeeCoordsFor(qty, "z") + X ,Y, Z = np.meshgrid(x,y,z,indexing="ij") + return np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape((len(x), len(y), len(z), ndim)) + def yeeCoordsFor(self, qty, direction, withGhosts=True, **kwargs): """ from a qty and a direction, returns a 1d array containing diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 2537259fb..7e4733ca3 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -194,9 +194,13 @@ def check_boundaries(ndim, **kwargs): 3: [4, 5, 8, 9, 25] }, 3: { - 1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27], - 2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64], - 3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125] + 1: [6, 12], + 2: [6, 12], + 3: [6, 12] + # full below + # 1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27], + # 2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64], + # 3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125] } } # Default refined_particle_nbr per dim/interp is considered index 0 of list def check_refined_particle_nbr(ndim, **kwargs): diff --git a/pyphare/pyphare/pharesee/geometry.py b/pyphare/pyphare/pharesee/geometry.py index d259fa8a6..b64b02f06 100644 --- a/pyphare/pyphare/pharesee/geometry.py +++ b/pyphare/pyphare/pharesee/geometry.py @@ -49,13 +49,25 @@ 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): if border == "upper": @@ -68,15 +80,22 @@ 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)], @@ -104,13 +123,33 @@ def periodicity_shifts(domain_box): "topleftright" : [ *shifts["leftright"], *shifts["top"], shifts["topleft"][-1], shifts["topright"][-1]], - "bottomtopleftright" : [ # one patch covers domain + "leftrightbottomtop" : [ # one patch covers domain *shifts["bottomleft"], *shifts["topright"], shifts["bottomright"][-1], shifts["topleft"][-1]] }) 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"], + ] + }) + + shifts = {"".join(sorted(k)) : l for k,l in shifts.items()} return shifts @@ -206,7 +245,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(): @@ -249,6 +288,7 @@ def hierarchy_overlaps(hierarchy, time=0): + def get_periodic_list(patches, domain_box, n_ghosts): """ given a list of patches and a domain box the function @@ -278,38 +318,50 @@ def get_periodic_list(patches, domain_box, n_ghosts): shift_patch(first_patch, domain_box.shape) sorted_patches.append(first_patch) - if dim == 2: + 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: + def borders_per(box): + return "".join([key for key, side in sides.items() if box * side is not None]) - in_sides = borders_per(boxm.grow(patch.box, n_ghosts)) + for patch in patches: - 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) + 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 + + def ghost_area_boxes(hierarchy, quantities, levelNbrs=[], time=0): """ this function returns boxes representing ghost cell boxes for all levels @@ -412,18 +464,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None): 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.py b/pyphare/pyphare/pharesee/hierarchy.py index 4046e6fe3..2c7c2ec03 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -100,13 +100,7 @@ def select(self, box): overlap = box * gbox if overlap is not None: - lower = self.layout.AMRToLocal(overlap.lower) - upper = self.layout.AMRToLocal(overlap.upper) - - if box.ndim == 1: - return self.dataset[lower[0] : upper[0] + 1] - if box.ndim == 2: - return self.dataset[lower[0] : upper[0] + 1, lower[1] : upper[1] + 1] + return boxm.select(self.dataset, self.layout.AMRBoxToLocal(overlap)) return np.array([]) def __getitem__(self, box): @@ -182,6 +176,8 @@ def grid(): return mesh + + class ParticleData(PatchData): """ Concrete type of PatchData representing particles in a region @@ -1615,6 +1611,19 @@ def hierarchy_from_sim(simulator, qty, pop=""): return PatchHierarchy(patch_levels, domain_box, time=simulator.currentTime()) +def hierarchy_from_box(domain_box, ghosts_nbr): + """ + constructs a basic hierarchy with one patch for level 0 as the entire domain + """ + layout = GridLayout(domain_box, np.asarray([0] * domain_box.ndim), [.1] * domain_box.ndim, 1) + pdata = PatchData(layout, "qty") + object.__setattr__(pdata, "ghosts_nbr", np.asarray(ghosts_nbr)) + object.__setattr__(pdata, "ghost_box", boxm.grow(layout.box, ghosts_nbr)) + return PatchHierarchy( + {0 : PatchLevel(0, [Patch({"qty": pdata})])}, domain_box + ) + + def hierarchy_from( simulator=None, qty=None, pop="", h5_filename=None, time=None, hier=None ): 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/def.cmake b/res/cmake/def.cmake index 911d21737..e0fcf8b9f 100644 --- a/res/cmake/def.cmake +++ b/res/cmake/def.cmake @@ -269,3 +269,19 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests endif() +# -DwithPGO_GEN +if (withPGO_GEN) + set (PHARE_PGO_FLAGS -fprofile-generate) +endif(withPGO_GEN) + +# -DwithPGO_USE +if (withPGO_USE) + set (PHARE_PGO_FLAGS -fprofile-use) +endif(withPGO_USE) + +if (withPGO_GEN OR withPGO_USE) + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${PHARE_PGO_FLAGS}" ) + set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} ${PHARE_PGO_FLAGS}" ) + set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} ${PHARE_PGO_FLAGS}" ) +endif() # (withPGO_GEN OR withPGO_USE) + diff --git a/res/cmake/options.cmake b/res/cmake/options.cmake index a335b90d7..068196de0 100644 --- a/res/cmake/options.cmake +++ b/res/cmake/options.cmake @@ -66,10 +66,21 @@ option(withCaliper "Use LLNL Caliper" 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" ON) + # -DtestDuringBuild=ON enabled if devMode=ON, disabled if asan=ON (needs LD_PRELOAD) option(testDuringBuild "Runs C++ unit tests after they are built" OFF) + + +# -DwithPGO_GEN +option(withPGO_GEN "Activate generate step of PGO") +# -DwithPGO_A +option(withPGO_USE "Activate usage step of PGO") + + # Controlling the activation of tests if (NOT DEFINED PHARE_EXEC_LEVEL_MIN) set(PHARE_EXEC_LEVEL_MIN 1) 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..d2bd416a8 100644 --- a/src/amr/data/particles/refine/splitter.hpp +++ b/src/amr/data/particles/refine/splitter.hpp @@ -67,6 +67,7 @@ class PatternDispatcher using FineParticle = decltype(particles[0]); // may be a reference core::apply(patterns, [&](auto const& pattern) { + auto weight = static_cast(pattern.weight_); for (size_t rpIndex = 0; rpIndex < pattern.deltas_.size(); rpIndex++) { FineParticle fineParticle = particles[idx++]; @@ -114,25 +115,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 BrownDispatcher +struct PurplePattern { }; template -struct PinkDispatcher +struct BrownPattern +{ +}; +template +struct LimePattern +{ +}; +template +struct WhitePattern { }; diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp index 5e790bd25..fe9ee7932 100644 --- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp +++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp @@ -174,6 +174,8 @@ namespace amr void registerLevel(std::shared_ptr const& hierarchy, int const levelNumber) override { + PHARE_LOG_SCOPE("HybridHybridMessengerStrategy::registerLevel"); + auto const level = hierarchy->getPatchLevel(levelNumber); magSharedNodesRefiners_.registerLevel(hierarchy, level); @@ -282,6 +284,8 @@ namespace amr void initLevel(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, double const initDataTime) override { + PHARE_LOG_SCOPE("HybridHybridMessengerStrategy::initLevel"); + auto levelNumber = level.getLevelNumber(); magneticInitRefiners_.fill(levelNumber, initDataTime); diff --git a/src/core/data/grid/gridlayoutimplyee.hpp b/src/core/data/grid/gridlayoutimplyee.hpp index 98987fe49..34e35ec19 100644 --- a/src/core/data/grid/gridlayoutimplyee.hpp +++ b/src/core/data/grid/gridlayoutimplyee.hpp @@ -311,7 +311,7 @@ namespace core * the following is only valid when dual and primal do not have the same number of ghosts * and that depends on the interp order - * It is commented out because ghosts are hard coded to 5 for now. + * It is commented out because ghosts are hard coded to 3 for now. * if constexpr (interp_order == 1 || interp_order == 2 || interp_order == 4) return -1; @@ -331,7 +331,7 @@ namespace core * the following is only valid when dual and primal do not have the same number of ghosts * and that depends on the interp order - * It is commented out because ghosts are hard coded to 5 for now. + * It is commented out because ghosts are hard coded to 3 for now. * if constexpr (interp_order == 1 || interp_order == 2 || interp_order == 4) return 1; diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index c9bef9837..0c177b0bf 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 @@ -120,6 +121,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_; @@ -310,10 +315,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; } @@ -323,24 +328,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; } } @@ -348,7 +353,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 @@ -358,17 +400,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; } @@ -388,20 +436,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 @@ -418,7 +467,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()); @@ -447,7 +496,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 fd3e1fe09..2ab28e6f6 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -7,6 +7,7 @@ #include "core/utilities/meta/meta_utilities.hpp" #include "core/def.hpp" +#include #include #include #include @@ -145,6 +146,26 @@ struct Box else return 6; } + + auto as_points() const + { + std::vector> points; + if constexpr (dim == 1) + for (Type i = lower[0]; i <= upper[0]; ++i) + points.emplace_back(Point{i}); + + else if constexpr (dim == 2) + for (Type i = lower[0]; i <= upper[0]; ++i) + for (Type j = lower[1]; j <= upper[1]; ++j) + points.emplace_back(Point{i, j}); + else + for (Type i = lower[0]; i <= upper[0]; ++i) + for (Type j = lower[1]; j <= upper[1]; ++j) + for (Type k = lower[2]; k <= upper[2]; ++k) + points.emplace_back(Point{i, j, k}); + + return points; + } }; template @@ -199,7 +220,7 @@ class box_iterator template -Box(Point lower, Point upper) -> Box; +Box(Point lower, Point upper)->Box; 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/src/core/utilities/point/point.hpp b/src/core/utilities/point/point.hpp index 754b0a12e..c5c6c6097 100644 --- a/src/core/utilities/point/point.hpp +++ b/src/core/utilities/point/point.hpp @@ -178,7 +178,7 @@ namespace core template Point(Indexes... indexes) - -> Point>::type, sizeof...(indexes)>; + ->Point>::type, sizeof...(indexes)>; template auto& operator<<(std::ostream& os, Point const& p) @@ -204,6 +204,18 @@ NO_DISCARD PHARE::core::Point abs(PHARE::core::Point const return postive; } +template +auto to_string(PHARE::core::Point const& point) +{ + return point.str(); +} + +template +auto to_string(PHARE::core::Point&& point) +{ + return point.str(); +} + } // namespace std diff --git a/src/phare/phare.cpp b/src/phare/phare.cpp index fec4d4c06..16dab12ba 100644 --- a/src/phare/phare.cpp +++ b/src/phare/phare.cpp @@ -32,7 +32,8 @@ std::unique_ptr fromCommandLine(int argc, char return nullptr; } -int main(int argc, char** argv) + +int naim(int argc, char** argv) { std::string const welcome = R"~( _____ _ _ _____ ______ @@ -76,4 +77,14 @@ int main(int argc, char** argv) std::cout << simulator->currentTime() << "\n"; // time += simulator.timeStep(); } + + return 0; +} + + +#if !defined(PHARE_EXE_NO_MAIN) +int main(int argc, char** argv) +{ + return naim(argc, argv); } +#endif diff --git a/src/phare/phare_init.py b/src/phare/phare_init.py index 41aa24932..cb863be3f 100644 --- a/src/phare/phare_init.py +++ b/src/phare/phare_init.py @@ -20,7 +20,7 @@ max_nbr_levels=2, # (default=1) max nbr of levels in the AMR hierarchy refinement = "tagging", #refinement_boxes={"L0": {"B0": [(10, ), (20, )]}}, - diag_options={"format": "phareh5", "options": {"dir": "phare_outputs"}} + diag_options={"format": "phareh5", "options": {"dir": "phare_outputs", "mode":"overwrite"}} ) diff --git a/src/python3/CMakeLists.txt b/src/python3/CMakeLists.txt index d78a2d59d..298656ed4 100644 --- a/src/python3/CMakeLists.txt +++ b/src/python3/CMakeLists.txt @@ -2,23 +2,37 @@ cmake_minimum_required (VERSION 3.9) project(phare_python3) +set(PHARE_PYBIND_CXX_FLAGS -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) +set(PHARE_PYBIND_LD_FLAGS ) + +# set(PHARE_PYBIND_CXX_FLAGS "-DPHARE_EXEC_PYBIND=1 -nostartfiles -Wl,--entry=naim") +# set(PHARE_PYBIND_LD_FLAGS "${PHARE_PYBIND_CXX_FLAGS}") +# set(PHARE_PYBIND_CXX_FLAGS "${PHARE_PYBIND_CXX_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}") + pybind11_add_module(cpp cpp_simulator.cpp) target_link_libraries(cpp PUBLIC phare_simulator) -target_compile_options(cpp PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) +target_compile_options( + cpp + PRIVATE ${PHARE_WERROR_FLAGS} ${PHARE_PYBIND_CXX_FLAGS}) set_target_properties(cpp PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs" + LINK_FLAGS "${PHARE_PYBIND_LD_FLAGS}" ) + # this is on by default "pybind11_add_module" but can interfere with coverage so we disable it if coverage is enabled set_property(TARGET cpp PROPERTY INTERPROCEDURAL_OPTIMIZATION ${PHARE_INTERPROCEDURAL_OPTIMIZATION}) if (CMAKE_BUILD_TYPE STREQUAL "Debug") pybind11_add_module(cpp_dbg cpp_simulator.cpp) target_link_libraries(cpp_dbg PUBLIC phare_simulator) - target_compile_options(cpp_dbg PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg) + target_compile_options( + cpp_dbg + PRIVATE ${PHARE_WERROR_FLAGS} ${PHARE_PYBIND_CXX_FLAGS} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg) set_target_properties(cpp_dbg PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs" + LINK_FLAGS "${PHARE_PYBIND_LD_FLAGS}" ) set_property(TARGET cpp_dbg PROPERTY INTERPROCEDURAL_OPTIMIZATION ${PHARE_INTERPROCEDURAL_OPTIMIZATION}) endif (CMAKE_BUILD_TYPE STREQUAL "Debug") diff --git a/src/python3/cpp_simulator.cpp b/src/python3/cpp_simulator.cpp index a5c0a5d1c..03e55b582 100644 --- a/src/python3/cpp_simulator.cpp +++ b/src/python3/cpp_simulator.cpp @@ -24,3 +24,9 @@ PYBIND11_MODULE(PHARE_CPP_MOD_NAME, m) declarePatchData, 3>(m, "PatchDataVectorDouble_3D"); } } // namespace PHARE::pydata + +#if defined(PHARE_EXEC_PYBIND) +#define PHARE_EXE_NO_MAIN 1 +#include "phare/phare.cpp" +#undef PHARE_EXE_NO_MAIN +#endif diff --git a/src/simulator/simulator.hpp b/src/simulator/simulator.hpp index 5475925fa..401ad59b4 100644 --- a/src/simulator/simulator.hpp +++ b/src/simulator/simulator.hpp @@ -78,6 +78,7 @@ class Simulator : public ISimulator Simulator(PHARE::initializer::PHAREDict const& dict, std::shared_ptr const& hierarchy); + ~Simulator() { if (coutbuf != nullptr) diff --git a/tests/amr/data/field/refine/test_refine_field.py b/tests/amr/data/field/refine/test_refine_field.py index 88f0bbd77..d644fe694 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/refine/input/input_1d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt index 86e8b0c0c..496be4bf0 100644 --- a/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt +++ b/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt @@ -1,6 +1,6 @@ CartesianGridGeometry { - domain_boxes = [ (0) , (64) ] + domain_boxes = [ (0) , (20) ] x_lo = 0.0 x_up = 1.0 periodic_dimension = 1 diff --git a/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt index 1b0d9ef37..7d7d5d6bb 100644 --- a/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt +++ b/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt @@ -1,6 +1,6 @@ CartesianGridGeometry { - domain_boxes = [ (0, 0) , (64, 64) ] + domain_boxes = [ (0, 0) , (20, 20) ] x_lo = 0.0, 0.0 x_up = 1.0, 1.0 periodic_dimension = 1, 1 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/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 dc0f01635..d60e5e849 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -796,6 +796,80 @@ 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; + + using PHARE_TYPES = PHARE::core::PHARE_Types; + using NdArray_t = typename PHARE_TYPES::Array_t; + using ParticleArray_t = typename PHARE_TYPES::ParticleArray_t; + using GridLayout_t = typename PHARE_TYPES::GridLayout_t; + constexpr static auto safeLayer = static_cast(1 + ghostWidthForParticles()); + + ACollectionOfParticles_3d() + : particles{grow(layout.AMRBox(), safeLayer)} + , rho{"field", HybridQuantity::Scalar::rho, nx, ny, nz} + , vx{"v_x", HybridQuantity::Scalar::Vx, nx, ny, nz} + , vy{"v_y", HybridQuantity::Scalar::Vy, nx, ny, nz} + , vz{"v_z", HybridQuantity::Scalar::Vz, nx, ny, nz} + , v{"v", HybridQuantity::Vector::V} + { + v.setBuffer("v_x", &vx); + v.setBuffer("v_y", &vy); + v.setBuffer("v_z", &vz); + + 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); + } + + GridLayout> layout{ + ConstArray(.1), {nx, ny}, ConstArray(0)}; + + ParticleArray_t particles; + Field rho, vx, vy, vz; + VecField v; + Interpolator interpolator; +}; +TYPED_TEST_SUITE_P(ACollectionOfParticles_3d); + + +TYPED_TEST_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d) +{ + // EXPECT_DOUBLE_EQ(this->rho(7, 7, 7), 1.0); + // EXPECT_DOUBLE_EQ(this->vx(7, 7, 7), 2.0); + // EXPECT_DOUBLE_EQ(this->vy(7, 7, 7), -1.0); + // EXPECT_DOUBLE_EQ(this->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 3a5964a5f..b7b5378c1 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..5df839152 --- /dev/null +++ b/tests/diagnostic/test-diagnostics_3d.cpp @@ -0,0 +1,35 @@ + +#include + +#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/diagnostic/test_diagnostics.hpp b/tests/diagnostic/test_diagnostics.hpp index c39df301b..3033fb99e 100644 --- a/tests/diagnostic/test_diagnostics.hpp +++ b/tests/diagnostic/test_diagnostics.hpp @@ -20,6 +20,7 @@ constexpr unsigned NEW_HI5_FILE = HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate; + template auto checkField(HighFiveFile const& hifile, GridLayout const& layout, Field const& field, std::string const path, FieldFilter const ff = FieldFilter{}) diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index 56ca9e9d9..45f8e112b 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -93,12 +93,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): @@ -260,3 +265,33 @@ def clean_up_diags_dirs(self): for diag_dir in self.diag_dirs: if os.path.exists(diag_dir): shutil.rmtree(diag_dir) + + +def debug_tracer(): + """ + print live stack trace during execution + """ + + import sys, os + + def tracefunc(frame, event, arg, indent=[0]): + filename = os.path.basename(frame.f_code.co_filename) + line_number = frame.f_lineno + if event == "call": + indent[0] += 2 + print( + "-" * indent[0] + "> enter function", + frame.f_code.co_name, + f"{filename} {line_number}", + ) + elif event == "return": + print( + "<" + "-" * indent[0], + "exit function", + frame.f_code.co_name, + f"{filename} {line_number}", + ) + indent[0] -= 2 + return tracefunc + + sys.setprofile(tracefunc) diff --git a/tests/simulator/advance/CMakeLists.txt b/tests/simulator/advance/CMakeLists.txt index eb51cca8c..6c7e0caff 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) # not ready + # 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_2d.py b/tests/simulator/advance/test_fields_advance_2d.py index 0a5540897..bc822c0cd 100644 --- a/tests/simulator/advance/test_fields_advance_2d.py +++ b/tests/simulator/advance/test_fields_advance_2d.py @@ -62,6 +62,11 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts( largest_patch_size, smallest_patch_size = check_patch_size( ndim, interp_order=interp_order, cells=[60] * ndim ) + print( + "largest_patch_size, smallest_patch_size", + largest_patch_size, + smallest_patch_size, + ) datahier = self.getHierarchy( ndim, interp_order, @@ -85,8 +90,8 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts( ( { "L0": {"B0": Box2D(5, 20)}, - "L1": {"B0": Box2D(12, 38)}, - "L2": {"B0": Box2D(30, 52)}, + "L1": {"B0": Box2D(12, 37)}, + "L2": {"B0": Box2D(30, 51)}, } ) ), 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 51808cbdc..38ee5de50 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 91d428d09..d06d4f43e 100644 --- a/tests/simulator/per_test.hpp +++ b/tests/simulator/per_test.hpp @@ -122,6 +122,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/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py index cb96cad2d..2dcd42e7b 100644 --- a/tests/simulator/refinement/test_2d_10_core.py +++ b/tests/simulator/refinement/test_2d_10_core.py @@ -28,7 +28,7 @@ def config(diag_outputs, model_init={}, refinement_boxes=None): ph.global_vars.sim = None - Simulation( + sim = Simulation( smallest_patch_size=10, largest_patch_size=20, time_step_nbr=1, @@ -52,12 +52,15 @@ def config(diag_outputs, model_init={}, refinement_boxes=None): ) def density(x, y): + Lx = sim.simulation_domain()[0] return 1.0 def bx(x, y): return 0.1 def by(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] return 0.2 def bz(x, y): @@ -140,6 +143,7 @@ def get_hier(path): from tests.simulator.test_advance import AdvanceTestBase + cpp = cpp_lib() test = AdvanceTestBase(rethrow=True) # change to False for debugging images L0_diags = "phare_outputs/test_x_homo_0" diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 56c579308..769dd82d9 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -31,7 +31,7 @@ def _density(*xyz): hL = np.array(sim.simulation_domain()) / 2 _ = lambda i: -((xyz[i] - hL[i]) ** 2) - return 0.3 + np.exp(sum([_(i) for i in range(len(xyz))])) + return 0.3 + np.exp(sum([_(i) for i, v in enumerate(xyz)])) def getHierarchy( self, @@ -84,6 +84,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 @@ -264,29 +267,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: assert slice1.dtype == np.float64 @@ -305,12 +294,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) @@ -388,8 +374,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 @@ -416,7 +403,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( @@ -430,12 +417,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, @@ -508,16 +497,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 66bab88a8..427fa6ff1 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -52,7 +52,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 @@ -256,19 +255,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 @@ -324,16 +330,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[:] + + 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 + ) - def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order): + 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 @@ -351,117 +385,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 @@ -483,34 +458,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) @@ -523,11 +480,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 @@ -540,12 +497,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), @@ -563,7 +520,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) @@ -572,14 +529,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() @@ -591,7 +548,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( @@ -605,25 +562,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] @@ -658,6 +622,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 @@ -695,15 +662,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") diff --git a/tests/simulator/test_python_concurrent.py b/tests/simulator/test_python_concurrent.py deleted file mode 100644 index bde9dd6b4..000000000 --- a/tests/simulator/test_python_concurrent.py +++ /dev/null @@ -1,77 +0,0 @@ -""" - This script exists to minimize testing time by running all simulation/tests - concurrently without needing to wait for any particular file or set of tests -""" - -import os -import unittest -import multiprocessing - -from tests.simulator.test_validation import SimulatorValidation - -from tests.simulator.initialize.test_fields_init_1d import ( - InitializationTest as InitField1d, -) -from tests.simulator.initialize.test_particles_init_1d import ( - InitializationTest as InitParticles1d, -) - -from tests.simulator.advance.test_fields_advance_1d import AdvanceTest as AdvanceField1d -from tests.simulator.advance.test_particles_advance_1d import ( - AdvanceTest as AdvanceParticles1d, -) - -from tests.simulator.initialize.test_fields_init_2d import ( - InitializationTest as InitField2d, -) -from tests.simulator.initialize.test_particles_init_2d import ( - InitializationTest as InitParticles2d, -) - -from tests.simulator.advance.test_fields_advance_2d import AdvanceTest as AdvanceField2d -from tests.simulator.advance.test_particles_advance_2d import ( - AdvanceTest as AdvanceParticles2d, -) - - -N_CORES = ( - int(os.environ["N_CORES"]) - if "N_CORES" in os.environ - else multiprocessing.cpu_count() -) -MPI_RUN = int(os.environ["MPI_RUN"]) if "MPI_RUN" in os.environ else 1 -PRINT = int(os.environ["PRINT"]) if "PRINT" in os.environ else 0 - - -def test_cmd(clazz, test_id): - return ( - f"mpirun -n {MPI_RUN} python3 -m {clazz.__module__} {clazz.__name__}.{test_id}" - ) - - -if __name__ == "__main__": - test_classes_to_run = [ - SimulatorValidation, - InitField1d, - InitParticles1d, - AdvanceField1d, - AdvanceParticles1d, - InitField2d, - InitParticles2d, - AdvanceField2d, - AdvanceParticles2d, - ] - - tests = [] - loader = unittest.TestLoader() - for test_class in test_classes_to_run: - for suite in loader.loadTestsFromTestCase(test_class): - tests += [test_cmd(type(suite), suite._testMethodName)] - - from tools.python3 import run_mp - - if PRINT: - for test in tests: - print(test) - else: - run_mp(tests, N_CORES, check=True) diff --git a/tools/python3/__init__.py b/tools/python3/__init__.py index 16db65902..6e8ecf55d 100644 --- a/tools/python3/__init__.py +++ b/tools/python3/__init__.py @@ -1,6 +1,18 @@ def decode_bytes(input): return input.decode("ascii", errors="ignore") +class RunTimer: + def __init__(self, cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwargs): + import time + import subprocess + self.cmd = cmd + start = time.time() + self.run = subprocess.run(self.cmd, shell=shell, capture_output=capture_output, check=check, **kwargs) + self.t = time.time() - start + self.stdout = self.run.stdout + self.stderr = self.run.stderr + + def run(cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwargs): """ @@ -11,10 +23,13 @@ def run(cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwa if print_cmd: print(f"running: {cmd}") try: - return subprocess.run(cmd, shell=shell, capture_output=capture_output, check=check, **kwargs) + return RunTimer(cmd, shell=shell, capture_output=capture_output, check=check, **kwargs) except subprocess.CalledProcessError as e: # only triggers on failure if check=True - print(f"run failed with error: {e}\n\t{e.stdout}\n\t{e.stderr} ") - raise RuntimeError(decode_bytes(e.stderr)) + what = f"run failed with error: {e}" + print(what) + if capture_output: + raise RuntimeError(decode_bytes(e.stderr)) + raise RuntimeError(what) def run_mp(cmds, N_CORES=None, **kwargs): """ @@ -31,9 +46,11 @@ def run_mp(cmds, N_CORES=None, **kwargs): results = [] for future in concurrent.futures.as_completed(jobs): try: - results += [future.result()] + proc = future.result() + results += [proc] if future.exception() is not None: raise future.exception() + print(proc.cmd, f"finished in {proc.t:.2f} seconds") except Exception as exc: if kwargs.get("check", False): executor.shutdown(wait=False, cancel_futures=True) @@ -43,11 +60,20 @@ def run_mp(cmds, N_CORES=None, **kwargs): return results +def find_on_path(file): + import os + for dir in os.environ["PATH"].split(os.pathsep): + full = os.path.join(dir, file) + if os.path.exists(full): + return full + return "" + + def binary_exists_on_path(bin): """ https://linux.die.net/man/1/which """ - return run(f"which {bin}").returncode == 0 + return len(find_on_path(bin)) def scan_dir(path, files_only=False, dirs_only=False, drop=[]): @@ -75,3 +101,4 @@ def pushd(new_cwd): yield finally: os.chdir(cwd) +