From eb507bba93e3d064ca88913b06fa1fade3af5dcd Mon Sep 17 00:00:00 2001 From: dekken Date: Tue, 5 Oct 2021 17:40:16 +0200 Subject: [PATCH] 3d-ish --- pyphare/pyphare/core/box.py | 44 +++ pyphare/pyphare/core/gridlayout.py | 16 + .../pyphare/pharein/maxwellian_fluid_model.py | 98 ++---- pyphare/pyphare/pharein/simulation.py | 10 +- pyphare/pyphare/pharesee/geometry.py | 111 ++++--- pyphare/pyphare/pharesee/hierarchy.py | 33 +- res/amr/splitting.yml | 40 +++ res/cmake/def.cmake | 16 + res/cmake/options.cmake | 6 + src/amr/data/particles/refine/split.h | 1 + src/amr/data/particles/refine/split_1d.h | 24 +- src/amr/data/particles/refine/split_2d.h | 52 +-- src/amr/data/particles/refine/split_3d.h | 302 ++++++++++++++++++ src/amr/data/particles/refine/splitter.h | 18 +- .../hybrid_hybrid_messenger_strategy.h | 4 + src/core/data/grid/gridlayoutimplyee.h | 4 +- src/core/data/ndarray/ndarray_vector.h | 164 ++++++++-- src/core/utilities/meta/meta_utilities.h | 6 +- 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.h | 1 + .../data/field/refine/test_refine_field.py | 4 + .../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 | 66 ++++ tests/diagnostic/CMakeLists.txt | 2 + tests/diagnostic/job_3d.py.in | 14 + tests/diagnostic/test-diagnostics_3d.cpp | 35 ++ tests/diagnostic/test_diagnostics.h | 3 +- tests/simulator/__init__.py | 35 +- tests/simulator/advance/CMakeLists.txt | 3 + .../advance/test_fields_advance_1d.py | 2 +- .../advance/test_fields_advance_2d.py | 3 +- .../advance/test_fields_advance_3d.py | 88 +++++ .../advance/test_particles_advance_1d.py | 2 +- .../advance/test_particles_advance_2d.py | 3 +- .../advance/test_particles_advance_3d.py | 44 +++ tests/simulator/initialize/CMakeLists.txt | 3 + .../initialize/test_fields_init_1d.py | 2 +- .../initialize/test_fields_init_2d.py | 4 +- .../initialize/test_fields_init_3d.py | 44 +++ .../initialize/test_particles_init_1d.py | 2 +- .../initialize/test_particles_init_2d.py | 2 +- .../initialize/test_particles_init_3d.py | 94 ++++++ tests/simulator/per_test.h | 8 + tests/simulator/refined_particle_nbr.py | 124 +++---- tests/simulator/test_advance.py | 39 +-- tests/simulator/test_initialization.py | 178 +++++------ tests/simulator/test_python_concurrent.py | 59 +++- 54 files changed, 1555 insertions(+), 408 deletions(-) create mode 100644 src/amr/data/particles/refine/split_3d.h 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 diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py index 15d5e571d..46ae83ec9 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,46 @@ 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): + print("self.data.shape", self.data.shape) + data = self.__getitem__(box_or_slice) + print("data.shape", data.shape) + print("data[:].shape", data[:].shape) + data[:] = val diff --git a/pyphare/pyphare/core/gridlayout.py b/pyphare/pyphare/core/gridlayout.py index dc17663e0..5ae5e976b 100644 --- a/pyphare/pyphare/core/gridlayout.py +++ b/pyphare/pyphare/core/gridlayout.py @@ -275,6 +275,22 @@ 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): """ from a qty and a direction, returns a 1d array containing diff --git a/pyphare/pyphare/pharein/maxwellian_fluid_model.py b/pyphare/pyphare/pharein/maxwellian_fluid_model.py index 0a673f386..053709f39 100644 --- a/pyphare/pyphare/pharein/maxwellian_fluid_model.py +++ b/pyphare/pyphare/pharein/maxwellian_fluid_model.py @@ -142,104 +142,54 @@ def to_dict(self): #------------------------------------------------------------------------------ def validate(self, sim): - import itertools import numpy as np from ..core.gridlayout import GridLayout, directions as all_directions from ..core.box import Box from pyphare.pharein import py_fn_wrapper - dl = np.array(sim.dl) directions = all_directions[:sim.ndim] - domain_box = Box([0] * sim.ndim, np.asarray(sim.cells) - 1) assert len(sim.origin) == domain_box.ndim + ndim = domain_box.ndim def _primal_directions(qty): return [1 if yee_element_is_primal(qty, direction) else 0 for direction in directions] - def _nbrGhosts(qty, layout, primal_directions): + def _nbrGhosts(layout, qty, primal_directions): return [ - layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(sim.ndim) + layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(ndim) ] - def qty_shifts(qty, nbrGhosts): - shifts = [] - for dim in range(sim.ndim): - if yee_element_is_primal(qty, directions[dim]): - shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim] + 1) * dl[dim]] - else: - shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim]) * dl[dim] + (.5 * dl[dim])] - return shifts + def split_to_columns(ar): + ar = ar.reshape(-1, ndim) # drop outer arrays if any + return [ar[:,dim] for dim in range(ndim)] - def _1d_check(layout, nbrGhosts, primal_directions, qty, fn): - fnX = py_fn_wrapper(fn)(layout.yeeCoordsFor(qty, "x")) - include_border = primal_directions[0] - np.testing.assert_allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15) - return ( - np.allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15) - ) + def _nd_check(layout, qty, nbrGhosts, fn): + from pyphare.pharesee.geometry import hierarchy_overlaps + from pyphare.pharesee.hierarchy import hierarchy_from_box + from ..core import box as boxm + + hier = hierarchy_from_box(domain_box, nbrGhosts) + coords = layout.meshCoords(qty) + + for ilvl, overlaps in hierarchy_overlaps(hier).items(): + for overlap in overlaps: + (pd1, pd2), box, offsets = overlap["pdatas"], overlap["box"], overlap["offset"] + slice1 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0]))) + slice2 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1]))) + if not np.allclose(fn(*split_to_columns(slice1)) , fn(*split_to_columns(slice2)), atol=1e-15): + return False + return True - def split_to_columns(ar): - ar = ar.reshape(-1, sim.ndim) # drop outer arrays if any - return [ar[:,dim] for dim in range(sim.ndim)] - - - def _2d_check(nbrGhosts, primal_directions, qty, fn, shift): - layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order) - x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])] - y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])] - X,Y = np.meshgrid(x,y,indexing="ij") - coords = np.array([X.flatten(),Y.flatten()]).T - top = coords[:len(x)] - bot = coords[-len(x):] - left = coords[0 :: len(x)] - right = coords[len(x) - 1 :: len(x)] - return ( - np.allclose(fn(*split_to_columns(left)) , fn(*split_to_columns(right)), atol=1e-15) and - np.allclose(fn(*split_to_columns(top)) , fn(*split_to_columns(bot)), atol=1e-15) - ) - - def _3d_check(nbrGhosts, primal_directions, qty, fn, shift): - def get_3d_faces(M): - def get_face(M, dim, front_side): - return M[tuple((0 if front_side else -1) if i == dim else slice(None) for i in range(M.ndim))] - return [get_face(M, dim, front_side) for dim in range(sim.ndim) for front_side in (True, False)] - - layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order) - x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])] - y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])] - z = layout.yeeCoordsFor(qty, "z")[nbrGhosts[2]:-(nbrGhosts[2]-1)-(primal_directions[2])] - coords = np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape((len(x), len(y), len(z), sim.ndim)) - left, right, top, bot, front, back = get_3d_faces(coords) - return ( - np.allclose(fn(split_to_columns(left)) , fn(split_to_columns(right)), atol=1e-15) and - np.allclose(fn(split_to_columns(top)) , fn(split_to_columns(bot)), atol=1e-15) and - np.allclose(fn(split_to_columns(front)), fn(split_to_columns(back)), atol=1e-15) - ) def periodic_function_check(vec_field, dic): valid = True for xyz in ["x", "y", "z"]: qty = vec_field + xyz - fn = dic[qty] - layout = GridLayout(domain_box, sim.origin, sim.dl, sim.interp_order) - primal_directions = _primal_directions(qty) - nbrGhosts = _nbrGhosts(qty, layout, primal_directions) - - if sim.ndim == 1: - valid &= _1d_check(layout, nbrGhosts[0], primal_directions, qty, fn) - - if sim.ndim == 2: - permutations = itertools.product(*qty_shifts(qty, nbrGhosts)) - valid &= all([_2d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations]) - - if sim.ndim == 3: - # untested block - add in during 3d/splitting PR https://github.com/PHAREHUB/PHARE/pull/314 - permutations = itertools.product(*qty_shifts(qty, nbrGhosts)) - valid &= all([_3d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations]) - + nbrGhosts = _nbrGhosts(layout, qty, _primal_directions(qty)) + valid &= _nd_check(layout, qty, nbrGhosts, fn=dic[qty]) return valid if sim.boundary_types[0] == "periodic": diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 70dd0f100..872380187 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..fa75c7745 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)), + } + + 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)), } - raise ValueError("Unhandeled dimension") + + def touch_domain_border(box, domain_box, border): if border == "upper": @@ -70,13 +82,18 @@ def periodicity_shifts(domain_box): 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 +121,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 +243,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 +286,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 +316,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 +462,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 9e93f1325..a143374ca 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -84,13 +84,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): @@ -147,6 +141,19 @@ def __init__(self, layout, field_name, data, **kwargs): + def meshgrid(self, select=None): + def grid(): + if self.ndim == 1: + return [self.x] + if self.ndim == 2: + return np.meshgrid(self.x, self.y, indexing="ij") + return np.meshgrid(self.x, self.y, self.z, indexing="ij") + mesh = grid() + if select is not None: + return tuple(g[select] for g in mesh) + return mesh + + class ParticleData(PatchData): """ @@ -1388,6 +1395,18 @@ 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 ed218302f..2787e17ef 100644 --- a/res/cmake/def.cmake +++ b/res/cmake/def.cmake @@ -62,3 +62,19 @@ endif(ubsan) # msan is not supported - it's not practical to configure - use valgrind +# -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 2f2ae4280..947e35999 100644 --- a/res/cmake/options.cmake +++ b/res/cmake/options.cmake @@ -70,6 +70,12 @@ option(withCaliper "Use LLNL Caliper" OFF) option(lowResourceTests "Disable heavy tests for CI (2d/3d/etc" 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.h b/src/amr/data/particles/refine/split.h index f9b80a17f..021d91504 100644 --- a/src/amr/data/particles/refine/split.h +++ b/src/amr/data/particles/refine/split.h @@ -3,6 +3,7 @@ #include "amr/data/particles/refine/split_1d.h" #include "amr/data/particles/refine/split_2d.h" +#include "amr/data/particles/refine/split_3d.h" #endif // endif PHARE_SPLIT_H diff --git a/src/amr/data/particles/refine/split_1d.h b/src/amr/data/particles/refine/split_1d.h index a8553df97..5506c1f68 100644 --- a/src/amr/data/particles/refine/split_1d.h +++ b/src/amr/data/particles/refine/split_1d.h @@ -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.h b/src/amr/data/particles/refine/split_2d.h index e68c4d6b2..f229471d4 100644 --- a/src/amr/data/particles/refine/split_2d.h +++ b/src/amr/data/particles/refine/split_2d.h @@ -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.h b/src/amr/data/particles/refine/split_3d.h new file mode 100644 index 000000000..b59d5a688 --- /dev/null +++ b/src/amr/data/particles/refine/split_3d.h @@ -0,0 +1,302 @@ +/* +Splitting reference material can be found @ + https://github.com/PHAREHUB/PHARE/wiki/SplitPattern +*/ + +#ifndef PHARE_SPLIT_3D_H +#define PHARE_SPLIT_3D_H + +#include +#include +#include "core/utilities/point/point.h" +#include "core/utilities/types.h" +#include "splitter.h" + +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_H*/ diff --git a/src/amr/data/particles/refine/splitter.h b/src/amr/data/particles/refine/splitter.h index ee996b8e6..76deb864e 100644 --- a/src/amr/data/particles/refine/splitter.h +++ b/src/amr/data/particles/refine/splitter.h @@ -109,25 +109,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.h b/src/amr/messengers/hybrid_hybrid_messenger_strategy.h index cc992bc2c..a46f808eb 100644 --- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.h +++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.h @@ -138,6 +138,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); magneticSharedNodes_.registerLevel(hierarchy, level); @@ -226,6 +228,8 @@ namespace amr void initLevel(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, double const initDataTime) override { + PHARE_LOG_SCOPE("HybridHybridMessengerStrategy::initLevel"); + auto levelNumber = level.getLevelNumber(); magneticInit_.fill(levelNumber, initDataTime); diff --git a/src/core/data/grid/gridlayoutimplyee.h b/src/core/data/grid/gridlayoutimplyee.h index 630098b44..94be4a197 100644 --- a/src/core/data/grid/gridlayoutimplyee.h +++ b/src/core/data/grid/gridlayoutimplyee.h @@ -310,7 +310,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; @@ -330,7 +330,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.h b/src/core/data/ndarray/ndarray_vector.h index 2a3f4aebe..a3eee2000 100644 --- a/src/core/data/ndarray/ndarray_vector.h +++ b/src/core/data/ndarray/ndarray_vector.h @@ -2,6 +2,7 @@ #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_H #include +#include #include #include #include @@ -102,15 +103,17 @@ class MaskedView auto operator=(data_type value) { mask_.fill(array_, value); } auto xstart() const { return mask_.min(); } - auto xend() const { return shape_[0] - 1 - mask_.max(); } auto ystart() const { return mask_.min(); } - 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_; @@ -313,10 +316,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; } @@ -326,24 +329,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; } } @@ -351,7 +354,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 @@ -361,17 +401,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; } @@ -391,20 +437,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 @@ -421,7 +468,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()); @@ -450,7 +497,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/meta/meta_utilities.h b/src/core/utilities/meta/meta_utilities.h index 7b938bd49..1f25e8a0a 100644 --- a/src/core/utilities/meta/meta_utilities.h +++ b/src/core/utilities/meta/meta_utilities.h @@ -78,7 +78,11 @@ namespace core 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>, + + SimulatorOption, InterpConst<1>, 6, 12 /*, 27*/>, + SimulatorOption, InterpConst<2>, 6, 12>, + SimulatorOption, InterpConst<3>, 6, 12>>{}; } diff --git a/src/phare/phare.cpp b/src/phare/phare.cpp index 2283c227e..b19c62d05 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 8b4565912..974a822b0 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_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) # pybind fails with Werror +target_compile_options( + cpp + PRIVATE ${PHARE_FLAGS} ${PHARE_PYBIND_CXX_FLAGS}) # pybind fails with Werror 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_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg) + target_compile_options( + cpp_dbg + PRIVATE ${PHARE_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 3b914914f..7c370f592 100644 --- a/src/python3/cpp_simulator.cpp +++ b/src/python3/cpp_simulator.cpp @@ -234,3 +234,9 @@ PYBIND11_MODULE(PHARE_CPP_MOD_NAME, m) } // 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.h b/src/simulator/simulator.h index a4405a2cd..910a55790 100644 --- a/src/simulator/simulator.h +++ b/src/simulator/simulator.h @@ -68,6 +68,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 fdd8fbbac..249abd411 100644 --- a/tests/amr/data/field/refine/test_refine_field.py +++ b/tests/amr/data/field/refine/test_refine_field.py @@ -152,6 +152,10 @@ def refine(field, **kwargs): ) + if fine_box.ndim == 3: + assert False + + # the fine data is ghost_nbr larger than expected to include ghost values from the coarser, so we drop the outer values if field.box.ndim == 1: fine_data = fine_data[field.ghosts_nbr[0] : -field.ghosts_nbr[0]] 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 471cac85f..29da4c0b8 100644 --- a/tests/amr/data/particles/refine/test_split.cpp +++ b/tests/amr/data/particles/refine/test_split.cpp @@ -20,7 +20,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 665ecf94f..6505ac0aa 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 a14102307..44a21088e 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -747,6 +747,72 @@ 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 std::size_t dim = 3; + static constexpr std::uint32_t nx = 15, ny = 15, nz = 15; + static constexpr int start = 0, end = 5; + + ACollectionOfParticles_3d() + : 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(std::begin(particles), std::end(particles), rho, v, layout); + } + + GridLayout> layout{ + ConstArray(.1), {nx, ny}, ConstArray(0)}; + + ParticleArray particles; + Field, typename HybridQuantity::Scalar> rho, vx, vy, vz; + VecField, HybridQuantity> 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 32aad9384..77fd27623 100644 --- a/tests/diagnostic/CMakeLists.txt +++ b/tests/diagnostic/CMakeLists.txt @@ -32,9 +32,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.h b/tests/diagnostic/test_diagnostics.h index 971877b31..57870512a 100644 --- a/tests/diagnostic/test_diagnostics.h +++ b/tests/diagnostic/test_diagnostics.h @@ -65,7 +65,7 @@ void validateFluidGhosts(Data const& data, GridLayout const& layout, Field const std::size_t nans = 0; for (auto const& v : data) if (std::isnan(v)) - nans++; + ++nans; auto phyStartIdx = layout.physicalStartIndex(field.physicalQuantity(), core::Direction::X); core::NdArrayMask mask{0, phyStartIdx - 2}; @@ -105,7 +105,6 @@ void validateFluidGhosts(Data const& data, GridLayout const& layout, Field const } else if constexpr (dim == 3) { - throw std::runtime_error("not implemented"); } } diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index dad4fe1d7..c00828517 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -11,6 +11,8 @@ def parse_cli_args(pop_from_sys = True): sys.argv = [sys.argv[0]] return r + + # Block accidental dictionary key rewrites class NoOverwriteDict(dict): def __init__(self, dict): @@ -81,12 +83,13 @@ def density_2d_periodic(sim, x, y): xx, yy = meshify(x, y) return np.exp(-(xx-0.5*xmax)**2)*np.exp(-(yy-ymax/2.)**2) + background_particles -# 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.)**2)*np.exp(-(zz-zmax/2.)**2) + background_particles + return r def defaultPopulationSettings(sim, density_fn, vbulk_fn): @@ -178,3 +181,23 @@ def _diff(slice0): boxes += [Box([x, y, z], [x, y, z])] return boxes + + + +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..a0b5cf289 100644 --- a/tests/simulator/advance/CMakeLists.txt +++ b/tests/simulator/advance/CMakeLists.txt @@ -18,6 +18,9 @@ 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}) + 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() endif() diff --git a/tests/simulator/advance/test_fields_advance_1d.py b/tests/simulator/advance/test_fields_advance_1d.py index fc9579ea4..474fd8691 100644 --- a/tests/simulator/advance/test_fields_advance_1d.py +++ b/tests/simulator/advance/test_fields_advance_1d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box1D, nDBox +from pyphare.core.box import Box1D from tests.simulator.test_advance import AdvanceTestBase import matplotlib diff --git a/tests/simulator/advance/test_fields_advance_2d.py b/tests/simulator/advance/test_fields_advance_2d.py index 7ccfbc387..ea46f829f 100644 --- a/tests/simulator/advance/test_fields_advance_2d.py +++ b/tests/simulator/advance/test_fields_advance_2d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box2D, nDBox +from pyphare.core.box import Box2D from tests.simulator.test_advance import AdvanceTestBase import matplotlib @@ -50,6 +50,7 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self, 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}/{interp_order}/{self.ddt_test_id()}" 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(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size, time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, nbr_part_per_cell=ppc) 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..81f0c238e --- /dev/null +++ b/tests/simulator/advance/test_fields_advance_3d.py @@ -0,0 +1,88 @@ +""" + 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 +from ddt import ddt, data, unpack +from pyphare.core.box import Box3D +from tests.simulator.test_advance import AdvanceTestBase + +import matplotlib + +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 + diag_outputs=f"phare_overlaped_fields_are_equal_{ndim}_{self.ddt_test_id()}" + datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, cells=cells, + time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, 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 + diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts/{ndim}/{interp_order}/{self.ddt_test_id()}" + largest_patch_size, smallest_patch_size = check_patch_size(ndim, interp_order=interp_order, cells=[cells] * ndim) + datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, cells=cells, + smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size, + time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, 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) + + + # @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_1d.py b/tests/simulator/advance/test_particles_advance_1d.py index 145bf184d..a650f4505 100644 --- a/tests/simulator/advance/test_particles_advance_1d.py +++ b/tests/simulator/advance/test_particles_advance_1d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box1D, nDBox +from pyphare.core.box import Box1D from tests.simulator.test_advance import AdvanceTestBase import matplotlib diff --git a/tests/simulator/advance/test_particles_advance_2d.py b/tests/simulator/advance/test_particles_advance_2d.py index 69b2e7752..cbb229e52 100644 --- a/tests/simulator/advance/test_particles_advance_2d.py +++ b/tests/simulator/advance/test_particles_advance_2d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box2D, nDBox +from pyphare.core.box import Box2D from tests.simulator.test_advance import AdvanceTestBase import matplotlib @@ -33,6 +33,7 @@ def test_overlapped_particledatas_have_identical_particles(self, interp_order, r 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) 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..dc3139148 --- /dev/null +++ b/tests/simulator/advance/test_particles_advance_3d.py @@ -0,0 +1,44 @@ +""" + 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 +from ddt import ddt, data, unpack +from pyphare.core.box import Box3D +from tests.simulator.test_advance import AdvanceTestBase + +import matplotlib + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc = 10 + +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..52a751df6 100644 --- a/tests/simulator/initialize/CMakeLists.txt +++ b/tests/simulator/initialize/CMakeLists.txt @@ -18,6 +18,9 @@ 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}) + 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() diff --git a/tests/simulator/initialize/test_fields_init_1d.py b/tests/simulator/initialize/test_fields_init_1d.py index e2f4f5311..2380e1702 100644 --- a/tests/simulator/initialize/test_fields_init_1d.py +++ b/tests/simulator/initialize/test_fields_init_1d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box1D, nDBox + from tests.simulator.test_initialization import InitializationTest import matplotlib diff --git a/tests/simulator/initialize/test_fields_init_2d.py b/tests/simulator/initialize/test_fields_init_2d.py index a37c47599..a5b72fc9d 100644 --- a/tests/simulator/initialize/test_fields_init_2d.py +++ b/tests/simulator/initialize/test_fields_init_2d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box2D, nDBox + from tests.simulator.test_initialization import InitializationTest import matplotlib @@ -22,7 +22,7 @@ class InitializationTest(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..b4059c8c6 --- /dev/null +++ b/tests/simulator/initialize/test_fields_init_3d.py @@ -0,0 +1,44 @@ +""" + 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 numpy as np +import unittest +from ddt import ddt, data, unpack + +from tests.simulator.test_initialization import InitializationTest + +import matplotlib + +matplotlib.use("Agg") # for systems without GUI + +ndim = 3 +interp_orders = [1, 2, 3] +ppc, cells = 10, 20 + +@ddt +class InitializationTest(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 5ea562f4d..d79ddc553 100644 --- a/tests/simulator/initialize/test_particles_init_1d.py +++ b/tests/simulator/initialize/test_particles_init_1d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box1D, nDBox +from pyphare.core.box import Box1D from tests.simulator.test_initialization import InitializationTest import matplotlib diff --git a/tests/simulator/initialize/test_particles_init_2d.py b/tests/simulator/initialize/test_particles_init_2d.py index 023de7eea..2fbd0f133 100644 --- a/tests/simulator/initialize/test_particles_init_2d.py +++ b/tests/simulator/initialize/test_particles_init_2d.py @@ -5,7 +5,7 @@ import unittest from ddt import ddt, data, unpack -from pyphare.core.box import Box, Box2D, nDBox +from pyphare.core.box import Box2D from tests.simulator.test_initialization import InitializationTest import matplotlib 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..7f45f60a0 --- /dev/null +++ b/tests/simulator/initialize/test_particles_init_3d.py @@ -0,0 +1,94 @@ +""" + 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 +from ddt import ddt, data, unpack +from pyphare.core.box import Box, Box3D, nDBox +from tests.simulator.test_initialization import InitializationTest + +import matplotlib + +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 InitializationTest(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( + interp_order, + refinement_boxes, + "particles", + ndim=ndim, + cells=cells, nbr_part_per_cell=ppc, + diag_outputs=f"phare_outputs/test_levelghost/{ndim}/{interp_order}/{self.ddt_test_id()}", + ) + ) + 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.h b/tests/simulator/per_test.h index 925160641..9c59e49e1 100644 --- a/tests/simulator/per_test.h +++ b/tests/simulator/per_test.h @@ -121,6 +121,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); + namespace PHARE diff --git a/tests/simulator/refined_particle_nbr.py b/tests/simulator/refined_particle_nbr.py index ff1078a09..4f1bc047f 100644 --- a/tests/simulator/refined_particle_nbr.py +++ b/tests/simulator/refined_particle_nbr.py @@ -9,7 +9,6 @@ import numpy as np import pyphare.pharein as ph from pyphare.simulator.simulator import Simulator -from tests.simulator import NoOverwriteDict from tests.simulator import populate_simulation from pyphare.cpp import splitter_type @@ -28,41 +27,32 @@ 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 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() @@ -76,54 +66,76 @@ 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) - ) + 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 self.assertTrue(max_per_pop < prev_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) + + # 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 tearDown(self): - # needed in case exception is raised in test and Simulator - # not reset properly - if self.simulator is not None: - self.simulator.reset() + def _check_deltas_and_weights(self, dim, interp, refined_particle_nbr): + yaml_dim = self.yaml_root["dimension_" + str(dim)] + yaml_interp = yaml_dim["interp_" + str(interp)] + yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)] + + yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")] + yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")] + + splitter_t = splitter_type(dim, interp, refined_particle_nbr) + np.testing.assert_allclose(yaml_delta, splitter_t.delta) + np.testing.assert_allclose(yaml_weight, splitter_t.weight) + + + def test_values(self): + from pyphare.pharein.simulation import valid_refined_particle_nbr + for dim in range(1, 4): + for interp in range(1, 4): + for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]: + self._check_deltas_and_weights(dim, interp, refined_particle_nbr) if __name__ == "__main__": - unittest.main() + # the following ensures the order of execution so delta/weights are verified first + suite = unittest.TestSuite() + suite.addTest(SimulatorRefinedParticleNbr('test_values')) + for dim in range(1, 4): + suite.addTest(SimulatorRefinedParticleNbr('test_' + str(dim) + 'd')) + unittest.TextTestRunner(failfast=True).run(suite) diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 7f1c42293..48448c855 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -207,13 +207,14 @@ def vthz(*xyz): return mom_hier + + + def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time): checks = 0 for ilvl, overlaps in hierarchy_overlaps(datahier, coarsest_time).items(): for overlap in overlaps: - pd1, pd2 = overlap["pdatas"] - box = overlap["box"] - offsets = overlap["offset"] + (pd1, pd2), box, offsets = overlap["pdatas"], overlap["box"], overlap["offset"] self.assertEqual(pd1.quantity, pd2.quantity) @@ -229,22 +230,11 @@ 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])) - loc_b2 = 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] + slice1 = boxm.select(pd1.dataset, boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0]))) + slice2 = boxm.select(pd2.dataset, boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1]))) + assert slice1.dtype == np.float64 try: - assert slice1.dtype == np.float64 np.testing.assert_equal(slice1, slice2) checks += 1 except AssertionError as e: @@ -315,8 +305,8 @@ def _test_overlapped_particledatas_have_identical_particles(self, ndim, interp_o 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 @@ -339,7 +329,7 @@ def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100): - def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_boxes, **kwargs): + def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_boxes, cells=60, **kwargs): print("test_field_coarsening_via_subcycles for dim/interp : {}/{}".format(dim, interp_order)) from tests.amr.data.field.coarsening.test_coarsen_field import coarsen @@ -348,7 +338,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box time_step_nbr=3 diag_outputs=f"subcycle_coarsening/{dim}/{interp_order}/{self.ddt_test_id()}" - datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=60, + datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=cells, diag_outputs=diag_outputs, time_step=0.001, extra_diag_options={"fine_dump_lvl_max": 10}, time_step_nbr=time_step_nbr, @@ -411,12 +401,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box 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, coarse_pd, fine_pd, coarseBox, fine_pdDataset, afterCoarse) np.testing.assert_allclose(coarse_pdDataset, afterCoarse, atol=1e-16, rtol=0) diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 81d043489..cc75c890b 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -8,7 +8,7 @@ from pyphare.pharein.diagnostics import ParticleDiagnostics, FluidDiagnostics, ElectromagDiagnostics from pyphare.pharein import ElectronModel from pyphare.pharein.simulation import Simulation -from pyphare.pharesee.geometry import level_ghost_boxes, hierarchy_overlaps, touch_domain_border +from pyphare.pharesee.geometry import level_ghost_boxes from pyphare.pharesee.particles import aggregate as aggregate_particles import pyphare.core.box as boxm from pyphare.core.box import nDBox @@ -219,11 +219,14 @@ def vthz(*xyz): - 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)) - hier = self.getHierarchy(interp_order, refinement_boxes=None, qty="b", ndim=dim, + now = self.datetime_now() + hier = self.getHierarchy(interp_order, refinement_boxes=None, qty="b", ndim=dim, 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 model = global_vars.sim.model @@ -269,16 +272,28 @@ def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs): np.testing.assert_allclose(bz, bz_fn(xbz, ybz).reshape(bz.shape), atol=1e-16, rtol=0) 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) + 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): + + + + def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order, ppc=100, **kwargs): hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 19)}}, - "moments", nbr_part_per_cell=100, beam=True, ndim=dim, # ppc needs to be 10000? - diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}") + "moments", nbr_part_per_cell=ppc, beam=True, ndim=dim, # ppc needs to be 10000? + diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}", **kwargs) from pyphare.pharein import global_vars model = global_vars.sim.model @@ -294,86 +309,39 @@ 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): - 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] + 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)]) - fbx = beam_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fby = beam_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - fbz = beam_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] + def domain(patch_data): + if dim == 1: + return patch_data.dataset[select] + return patch_data.dataset[:].reshape(patch.box.shape + (nbrGhosts * 2) + 1)[select] - ni = density[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] + 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 - vxact = (fpx + fbx)/ni - vyact = (fpy + fby)/ni - vzact = (fpz + fbz)/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) - 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) - - - def _test_density_is_as_provided_by_user(self, dim, interp_order): - nbParts = {1 : 10000, 2: 1000} + def _test_density_is_as_provided_by_user(self, dim, interp_order, **kwargs): + nbParts = {1 : 10000, 2: 1000, 3: 20} print("test_density_is_as_provided_by_user : interp_order : {}".format(interp_order)) - hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 20)}}, + hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 5, 14)}}, qty="moments", nbr_part_per_cell=nbParts[dim], beam=True, ndim=dim, - diag_outputs=f"test_density/{dim}/{interp_order}/{self.ddt_test_id()}") + diag_outputs=f"test_density/{dim}/{interp_order}/{self.ddt_test_id()}", **kwargs) from pyphare.pharein import global_vars model = global_vars.sim.model @@ -393,6 +361,7 @@ 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(dim)]) if dim == 1: protons_expected = proton_density_fn(x[nbrGhosts:-nbrGhosts]) @@ -413,19 +382,18 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order): self.assertTrue(dev < 6e-3, '{} has dev = {}'.format(name, dev)) # empirical value obtained from test prints if dim == 2: - y = patch.patch_datas["rho"].y - xx, yy = np.meshgrid(x, y, indexing="ij") + xx, yy = patch.patch_datas["rho"].meshgrid() - x0 = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] - y0 = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts] + x0 = xx[select] + y0 = yy[select] 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] + ion_actual = ion_density[select] + beam_actual = beam_density[select] + protons_actual = proton_density[select] names = ("ions", "protons", "beam") expected = (ion_expected, protons_expected, beam_expected) @@ -436,15 +404,39 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order): print("sigma(user density - {} density) = {}".format(name, dev)) self.assertLess(dev, 3e-2, '{} has dev = {}'.format(name, dev)) # empirical value obtained from test prints + if dim == 3: + xx, yy, zz = patch.patch_datas["rho"].meshgrid() + + x0 = xx[select] + y0 = yy[select] + + protons_expected = proton_density_fn(x0, y0) + beam_expected = beam_density_fn(x0, y0) + ion_expected = protons_expected + beam_expected + + ion_actual = ion_density[select] + beam_actual = beam_density[select] + protons_actual = proton_density[select] + + names = ("ions", "protons", "beam") + expected = (ion_expected, protons_expected, beam_expected) + actual = (ion_actual, protons_actual, beam_actual) + devs = {name:np.std(expected-actual) for name, expected, actual in zip(names, expected, actual)} + + for name,dev in devs.items(): + print("sigma(user density - {} density) = {}".format(name, dev)) + self.assertLess(dev, 35e-2, '{} has dev = {}'.format(name, dev)) # empirical value obtained from test prints - def _test_density_decreases_as_1overSqrtN(self, dim, interp_order): + def _test_density_decreases_as_1overSqrtN(self, dim, interp_order, nbr_particles=None, cells=960): import matplotlib.pyplot as plt print(f"test_density_decreases_as_1overSqrtN, interp_order = {interp_order}") - nbr_particles = np.asarray([100, 1000, 5000, 10000]) + if nbr_particles is None: + nbr_particles = np.asarray([100, 1000, 5000, 10000]) + noise = np.zeros(len(nbr_particles)) for inbr,nbrpart in enumerate(nbr_particles): @@ -453,9 +445,9 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order): nbr_part_per_cell=nbrpart, diag_outputs=f"1overSqrtN/{dim}/{interp_order}/{nbrpart}", density=lambda x:np.zeros_like(x)+1., - smallest_patch_size=480, - largest_patch_size=480, - cells=960, + smallest_patch_size=int(cells/2), + largest_patch_size=int(cells/2), + cells=cells, dl=0.0125) from pyphare.pharein import global_vars @@ -485,7 +477,6 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order): plt.close("all") - plt.figure() plt.plot(nbr_particles, noise/noise[0], label=r"$\sigma/\sigma_0$") plt.plot(nbr_particles, 1/np.sqrt(nbr_particles/nbr_particles[0]), label=r"$1/sqrt(nppc/nppc0)$") @@ -508,17 +499,16 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order): - def _test_nbr_particles_per_cell_is_as_provided(self, dim, interp_order, default_ppc=100): + def _test_nbr_particles_per_cell_is_as_provided(self, dim, interp_order, ppc=100, **kwargs): ddt_test_id = self.ddt_test_id() datahier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 20)}}, "particles", ndim=dim, - diag_outputs=f"ppc/{dim}/{interp_order}/{ddt_test_id}") + diag_outputs=f"ppc/{dim}/{interp_order}/{ddt_test_id}", nbr_part_per_cell=ppc, **kwargs) for patch in datahier.level(0).patches: pd = patch.patch_datas["protons_particles"] icells = pd.dataset[patch.box].iCells H, edges = np.histogramdd(icells) - self.assertTrue((H == default_ppc).all()) - + self.assertTrue((H == ppc).all()) diff --git a/tests/simulator/test_python_concurrent.py b/tests/simulator/test_python_concurrent.py index 033df68ab..72a6a8b5b 100644 --- a/tests/simulator/test_python_concurrent.py +++ b/tests/simulator/test_python_concurrent.py @@ -11,36 +11,57 @@ 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 +from tests.simulator.initialize.test_fields_init_3d import InitializationTest as InitField3d +from tests.simulator.initialize.test_particles_init_3d import InitializationTest as InitParticles3d +from tests.simulator.advance.test_fields_advance_3d import AdvanceTest as AdvanceField3d +from tests.simulator.advance.test_particles_advance_3d import AdvanceTest as AdvanceParticles3d + N_CORES = int(os.environ["N_CORES"]) if "N_CORES" in os.environ else multiprocessing.cpu_count() +DO_DIM = int(os.environ["DO_DIM"]) if "DO_DIM" in os.environ else 0 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 +ALL_DIMS = DO_DIM == 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 - ] + test_classes_to_run = [] + + if ALL_DIMS: + test_classes_to_run += [SimulatorValidation] + + if ALL_DIMS or DO_DIM == 1: + test_classes_to_run += [ + InitField1d, + InitParticles1d, + AdvanceField1d, + AdvanceParticles1d, + ] + if ALL_DIMS or DO_DIM == 2: + test_classes_to_run += [ + InitField2d, + InitParticles2d, + AdvanceField2d, + AdvanceParticles2d, + ] + if ALL_DIMS or DO_DIM == 3: + test_classes_to_run += [ + InitField3d, + InitParticles3d, + AdvanceField3d, + AdvanceParticles3d + ] tests = [] loader = unittest.TestLoader() @@ -48,5 +69,13 @@ def test_cmd(clazz, test_id): for suite in loader.loadTestsFromTestCase(test_class): tests += [test_cmd(type(suite), suite._testMethodName)] - from tools.python3 import run_mp - run_mp(tests, N_CORES, check=True) + if PRINT: + for test in tests: + print(test) + else: + from tools.python3 import run_mp + try: + run_mp(tests, N_CORES, check=True) + except: + import sys + sys.exit(1) \ No newline at end of file