From 0af624b7c5bf460793155d65d9c4d86867de364a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 20 Oct 2022 21:39:22 -0700 Subject: [PATCH] Helpers: to_numpy/cupy --- MANIFEST.in | 4 ++ src/amrex/Array4.py | 88 +++++++++++++++++++++++++++++++++++ src/amrex/space1d/__init__.py | 17 +++++++ src/amrex/space2d/__init__.py | 17 +++++++ src/amrex/space3d/__init__.py | 17 +++++++ tests/test_multifab.py | 26 ++++++----- 6 files changed, 157 insertions(+), 12 deletions(-) create mode 100644 src/amrex/Array4.py diff --git a/MANIFEST.in b/MANIFEST.in index 0b0abb03..9a61a6d7 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -6,6 +6,10 @@ recursive-include cmake * recursive-include src * recursive-include tests * +# avoid accidentially copying compiled Python files +global-exclude */__pycache__/* +global-exclude *.pyc + # see .gitignore prune cmake-build* prune .spack-env* diff --git a/src/amrex/Array4.py b/src/amrex/Array4.py new file mode 100644 index 00000000..27299ed1 --- /dev/null +++ b/src/amrex/Array4.py @@ -0,0 +1,88 @@ +""" +This file is part of pyAMReX + +Copyright 2023 AMReX community +Authors: Axel Huebl +License: BSD-3-Clause-LBNL +""" + + +def array4_to_numpy(self, copy=False, order="F"): + """ + Provide a Numpy view into an Array4. + + Note on the order of indices: + By default, this is as in AMReX in Fortran contiguous order, indexing as + x,y,z. This has performance implications for use in external libraries such + as cupy. + The order="C" option will index as z,y,x and perform better with cupy. + https://github.com/AMReX-Codes/pyamrex/issues/55#issuecomment-1579610074 + + Parameters + ---------- + self : amrex.Array4_* + An Array4 class in pyAMReX + copy : bool, optional + Copy the data if true, otherwise create a view (default). + order : string, optional + F order (default) or C. C is faster with external libraries. + + Returns + ------- + np.array + A numpy n-dimensional array. + """ + import numpy as np + + if order == "F": + return np.array(self, copy=copy).T + elif order == "C": + return np.array(self, copy=copy) + else: + raise ValueError("The order argument must be F or C.") + + +def array4_to_cupy(self, copy=False, order="F"): + """ + Provide a Cupy view into an Array4. + + Note on the order of indices: + By default, this is as in AMReX in Fortran contiguous order, indexing as + x,y,z. This has performance implications for use in external libraries such + as cupy. + The order="C" option will index as z,y,x and perform better with cupy. + https://github.com/AMReX-Codes/pyamrex/issues/55#issuecomment-1579610074 + + Parameters + ---------- + self : amrex.Array4_* + An Array4 class in pyAMReX + copy : bool, optional + Copy the data if true, otherwise create a view (default). + order : string, optional + F order (default) or C. C is faster with external libraries. + + Returns + ------- + cupy.array + A numpy n-dimensional array. + + Raises + ------ + ImportError + Raises an exception if cupy is not installed + """ + import cupy as cp + + if order == "F": + return cp.array(self, copy=copy).T + elif order == "C": + return cp.array(self, copy=copy) + else: + raise ValueError("The order argument must be F or C.") + + +def register_Array4_extension(Array4_type): + """Array4 helper methods""" + Array4_type.to_numpy = array4_to_numpy + Array4_type.to_cupy = array4_to_cupy diff --git a/src/amrex/space1d/__init__.py b/src/amrex/space1d/__init__.py index 6e042339..fbdfbf6e 100644 --- a/src/amrex/space1d/__init__.py +++ b/src/amrex/space1d/__init__.py @@ -41,3 +41,20 @@ def Print(*args, **kwargs): print(*args, **kwargs) elif ParallelDescriptor.IOProcessor(): print(*args, **kwargs) + + +from ..Array4 import register_Array4_extension + +register_Array4_extension(Array4_float) +register_Array4_extension(Array4_double) +register_Array4_extension(Array4_longdouble) + +register_Array4_extension(Array4_short) +register_Array4_extension(Array4_int) +register_Array4_extension(Array4_long) +register_Array4_extension(Array4_longlong) + +register_Array4_extension(Array4_ushort) +register_Array4_extension(Array4_uint) +register_Array4_extension(Array4_ulong) +register_Array4_extension(Array4_ulonglong) diff --git a/src/amrex/space2d/__init__.py b/src/amrex/space2d/__init__.py index 235189dd..91d10b03 100644 --- a/src/amrex/space2d/__init__.py +++ b/src/amrex/space2d/__init__.py @@ -41,3 +41,20 @@ def Print(*args, **kwargs): print(*args, **kwargs) elif ParallelDescriptor.IOProcessor(): print(*args, **kwargs) + + +from ..Array4 import register_Array4_extension + +register_Array4_extension(Array4_float) +register_Array4_extension(Array4_double) +register_Array4_extension(Array4_longdouble) + +register_Array4_extension(Array4_short) +register_Array4_extension(Array4_int) +register_Array4_extension(Array4_long) +register_Array4_extension(Array4_longlong) + +register_Array4_extension(Array4_ushort) +register_Array4_extension(Array4_uint) +register_Array4_extension(Array4_ulong) +register_Array4_extension(Array4_ulonglong) diff --git a/src/amrex/space3d/__init__.py b/src/amrex/space3d/__init__.py index 50e7897c..63cdbd61 100644 --- a/src/amrex/space3d/__init__.py +++ b/src/amrex/space3d/__init__.py @@ -41,3 +41,20 @@ def Print(*args, **kwargs): print(*args, **kwargs) elif ParallelDescriptor.IOProcessor(): print(*args, **kwargs) + + +from ..Array4 import register_Array4_extension + +register_Array4_extension(Array4_float) +register_Array4_extension(Array4_double) +register_Array4_extension(Array4_longdouble) + +register_Array4_extension(Array4_short) +register_Array4_extension(Array4_int) +register_Array4_extension(Array4_long) +register_Array4_extension(Array4_longlong) + +register_Array4_extension(Array4_ushort) +register_Array4_extension(Array4_uint) +register_Array4_extension(Array4_ulong) +register_Array4_extension(Array4_ulonglong) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 7ad60c0d..80727615 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -45,18 +45,17 @@ def test_mfab_loop(make_mfab): # numpy representation: non-copying view, including the # guard/ghost region - # note: in numpy, indices are in C-order! - marr_np = np.array(marr, copy=False) + marr_np = marr.to_numpy() # check the values at start/end are the same: first component assert marr_np[0, 0, 0, 0] == marr[bx.small_end] - assert marr_np[0, -1, -1, -1] == marr[bx.big_end] + assert marr_np[-1, -1, -1, 0] == marr[bx.big_end] # same check, but for all components for n in range(mfab.num_comp): small_end_comp = list(bx.small_end) + [n] big_end_comp = list(bx.big_end) + [n] - assert marr_np[n, 0, 0, 0] == marr[small_end_comp] - assert marr_np[n, -1, -1, -1] == marr[big_end_comp] + assert marr_np[0, 0, 0, n] == marr[small_end_comp] + assert marr_np[-1, -1, -1, n] == marr[big_end_comp] # now we do some faster assignments, using range based access # this should fail as out-of-bounds, but does not @@ -64,7 +63,7 @@ def test_mfab_loop(make_mfab): # marr_np[24:200, :, :, :] = 42. # all components and all indices set at once to 42 - marr_np[:, :, :, :] = 42.0 + marr_np[()] = 42.0 # values in start & end still match? assert marr_np[0, 0, 0, 0] == marr[bx.small_end] @@ -210,10 +209,11 @@ def test_mfab_ops_cuda_cupy(make_mfab_device): with cupy.profiler.time_range("assign 3 [()]", color_id=0): for mfi in mfab_device: bx = mfi.tilebox().grow(ngv) - marr = mfab_device.array(mfi) - marr_cupy = cp.array(marr, copy=False) + marr_cupy = mfab_device.array(mfi).to_cupy(order="C") # print(marr_cupy.shape) # 1, 32, 32, 32 # print(marr_cupy.dtype) # float64 + # performance: + # https://github.com/AMReX-Codes/pyamrex/issues/55#issuecomment-1579610074 # write and read into the marr_cupy marr_cupy[()] = 3.0 @@ -244,8 +244,11 @@ def set_to_five(mm): for mfi in mfab_device: bx = mfi.tilebox().grow(ngv) - marr = mfab_device.array(mfi) - marr_cupy = cp.array(marr, copy=False) + marr_cupy = mfab_device.array(mfi).to_cupy(order="F") + # print(marr_cupy.shape) # 32, 32, 32, 1 + # print(marr_cupy.dtype) # float64 + # performance: + # https://github.com/AMReX-Codes/pyamrex/issues/55#issuecomment-1579610074 # write and read into the marr_cupy fives_cp = set_to_five(marr_cupy) @@ -266,8 +269,7 @@ def set_to_seven(x): for mfi in mfab_device: bx = mfi.tilebox().grow(ngv) - marr = mfab_device.array(mfi) - marr_cupy = cp.array(marr, copy=False) + marr_cupy = mfab_device.array(mfi).to_cupy(order="C") # write and read into the marr_cupy set_to_seven(marr_cupy)