From 3e5266da4b469ff13ded4c9c5527a0117412954a Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 26 Sep 2024 16:55:25 -0700 Subject: [PATCH 01/27] Implements global indexing --- src/amrex/extensions/MultiFab.py | 479 ++++++++++++++++++++++++++++++- 1 file changed, 477 insertions(+), 2 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index ab0f1379..a7e1ff8f 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -1,13 +1,22 @@ """ This file is part of pyAMReX -Copyright 2023 AMReX community -Authors: Axel Huebl +Copyright 2024 AMReX community +Authors: Axel Huebl, David Grote License: BSD-3-Clause-LBNL """ from .Iterator import next +import numpy as np + +try: + from mpi4py import MPI as mpi + comm_world = mpi.COMM_WORLD + npes = comm_world.Get_size() +except ImportError: + npes = 1 + def mf_to_numpy(self, copy=False, order="F"): """ @@ -173,6 +182,467 @@ def copy_multifab(amr, self): return mf +def imesh(self, idir, include_ghosts=False): + """Returns the integer mesh along the specified direction with the appropriate centering. + This is the location of the data points in grid cell units. + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + direction : integer + Zero based direction number. + In a typical Cartesian case, 0 would be 'x' direction. + include_ghosts : bool, default=False + Whether or not ghost cells are included in the mesh. + """ + + min_box = self.box_array().minimal_box() + ilo = min_box.small_end[idir] + ihi = min_box.big_end[idir] + + if include_ghosts: + # The ghost cells are added to the upper and lower end of the global domain. + nghosts = self.n_grow_vect + ilo -= nghosts[idir] + ihi += nghosts[idir] + + # The centering shift + ix_type = self.box_array().ix_type() + if ix_type.node_centered(idir): + # node centered + shift = 0.0 + else: + # cell centered + shift = 0.5 + + return np.arange(ilo, ihi + 1) + shift + + +def shape(self, include_ghosts=False): + """Returns the shape of the global array + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + include_ghosts : bool, default=False + Whether or not ghost cells are included + """ + min_box = self.box_array().minimal_box() + result = min_box.size + if include_ghosts: + result = result + 2*self.n_grow_vect + result = list(result) + [self.nComp] + return tuple(result) + + +def _get_indices(index, missing): + """Expand the index list to length three. + + Parameters + ---------- + index: sequence of length dims + The indices for each dim + + missing: + The value used to fill in the extra dimensions added + """ + return list(index) + (3 - len(index))*[missing] + + +def _get_min_indices(self, include_ghosts): + """Returns the minimum indices, expanded to length 3 + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + include_ghosts : bool, default=False + Whether or not ghost cells are included + """ + min_box = self.box_array().minimal_box() + if include_ghosts: + min_box.grow(self.n_grow_vect) + return _get_indices(min_box.small_end, 0) + + +def _get_max_indices(self, include_ghosts): + """Returns the maximum indices, expanded to length 3. + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + include_ghosts : bool, default=False + Whether or not ghost cells are included + """ + min_box = self.box_array().minimal_box() + if include_ghosts: + min_box.grow(self.n_grow_vect) + return _get_indices(min_box.big_end, 0) + + +def _fix_index(self, ii, imax, d, include_ghosts): + """Handle negative index, wrapping them as needed. + When ghost cells are included, the indices are + shifted by the number of ghost cells before being wrapped. + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + ii : integer + The index to be wrapped + imax : integer + The maximum value that the index could have + d : integer + The direction of the index + include_ghosts: bool + Whether or not ghost cells are included + """ + nghosts = list(_get_indices(self.n_grow_vect, 0)) + [0] + if include_ghosts: + ii += nghosts[d] + if ii < 0: + ii += imax + if include_ghosts: + ii -= nghosts[d] + return ii + + +def _find_start_stop(self, ii, imin, imax, d, include_ghosts): + """Given the input index, calculate the start and stop range of the indices. + + Parameters + ---------- + ii : None, slice, or integer + Input index, either None, a slice object, or an integer. + Note that ii can be negative. + imin : integer + The global lowest lower bound in the specified direction. + This can include the ghost cells. + imax : integer + The global highest upper bound in the specified direction. + This can include the ghost cells. + This should be the max index + 1. + d : integer + The dimension number, 0, 1, 2, or 3 (3 being the components) + include_ghosts : bool + Whether or not ghost cells are included + + If ii is a slice, the start and stop values are used directly, + unless they are None, then the lower or upper bound is used. + An assertion checks if the indices are within the bounds. + """ + if ii is None: + iistart = imin + iistop = imax + elif isinstance(ii, slice): + if ii.start is None: + iistart = imin + else: + iistart = _fix_index(self, ii.start, imax, d, include_ghosts) + if ii.stop is None: + iistop = imax + else: + iistop = _fix_index(self, ii.stop, imax, d, include_ghosts) + else: + ii = _fix_index(self, ii, imax, d, include_ghosts) + iistart = ii + iistop = ii + 1 + assert imin <= iistart <= imax, Exception( + f"Dimension {d+1} lower index is out of bounds" + ) + assert imin <= iistop <= imax, Exception( + f"Dimension {d+1} upper index is out of bounds" + ) + return iistart, iistop + +def _get_field(self, mfi, include_ghosts): + """Return the field at the given mfi. + If include ghosts is true, return the whole array, otherwise + return the interior slice that does not include the ghosts. + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + mfi : amrex.MFIiter + Index to the FAB of the MultiFab + include_ghosts : bool, default=False + Whether or not ghost cells are included + """ + # Note that the array will always have 4 dimensions. + # even when dims < 3. + # The transpose is taken since the Python array interface to Array4 in + # self.array(mfi) is in C ordering. + # Note: transposing creates a view and not a copy. + import inspect + amr = inspect.getmodule(self) + if amr.Config.have_gpu: + device_arr = self.array(mfi).to_cupy(copy=False, order="F") + else: + device_arr = self.array(mfi).to_numpy(copy=False, order="F") + if not include_ghosts: + device_arr = device_arr[ + tuple([slice(ng, -ng) for ng in self.n_grow_vect]) + ] + return device_arr + + +def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, include_ghosts): + """Return the slices where the block intersects with the global slice. + If the block does not intersect, return None. + This also shifts the block slices by the number of ghost cells in the + MultiFab arrays since the arrays include the ghost cells. + + Parameters + ---------- + mfi : MFIter + The MFIter instance for the current block, + starts : sequence + The minimum indices of the global slice. + These can be negative. + stops : sequence + The maximum indices of the global slice. + These can be negative. + icstart : integer + The minimum component index of the global slice. + These can be negative. + icstops : integer + The maximum component index of the global slice. + These can be negative. + include_ghosts : bool, default=False + Whether or not ghost cells are included + + Returns + ------- + block_slices : tuple or None + The slices of the intersections relative to the block + global_slices : tuple or None + The slices of the intersections relative to the global array where the data from individual block will go + """ + box = mfi.tilebox() + if include_ghosts: + box.grow(self.n_grow_vect) + + ilo = _get_indices(box.small_end, 0) + ihi = _get_indices(box.big_end, 0) + + # Add 1 to the upper end to be consistent with the slicing notation + ihi_p1 = [i + 1 for i in ihi] + i1 = np.maximum(starts, ilo) + i2 = np.minimum(stops, ihi_p1) + + if np.all(i1 < i2): + block_slices = [] + global_slices = [] + for i in range(3): + block_slices.append(slice(i1[i] - ilo[i], i2[i] - ilo[i])) + global_slices.append(slice(i1[i] - starts[i], i2[i] - starts[i])) + + block_slices.append(slice(icstart, icstop)) + global_slices.append(slice(0, icstop - icstart)) + + return tuple(block_slices), tuple(global_slices) + else: + return None, None + + +def __getitem__(self, index): + """Returns slice of the MultiFab using global indexing. + The shape of the object returned depends on the number of ix, iy and iz specified, which + can be from none to all three. Note that the values of ix, iy and iz are + relative to the fortran indexing, meaning that 0 is the lower boundary + of the whole domain, and in fortran ordering, i.e. [ix,iy,iz]. + This allows negative indexing, though with ghosts cells included, the first n-ghost negative + indices will refer to the lower guard cells. + + Parameters + ---------- + index : integer, or sequence of integers or slices, or Ellipsis + Index of the slice to return + """ + # Temporary value until fixed + include_ghosts = False + + # Get the number of dimensions. Is there a cleaner way to do this? + dims = len(self.n_grow_vect) + + # Note that the index can have negative values (which wrap around) and has 1 added to the upper + # limit using python style slicing + if index == Ellipsis: + index = dims * [slice(None)] + elif isinstance(index, slice): + # If only one slice passed in, it was not wrapped in a list + index = [index] + + if len(index) < dims + 1: + # Add extra dims to index, including for the component. + # These are the dims left out and assumed to extend over the full size of the dim + index = list(index) + while len(index) < dims + 1: + index.append(slice(None)) + elif len(index) > dims + 1: + raise Exception("Too many indices given") + + # Expand the indices to length 3 + ii = _get_indices(index, None) + ic = index[-1] + + # Global extent. These include the ghost cells when include_ghosts is True + ixmin, iymin, izmin = _get_min_indices(self, include_ghosts) + ixmax, iymax, izmax = _get_max_indices(self, include_ghosts) + + # Setup the size of the array to be returned + ixstart, ixstop = _find_start_stop(self, ii[0], ixmin, ixmax + 1, 0, include_ghosts) + iystart, iystop = _find_start_stop(self, ii[1], iymin, iymax + 1, 1, include_ghosts) + izstart, izstop = _find_start_stop(self, ii[2], izmin, izmax + 1, 2, include_ghosts) + icstart, icstop = _find_start_stop(self, ic, 0, self.n_comp, 3, include_ghosts) + + # Gather the data to be included in a list to be sent to other processes + starts = [ixstart, iystart, izstart] + stops = [ixstop, iystop, izstop] + datalist = [] + for mfi in self: + block_slices, global_slices = _get_intersect_slice( + self, mfi, starts, stops, icstart, icstop, include_ghosts + ) + if global_slices is not None: + # Note that the array will always have 4 dimensions. + device_arr = _get_field(self, mfi, include_ghosts) + slice_arr = device_arr[block_slices] + try: + # Copy data from host to device using cupy syntax + slice_arr = slice_arr.get() + except AttributeError: + # Array is already a numpy array on the host + pass + datalist.append((global_slices, slice_arr)) + + # Gather the data from all processors + if npes == 1: + all_datalist = [datalist] + else: + all_datalist = comm_world.allgather(datalist) + + # Create the array to be returned + result_shape = ( + max(0, ixstop - ixstart), + max(0, iystop - iystart), + max(0, izstop - izstart), + max(0, icstop - icstart), + ) + + # Now, copy the data into the result array + result_global = None + for datalist in all_datalist: + for global_slices, f_arr in datalist: + if result_global is None: + # Delay allocation to here so that the type can be obtained + result_global = np.zeros(result_shape, dtype=f_arr.dtype) + result_global[global_slices] = f_arr + + if result_global is None: + # Something went wrong with the index and no data was found. Return an empty array. + result_global = np.zeros(0) + + # Remove dimensions of length 1, and if all dimensions + # are removed, return a scalar (that's what the [()] does) + return result_global.squeeze()[()] + + +def __setitem__(self, index, value): + """Sets slices of a decomposed array. + The shape of the input object depends on the number of arguments specified, which can + be from none to all three. + This allows negative indexing, though with ghosts cells included, the first n-ghost negative + indices will refer to the lower guard cells. + + Parameters + ---------- + index : integer, or sequence of integers or slices, or Ellipsis + The slice to set + value : scalar or array + Input value to assign to the specified slice of the MultiFab + """ + # Temporary value until fixed + include_ghosts = False + # Get the number of dimensions. Is there a cleaner way to do this? + dims = len(self.n_grow_vect) + + # Note that the index can have negative values (which wrap around) and has 1 added to the upper + # limit using python style slicing + if index == Ellipsis: + index = tuple(dims * [slice(None)]) + elif isinstance(index, slice): + # If only one slice passed in, it was not wrapped in a list + index = [index] + + if len(index) < dims + 1: + # Add extra dims to index, including for the component. + # These are the dims left out and assumed to extend over the full size of the dim. + index = list(index) + while len(index) < dims + 1: + index.append(slice(None)) + elif len(index) > dims + 1: + raise Exception("Too many indices given") + + # Expand the indices to length 3 + ii = _get_indices(index, None) + ic = index[-1] + + # Global extent. These include the ghost cells when include_ghosts is True + ixmin, iymin, izmin = _get_min_indices(self, include_ghosts) + ixmax, iymax, izmax = _get_max_indices(self, include_ghosts) + + # Setup the size of the global array to be set + ixstart, ixstop = _find_start_stop(self, ii[0], ixmin, ixmax + 1, 0, include_ghosts) + iystart, iystop = _find_start_stop(self, ii[1], iymin, iymax + 1, 1, include_ghosts) + izstart, izstop = _find_start_stop(self, ii[2], izmin, izmax + 1, 2, include_ghosts) + icstart, icstop = _find_start_stop(self, ic, 0, self.n_comp, 3, include_ghosts) + + if isinstance(value, np.ndarray): + # Expand the shape of the input array to match the shape of the global array + # (it needs to be 4-D). + # This converts value to an array if needed, and the [...] grabs a view so + # that the shape change below doesn't affect value. + value3d = np.array(value)[...] + global_shape = list(value3d.shape) + # The shape of 1 is added for the extra dimensions and when index is an integer + # (in which case the dimension was not in the input array). + if not isinstance(ii[0], slice): + global_shape[0:0] = [1] + if not isinstance(ii[1], slice): + global_shape[1:1] = [1] + if not isinstance(ii[2], slice): + global_shape[2:2] = [1] + if not isinstance(ic, slice) or len(global_shape) < 4: + global_shape[3:3] = [1] + value3d.shape = global_shape + + if libwarpx.libwarpx_so.Config.have_gpu: + # check if cupy is available for use + xp, cupy_status = load_cupy() + if cupy_status is not None: + libwarpx.amr.Print(cupy_status) + + starts = [ixstart, iystart, izstart] + stops = [ixstop, iystop, izstop] + for mfi in self: + block_slices, global_slices = _get_intersect_slice( + self, mfi, starts, stops, icstart, icstop, include_ghosts + ) + if global_slices is not None: + mf_arr = _get_field(self, mfi, include_ghosts) + if isinstance(value, np.ndarray): + # The data is copied from host to device automatically if needed + mf_arr[block_slices] = value3d[global_slices] + else: + mf_arr[block_slices] = value + + def register_MultiFab_extension(amr): """MultiFab helper methods""" @@ -191,3 +661,8 @@ def register_MultiFab_extension(amr): amr.MultiFab.copy = lambda self: copy_multifab(amr, self) amr.MultiFab.copy.__doc__ = copy_multifab.__doc__ + + amr.MultiFab.imesh = imesh + amr.MultiFab.shape = shape + amr.MultiFab.__getitem__ = __getitem__ + amr.MultiFab.__setitem__ = __setitem__ From b454490ac263a49b1003499f5d6c05fd2004f94c Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 26 Sep 2024 18:17:12 -0700 Subject: [PATCH 02/27] Implemented fix from WarpX PR #4370 --- src/amrex/extensions/MultiFab.py | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index a7e1ff8f..c872eacf 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -391,7 +391,8 @@ def _get_field(self, mfi, include_ghosts): return device_arr -def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, include_ghosts): +def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, + include_ghosts, with_internal_ghosts): """Return the slices where the block intersects with the global slice. If the block does not intersect, return None. This also shifts the block slices by the number of ghost cells in the @@ -415,6 +416,8 @@ def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, include_ghos These can be negative. include_ghosts : bool, default=False Whether or not ghost cells are included + with_internal_ghosts: bool + Whether the internal ghosts are included in the slices Returns ------- @@ -424,11 +427,25 @@ def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, include_ghos The slices of the intersections relative to the global array where the data from individual block will go """ box = mfi.tilebox() + box_small_end = box.small_end + box_big_end = box.big_end if include_ghosts: - box.grow(self.n_grow_vect) + nghosts = self.mf.n_grow_vect + box.grow(nghosts) + if with_internal_ghosts: + box_small_end = box.small_end + box_big_end = box.big_end + else: + min_box = self.box_array().minimal_box() + for i in range(len(nghosts)): + if box_small_end[i] == min_box.small_end[i]: + box_small_end[i] -= nghosts[i] + if box_big_end[i] == min_box.big_end[i]: + box_big_end[i] += nghosts[i] - ilo = _get_indices(box.small_end, 0) - ihi = _get_indices(box.big_end, 0) + boxlo = _get_indices(box.small_end, 0) + ilo = _get_indices(box_small_end, 0) + ihi = _get_indices(box_big_end, 0) # Add 1 to the upper end to be consistent with the slicing notation ihi_p1 = [i + 1 for i in ihi] @@ -439,7 +456,7 @@ def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, include_ghos block_slices = [] global_slices = [] for i in range(3): - block_slices.append(slice(i1[i] - ilo[i], i2[i] - ilo[i])) + block_slices.append(slice(i1[i] - boxlo[i], i2[i] - boxlo[i])) global_slices.append(slice(i1[i] - starts[i], i2[i] - starts[i])) block_slices.append(slice(icstart, icstop)) @@ -507,7 +524,7 @@ def __getitem__(self, index): datalist = [] for mfi in self: block_slices, global_slices = _get_intersect_slice( - self, mfi, starts, stops, icstart, icstop, include_ghosts + self, mfi, starts, stops, icstart, icstop, include_ghosts, False ) if global_slices is not None: # Note that the array will always have 4 dimensions. @@ -632,7 +649,7 @@ def __setitem__(self, index, value): stops = [ixstop, iystop, izstop] for mfi in self: block_slices, global_slices = _get_intersect_slice( - self, mfi, starts, stops, icstart, icstop, include_ghosts + self, mfi, starts, stops, icstart, icstop, include_ghosts, True ) if global_slices is not None: mf_arr = _get_field(self, mfi, include_ghosts) From 5c064957c6c13688f2d1ed70313f07c9e4d317ec Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 27 Sep 2024 09:14:59 -0700 Subject: [PATCH 03/27] Made shape a propery and added shape_with_ghosts --- src/amrex/extensions/MultiFab.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index c872eacf..3236c856 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -232,11 +232,22 @@ def shape(self, include_ghosts=False): min_box = self.box_array().minimal_box() result = min_box.size if include_ghosts: - result = result + 2*self.n_grow_vect + result = result + self.n_grow_vect*2 result = list(result) + [self.nComp] return tuple(result) +def shape_with_ghost(self): + """Returns the shape of the global array including ghost cells + + Parameters + ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX + """ + return shape(self, include_ghosts=True) + + def _get_indices(index, missing): """Expand the index list to length three. @@ -680,6 +691,7 @@ def register_MultiFab_extension(amr): amr.MultiFab.copy.__doc__ = copy_multifab.__doc__ amr.MultiFab.imesh = imesh - amr.MultiFab.shape = shape + amr.MultiFab.shape = property(shape) + amr.MultiFab.shape_with_ghost = property(shape_with_ghost) amr.MultiFab.__getitem__ = __getitem__ amr.MultiFab.__setitem__ = __setitem__ From 4c2f800b204ad864c7980704e61fb90ed7684536 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 27 Sep 2024 09:15:47 -0700 Subject: [PATCH 04/27] Fix typo --- src/amrex/extensions/MultiFab.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 3236c856..8ac2883d 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -237,7 +237,7 @@ def shape(self, include_ghosts=False): return tuple(result) -def shape_with_ghost(self): +def shape_with_ghosts(self): """Returns the shape of the global array including ghost cells Parameters @@ -692,6 +692,6 @@ def register_MultiFab_extension(amr): amr.MultiFab.imesh = imesh amr.MultiFab.shape = property(shape) - amr.MultiFab.shape_with_ghost = property(shape_with_ghost) + amr.MultiFab.shape_with_ghosts = property(shape_with_ghosts) amr.MultiFab.__getitem__ = __getitem__ amr.MultiFab.__setitem__ = __setitem__ From c704851f2d967f2471fac5a83665de6554001dc3 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 27 Sep 2024 09:34:09 -0700 Subject: [PATCH 05/27] Update documentation --- src/amrex/extensions/MultiFab.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 8ac2883d..b102be1d 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -582,7 +582,7 @@ def __getitem__(self, index): def __setitem__(self, index, value): - """Sets slices of a decomposed array. + """Sets slices of a decomposed array using global indexing. The shape of the input object depends on the number of arguments specified, which can be from none to all three. This allows negative indexing, though with ghosts cells included, the first n-ghost negative From 162fc85e7c6bf13581cc956dbde3e8903aa9756e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 27 Sep 2024 16:56:10 +0000 Subject: [PATCH 06/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/amrex/extensions/MultiFab.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index b102be1d..1c9660ae 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -6,12 +6,13 @@ License: BSD-3-Clause-LBNL """ -from .Iterator import next - import numpy as np +from .Iterator import next + try: from mpi4py import MPI as mpi + comm_world = mpi.COMM_WORLD npes = comm_world.Get_size() except ImportError: @@ -232,7 +233,7 @@ def shape(self, include_ghosts=False): min_box = self.box_array().minimal_box() result = min_box.size if include_ghosts: - result = result + self.n_grow_vect*2 + result = result + self.n_grow_vect * 2 result = list(result) + [self.nComp] return tuple(result) @@ -259,7 +260,7 @@ def _get_indices(index, missing): missing: The value used to fill in the extra dimensions added """ - return list(index) + (3 - len(index))*[missing] + return list(index) + (3 - len(index)) * [missing] def _get_min_indices(self, include_ghosts): @@ -370,6 +371,7 @@ def _find_start_stop(self, ii, imin, imax, d, include_ghosts): ) return iistart, iistop + def _get_field(self, mfi, include_ghosts): """Return the field at the given mfi. If include ghosts is true, return the whole array, otherwise @@ -390,20 +392,20 @@ def _get_field(self, mfi, include_ghosts): # self.array(mfi) is in C ordering. # Note: transposing creates a view and not a copy. import inspect + amr = inspect.getmodule(self) if amr.Config.have_gpu: device_arr = self.array(mfi).to_cupy(copy=False, order="F") else: device_arr = self.array(mfi).to_numpy(copy=False, order="F") if not include_ghosts: - device_arr = device_arr[ - tuple([slice(ng, -ng) for ng in self.n_grow_vect]) - ] + device_arr = device_arr[tuple([slice(ng, -ng) for ng in self.n_grow_vect])] return device_arr -def _get_intersect_slice(self, mfi, starts, stops, icstart, icstop, - include_ghosts, with_internal_ghosts): +def _get_intersect_slice( + self, mfi, starts, stops, icstart, icstop, include_ghosts, with_internal_ghosts +): """Return the slices where the block intersects with the global slice. If the block does not intersect, return None. This also shifts the block slices by the number of ghost cells in the From 682dba54d1efa96539b71f432b87782f79c2efdc Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 30 Sep 2024 09:35:17 -0700 Subject: [PATCH 07/27] Remove references to warpx --- src/amrex/extensions/MultiFab.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 1c9660ae..483b7e44 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -597,6 +597,10 @@ def __setitem__(self, index, value): value : scalar or array Input value to assign to the specified slice of the MultiFab """ + import inspect + + amr = inspect.getmodule(self) + # Temporary value until fixed include_ghosts = False # Get the number of dimensions. Is there a cleaner way to do this? @@ -652,11 +656,11 @@ def __setitem__(self, index, value): global_shape[3:3] = [1] value3d.shape = global_shape - if libwarpx.libwarpx_so.Config.have_gpu: + if amr.Config.have_gpu: # check if cupy is available for use xp, cupy_status = load_cupy() if cupy_status is not None: - libwarpx.amr.Print(cupy_status) + amr.Print(cupy_status) starts = [ixstart, iystart, izstart] stops = [ixstop, iystop, izstop] From 1d0b8ca6a2df0d3ddb99dd912c6469fd0fc7348e Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 30 Sep 2024 10:34:39 -0700 Subject: [PATCH 08/27] Updated the documentation --- src/amrex/extensions/MultiFab.py | 43 +++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 483b7e44..deb83e73 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -481,18 +481,24 @@ def _get_intersect_slice( def __getitem__(self, index): - """Returns slice of the MultiFab using global indexing. + """Returns slice of the MultiFab using global indexing, as a numpy array. + This uses numpy array indexing, with the indexing relative to the global array. + The slice ranges can cross multiple blocks and the result will be gathered into a single + array. + + In an MPI context, this is a global operation. An "allgather" is performed so that the full + result is returned on all processors. + The shape of the object returned depends on the number of ix, iy and iz specified, which - can be from none to all three. Note that the values of ix, iy and iz are - relative to the fortran indexing, meaning that 0 is the lower boundary - of the whole domain, and in fortran ordering, i.e. [ix,iy,iz]. - This allows negative indexing, though with ghosts cells included, the first n-ghost negative - indices will refer to the lower guard cells. + can be from none to all three. Note that the values of ix, iy and iz are in fortran ordering, + i.e. [ix,iy,iz], and that 0 is the lower boundary of the whole domain. Parameters ---------- - index : integer, or sequence of integers or slices, or Ellipsis - Index of the slice to return + index : the index using numpy style indexing + Index of the slice to return. + The slice for each dimension can be ":" for the whole range, a slice instance, or an integer. + An Ellipsis can be used to represent the full range of multiple dimensions, as with numpy. """ # Temporary value until fixed include_ghosts = False @@ -584,16 +590,23 @@ def __getitem__(self, index): def __setitem__(self, index, value): - """Sets slices of a decomposed array using global indexing. - The shape of the input object depends on the number of arguments specified, which can - be from none to all three. - This allows negative indexing, though with ghosts cells included, the first n-ghost negative - indices will refer to the lower guard cells. + """Sets the slice of the MultiFab using global indexing. + This uses numpy array indexing, with the indexing relative to the global array. + The slice ranges can cross multiple blocks and the value will be distributed accordingly. + + In an MPI context, this is a local operation. On each processor, the blocks within the slice + range will be set to the value. + + The shape of the value depends on the number of ix, iy and iz specified, which + can be from none to all three. Note that the values of ix, iy and iz are in fortran ordering, + i.e. [ix,iy,iz], and that 0 is the lower boundary of the whole domain. Parameters ---------- - index : integer, or sequence of integers or slices, or Ellipsis - The slice to set + index : the index using numpy style indexing + Index of the slice to return. + The slice for each dimension can be ":" for the whole range, a slice instance, or an integer. + An Ellipsis can be used to represent the full range of multiple dimensions, as with numpy. value : scalar or array Input value to assign to the specified slice of the MultiFab """ From 54b69d44f4c0959418273b7aa5ab0e55753a10f4 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 30 Sep 2024 13:48:33 -0700 Subject: [PATCH 09/27] In __setitem__, always apply the value to both valid and ghost cells --- src/amrex/extensions/MultiFab.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index deb83e73..af852a1b 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -443,18 +443,29 @@ def _get_intersect_slice( box_small_end = box.small_end box_big_end = box.big_end if include_ghosts: - nghosts = self.mf.n_grow_vect + nghosts = self.n_grow_vect box.grow(nghosts) if with_internal_ghosts: box_small_end = box.small_end box_big_end = box.big_end else: + # Only expand the box to include the ghost cells at the edge of the domain min_box = self.box_array().minimal_box() for i in range(len(nghosts)): if box_small_end[i] == min_box.small_end[i]: box_small_end[i] -= nghosts[i] if box_big_end[i] == min_box.big_end[i]: box_big_end[i] += nghosts[i] + else: + if with_internal_ghosts: + # Expand the box to include the ghost cells within the domain + nghosts = self.n_grow_vect + min_box = self.box_array().minimal_box() + for i in range(len(nghosts)): + if box_small_end[i] > min_box.small_end[i]: + box_small_end[i] -= nghosts[i] + if box_big_end[i] < min_box.big_end[i]: + box_big_end[i] += nghosts[i] boxlo = _get_indices(box.small_end, 0) ilo = _get_indices(box_small_end, 0) @@ -593,6 +604,7 @@ def __setitem__(self, index, value): """Sets the slice of the MultiFab using global indexing. This uses numpy array indexing, with the indexing relative to the global array. The slice ranges can cross multiple blocks and the value will be distributed accordingly. + Note that this will apply the value to both valid and ghost cells. In an MPI context, this is a local operation. On each processor, the blocks within the slice range will be set to the value. From eb83663fa63d0279a3f00ec64659a27149e9193d Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 30 Sep 2024 13:51:54 -0700 Subject: [PATCH 10/27] Remove unneeded code --- src/amrex/extensions/MultiFab.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index af852a1b..36f2e3e0 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -622,10 +622,6 @@ def __setitem__(self, index, value): value : scalar or array Input value to assign to the specified slice of the MultiFab """ - import inspect - - amr = inspect.getmodule(self) - # Temporary value until fixed include_ghosts = False # Get the number of dimensions. Is there a cleaner way to do this? @@ -681,12 +677,6 @@ def __setitem__(self, index, value): global_shape[3:3] = [1] value3d.shape = global_shape - if amr.Config.have_gpu: - # check if cupy is available for use - xp, cupy_status = load_cupy() - if cupy_status is not None: - amr.Print(cupy_status) - starts = [ixstart, iystart, izstart] stops = [ixstop, iystop, izstop] for mfi in self: From 88ee98618f1c5ef75add7e669665ecca3a0c0334 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 1 Oct 2024 11:37:41 -0700 Subject: [PATCH 11/27] Fix MPI usage --- src/amrex/extensions/MultiFab.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 36f2e3e0..7c194320 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -10,14 +10,6 @@ from .Iterator import next -try: - from mpi4py import MPI as mpi - - comm_world = mpi.COMM_WORLD - npes = comm_world.Get_size() -except ImportError: - npes = 1 - def mf_to_numpy(self, copy=False, order="F"): """ @@ -569,9 +561,21 @@ def __getitem__(self, index): datalist.append((global_slices, slice_arr)) # Gather the data from all processors + import inspect + + amr = inspect.getmodule(self) + if amr.Config.have_mpi: + npes = amr.ParallelDescriptor.NProcs() + else: + npes = 1 if npes == 1: all_datalist = [datalist] else: + try: + from mpi4py import MPI as mpi + comm_world = mpi.COMM_WORLD + except ImportError: + raise Exception("MultiFab.__getitem__ requires mpi4py") all_datalist = comm_world.allgather(datalist) # Create the array to be returned From c0c5295ed1a4aa8cba7e9d6934a97892add10bfe Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 1 Oct 2024 18:53:17 +0000 Subject: [PATCH 12/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/amrex/extensions/MultiFab.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 7c194320..ddf3d707 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -573,6 +573,7 @@ def __getitem__(self, index): else: try: from mpi4py import MPI as mpi + comm_world = mpi.COMM_WORLD except ImportError: raise Exception("MultiFab.__getitem__ requires mpi4py") From a24e8c41d5a8146dc03e1cf17968f799d848117a Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 14:11:54 -0700 Subject: [PATCH 13/27] Reworked the indexing and how to handle ghost cells --- src/amrex/extensions/MultiFab.py | 356 +++++++++++-------------------- 1 file changed, 124 insertions(+), 232 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index ddf3d707..4ed4bfde 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -256,7 +256,7 @@ def _get_indices(index, missing): def _get_min_indices(self, include_ghosts): - """Returns the minimum indices, expanded to length 3 + """Returns the minimum indices, expanded to length 3+1 Parameters ---------- @@ -268,11 +268,11 @@ def _get_min_indices(self, include_ghosts): min_box = self.box_array().minimal_box() if include_ghosts: min_box.grow(self.n_grow_vect) - return _get_indices(min_box.small_end, 0) + return _get_indices(min_box.small_end, 0) + [0] def _get_max_indices(self, include_ghosts): - """Returns the maximum indices, expanded to length 3. + """Returns the maximum indices, expanded to length 3+1 Parameters ---------- @@ -284,90 +284,88 @@ def _get_max_indices(self, include_ghosts): min_box = self.box_array().minimal_box() if include_ghosts: min_box.grow(self.n_grow_vect) - return _get_indices(min_box.big_end, 0) + return _get_indices(min_box.big_end, 0) + [self.n_comp - 1] -def _fix_index(self, ii, imax, d, include_ghosts): - """Handle negative index, wrapping them as needed. - When ghost cells are included, the indices are - shifted by the number of ghost cells before being wrapped. +def _handle_imaginary_negative_index(index, imin, imax): + """This convects imaginary and negative indices to the actual value""" + if isinstance(index, complex): + assert (index.real == 0.), "Ghost indices must be purely imaginary" + ii = int(index.imag) + if ii <= 0: + result = imin + ii + else: + result = imax + ii + elif index < 0: + result = index + (imax - imin) + 1 + else: + result = index + return result - Parameters - ---------- - self : amrex.MultiFab - A MultiFab class in pyAMReX - ii : integer - The index to be wrapped - imax : integer - The maximum value that the index could have - d : integer - The direction of the index - include_ghosts: bool - Whether or not ghost cells are included - """ - nghosts = list(_get_indices(self.n_grow_vect, 0)) + [0] - if include_ghosts: - ii += nghosts[d] - if ii < 0: - ii += imax - if include_ghosts: - ii -= nghosts[d] - return ii +def _process_index(self, index): + """Convert the input index into a list of slices""" + # Get the number of dimensions. Is there a cleaner way to do this? + dims = len(self.n_grow_vect) -def _find_start_stop(self, ii, imin, imax, d, include_ghosts): - """Given the input index, calculate the start and stop range of the indices. + if index == Ellipsis: + index = [] # This will be filled in below to cover the valid cells + elif isinstance(index, slice) or isinstance(index, int): + # If only one slice or integer passed in, it was not wrapped in a tuple + index = [index] + elif isinstance(index, tuple): + index = list(index) + for i in range(len(index)): + if index[i] == Ellipsis: + index = index[:i] + (dims + 2 - len(index))*[slice(None)] + index[i+1:] + break + else: + raise Exception("MultiFab.__getitem__: unexpected index type") - Parameters - ---------- - ii : None, slice, or integer - Input index, either None, a slice object, or an integer. - Note that ii can be negative. - imin : integer - The global lowest lower bound in the specified direction. - This can include the ghost cells. - imax : integer - The global highest upper bound in the specified direction. - This can include the ghost cells. - This should be the max index + 1. - d : integer - The dimension number, 0, 1, 2, or 3 (3 being the components) - include_ghosts : bool - Whether or not ghost cells are included + if len(index) < dims + 1: + # Add extra dims to index, including for the component. + # These are the dims left out and assumed to extend over the valid cells + index = index + (dims + 1 - len(index))*[slice(None)] + elif len(index) > dims + 1: + raise Exception("Too many indices given") - If ii is a slice, the start and stop values are used directly, - unless they are None, then the lower or upper bound is used. - An assertion checks if the indices are within the bounds. - """ - if ii is None: - iistart = imin - iistop = imax - elif isinstance(ii, slice): - if ii.start is None: - iistart = imin - else: - iistart = _fix_index(self, ii.start, imax, d, include_ghosts) - if ii.stop is None: - iistop = imax + # Expand index to length 3 + 1 + index = _get_indices(index[:-1], 0) + [index[-1]] + + mins = _get_min_indices(self, False) + maxs = _get_max_indices(self, False) + mins_with_ghost = _get_min_indices(self, True) + maxs_with_ghost = _get_max_indices(self, True) + + # Replace all None's in the slices with the bounds of the valid cells, + # handle imaginary indices that specify ghost cells, and adjust negative indices + for i in range(4): + if isinstance(index[i], slice): + if index[i].start is None: + start = mins[i] + else: + start = _handle_imaginary_negative_index(index[i].start, mins[i], maxs[i]) + if index[i].stop is None: + stop = maxs[i] + 1 + else: + stop = _handle_imaginary_negative_index(index[i].stop, mins[i], maxs[i]) + index[i] = slice(start, stop, index[i].step) + elif isinstance(index[i], tuple): + # The slice includes ghosts + assert len(index[i]) == 0, "Indicator to include all ghost cells must be an empty tuple" + index[i] = slice(mins_with_ghost[i], maxs_with_ghost[i] + 1) else: - iistop = _fix_index(self, ii.stop, imax, d, include_ghosts) - else: - ii = _fix_index(self, ii, imax, d, include_ghosts) - iistart = ii - iistop = ii + 1 - assert imin <= iistart <= imax, Exception( - f"Dimension {d+1} lower index is out of bounds" - ) - assert imin <= iistop <= imax, Exception( - f"Dimension {d+1} upper index is out of bounds" - ) - return iistart, iistop + ii = _handle_imaginary_negative_index(index[i], mins[i], maxs[i]) + assert (mins_with_ghost[i] <= ii and ii <= maxs_with_ghost[i]), "Index out of range" + index[i] = slice(ii, ii + 1) + + return index -def _get_field(self, mfi, include_ghosts): +def _get_field(self, mfi): """Return the field at the given mfi. - If include ghosts is true, return the whole array, otherwise - return the interior slice that does not include the ghosts. + If the field is on a GPU, a cupy reference to it is returned, + otherwise a numpy reference. Parameters ---------- @@ -375,8 +373,6 @@ def _get_field(self, mfi, include_ghosts): A MultiFab class in pyAMReX mfi : amrex.MFIiter Index to the FAB of the MultiFab - include_ghosts : bool, default=False - Whether or not ghost cells are included """ # Note that the array will always have 4 dimensions. # even when dims < 3. @@ -390,14 +386,10 @@ def _get_field(self, mfi, include_ghosts): device_arr = self.array(mfi).to_cupy(copy=False, order="F") else: device_arr = self.array(mfi).to_numpy(copy=False, order="F") - if not include_ghosts: - device_arr = device_arr[tuple([slice(ng, -ng) for ng in self.n_grow_vect])] return device_arr -def _get_intersect_slice( - self, mfi, starts, stops, icstart, icstop, include_ghosts, with_internal_ghosts -): +def _get_intersect_slice(self, mfi, index, with_internal_ghosts): """Return the slices where the block intersects with the global slice. If the block does not intersect, return None. This also shifts the block slices by the number of ghost cells in the @@ -405,22 +397,12 @@ def _get_intersect_slice( Parameters ---------- + self : amrex.MultiFab + A MultiFab class in pyAMReX mfi : MFIter The MFIter instance for the current block, - starts : sequence - The minimum indices of the global slice. - These can be negative. - stops : sequence - The maximum indices of the global slice. - These can be negative. - icstart : integer - The minimum component index of the global slice. - These can be negative. - icstops : integer - The maximum component index of the global slice. - These can be negative. - include_ghosts : bool, default=False - Whether or not ghost cells are included + index : sequence + The list indices of the global slice. with_internal_ghosts: bool Whether the internal ghosts are included in the slices @@ -434,30 +416,21 @@ def _get_intersect_slice( box = mfi.tilebox() box_small_end = box.small_end box_big_end = box.big_end - if include_ghosts: - nghosts = self.n_grow_vect - box.grow(nghosts) - if with_internal_ghosts: - box_small_end = box.small_end - box_big_end = box.big_end - else: - # Only expand the box to include the ghost cells at the edge of the domain - min_box = self.box_array().minimal_box() - for i in range(len(nghosts)): - if box_small_end[i] == min_box.small_end[i]: - box_small_end[i] -= nghosts[i] - if box_big_end[i] == min_box.big_end[i]: - box_big_end[i] += nghosts[i] + + # This always include ghost cells since they are included in the MF arrays + nghosts = self.n_grow_vect + box.grow(nghosts) + if with_internal_ghosts: + box_small_end = box.small_end + box_big_end = box.big_end else: - if with_internal_ghosts: - # Expand the box to include the ghost cells within the domain - nghosts = self.n_grow_vect - min_box = self.box_array().minimal_box() - for i in range(len(nghosts)): - if box_small_end[i] > min_box.small_end[i]: - box_small_end[i] -= nghosts[i] - if box_big_end[i] < min_box.big_end[i]: - box_big_end[i] += nghosts[i] + # Only expand the box to include the ghost cells at the edge of the domain + min_box = self.box_array().minimal_box() + for i in range(len(nghosts)): + if box_small_end[i] == min_box.small_end[i]: + box_small_end[i] -= nghosts[i] + if box_big_end[i] == min_box.big_end[i]: + box_big_end[i] += nghosts[i] boxlo = _get_indices(box.small_end, 0) ilo = _get_indices(box_small_end, 0) @@ -465,18 +438,15 @@ def _get_intersect_slice( # Add 1 to the upper end to be consistent with the slicing notation ihi_p1 = [i + 1 for i in ihi] - i1 = np.maximum(starts, ilo) - i2 = np.minimum(stops, ihi_p1) + i1 = np.maximum([index[0].start, index[1].start, index[2].start], ilo) + i2 = np.minimum([index[0].stop, index[1].stop, index[2].stop], ihi_p1) if np.all(i1 < i2): - block_slices = [] - global_slices = [] - for i in range(3): - block_slices.append(slice(i1[i] - boxlo[i], i2[i] - boxlo[i])) - global_slices.append(slice(i1[i] - starts[i], i2[i] - starts[i])) + block_slices = [slice(i1[i] - boxlo[i], i2[i] - boxlo[i]) for i in range(3)] + global_slices = [slice(i1[i] - index[i].start, i2[i] - index[i].start) for i in range(3)] - block_slices.append(slice(icstart, icstop)) - global_slices.append(slice(0, icstop - icstart)) + block_slices.append(index[3]) + global_slices.append(slice(0, index[3].stop - index[3].start)) return tuple(block_slices), tuple(global_slices) else: @@ -492,65 +462,28 @@ def __getitem__(self, index): In an MPI context, this is a global operation. An "allgather" is performed so that the full result is returned on all processors. - The shape of the object returned depends on the number of ix, iy and iz specified, which - can be from none to all three. Note that the values of ix, iy and iz are in fortran ordering, - i.e. [ix,iy,iz], and that 0 is the lower boundary of the whole domain. + Note that the index is in fortran ordering and that 0 is the lower boundary of the whole domain. + + The default range of the indices includes only the valid cells. The ":" index will include all of + the valid cels and no ghost cells. The ghost cells can be accessed using imaginary numbers, with + negative imaginary numbers for the lower ghost cells, and positive for the upper ghost cells. + The index "[-1j]" for example refers to the first lower ghost cell, and "[1j]" to the first upper + ghost cell. To access all cells, ghosts and valid cells, use an empty tuple for the index, i.e. "[()]". Parameters ---------- index : the index using numpy style indexing Index of the slice to return. - The slice for each dimension can be ":" for the whole range, a slice instance, or an integer. - An Ellipsis can be used to represent the full range of multiple dimensions, as with numpy. """ - # Temporary value until fixed - include_ghosts = False - - # Get the number of dimensions. Is there a cleaner way to do this? - dims = len(self.n_grow_vect) - - # Note that the index can have negative values (which wrap around) and has 1 added to the upper - # limit using python style slicing - if index == Ellipsis: - index = dims * [slice(None)] - elif isinstance(index, slice): - # If only one slice passed in, it was not wrapped in a list - index = [index] - - if len(index) < dims + 1: - # Add extra dims to index, including for the component. - # These are the dims left out and assumed to extend over the full size of the dim - index = list(index) - while len(index) < dims + 1: - index.append(slice(None)) - elif len(index) > dims + 1: - raise Exception("Too many indices given") - - # Expand the indices to length 3 - ii = _get_indices(index, None) - ic = index[-1] - - # Global extent. These include the ghost cells when include_ghosts is True - ixmin, iymin, izmin = _get_min_indices(self, include_ghosts) - ixmax, iymax, izmax = _get_max_indices(self, include_ghosts) - - # Setup the size of the array to be returned - ixstart, ixstop = _find_start_stop(self, ii[0], ixmin, ixmax + 1, 0, include_ghosts) - iystart, iystop = _find_start_stop(self, ii[1], iymin, iymax + 1, 1, include_ghosts) - izstart, izstop = _find_start_stop(self, ii[2], izmin, izmax + 1, 2, include_ghosts) - icstart, icstop = _find_start_stop(self, ic, 0, self.n_comp, 3, include_ghosts) + index = _process_index(self, index) # Gather the data to be included in a list to be sent to other processes - starts = [ixstart, iystart, izstart] - stops = [ixstop, iystop, izstop] datalist = [] for mfi in self: - block_slices, global_slices = _get_intersect_slice( - self, mfi, starts, stops, icstart, icstop, include_ghosts, False - ) + block_slices, global_slices = _get_intersect_slice(self, mfi, index, False) if global_slices is not None: # Note that the array will always have 4 dimensions. - device_arr = _get_field(self, mfi, include_ghosts) + device_arr = _get_field(self, mfi) slice_arr = device_arr[block_slices] try: # Copy data from host to device using cupy syntax @@ -580,12 +513,7 @@ def __getitem__(self, index): all_datalist = comm_world.allgather(datalist) # Create the array to be returned - result_shape = ( - max(0, ixstop - ixstart), - max(0, iystop - iystart), - max(0, izstop - izstart), - max(0, icstop - icstart), - ) + result_shape = tuple([max(0, ii.stop - ii.start) for ii in index]) # Now, copy the data into the result array result_global = None @@ -614,54 +542,22 @@ def __setitem__(self, index, value): In an MPI context, this is a local operation. On each processor, the blocks within the slice range will be set to the value. - The shape of the value depends on the number of ix, iy and iz specified, which - can be from none to all three. Note that the values of ix, iy and iz are in fortran ordering, - i.e. [ix,iy,iz], and that 0 is the lower boundary of the whole domain. + Note that the index is in fortran ordering and that 0 is the lower boundary of the whole domain. + + The default range of the indices includes only the valid cells. The ":" index will include all of + the valid cels and no ghost cells. The ghost cells can be accessed using imaginary numbers, with + negative imaginary numbers for the lower ghost cells, and positive for the upper ghost cells. + The index "[-1j]" for example refers to the first lower ghost cell, and "[1j]" to the first upper + ghost cell. To access all cells, ghosts and valid cells, use an empty tuple for the index, i.e. "[()]". Parameters ---------- index : the index using numpy style indexing Index of the slice to return. - The slice for each dimension can be ":" for the whole range, a slice instance, or an integer. - An Ellipsis can be used to represent the full range of multiple dimensions, as with numpy. value : scalar or array Input value to assign to the specified slice of the MultiFab """ - # Temporary value until fixed - include_ghosts = False - # Get the number of dimensions. Is there a cleaner way to do this? - dims = len(self.n_grow_vect) - - # Note that the index can have negative values (which wrap around) and has 1 added to the upper - # limit using python style slicing - if index == Ellipsis: - index = tuple(dims * [slice(None)]) - elif isinstance(index, slice): - # If only one slice passed in, it was not wrapped in a list - index = [index] - - if len(index) < dims + 1: - # Add extra dims to index, including for the component. - # These are the dims left out and assumed to extend over the full size of the dim. - index = list(index) - while len(index) < dims + 1: - index.append(slice(None)) - elif len(index) > dims + 1: - raise Exception("Too many indices given") - - # Expand the indices to length 3 - ii = _get_indices(index, None) - ic = index[-1] - - # Global extent. These include the ghost cells when include_ghosts is True - ixmin, iymin, izmin = _get_min_indices(self, include_ghosts) - ixmax, iymax, izmax = _get_max_indices(self, include_ghosts) - - # Setup the size of the global array to be set - ixstart, ixstop = _find_start_stop(self, ii[0], ixmin, ixmax + 1, 0, include_ghosts) - iystart, iystop = _find_start_stop(self, ii[1], iymin, iymax + 1, 1, include_ghosts) - izstart, izstop = _find_start_stop(self, ii[2], izmin, izmax + 1, 2, include_ghosts) - icstart, icstop = _find_start_stop(self, ic, 0, self.n_comp, 3, include_ghosts) + index = _process_index(self, index) if isinstance(value, np.ndarray): # Expand the shape of the input array to match the shape of the global array @@ -672,24 +568,20 @@ def __setitem__(self, index, value): global_shape = list(value3d.shape) # The shape of 1 is added for the extra dimensions and when index is an integer # (in which case the dimension was not in the input array). - if not isinstance(ii[0], slice): + if (index[0].stop - index[0].start) == 1: global_shape[0:0] = [1] - if not isinstance(ii[1], slice): + if (index[1].stop - index[1].start) == 1: global_shape[1:1] = [1] - if not isinstance(ii[2], slice): + if (index[2].stop - index[2].start) == 1: global_shape[2:2] = [1] - if not isinstance(ic, slice) or len(global_shape) < 4: + if (index[3].stop - index[3].start) == 1 or len(global_shape) < 4: global_shape[3:3] = [1] value3d.shape = global_shape - starts = [ixstart, iystart, izstart] - stops = [ixstop, iystop, izstop] for mfi in self: - block_slices, global_slices = _get_intersect_slice( - self, mfi, starts, stops, icstart, icstop, include_ghosts, True - ) + block_slices, global_slices = _get_intersect_slice(self, mfi, index, True) if global_slices is not None: - mf_arr = _get_field(self, mfi, include_ghosts) + mf_arr = _get_field(self, mfi) if isinstance(value, np.ndarray): # The data is copied from host to device automatically if needed mf_arr[block_slices] = value3d[global_slices] From 810847bf7b4fa45db63c28da7894588c50a928d8 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 16:31:58 -0700 Subject: [PATCH 14/27] Add example use to tests/test_multifab.py --- docs/source/usage/compute.rst | 8 ++++++++ tests/test_multifab.py | 18 ++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/docs/source/usage/compute.rst b/docs/source/usage/compute.rst index 24bfa5bd..5e7fced1 100644 --- a/docs/source/usage/compute.rst +++ b/docs/source/usage/compute.rst @@ -48,6 +48,14 @@ This is how to iterate and potentially compute for all blocks assigned to a loca :start-after: # Manual: Compute Mfab Detailed START :end-before: # Manual: Compute Mfab Detailed END + .. tab-item:: Global + + .. literalinclude:: ../../../tests/test_multifab.py + :language: python3 + :dedent: 4 + :start-after: # Manual: Compute Mfab Global START + :end-before: # Manual: Compute Mfab Global END + For a complete physics example that uses CPU/GPU agnostic Python code for computation on fields, see: diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 46b582bc..ed6e73d8 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -63,6 +63,24 @@ def test_mfab_numpy(mfab): field[()] = 42.0 # Manual: Compute Mfab Simple END + # Manual: Compute Mfab Global START + # finest active MR level, get from a + # simulation's AmrMesh object, e.g.: + # finest_level = sim.finest_level + finest_level = 0 # no MR + + # iterate over mesh-refinement levels + for lev in range(finest_level + 1): + # get an existing MultiFab, e.g., + # from a simulation: + # mfab = sim.get_field(lev=lev) + # Config = sim.extension.Config + + # Using global indexing + mfab[...] = 42.0 + + # Manual: Compute Mfab Global END + @pytest.mark.skipif(amr.Config.have_gpu, reason="This test only runs on CPU") def test_mfab_loop_slow(mfab): From 3e33c6b5de3b82de29a05a139025e44982f7db67 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 3 Oct 2024 23:32:25 +0000 Subject: [PATCH 15/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/amrex/extensions/MultiFab.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 4ed4bfde..e91f9b02 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -290,7 +290,7 @@ def _get_max_indices(self, include_ghosts): def _handle_imaginary_negative_index(index, imin, imax): """This convects imaginary and negative indices to the actual value""" if isinstance(index, complex): - assert (index.real == 0.), "Ghost indices must be purely imaginary" + assert index.real == 0.0, "Ghost indices must be purely imaginary" ii = int(index.imag) if ii <= 0: result = imin + ii @@ -317,7 +317,9 @@ def _process_index(self, index): index = list(index) for i in range(len(index)): if index[i] == Ellipsis: - index = index[:i] + (dims + 2 - len(index))*[slice(None)] + index[i+1:] + index = ( + index[:i] + (dims + 2 - len(index)) * [slice(None)] + index[i + 1 :] + ) break else: raise Exception("MultiFab.__getitem__: unexpected index type") @@ -325,7 +327,7 @@ def _process_index(self, index): if len(index) < dims + 1: # Add extra dims to index, including for the component. # These are the dims left out and assumed to extend over the valid cells - index = index + (dims + 1 - len(index))*[slice(None)] + index = index + (dims + 1 - len(index)) * [slice(None)] elif len(index) > dims + 1: raise Exception("Too many indices given") @@ -344,7 +346,9 @@ def _process_index(self, index): if index[i].start is None: start = mins[i] else: - start = _handle_imaginary_negative_index(index[i].start, mins[i], maxs[i]) + start = _handle_imaginary_negative_index( + index[i].start, mins[i], maxs[i] + ) if index[i].stop is None: stop = maxs[i] + 1 else: @@ -352,11 +356,15 @@ def _process_index(self, index): index[i] = slice(start, stop, index[i].step) elif isinstance(index[i], tuple): # The slice includes ghosts - assert len(index[i]) == 0, "Indicator to include all ghost cells must be an empty tuple" + assert ( + len(index[i]) == 0 + ), "Indicator to include all ghost cells must be an empty tuple" index[i] = slice(mins_with_ghost[i], maxs_with_ghost[i] + 1) else: ii = _handle_imaginary_negative_index(index[i], mins[i], maxs[i]) - assert (mins_with_ghost[i] <= ii and ii <= maxs_with_ghost[i]), "Index out of range" + assert ( + mins_with_ghost[i] <= ii and ii <= maxs_with_ghost[i] + ), "Index out of range" index[i] = slice(ii, ii + 1) return index @@ -443,7 +451,9 @@ def _get_intersect_slice(self, mfi, index, with_internal_ghosts): if np.all(i1 < i2): block_slices = [slice(i1[i] - boxlo[i], i2[i] - boxlo[i]) for i in range(3)] - global_slices = [slice(i1[i] - index[i].start, i2[i] - index[i].start) for i in range(3)] + global_slices = [ + slice(i1[i] - index[i].start, i2[i] - index[i].start) for i in range(3) + ] block_slices.append(index[3]) global_slices.append(slice(0, index[3].stop - index[3].start)) From f088d4e4a8fa7959dd3360c76c52d8d70469b671 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 16:49:27 -0700 Subject: [PATCH 16/27] Add further example to tests/test_multifab.py --- tests/test_multifab.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index ed6e73d8..3b4fd447 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -77,8 +77,18 @@ def test_mfab_numpy(mfab): # Config = sim.extension.Config # Using global indexing + # Set all valid cells mfab[...] = 42.0 + # Set a range of cells. Indices are in Fortran order. + # First dimension, sets from first lower guard cell to first upper guard cell. + # - Imaginary indices refer to the guard cells, negative lower, positive upper. + # Second dimension, sets all valid cells. + # Third dimension, sets all valid and ghost cells + # - The empty tuple is used to specify the range to include all valid and ghost cells. + # Components dimension, sets second component. + mfab[-1j:2j,:,(),2] = np.full((nx+2, ny, nz+2*nghosts), 42.) + # Manual: Compute Mfab Global END From bbe6db22a43f9f7136f7f6c8e2f4602a46356cc6 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 16:56:35 -0700 Subject: [PATCH 17/27] Add example of get in tests/test_multifab.py --- tests/test_multifab.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 3b4fd447..099ac82c 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -89,6 +89,14 @@ def test_mfab_numpy(mfab): # Components dimension, sets second component. mfab[-1j:2j,:,(),2] = np.full((nx+2, ny, nz+2*nghosts), 42.) + # Get a range of cells + # Get the data along the valid cells in the first dimension (gathering data across blocks + # and processors), at the first upper guard cell in the second dimensionn, and cell 16 of + # the third (with 16 being relative to 0 which is the lower end of the full domain). + # Note that in an MPI context, this is a global operation, so caution is required when + # scaling to large numbers of processors. + mfslice = mfab[:,1j,16] + # Manual: Compute Mfab Global END From eb5e6761a004bb8398da86e5f9247d17ac1f4ca9 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 3 Oct 2024 23:56:57 +0000 Subject: [PATCH 18/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_multifab.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 099ac82c..e236766d 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -87,7 +87,7 @@ def test_mfab_numpy(mfab): # Third dimension, sets all valid and ghost cells # - The empty tuple is used to specify the range to include all valid and ghost cells. # Components dimension, sets second component. - mfab[-1j:2j,:,(),2] = np.full((nx+2, ny, nz+2*nghosts), 42.) + mfab[-1j:2j, :, (), 2] = np.full((nx + 2, ny, nz + 2 * nghosts), 42.0) # Get a range of cells # Get the data along the valid cells in the first dimension (gathering data across blocks @@ -95,7 +95,7 @@ def test_mfab_numpy(mfab): # the third (with 16 being relative to 0 which is the lower end of the full domain). # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. - mfslice = mfab[:,1j,16] + mfslice = mfab[:, 1j, 16] # Manual: Compute Mfab Global END From c6e908b05afc28b818406d3f7d2a8dfa9f79b24e Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 17:19:19 -0700 Subject: [PATCH 19/27] Fix tests/test_multifab.py --- tests/test_multifab.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index e236766d..536f7133 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -87,7 +87,7 @@ def test_mfab_numpy(mfab): # Third dimension, sets all valid and ghost cells # - The empty tuple is used to specify the range to include all valid and ghost cells. # Components dimension, sets second component. - mfab[-1j:2j, :, (), 2] = np.full((nx + 2, ny, nz + 2 * nghosts), 42.0) + mfab[-1j:2j, :, (), 2] = 42. # Get a range of cells # Get the data along the valid cells in the first dimension (gathering data across blocks From d8667006bb30cd18b79f59641353245aed63edfa Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 4 Oct 2024 00:20:21 +0000 Subject: [PATCH 20/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_multifab.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 536f7133..0cee6349 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -87,7 +87,7 @@ def test_mfab_numpy(mfab): # Third dimension, sets all valid and ghost cells # - The empty tuple is used to specify the range to include all valid and ghost cells. # Components dimension, sets second component. - mfab[-1j:2j, :, (), 2] = 42. + mfab[-1j:2j, :, (), 2] = 42.0 # Get a range of cells # Get the data along the valid cells in the first dimension (gathering data across blocks From f59a6c11480b84e03a15ae2e55792cc1ca0d15d8 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 17:22:34 -0700 Subject: [PATCH 21/27] Further fix in tests/test_multifab.py --- tests/test_multifab.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 0cee6349..7e264439 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -96,6 +96,7 @@ def test_mfab_numpy(mfab): # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. mfslice = mfab[:, 1j, 16] + mfab[:, 1j, 16] = 2*mfslice # Manual: Compute Mfab Global END From a6d50f0a4f610ab70aaea703b0176e2660432ba7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 4 Oct 2024 00:23:16 +0000 Subject: [PATCH 22/27] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_multifab.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 7e264439..a76919b0 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -96,7 +96,7 @@ def test_mfab_numpy(mfab): # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. mfslice = mfab[:, 1j, 16] - mfab[:, 1j, 16] = 2*mfslice + mfab[:, 1j, 16] = 2 * mfslice # Manual: Compute Mfab Global END From 5d9efbd902dc4f12fad4305fe4c6fb9f6a272da2 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 17:52:20 -0700 Subject: [PATCH 23/27] Fix out of bounds in tests/test_multifab.py --- tests/test_multifab.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index a76919b0..6b29be93 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -86,8 +86,8 @@ def test_mfab_numpy(mfab): # Second dimension, sets all valid cells. # Third dimension, sets all valid and ghost cells # - The empty tuple is used to specify the range to include all valid and ghost cells. - # Components dimension, sets second component. - mfab[-1j:2j, :, (), 2] = 42.0 + # Components dimension, sets first component. + mfab[-1j:2j, :, (), 0] = 42.0 # Get a range of cells # Get the data along the valid cells in the first dimension (gathering data across blocks From 1b2d99346c039fbe480eacea7e7abf4f4ea11b7a Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 18:00:42 -0700 Subject: [PATCH 24/27] More fixes in tests/test_multifab.py --- tests/test_multifab.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 6b29be93..2501943b 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -95,8 +95,8 @@ def test_mfab_numpy(mfab): # the third (with 16 being relative to 0 which is the lower end of the full domain). # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. - mfslice = mfab[:, 1j, 16] - mfab[:, 1j, 16] = 2 * mfslice + mfslice = mfab[:, 1j, 2] + mfab[:, 1j, 2] = 2 * mfslice # Manual: Compute Mfab Global END From d4b569de18a2bb38ea205cc1066d038ffc37c399 Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 3 Oct 2024 18:38:47 -0700 Subject: [PATCH 25/27] Check for CI test --- src/amrex/extensions/MultiFab.py | 1 + tests/test_multifab.py | 3 +++ 2 files changed, 4 insertions(+) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index e91f9b02..3761ceb6 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -511,6 +511,7 @@ def __getitem__(self, index): npes = amr.ParallelDescriptor.NProcs() else: npes = 1 + print(f"\n\n\n\n npes = {npes}\n\n\n") if npes == 1: all_datalist = [datalist] else: diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 2501943b..34d27762 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -95,6 +95,9 @@ def test_mfab_numpy(mfab): # the third (with 16 being relative to 0 which is the lower end of the full domain). # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. + print(f"small end {mfab.box_array().minimal_box().small_end}") + print(f"big end {mfab.box_array().minimal_box().big_end}") + print(f"nghosts {mfab.n_grow_vect}") mfslice = mfab[:, 1j, 2] mfab[:, 1j, 2] = 2 * mfslice From 70c4692b8159843bcb612223aa02e8c44a6231c9 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 4 Oct 2024 09:08:12 -0700 Subject: [PATCH 26/27] More for CI test --- tests/test_multifab.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 34d27762..ce64c350 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -95,9 +95,9 @@ def test_mfab_numpy(mfab): # the third (with 16 being relative to 0 which is the lower end of the full domain). # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. - print(f"small end {mfab.box_array().minimal_box().small_end}") - print(f"big end {mfab.box_array().minimal_box().big_end}") - print(f"nghosts {mfab.n_grow_vect}") + print(f"small end {list(mfab.box_array().minimal_box().small_end)}") + print(f"big end {list(mfab.box_array().minimal_box().big_end)}") + print(f"nghosts {list(mfab.n_grow_vect)}") mfslice = mfab[:, 1j, 2] mfab[:, 1j, 2] = 2 * mfslice From ab25853729e4cc61e611db2b0098f6052007346a Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 4 Oct 2024 09:22:44 -0700 Subject: [PATCH 27/27] Another fix for the CI test --- src/amrex/extensions/MultiFab.py | 1 - tests/test_multifab.py | 9 ++++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/amrex/extensions/MultiFab.py b/src/amrex/extensions/MultiFab.py index 3761ceb6..e91f9b02 100644 --- a/src/amrex/extensions/MultiFab.py +++ b/src/amrex/extensions/MultiFab.py @@ -511,7 +511,6 @@ def __getitem__(self, index): npes = amr.ParallelDescriptor.NProcs() else: npes = 1 - print(f"\n\n\n\n npes = {npes}\n\n\n") if npes == 1: all_datalist = [datalist] else: diff --git a/tests/test_multifab.py b/tests/test_multifab.py index ce64c350..becd777a 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -95,11 +95,10 @@ def test_mfab_numpy(mfab): # the third (with 16 being relative to 0 which is the lower end of the full domain). # Note that in an MPI context, this is a global operation, so caution is required when # scaling to large numbers of processors. - print(f"small end {list(mfab.box_array().minimal_box().small_end)}") - print(f"big end {list(mfab.box_array().minimal_box().big_end)}") - print(f"nghosts {list(mfab.n_grow_vect)}") - mfslice = mfab[:, 1j, 2] - mfab[:, 1j, 2] = 2 * mfslice + if mfab.n_grow_vect.max > 0: + mfslice = mfab[:, 1j, 2] + # The assignment is to the last valid cell of the second dimension. + mfab[:, -1, 2] = 2 * mfslice # Manual: Compute Mfab Global END