From 44dea154311e9f5005db7ab4b8a05965046dbb4e Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 25 May 2023 10:48:07 -0400 Subject: [PATCH 01/11] First sketch of Python wrapper --- .gitignore | 3 ++ python/ffi.py | 112 ++++++++++++++++++++++++++++++++++++++++++++++++ python/pvfmm.py | 39 +++++++++++++++++ 3 files changed, 154 insertions(+) create mode 100644 python/ffi.py create mode 100644 python/pvfmm.py diff --git a/.gitignore b/.gitignore index f1d900f..3fe5ccf 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,6 @@ libtool .project .settings/ +build/ + +*.pyc diff --git a/python/ffi.py b/python/ffi.py new file mode 100644 index 0000000..f42ff4d --- /dev/null +++ b/python/ffi.py @@ -0,0 +1,112 @@ +import ctypes + +import numpy as np +from mpi4py import MPI + +# boilerplate from https://github.com/mpi4py/mpi4py/blob/master/demo/wrap-ctypes/helloworld.py +if MPI._sizeof(MPI.Comm) == ctypes.sizeof(ctypes.c_int): + MPI_Comm = ctypes.c_int +else: + MPI_Comm = ctypes.c_void_p + + +def get_MPI_COMM(comm) -> MPI_Comm: + comm_ptr = MPI._addressof(comm) + return MPI_Comm.from_address(comm_ptr) + +PVFMMKernel = ctypes.c_uint # enum + +# try to load lib +# hardcoded path for now +SHARED_LIB = ctypes.CDLL("./build/libpvfmm.so") + + +# somwhat automatically generated: +# clang2py --clang-args="-I/mnt/sw/nix/store/z5w5a7pr5cmdbds0pn9ajdgy0jg71sl6-gcc-11.3.0/lib/gcc/x86_64-pc-linux-gnu/11.3.0/include/" include/pvfmm.h -l ~/pvfmm/build/libpvfmm.so > python/auto.py + +PVFMMCreateVolumeFMMD = SHARED_LIB.PVFMMCreateVolumeFMMD +PVFMMCreateVolumeFMMD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeFMMD.argtypes = [ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] + +PVFMMCreateVolumeFMMF = SHARED_LIB.PVFMMCreateVolumeFMMF +PVFMMCreateVolumeFMMF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeFMMF.argtypes = [ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] + +PVFMMCreateVolumeTreeD = SHARED_LIB.PVFMMCreateVolumeTreeD +PVFMMCreateVolumeTreeD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeD.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.CFUNCTYPE(None, ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)), ctypes.POINTER(None), ctypes.POINTER(ctypes.c_double), ctypes.c_long, MPI_Comm, ctypes.c_double, ctypes.c_int, ctypes.c_bool, ctypes.c_int] +PVFMMCreateVolumeTreeF = SHARED_LIB.PVFMMCreateVolumeTreeF +PVFMMCreateVolumeTreeF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeF.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.CFUNCTYPE(None, ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)), ctypes.POINTER(None), ctypes.POINTER(ctypes.c_float), ctypes.c_long, MPI_Comm, ctypes.c_float, ctypes.c_int, ctypes.c_bool, ctypes.c_int] +PVFMMCreateVolumeTreeFromCoeffD = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffD +PVFMMCreateVolumeTreeFromCoeffD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeFromCoeffD.argtypes = [ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, MPI_Comm, ctypes.c_bool] +PVFMMCreateVolumeTreeFromCoeffF = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffF +PVFMMCreateVolumeTreeFromCoeffF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeFromCoeffF.argtypes = [ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, MPI_Comm, ctypes.c_bool] +PVFMMEvaluateVolumeFMMD = SHARED_LIB.PVFMMEvaluateVolumeFMMD +PVFMMEvaluateVolumeFMMD.restype = None +PVFMMEvaluateVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long] +PVFMMEvaluateVolumeFMMF = SHARED_LIB.PVFMMEvaluateVolumeFMMF +PVFMMEvaluateVolumeFMMF.restype = None +PVFMMEvaluateVolumeFMMF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long] +PVFMMDestroyVolumeFMMD = SHARED_LIB.PVFMMDestroyVolumeFMMD +PVFMMDestroyVolumeFMMD.restype = None +PVFMMDestroyVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeFMMF = SHARED_LIB.PVFMMDestroyVolumeFMMF +PVFMMDestroyVolumeFMMF.restype = None +PVFMMDestroyVolumeFMMF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeTreeD = SHARED_LIB.PVFMMDestroyVolumeTreeD +PVFMMDestroyVolumeTreeD.restype = None +PVFMMDestroyVolumeTreeD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeTreeF = SHARED_LIB.PVFMMDestroyVolumeTreeF +PVFMMDestroyVolumeTreeF.restype = None +PVFMMDestroyVolumeTreeF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMGetLeafCountD = SHARED_LIB.PVFMMGetLeafCountD +PVFMMGetLeafCountD.restype = ctypes.c_long +PVFMMGetLeafCountD.argtypes = [ctypes.POINTER(None)] +PVFMMGetLeafCountF = SHARED_LIB.PVFMMGetLeafCountF +PVFMMGetLeafCountF.restype = ctypes.c_long +PVFMMGetLeafCountF.argtypes = [ctypes.POINTER(None)] +PVFMMGetLeafCoordD = SHARED_LIB.PVFMMGetLeafCoordD +PVFMMGetLeafCoordD.restype = None +PVFMMGetLeafCoordD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)] +PVFMMGetLeafCoordF = SHARED_LIB.PVFMMGetLeafCoordF +PVFMMGetLeafCoordF.restype = None +PVFMMGetLeafCoordF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)] +PVFMMGetPotentialCoeffD = SHARED_LIB.PVFMMGetPotentialCoeffD +PVFMMGetPotentialCoeffD.restype = None +PVFMMGetPotentialCoeffD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)] +PVFMMGetPotentialCoeffF = SHARED_LIB.PVFMMGetPotentialCoeffF +PVFMMGetPotentialCoeffF.restype = None +PVFMMGetPotentialCoeffF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)] +PVFMMCoeff2NodesD = SHARED_LIB.PVFMMCoeff2NodesD +PVFMMCoeff2NodesD.restype = None +PVFMMCoeff2NodesD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)] +PVFMMCoeff2NodesF = SHARED_LIB.PVFMMCoeff2NodesF +PVFMMCoeff2NodesF.restype = None +PVFMMCoeff2NodesF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)] +PVFMMNodes2CoeffD = SHARED_LIB.PVFMMNodes2CoeffD +PVFMMNodes2CoeffD.restype = None +PVFMMNodes2CoeffD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)] +PVFMMNodes2CoeffF = SHARED_LIB.PVFMMNodes2CoeffF +PVFMMNodes2CoeffF.restype = None +PVFMMNodes2CoeffF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)] +PVFMMCreateContextD = SHARED_LIB.PVFMMCreateContextD +PVFMMCreateContextD.restype = ctypes.POINTER(None) +PVFMMCreateContextD.argtypes = [ctypes.c_double, ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] +PVFMMCreateContextF = SHARED_LIB.PVFMMCreateContextF +PVFMMCreateContextF.restype = ctypes.POINTER(None) +PVFMMCreateContextF.argtypes = [ctypes.c_float, ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] +PVFMMEvalD = SHARED_LIB.PVFMMEvalD +PVFMMEvalD.restype = None +PVFMMEvalD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(None), ctypes.c_int] +PVFMMEvalF = SHARED_LIB.PVFMMEvalF +PVFMMEvalF.restype = None +PVFMMEvalF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(None), ctypes.c_int] +PVFMMDestroyContextD = SHARED_LIB.PVFMMDestroyContextD +PVFMMDestroyContextD.restype = None +PVFMMDestroyContextD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyContextF = SHARED_LIB.PVFMMDestroyContextF +PVFMMDestroyContextF.restype = None +PVFMMDestroyContextF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] diff --git a/python/pvfmm.py b/python/pvfmm.py new file mode 100644 index 0000000..f39ccd9 --- /dev/null +++ b/python/pvfmm.py @@ -0,0 +1,39 @@ +# module load gcc/11.3.0 openmpi python/3.10.8 +# source ~/venvs/pvfmm/bin/activate + +import ctypes +import numpy as np +from mpi4py import MPI +from enum import Enum +import ffi + + +class FMMKernel(Enum): + """Mirroring PVFMMKernel in pvfmm.h""" + + PVFMMLaplacePotential = 0 + PVFMMLaplaceGradient = 1 + PVFMMStokesPressure = 2 + PVFMMStokesVelocity = 3 + PVFMMStokesVelocityGrad = 4 + PVFMMBiotSavartPotential = 5 + +class FMMDoubleVolumeContext: + def __init__( + self, + multipole_order: int, + chebyshev_degree: int, + kernel: FMMKernel, + comm: MPI.Comm, + ): + if multipole_order <= 0 or multipole_order % 2 != 0: + raise ValueError("multipole order must be even and postive") + self._ptr = ffi.PVFMMCreateVolumeFMMD( + multipole_order, chebyshev_degree, int(kernel.value), ffi.get_MPI_COMM(comm) + ) + + def __del__(self): + if hasattr(self, "_ptr"): + ffi.PVFMMDestroyVolumeFMMD(ctypes.pointer(self._ptr)) + +# fmm = FMMContext(10,14, FMMKernel.PVFMMStokesVelocity, MPI.COMM_WORLD) From 32b8f07e5f9b545f9ead258dd488989306847de0 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 25 May 2023 10:56:58 -0400 Subject: [PATCH 02/11] Formatting --- python/ffi.py | 163 +++++++++++++++++++++++++++++++++++++++++++----- python/pvfmm.py | 23 +++++++ 2 files changed, 169 insertions(+), 17 deletions(-) diff --git a/python/ffi.py b/python/ffi.py index f42ff4d..4b5910f 100644 --- a/python/ffi.py +++ b/python/ffi.py @@ -14,7 +14,8 @@ def get_MPI_COMM(comm) -> MPI_Comm: comm_ptr = MPI._addressof(comm) return MPI_Comm.from_address(comm_ptr) -PVFMMKernel = ctypes.c_uint # enum + +PVFMMKernel = ctypes.c_uint # enum # try to load lib # hardcoded path for now @@ -34,22 +35,88 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMCreateVolumeTreeD = SHARED_LIB.PVFMMCreateVolumeTreeD PVFMMCreateVolumeTreeD.restype = ctypes.POINTER(None) -PVFMMCreateVolumeTreeD.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.CFUNCTYPE(None, ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)), ctypes.POINTER(None), ctypes.POINTER(ctypes.c_double), ctypes.c_long, MPI_Comm, ctypes.c_double, ctypes.c_int, ctypes.c_bool, ctypes.c_int] +PVFMMCreateVolumeTreeD.argtypes = [ + ctypes.c_int, + ctypes.c_int, + ctypes.CFUNCTYPE( + None, + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(None), + ), + ctypes.POINTER(None), + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + MPI_Comm, + ctypes.c_double, + ctypes.c_int, + ctypes.c_bool, + ctypes.c_int, +] PVFMMCreateVolumeTreeF = SHARED_LIB.PVFMMCreateVolumeTreeF PVFMMCreateVolumeTreeF.restype = ctypes.POINTER(None) -PVFMMCreateVolumeTreeF.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.CFUNCTYPE(None, ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)), ctypes.POINTER(None), ctypes.POINTER(ctypes.c_float), ctypes.c_long, MPI_Comm, ctypes.c_float, ctypes.c_int, ctypes.c_bool, ctypes.c_int] +PVFMMCreateVolumeTreeF.argtypes = [ + ctypes.c_int, + ctypes.c_int, + ctypes.CFUNCTYPE( + None, + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(None), + ), + ctypes.POINTER(None), + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + MPI_Comm, + ctypes.c_float, + ctypes.c_int, + ctypes.c_bool, + ctypes.c_int, +] PVFMMCreateVolumeTreeFromCoeffD = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffD PVFMMCreateVolumeTreeFromCoeffD.restype = ctypes.POINTER(None) -PVFMMCreateVolumeTreeFromCoeffD.argtypes = [ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, MPI_Comm, ctypes.c_bool] +PVFMMCreateVolumeTreeFromCoeffD.argtypes = [ + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + MPI_Comm, + ctypes.c_bool, +] PVFMMCreateVolumeTreeFromCoeffF = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffF PVFMMCreateVolumeTreeFromCoeffF.restype = ctypes.POINTER(None) -PVFMMCreateVolumeTreeFromCoeffF.argtypes = [ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, MPI_Comm, ctypes.c_bool] +PVFMMCreateVolumeTreeFromCoeffF.argtypes = [ + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + MPI_Comm, + ctypes.c_bool, +] PVFMMEvaluateVolumeFMMD = SHARED_LIB.PVFMMEvaluateVolumeFMMD PVFMMEvaluateVolumeFMMD.restype = None -PVFMMEvaluateVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long] +PVFMMEvaluateVolumeFMMD.argtypes = [ + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(None), + ctypes.POINTER(None), + ctypes.c_long, +] PVFMMEvaluateVolumeFMMF = SHARED_LIB.PVFMMEvaluateVolumeFMMF PVFMMEvaluateVolumeFMMF.restype = None -PVFMMEvaluateVolumeFMMF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long] +PVFMMEvaluateVolumeFMMF.argtypes = [ + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(None), + ctypes.POINTER(None), + ctypes.c_long, +] PVFMMDestroyVolumeFMMD = SHARED_LIB.PVFMMDestroyVolumeFMMD PVFMMDestroyVolumeFMMD.restype = None PVFMMDestroyVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] @@ -76,34 +143,96 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMGetLeafCoordF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)] PVFMMGetPotentialCoeffD = SHARED_LIB.PVFMMGetPotentialCoeffD PVFMMGetPotentialCoeffD.restype = None -PVFMMGetPotentialCoeffD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)] +PVFMMGetPotentialCoeffD.argtypes = [ + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(None), +] PVFMMGetPotentialCoeffF = SHARED_LIB.PVFMMGetPotentialCoeffF PVFMMGetPotentialCoeffF.restype = None -PVFMMGetPotentialCoeffF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)] +PVFMMGetPotentialCoeffF.argtypes = [ + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(None), +] PVFMMCoeff2NodesD = SHARED_LIB.PVFMMCoeff2NodesD PVFMMCoeff2NodesD.restype = None -PVFMMCoeff2NodesD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)] +PVFMMCoeff2NodesD.argtypes = [ + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_double), +] PVFMMCoeff2NodesF = SHARED_LIB.PVFMMCoeff2NodesF PVFMMCoeff2NodesF.restype = None -PVFMMCoeff2NodesF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)] +PVFMMCoeff2NodesF.argtypes = [ + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_float), +] PVFMMNodes2CoeffD = SHARED_LIB.PVFMMNodes2CoeffD PVFMMNodes2CoeffD.restype = None -PVFMMNodes2CoeffD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_double)] +PVFMMNodes2CoeffD.argtypes = [ + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_double), +] PVFMMNodes2CoeffF = SHARED_LIB.PVFMMNodes2CoeffF PVFMMNodes2CoeffF.restype = None -PVFMMNodes2CoeffF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float)] +PVFMMNodes2CoeffF.argtypes = [ + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + ctypes.c_int, + ctypes.c_int, + ctypes.POINTER(ctypes.c_float), +] PVFMMCreateContextD = SHARED_LIB.PVFMMCreateContextD PVFMMCreateContextD.restype = ctypes.POINTER(None) -PVFMMCreateContextD.argtypes = [ctypes.c_double, ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] +PVFMMCreateContextD.argtypes = [ + ctypes.c_double, + ctypes.c_int, + ctypes.c_int, + PVFMMKernel, + MPI_Comm, +] PVFMMCreateContextF = SHARED_LIB.PVFMMCreateContextF PVFMMCreateContextF.restype = ctypes.POINTER(None) -PVFMMCreateContextF.argtypes = [ctypes.c_float, ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] +PVFMMCreateContextF.argtypes = [ + ctypes.c_float, + ctypes.c_int, + ctypes.c_int, + PVFMMKernel, + MPI_Comm, +] PVFMMEvalD = SHARED_LIB.PVFMMEvalD PVFMMEvalD.restype = None -PVFMMEvalD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_long, ctypes.POINTER(None), ctypes.c_int] +PVFMMEvalD.argtypes = [ + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), + ctypes.c_long, + ctypes.POINTER(None), + ctypes.c_int, +] PVFMMEvalF = SHARED_LIB.PVFMMEvalF PVFMMEvalF.restype = None -PVFMMEvalF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_long, ctypes.POINTER(None), ctypes.c_int] +PVFMMEvalF.argtypes = [ + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.c_long, + ctypes.POINTER(None), + ctypes.c_int, +] PVFMMDestroyContextD = SHARED_LIB.PVFMMDestroyContextD PVFMMDestroyContextD.restype = None PVFMMDestroyContextD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] diff --git a/python/pvfmm.py b/python/pvfmm.py index f39ccd9..f384cd7 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -18,6 +18,7 @@ class FMMKernel(Enum): PVFMMStokesVelocityGrad = 4 PVFMMBiotSavartPotential = 5 + class FMMDoubleVolumeContext: def __init__( self, @@ -36,4 +37,26 @@ def __del__(self): if hasattr(self, "_ptr"): ffi.PVFMMDestroyVolumeFMMD(ctypes.pointer(self._ptr)) +# SLOW # fmm = FMMContext(10,14, FMMKernel.PVFMMStokesVelocity, MPI.COMM_WORLD) + +class FFMDoubleParticleContext: + def __init__( + self, + box_size: int, + max_points: int, + multipole_order: int, + kernel: FMMKernel, + comm: MPI.Comm, + ): + if multipole_order <= 0 or multipole_order % 2 != 0: + raise ValueError("multipole order must be even and postive") + self._ptr = ffi.PVFMMCreateContextD(box_size, + multipole_order, max_points, int(kernel.value), ffi.get_MPI_COMM(comm) + ) + + def __del__(self): + if hasattr(self, "_ptr"): + ffi.PVFMMDestroyContextD(ctypes.pointer(self._ptr)) + +# fmm = FFMDoubleParticleContext(-1, 1000,10,FMMKernel.PVFMMBiotSavartPotential, MPI.COMM_WORLD) From a670e133d94956e82780ce9c5584a477abca2bf4 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 31 Jul 2023 17:21:32 -0400 Subject: [PATCH 03/11] Add example based on example-c --- python/example1.py | 92 ++++++++++++++++++++++++++++++++++++++++++++++ python/ffi.py | 83 +++++++++++++++++++++++------------------ python/notes.md | 7 ++++ python/pvfmm.py | 55 ++++++++++++++++++++++----- 4 files changed, 192 insertions(+), 45 deletions(-) create mode 100644 python/example1.py create mode 100644 python/notes.md diff --git a/python/example1.py b/python/example1.py new file mode 100644 index 0000000..6ab0602 --- /dev/null +++ b/python/example1.py @@ -0,0 +1,92 @@ +# module load gcc/11.3.0 openmpi python-mpi/3.10.8 fftw/mpi openblas +# export OMP_NUM_THREADS=16 +# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example1.py + +# based on examples/src/example-c.c + +import time + +from mpi4py import MPI +import numpy as np +import numba + +import pvfmm + +@numba.njit(parallel=True) +def BiotSavart(src_X, src_V, trg_X): + Ns = len(src_X) // 3 + Nt = len(trg_X) // 3 + trg_V = np.zeros(Nt * 3) + oofp = 1 / (4 * np.pi) + for t in numba.prange(Nt): + for s in range(Ns): + X0 = trg_X[t * 3 + 0] - src_X[s * 3 + 0] + X1 = trg_X[t * 3 + 1] - src_X[s * 3 + 1] + X2 = trg_X[t * 3 + 2] - src_X[s * 3 + 2] + r2 = X0 * X0 + X1 * X1 + X2 * X2 + rinv = 1 / np.sqrt(r2) if r2 > 0 else 0 + rinv3 = rinv * rinv * rinv + + trg_V[t * 3 + 0] += ( + (src_V[s * 3 + 1] * X2 - src_V[s * 3 + 2] * X1) * rinv3 * oofp + ) + trg_V[t * 3 + 1] += ( + (src_V[s * 3 + 2] * X0 - src_V[s * 3 + 0] * X2) * rinv3 * oofp + ) + trg_V[t * 3 + 2] += ( + (src_V[s * 3 + 0] * X1 - src_V[s * 3 + 1] * X0) * rinv3 * oofp + ) + return trg_V + + +def test_FMM(ctx): + kdim = (3, 3) + Ns = 20000 + Nt = 20000 + + # Initialize target coordinates + trg_X = np.random.rand(Nt * 3) + + # Initialize source coordinates and density + src_X = np.random.rand(Ns * 3) + src_V = np.random.rand(Ns * kdim[0]) - 0.5 + + # FMM evaluation + setup = 1 + tick = time.perf_counter_ns() + trg_V = ctx.evaluate(src_X, src_V, None, trg_X, setup) + tock = time.perf_counter_ns() + print("FMM evaluation time (with setup) :", (tock - tick) / 1e9) + + # FMM evaluation (without setup) + setup = 0 + tick = time.perf_counter_ns() + trg_V2 = ctx.evaluate(src_X, src_V, None, trg_X, setup) + tock = time.perf_counter_ns() + print("FMM evaluation time (without setup) :", (tock - tick) / 1e9) + + # Direct evaluation + tick = time.perf_counter_ns() + trg_V0 = BiotSavart(src_X, src_V, trg_X) + tock = time.perf_counter_ns() + print("Direct evaluation time :", (tock - tick) / 1e9) + + max_val = np.max(np.abs(trg_V0)) + max_err = np.max(np.abs(trg_V - trg_V0)) + print("Max relative error :", max_err / max_val) + +if __name__ == "__main__": + # MPI init handled by mpi4py import + box_size = -1 + points_per_box = 1000 + multipole_order = 10 + kernel = pvfmm.FMMKernel.PVFMMBiotSavartPotential + + print("Loaded.") + ctx = pvfmm.FFMDoubleParticleContext( + box_size, points_per_box, multipole_order, kernel, MPI.COMM_WORLD + ) + print("Running!") + test_FMM(ctx) + + # MPI finalize handled by mpi4py atexit handler diff --git a/python/ffi.py b/python/ffi.py index 4b5910f..5bc7cb1 100644 --- a/python/ffi.py +++ b/python/ffi.py @@ -1,6 +1,7 @@ import ctypes import numpy as np +from numpy.ctypeslib import ndpointer from mpi4py import MPI # boilerplate from https://github.com/mpi4py/mpi4py/blob/master/demo/wrap-ctypes/helloworld.py @@ -14,6 +15,16 @@ def get_MPI_COMM(comm) -> MPI_Comm: comm_ptr = MPI._addressof(comm) return MPI_Comm.from_address(comm_ptr) +def wrapped_ndptr(*args, **kwargs): + base = ndpointer(*args, **kwargs) + def from_param(cls, obj): + if obj is None: + return obj + return base.from_param(obj) + return type(base.__name__, (base,), {'from_param': classmethod(from_param)}) + +double_array = wrapped_ndptr(dtype=ctypes.c_double, flags=("C_CONTIGUOUS")) +float_array = wrapped_ndptr(dtype=ctypes.c_float, flags=("C_CONTIGUOUS")) PVFMMKernel = ctypes.c_uint # enum @@ -40,13 +51,13 @@ def get_MPI_COMM(comm) -> MPI_Comm: ctypes.c_int, ctypes.CFUNCTYPE( None, - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.c_long, - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.POINTER(None), ), ctypes.POINTER(None), - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.c_long, MPI_Comm, ctypes.c_double, @@ -61,13 +72,13 @@ def get_MPI_COMM(comm) -> MPI_Comm: ctypes.c_int, ctypes.CFUNCTYPE( None, - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.c_long, - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.POINTER(None), ), ctypes.POINTER(None), - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.c_long, MPI_Comm, ctypes.c_float, @@ -81,9 +92,9 @@ def get_MPI_COMM(comm) -> MPI_Comm: ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), + double_array, + double_array, + double_array, ctypes.c_long, MPI_Comm, ctypes.c_bool, @@ -94,9 +105,9 @@ def get_MPI_COMM(comm) -> MPI_Comm: ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), + float_array, + float_array, + float_array, ctypes.c_long, MPI_Comm, ctypes.c_bool, @@ -104,7 +115,7 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMEvaluateVolumeFMMD = SHARED_LIB.PVFMMEvaluateVolumeFMMD PVFMMEvaluateVolumeFMMD.restype = None PVFMMEvaluateVolumeFMMD.argtypes = [ - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long, @@ -112,7 +123,7 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMEvaluateVolumeFMMF = SHARED_LIB.PVFMMEvaluateVolumeFMMF PVFMMEvaluateVolumeFMMF.restype = None PVFMMEvaluateVolumeFMMF.argtypes = [ - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.POINTER(None), ctypes.POINTER(None), ctypes.c_long, @@ -137,57 +148,57 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMGetLeafCountF.argtypes = [ctypes.POINTER(None)] PVFMMGetLeafCoordD = SHARED_LIB.PVFMMGetLeafCoordD PVFMMGetLeafCoordD.restype = None -PVFMMGetLeafCoordD.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(None)] +PVFMMGetLeafCoordD.argtypes = [double_array, ctypes.POINTER(None)] PVFMMGetLeafCoordF = SHARED_LIB.PVFMMGetLeafCoordF PVFMMGetLeafCoordF.restype = None -PVFMMGetLeafCoordF.argtypes = [ctypes.POINTER(ctypes.c_float), ctypes.POINTER(None)] +PVFMMGetLeafCoordF.argtypes = [float_array, ctypes.POINTER(None)] PVFMMGetPotentialCoeffD = SHARED_LIB.PVFMMGetPotentialCoeffD PVFMMGetPotentialCoeffD.restype = None PVFMMGetPotentialCoeffD.argtypes = [ - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.POINTER(None), ] PVFMMGetPotentialCoeffF = SHARED_LIB.PVFMMGetPotentialCoeffF PVFMMGetPotentialCoeffF.restype = None PVFMMGetPotentialCoeffF.argtypes = [ - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.POINTER(None), ] PVFMMCoeff2NodesD = SHARED_LIB.PVFMMCoeff2NodesD PVFMMCoeff2NodesD.restype = None PVFMMCoeff2NodesD.argtypes = [ - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_double), + double_array, ] PVFMMCoeff2NodesF = SHARED_LIB.PVFMMCoeff2NodesF PVFMMCoeff2NodesF.restype = None PVFMMCoeff2NodesF.argtypes = [ - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_float), + float_array, ] PVFMMNodes2CoeffD = SHARED_LIB.PVFMMNodes2CoeffD PVFMMNodes2CoeffD.restype = None PVFMMNodes2CoeffD.argtypes = [ - ctypes.POINTER(ctypes.c_double), + double_array, ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_double), + double_array, ] PVFMMNodes2CoeffF = SHARED_LIB.PVFMMNodes2CoeffF PVFMMNodes2CoeffF.restype = None PVFMMNodes2CoeffF.argtypes = [ - ctypes.POINTER(ctypes.c_float), + float_array, ctypes.c_long, ctypes.c_int, ctypes.c_int, - ctypes.POINTER(ctypes.c_float), + float_array, ] PVFMMCreateContextD = SHARED_LIB.PVFMMCreateContextD PVFMMCreateContextD.restype = ctypes.POINTER(None) @@ -210,12 +221,12 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMEvalD = SHARED_LIB.PVFMMEvalD PVFMMEvalD.restype = None PVFMMEvalD.argtypes = [ - ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), + double_array, + double_array, + double_array, ctypes.c_long, - ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), + double_array, + double_array, ctypes.c_long, ctypes.POINTER(None), ctypes.c_int, @@ -223,12 +234,12 @@ def get_MPI_COMM(comm) -> MPI_Comm: PVFMMEvalF = SHARED_LIB.PVFMMEvalF PVFMMEvalF.restype = None PVFMMEvalF.argtypes = [ - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), + float_array, + float_array, + float_array, ctypes.c_long, - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), + float_array, + float_array, ctypes.c_long, ctypes.POINTER(None), ctypes.c_int, diff --git a/python/notes.md b/python/notes.md new file mode 100644 index 0000000..fe03305 --- /dev/null +++ b/python/notes.md @@ -0,0 +1,7 @@ +- no error handling +- examples (translate): + - [x] example-c.c + - [ ] example2-c.c test1 +- MPI init (not needed - see mpi4py doc: https://mpi4py.readthedocs.io/en/stable/overview.html#initialization-and-exit) +- function pointers for PVFMMCreateVolumeTree functions (numba for speed: https://numba.readthedocs.io/en/stable/user/cfunc.html) +- first run can take a long time, expected. diff --git a/python/pvfmm.py b/python/pvfmm.py index f384cd7..5a0f16b 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -1,10 +1,15 @@ -# module load gcc/11.3.0 openmpi python/3.10.8 +# module load gcc/11.3.0 openmpi python/3.10.8 openblas fftw3 # source ~/venvs/pvfmm/bin/activate import ctypes import numpy as np -from mpi4py import MPI from enum import Enum +from typing import Optional + +# calls MPI_Init_thread() and sets up MPI_Finalize() for you. +# see https://mpi4py.readthedocs.io/en/stable/mpi4py.run.html +from mpi4py import MPI + import ffi @@ -37,13 +42,11 @@ def __del__(self): if hasattr(self, "_ptr"): ffi.PVFMMDestroyVolumeFMMD(ctypes.pointer(self._ptr)) -# SLOW -# fmm = FMMContext(10,14, FMMKernel.PVFMMStokesVelocity, MPI.COMM_WORLD) class FFMDoubleParticleContext: def __init__( self, - box_size: int, + box_size: float, max_points: int, multipole_order: int, kernel: FMMKernel, @@ -51,12 +54,46 @@ def __init__( ): if multipole_order <= 0 or multipole_order % 2 != 0: raise ValueError("multipole order must be even and postive") - self._ptr = ffi.PVFMMCreateContextD(box_size, - multipole_order, max_points, int(kernel.value), ffi.get_MPI_COMM(comm) + self._ptr = ffi.PVFMMCreateContextD( + float(box_size), + max_points, + multipole_order, + int(kernel.value), + ffi.get_MPI_COMM(comm), ) def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyContextD(ctypes.pointer(self._ptr)) + ffi.PVFMMDestroyContextD(ctypes.pointer(ctypes.c_void_p(self._ptr))) -# fmm = FFMDoubleParticleContext(-1, 1000,10,FMMKernel.PVFMMBiotSavartPotential, MPI.COMM_WORLD) + def evaluate( + self, + src_pos: np.ndarray, + sl_den: np.ndarray, + dl_den: Optional[np.ndarray], + trg_pos: np.ndarray, + setup: bool, + ): + source_length = len(src_pos) + if source_length % 3 != 0: + raise ValueError( + "Source arrays must have a length which is a multiple of 3" + ) + if dl_den is not None and len(sl_den) != source_length: + raise ValueError("Source arrays must all be of the same length!") + if dl_den is not None and len(dl_den) != source_length * 2: + raise ValueError("Source arrays must all be of the same length!") + n_src = source_length // 3 + + target_length = len(trg_pos) + if target_length % 3 != 0: + raise ValueError( + "Target arrays must have a length which is a multiple of 3" + ) + n_trg = target_length // 3 + trg_val = np.empty(target_length) + + ffi.PVFMMEvalD( + src_pos, sl_den, dl_den, n_src, trg_pos, trg_val, n_trg, self._ptr, int(setup) + ) + return trg_val From c802659fa5894983e2cc9bfbae7bc190cffb75a7 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 1 Aug 2023 14:43:32 -0400 Subject: [PATCH 04/11] Add example-c to cmake --- CMakeLists.txt | 6 +++--- examples/CMakeLists.txt | 5 +++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b563ed4..68af91a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.10) project( PVFMM VERSION 1.2 - LANGUAGES CXX) + LANGUAGES CXX C) set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) set(CMAKE_CXX_STANDARD 14) @@ -28,7 +28,7 @@ set(MPI_CXX_SKIP_MPICXX true CACHE BOOL "If true, the MPI-2 C++ bindings are disabled using definitions.") - + # required compiler features find_package(MPI REQUIRED) find_package(OpenMP REQUIRED) @@ -44,7 +44,7 @@ else() find_package(FFTW REQUIRED) set(DEP_INC ${FFTW_INCLUDE_DIRS}) set(DEP_LIB ${FFTW_FLOAT_OPENMP_LIB} ${FFTW_FLOAT_LIB} - ${FFTW_DOUBLE_OPENMP_LIB} ${FFTW_DOUBLE_LIB} + ${FFTW_DOUBLE_OPENMP_LIB} ${FFTW_DOUBLE_LIB} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) endif() diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index d5c3626..1363647 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -5,6 +5,11 @@ else() unset(example_DEPS) endif() + +add_executable(example-c ./src/example-c.c) +target_link_libraries(example-c PRIVATE pvfmmStatic "${example_DEPS}") +target_include_directories(example-c PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) + add_executable(example1 ./src/example1.cpp) target_link_libraries(example1 PRIVATE pvfmmStatic "${example_DEPS}") target_include_directories(example1 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) From e70979709495ee326b3d069e9345b259ecbf9d4f Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 1 Aug 2023 14:43:53 -0400 Subject: [PATCH 05/11] Have setup be true by default --- python/example1.py | 6 ++---- python/pvfmm.py | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/python/example1.py b/python/example1.py index 6ab0602..edb491d 100644 --- a/python/example1.py +++ b/python/example1.py @@ -52,16 +52,14 @@ def test_FMM(ctx): src_V = np.random.rand(Ns * kdim[0]) - 0.5 # FMM evaluation - setup = 1 tick = time.perf_counter_ns() - trg_V = ctx.evaluate(src_X, src_V, None, trg_X, setup) + trg_V = ctx.evaluate(src_X, src_V, None, trg_X) tock = time.perf_counter_ns() print("FMM evaluation time (with setup) :", (tock - tick) / 1e9) # FMM evaluation (without setup) - setup = 0 tick = time.perf_counter_ns() - trg_V2 = ctx.evaluate(src_X, src_V, None, trg_X, setup) + trg_V2 = ctx.evaluate(src_X, src_V, None, trg_X, setup=False) tock = time.perf_counter_ns() print("FMM evaluation time (without setup) :", (tock - tick) / 1e9) diff --git a/python/pvfmm.py b/python/pvfmm.py index 5a0f16b..3dd00a4 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -72,7 +72,7 @@ def evaluate( sl_den: np.ndarray, dl_den: Optional[np.ndarray], trg_pos: np.ndarray, - setup: bool, + setup: bool = True, ): source_length = len(src_pos) if source_length % 3 != 0: From 1122089ba59f6cceacf3b7a4cf03995a141fd6ea Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 2 Aug 2023 15:57:27 -0400 Subject: [PATCH 06/11] Add example2. Lib needs work on output sizes in general case --- examples/CMakeLists.txt | 5 ++ examples/src/example2-c.c | 2 +- python/example1.py | 4 +- python/example2.py | 168 ++++++++++++++++++++++++++++++++++++++ python/ffi.py | 98 +++++++++++----------- python/pvfmm.py | 159 ++++++++++++++++++++++++++++++++++-- 6 files changed, 377 insertions(+), 59 deletions(-) create mode 100644 python/example2.py diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 1363647..81bdb20 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -10,6 +10,11 @@ add_executable(example-c ./src/example-c.c) target_link_libraries(example-c PRIVATE pvfmmStatic "${example_DEPS}") target_include_directories(example-c PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) +add_executable(example2-c ./src/example2-c.c) +target_link_libraries(example2-c PRIVATE pvfmmStatic "${example_DEPS}") +target_include_directories(example2-c PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) + + add_executable(example1 ./src/example1.cpp) target_link_libraries(example1 PRIVATE pvfmmStatic "${example_DEPS}") target_include_directories(example1 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) diff --git a/examples/src/example2-c.c b/examples/src/example2-c.c index e33e2ce..75dbfab 100644 --- a/examples/src/example2-c.c +++ b/examples/src/example2-c.c @@ -169,7 +169,7 @@ int main(int argc, char** argv) { // Build FMM translation operators void* fmm = PVFMMCreateVolumeFMMD(mult_order, cheb_deg, PVFMMStokesVelocity, comm); - //test1(fmm, kdim0, kdim1, cheb_deg, comm); + test1(fmm, kdim0, kdim1, cheb_deg, comm); test2(fmm, kdim0, kdim1, cheb_deg, comm); PVFMMDestroyVolumeFMMD(&fmm); diff --git a/python/example1.py b/python/example1.py index edb491d..a2bfc4e 100644 --- a/python/example1.py +++ b/python/example1.py @@ -12,6 +12,7 @@ import pvfmm + @numba.njit(parallel=True) def BiotSavart(src_X, src_V, trg_X): Ns = len(src_X) // 3 @@ -73,12 +74,13 @@ def test_FMM(ctx): max_err = np.max(np.abs(trg_V - trg_V0)) print("Max relative error :", max_err / max_val) + if __name__ == "__main__": # MPI init handled by mpi4py import box_size = -1 points_per_box = 1000 multipole_order = 10 - kernel = pvfmm.FMMKernel.PVFMMBiotSavartPotential + kernel = pvfmm.FMMKernel.BiotSavartPotential print("Loaded.") ctx = pvfmm.FFMDoubleParticleContext( diff --git a/python/example2.py b/python/example2.py new file mode 100644 index 0000000..5297cbe --- /dev/null +++ b/python/example2.py @@ -0,0 +1,168 @@ +# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example1.py + +# based on examples/src/example2-c.c + +import numpy as np +import numba +from numba import types +import ctypes +from mpi4py import MPI + +import pvfmm + +c_sig = types.void( + types.CPointer(types.double), + types.int64, + types.CPointer(types.double), + types.voidptr, +) + + +# see https://numba.readthedocs.io/en/stable/user/cfunc.html +@numba.cfunc(c_sig, nopython=True) +def fn_input(coord_, n, out_, _ctx): + coord = numba.carray(coord_, n * 3) + out = numba.carray(out_, n * 3) + L = 125 + for i in range(n): + idx = i * 3 + c = coord[idx : idx + 3] + r_2 = (c[0] - 0.5) ** 2 + (c[1] - 0.5) ** 2 + (c[2] - 0.5) ** 2 + out[idx + 0] = 0 + 2 * L * np.exp(-L * r_2) * (c[0] - 0.5) + out[idx + 1] = 4 * L**2 * (c[2] - 0.5) * (5 - 2 * L * r_2) * np.exp( + -L * r_2 + ) + 2 * L * np.exp(-L * r_2) * (c[1] - 0.5) + out[idx + 2] = -4 * L**2 * (c[1] - 0.5) * (5 - 2 * L * r_2) * np.exp( + -L * r_2 + ) + 2 * L * np.exp(-L * r_2) * (c[2] - 0.5) + + +@numba.njit +def fn_poten(coord): + n = len(coord) // 3 + dof = 3 + out = np.zeros(n * dof, dtype=np.float64) + L = 125 + + for i in range(n): + idx = i * dof + c = coord[idx : idx + dof] + r_2 = (c[0] - 0.5) ** 2 + (c[1] - 0.5) ** 2 + (c[2] - 0.5) ** 2 + out[idx + 0] = 0 + out[idx + 1] = 2 * L * (c[2] - 0.5) * np.exp(-L * r_2) + out[idx + 2] = -2 * L * (c[1] - 0.5) * np.exp(-L * r_2) + + return out + + +@numba.njit +def GetChebNodes(Nleaf, cheb_deg, depth, leaf_coord): + leaf_length = 1.0 / (1 << depth) + Ncheb = (cheb_deg + 1) ** 3 + cheb_coord = np.zeros(Nleaf * Ncheb * 3, dtype=np.float64) + + for leaf_idx in range(Nleaf): + for j2 in range(cheb_deg + 1): + for j1 in range(cheb_deg + 1): + for j0 in range(cheb_deg + 1): + node_idx = ( + leaf_idx * Ncheb + + (j2 * (cheb_deg + 1) + j1) * (cheb_deg + 1) + + j0 + ) + cheb_coord[node_idx * 3 + 0] = ( + leaf_coord[leaf_idx * 3 + 0] + + (1 - np.cos(np.pi * (j0 * 2 + 1) / (cheb_deg * 2 + 2))) + * leaf_length + * 0.5 + ) + cheb_coord[node_idx * 3 + 1] = ( + leaf_coord[leaf_idx * 3 + 1] + + (1 - np.cos(np.pi * (j1 * 2 + 1) / (cheb_deg * 2 + 2))) + * leaf_length + * 0.5 + ) + cheb_coord[node_idx * 3 + 2] = ( + leaf_coord[leaf_idx * 3 + 2] + + (1 - np.cos(np.pi * (j2 * 2 + 1) / (cheb_deg * 2 + 2))) + * leaf_length + * 0.5 + ) + + return cheb_coord + + +def test1(fmm, kdim0, kdim1, cheb_deg, comm): + Nt = 100 + trg_coord = np.random.rand(Nt * 3) + trg_value_ref = fn_poten(trg_coord) + tree = pvfmm.FMMDoubleVolumeTree.from_function( + cheb_deg, kdim0, fn_input.ctypes, None, trg_coord, comm, 1e-6, 100, False, 0 + ) + trg_value = tree.evaluate(fmm, Nt) + + max_val = np.max(np.abs(trg_value_ref)) + max_err = np.max(np.abs(trg_value - trg_value_ref)) + print("Maximum relative error = ", max_err / max_val) + + +def test2(fmm, kdim0, kdim1, cheb_deg, comm): + Ncheb = (cheb_deg + 1) * (cheb_deg + 1) * (cheb_deg + 1) + Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 + + depth = 3 + # skipping MPI (for now?) + + Nleaf = 1 << (3 * depth) + leaf_length = 1 / (1 << depth) + + leaf_coord = np.empty(Nleaf * 3, dtype=np.float64) + dens_value = np.empty(Nleaf * Ncheb * kdim0, dtype=np.float64) + + for leaf_idx in range(Nleaf): + leaf_coord[leaf_idx * 3 + 0] = ( + (leaf_idx // (1 << (depth * 0))) % (1 << depth) * leaf_length + ) + leaf_coord[leaf_idx * 3 + 1] = ( + (leaf_idx // (1 << (depth * 1))) % (1 << depth) * leaf_length + ) + leaf_coord[leaf_idx * 3 + 2] = ( + (leaf_idx // (1 << (depth * 2))) % (1 << depth) * leaf_length + ) + print("Getting nodes") + cheb_coord = GetChebNodes(Nleaf, cheb_deg, depth, leaf_coord) + fn_input( # ugly just because we used numba above + cheb_coord.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), + Nleaf * Ncheb, + dens_value.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), + None, + ) + dense_coeff = pvfmm.nodes_to_coeff(Nleaf, cheb_deg, kdim0, dens_value) + tree = pvfmm.FMMDoubleVolumeTree.from_coefficients( + cheb_deg, kdim0, leaf_coord, dense_coeff, None, comm, False + ) + tree.evaluate(fmm, 0) + potn_value = tree.get_values() + + Nleaf = tree.leaf_count() # does this change from above? + leaf_coord = tree.get_leaf_coordinates() + cheb_coord = GetChebNodes(Nleaf, cheb_deg, depth, leaf_coord) + potn_value_ref = fn_poten(cheb_coord) + + max_val = np.max(np.abs(potn_value_ref)) + max_err = np.max(np.abs(potn_value - potn_value_ref)) + print("Maximum relative error = ", max_err / max_val) + + +if __name__ == "__main__": + mult_order = 10 + cheb_deg = 14 + kdim0 = 3 + kdim1 = 3 + + comm = MPI.COMM_WORLD + fmm = pvfmm.FMMDoubleVolumeContext( + mult_order, cheb_deg, pvfmm.FMMKernel.StokesVelocity, comm + ) + test1(fmm, kdim0, kdim1, cheb_deg, comm) + test2(fmm, kdim0, kdim1, cheb_deg, comm) diff --git a/python/ffi.py b/python/ffi.py index 5bc7cb1..db69826 100644 --- a/python/ffi.py +++ b/python/ffi.py @@ -15,17 +15,31 @@ def get_MPI_COMM(comm) -> MPI_Comm: comm_ptr = MPI._addressof(comm) return MPI_Comm.from_address(comm_ptr) + def wrapped_ndptr(*args, **kwargs): - base = ndpointer(*args, **kwargs) - def from_param(cls, obj): - if obj is None: - return obj - return base.from_param(obj) - return type(base.__name__, (base,), {'from_param': classmethod(from_param)}) + base = ndpointer(*args, **kwargs) + + def from_param(cls, obj): + if obj is None: + return obj + return base.from_param(obj) + + return type(base.__name__, (base,), {"from_param": classmethod(from_param)}) + double_array = wrapped_ndptr(dtype=ctypes.c_double, flags=("C_CONTIGUOUS")) float_array = wrapped_ndptr(dtype=ctypes.c_float, flags=("C_CONTIGUOUS")) + +def volume_callback(type): + return ctypes.CFUNCTYPE( + None, ctypes.POINTER(type), ctypes.c_long, ctypes.POINTER(type), ctypes.c_void_p + ) + + +double_volume_callback = volume_callback(ctypes.c_double) +float_volume_callback = volume_callback(ctypes.c_float) + PVFMMKernel = ctypes.c_uint # enum # try to load lib @@ -37,26 +51,20 @@ def from_param(cls, obj): # clang2py --clang-args="-I/mnt/sw/nix/store/z5w5a7pr5cmdbds0pn9ajdgy0jg71sl6-gcc-11.3.0/lib/gcc/x86_64-pc-linux-gnu/11.3.0/include/" include/pvfmm.h -l ~/pvfmm/build/libpvfmm.so > python/auto.py PVFMMCreateVolumeFMMD = SHARED_LIB.PVFMMCreateVolumeFMMD -PVFMMCreateVolumeFMMD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeFMMD.restype = ctypes.c_void_p PVFMMCreateVolumeFMMD.argtypes = [ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] PVFMMCreateVolumeFMMF = SHARED_LIB.PVFMMCreateVolumeFMMF -PVFMMCreateVolumeFMMF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeFMMF.restype = ctypes.c_void_p PVFMMCreateVolumeFMMF.argtypes = [ctypes.c_int, ctypes.c_int, PVFMMKernel, MPI_Comm] PVFMMCreateVolumeTreeD = SHARED_LIB.PVFMMCreateVolumeTreeD -PVFMMCreateVolumeTreeD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeD.restype = ctypes.c_void_p PVFMMCreateVolumeTreeD.argtypes = [ ctypes.c_int, ctypes.c_int, - ctypes.CFUNCTYPE( - None, - double_array, - ctypes.c_long, - double_array, - ctypes.POINTER(None), - ), - ctypes.POINTER(None), + double_volume_callback, + ctypes.c_void_p, double_array, ctypes.c_long, MPI_Comm, @@ -66,18 +74,12 @@ def from_param(cls, obj): ctypes.c_int, ] PVFMMCreateVolumeTreeF = SHARED_LIB.PVFMMCreateVolumeTreeF -PVFMMCreateVolumeTreeF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeF.restype = ctypes.c_void_p PVFMMCreateVolumeTreeF.argtypes = [ ctypes.c_int, ctypes.c_int, - ctypes.CFUNCTYPE( - None, - float_array, - ctypes.c_long, - float_array, - ctypes.POINTER(None), - ), - ctypes.POINTER(None), + float_volume_callback, + ctypes.c_void_p, float_array, ctypes.c_long, MPI_Comm, @@ -87,7 +89,7 @@ def from_param(cls, obj): ctypes.c_int, ] PVFMMCreateVolumeTreeFromCoeffD = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffD -PVFMMCreateVolumeTreeFromCoeffD.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeFromCoeffD.restype = ctypes.c_void_p PVFMMCreateVolumeTreeFromCoeffD.argtypes = [ ctypes.c_long, ctypes.c_int, @@ -100,7 +102,7 @@ def from_param(cls, obj): ctypes.c_bool, ] PVFMMCreateVolumeTreeFromCoeffF = SHARED_LIB.PVFMMCreateVolumeTreeFromCoeffF -PVFMMCreateVolumeTreeFromCoeffF.restype = ctypes.POINTER(None) +PVFMMCreateVolumeTreeFromCoeffF.restype = ctypes.c_void_p PVFMMCreateVolumeTreeFromCoeffF.argtypes = [ ctypes.c_long, ctypes.c_int, @@ -116,53 +118,53 @@ def from_param(cls, obj): PVFMMEvaluateVolumeFMMD.restype = None PVFMMEvaluateVolumeFMMD.argtypes = [ double_array, - ctypes.POINTER(None), - ctypes.POINTER(None), + ctypes.c_void_p, + ctypes.c_void_p, ctypes.c_long, ] PVFMMEvaluateVolumeFMMF = SHARED_LIB.PVFMMEvaluateVolumeFMMF PVFMMEvaluateVolumeFMMF.restype = None PVFMMEvaluateVolumeFMMF.argtypes = [ float_array, - ctypes.POINTER(None), - ctypes.POINTER(None), + ctypes.c_void_p, + ctypes.c_void_p, ctypes.c_long, ] PVFMMDestroyVolumeFMMD = SHARED_LIB.PVFMMDestroyVolumeFMMD PVFMMDestroyVolumeFMMD.restype = None -PVFMMDestroyVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeFMMD.argtypes = [ctypes.POINTER(ctypes.c_void_p)] PVFMMDestroyVolumeFMMF = SHARED_LIB.PVFMMDestroyVolumeFMMF PVFMMDestroyVolumeFMMF.restype = None -PVFMMDestroyVolumeFMMF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeFMMF.argtypes = [ctypes.POINTER(ctypes.c_void_p)] PVFMMDestroyVolumeTreeD = SHARED_LIB.PVFMMDestroyVolumeTreeD PVFMMDestroyVolumeTreeD.restype = None -PVFMMDestroyVolumeTreeD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeTreeD.argtypes = [ctypes.POINTER(ctypes.c_void_p)] PVFMMDestroyVolumeTreeF = SHARED_LIB.PVFMMDestroyVolumeTreeF PVFMMDestroyVolumeTreeF.restype = None -PVFMMDestroyVolumeTreeF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyVolumeTreeF.argtypes = [ctypes.POINTER(ctypes.c_void_p)] PVFMMGetLeafCountD = SHARED_LIB.PVFMMGetLeafCountD PVFMMGetLeafCountD.restype = ctypes.c_long -PVFMMGetLeafCountD.argtypes = [ctypes.POINTER(None)] +PVFMMGetLeafCountD.argtypes = [ctypes.c_void_p] PVFMMGetLeafCountF = SHARED_LIB.PVFMMGetLeafCountF PVFMMGetLeafCountF.restype = ctypes.c_long -PVFMMGetLeafCountF.argtypes = [ctypes.POINTER(None)] +PVFMMGetLeafCountF.argtypes = [ctypes.c_void_p] PVFMMGetLeafCoordD = SHARED_LIB.PVFMMGetLeafCoordD PVFMMGetLeafCoordD.restype = None -PVFMMGetLeafCoordD.argtypes = [double_array, ctypes.POINTER(None)] +PVFMMGetLeafCoordD.argtypes = [double_array, ctypes.c_void_p] PVFMMGetLeafCoordF = SHARED_LIB.PVFMMGetLeafCoordF PVFMMGetLeafCoordF.restype = None -PVFMMGetLeafCoordF.argtypes = [float_array, ctypes.POINTER(None)] +PVFMMGetLeafCoordF.argtypes = [float_array, ctypes.c_void_p] PVFMMGetPotentialCoeffD = SHARED_LIB.PVFMMGetPotentialCoeffD PVFMMGetPotentialCoeffD.restype = None PVFMMGetPotentialCoeffD.argtypes = [ double_array, - ctypes.POINTER(None), + ctypes.c_void_p, ] PVFMMGetPotentialCoeffF = SHARED_LIB.PVFMMGetPotentialCoeffF PVFMMGetPotentialCoeffF.restype = None PVFMMGetPotentialCoeffF.argtypes = [ float_array, - ctypes.POINTER(None), + ctypes.c_void_p, ] PVFMMCoeff2NodesD = SHARED_LIB.PVFMMCoeff2NodesD PVFMMCoeff2NodesD.restype = None @@ -201,7 +203,7 @@ def from_param(cls, obj): float_array, ] PVFMMCreateContextD = SHARED_LIB.PVFMMCreateContextD -PVFMMCreateContextD.restype = ctypes.POINTER(None) +PVFMMCreateContextD.restype = ctypes.c_void_p PVFMMCreateContextD.argtypes = [ ctypes.c_double, ctypes.c_int, @@ -210,7 +212,7 @@ def from_param(cls, obj): MPI_Comm, ] PVFMMCreateContextF = SHARED_LIB.PVFMMCreateContextF -PVFMMCreateContextF.restype = ctypes.POINTER(None) +PVFMMCreateContextF.restype = ctypes.c_void_p PVFMMCreateContextF.argtypes = [ ctypes.c_float, ctypes.c_int, @@ -228,7 +230,7 @@ def from_param(cls, obj): double_array, double_array, ctypes.c_long, - ctypes.POINTER(None), + ctypes.c_void_p, ctypes.c_int, ] PVFMMEvalF = SHARED_LIB.PVFMMEvalF @@ -241,12 +243,12 @@ def from_param(cls, obj): float_array, float_array, ctypes.c_long, - ctypes.POINTER(None), + ctypes.c_void_p, ctypes.c_int, ] PVFMMDestroyContextD = SHARED_LIB.PVFMMDestroyContextD PVFMMDestroyContextD.restype = None -PVFMMDestroyContextD.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyContextD.argtypes = [ctypes.POINTER(ctypes.c_void_p)] PVFMMDestroyContextF = SHARED_LIB.PVFMMDestroyContextF PVFMMDestroyContextF.restype = None -PVFMMDestroyContextF.argtypes = [ctypes.POINTER(ctypes.POINTER(None))] +PVFMMDestroyContextF.argtypes = [ctypes.POINTER(ctypes.c_void_p)] diff --git a/python/pvfmm.py b/python/pvfmm.py index 3dd00a4..857200e 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -13,15 +13,28 @@ import ffi +def nodes_to_coeff(N_leaf: int, cheb_deg: int, dof: int, node_val: np.ndarray): + is_double = node_val.dtype == np.float64 + Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 + # XXX: is this valid? + coeff = np.empty(Ncoef * N_leaf * dof, dtype=node_val.dtype) + if is_double: + ffi.PVFMMNodes2CoeffD(coeff, N_leaf, cheb_deg, dof, node_val) + else: + ffi.PVFMMNodes2CoeffF(coeff, N_leaf, cheb_deg, dof, node_val) + + return coeff + + class FMMKernel(Enum): """Mirroring PVFMMKernel in pvfmm.h""" - PVFMMLaplacePotential = 0 - PVFMMLaplaceGradient = 1 - PVFMMStokesPressure = 2 - PVFMMStokesVelocity = 3 - PVFMMStokesVelocityGrad = 4 - PVFMMBiotSavartPotential = 5 + LaplacePotential = 0 + LaplaceGradient = 1 + StokesPressure = 2 + StokesVelocity = 3 + StokesVelocityGrad = 4 + BiotSavartPotential = 5 class FMMDoubleVolumeContext: @@ -40,7 +53,7 @@ def __init__( def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyVolumeFMMD(ctypes.pointer(self._ptr)) + ffi.PVFMMDestroyVolumeFMMD(ctypes.byref(ctypes.c_void_p(self._ptr))) class FFMDoubleParticleContext: @@ -64,7 +77,7 @@ def __init__( def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyContextD(ctypes.pointer(ctypes.c_void_p(self._ptr))) + ffi.PVFMMDestroyContextD(ctypes.byref(ctypes.c_void_p(self._ptr))) def evaluate( self, @@ -94,6 +107,134 @@ def evaluate( trg_val = np.empty(target_length) ffi.PVFMMEvalD( - src_pos, sl_den, dl_den, n_src, trg_pos, trg_val, n_trg, self._ptr, int(setup) + src_pos, + sl_den, + dl_den, + n_src, + trg_pos, + trg_val, + n_trg, + self._ptr, + int(setup), + ) + return trg_val + + +class FMMDoubleVolumeTree: + def __init__(self, ptr: ctypes.c_void_p, cheb_deg: int, data_dim: int, n_trg: int): + self._ptr = ptr + self.cheb_deg = cheb_deg + self.n_cheb = (cheb_deg + 1) ** 3 + self.n_coeff = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 + self.data_dim = data_dim + self.n_trg = n_trg + + @classmethod + def from_function( + cls, + cheb_deg: int, + data_dim: int, + fn: ffi.double_volume_callback, + context: ctypes.c_void_p, # TODO: investigate + trg_coord: np.ndarray, + comm: MPI.Comm, + tol: float, + max_pts: int, + periodic: bool, + init_depth: int, + ): + n_trg = len(trg_coord) // 3 + ptr = ffi.PVFMMCreateVolumeTreeD( + cheb_deg, + data_dim, + fn, + context, + trg_coord, + n_trg, + ffi.get_MPI_COMM(comm), + tol, + max_pts, + periodic, + init_depth, + ) + return cls(ptr, cheb_deg, data_dim, n_trg) + + @classmethod + def from_coefficients( + cls, + cheb_deg: int, + data_dim: int, + leaf_coord: np.ndarray, + fn_coeff: np.ndarray, + trg_coord: Optional[np.ndarray], + comm: MPI.Comm, + periodic: bool, + ): + if len(leaf_coord) % 3 != 0: + raise ValueError( + "Leaf coordinates must have a length which is a multiple of 3" + ) + N_leaf = len(leaf_coord) // 3 + fn_coeff_size = ( + N_leaf * data_dim * (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 + ) + if len(fn_coeff) != fn_coeff_size: + raise ValueError( + "Function coefficients array has the wrong length, required length " + + fn_coeff_size + ) + n_trg = len(trg_coord) // 3 if trg_coord is not None else 0 + + ptr = ffi.PVFMMCreateVolumeTreeFromCoeffD( + N_leaf, + cheb_deg, + data_dim, + leaf_coord, + fn_coeff, + trg_coord, + n_trg, + ffi.get_MPI_COMM(comm), + periodic, ) + return cls(ptr, cheb_deg, data_dim, n_trg) + + def __del__(self): + if hasattr(self, "_ptr"): + ffi.PVFMMDestroyVolumeTreeD(ctypes.byref(ctypes.c_void_p(self._ptr))) + + def evaluate(self, fmm: FMMDoubleVolumeContext, loc_size: int) -> np.ndarray: + # TODO: XXX: in C code this uses kdim1, but data_dim was passed kdim0 + trg_val = np.empty(self.n_trg * self.data_dim) + ffi.PVFMMEvaluateVolumeFMMD(trg_val, self._ptr, fmm._ptr, loc_size) + self._evaluated = True return trg_val + + def leaf_count(self): + return int(ffi.PVFMMGetLeafCountD(self._ptr)) + + def get_leaf_coordinates(self): + Nleaf = self.leaf_count() + leaf_coord = np.empty(Nleaf * 3) + ffi.PVFMMGetLeafCoordD(leaf_coord, self._ptr) + return leaf_coord + + # WHERE DO I GET KDIM0? + + def get_coefficients(self): + if not self._evaluated: + raise ValueError( + "Cannot get coefficients of an un-evaluated tree" + ) # TODO: true? + n_leaf = self.leaf_count() + # TODO: XXX: in C code this uses kdim1, but data_dim was passed kdim0 + coeff = np.empty(n_leaf * self.n_coeff * self.data_dim) + ffi.PVFMMGetPotentialCoeffD(coeff, self._ptr) + return coeff + + def get_values(self): + coeff = self.get_coefficients() + n_leaf = self.leaf_count() + # TODO: XXX: in C code these two lines use kdim1, but data_dim was passed kdim0 + value = np.empty(n_leaf * self.n_cheb * self.data_dim) + ffi.PVFMMCoeff2NodesD(value, n_leaf, self.cheb_deg, self.data_dim, coeff) + return value From 1dd9e8c1f956d56dbaee0ae44c9593112c5f2d69 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 2 Aug 2023 15:58:53 -0400 Subject: [PATCH 07/11] Typo --- python/example2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/example2.py b/python/example2.py index 5297cbe..3068db9 100644 --- a/python/example2.py +++ b/python/example2.py @@ -1,4 +1,4 @@ -# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example1.py +# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example2.py # based on examples/src/example2-c.c From 54b12f3543eac40a54990ff8a6532156e43ad007 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 29 Aug 2023 12:19:54 -0400 Subject: [PATCH 08/11] Generalize use of kernel dimension Can we even get rid of data_dim input? --- python/notes.md | 2 +- python/pvfmm.py | 68 ++++++++++++++++++++++++++++++++----------------- 2 files changed, 46 insertions(+), 24 deletions(-) diff --git a/python/notes.md b/python/notes.md index fe03305..1db6a61 100644 --- a/python/notes.md +++ b/python/notes.md @@ -1,7 +1,7 @@ - no error handling - examples (translate): - [x] example-c.c - - [ ] example2-c.c test1 + - [x] example2-c.c test1 - MPI init (not needed - see mpi4py doc: https://mpi4py.readthedocs.io/en/stable/overview.html#initialization-and-exit) - function pointers for PVFMMCreateVolumeTree functions (numba for speed: https://numba.readthedocs.io/en/stable/user/cfunc.html) - first run can take a long time, expected. diff --git a/python/pvfmm.py b/python/pvfmm.py index 857200e..b1f8179 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -1,6 +1,3 @@ -# module load gcc/11.3.0 openmpi python/3.10.8 openblas fftw3 -# source ~/venvs/pvfmm/bin/activate - import ctypes import numpy as np from enum import Enum @@ -12,8 +9,18 @@ import ffi +__all__ = [ + "FMMKernel", + "FMMDoubleVolumeContext", + "FFMDoubleParticleContext", + "FMMDoubleVolumeTree", + "nodes_to_coeff", +] + -def nodes_to_coeff(N_leaf: int, cheb_deg: int, dof: int, node_val: np.ndarray): +def nodes_to_coeff( + N_leaf: int, cheb_deg: int, dof: int, node_val: np.ndarray +) -> np.ndarray: is_double = node_val.dtype == np.float64 Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 # XXX: is this valid? @@ -37,6 +44,17 @@ class FMMKernel(Enum): BiotSavartPotential = 5 +# read out of calls to BuildKernel in pvfmm/include/kernel.txx +KERNEL_DIMS = { + FMMKernel.LaplacePotential: (1, 1), + FMMKernel.LaplaceGradient: (1, 3), + FMMKernel.StokesPressure: (3, 1), + FMMKernel.StokesVelocity: (3, 3), + FMMKernel.StokesVelocityGrad: (3, 9), + FMMKernel.BiotSavartPotential: (3, 3), +} + + class FMMDoubleVolumeContext: def __init__( self, @@ -45,10 +63,14 @@ def __init__( kernel: FMMKernel, comm: MPI.Comm, ): + self.kernel = kernel if multipole_order <= 0 or multipole_order % 2 != 0: raise ValueError("multipole order must be even and postive") self._ptr = ffi.PVFMMCreateVolumeFMMD( - multipole_order, chebyshev_degree, int(kernel.value), ffi.get_MPI_COMM(comm) + multipole_order, + chebyshev_degree, + int(self.kernel.value), + ffi.get_MPI_COMM(comm), ) def __del__(self): @@ -65,13 +87,14 @@ def __init__( kernel: FMMKernel, comm: MPI.Comm, ): + self.kernel = kernel if multipole_order <= 0 or multipole_order % 2 != 0: raise ValueError("multipole order must be even and postive") self._ptr = ffi.PVFMMCreateContextD( float(box_size), max_points, multipole_order, - int(kernel.value), + int(self.kernel.value), ffi.get_MPI_COMM(comm), ) @@ -86,7 +109,7 @@ def evaluate( dl_den: Optional[np.ndarray], trg_pos: np.ndarray, setup: bool = True, - ): + ) -> np.ndarray: source_length = len(src_pos) if source_length % 3 != 0: raise ValueError( @@ -128,6 +151,7 @@ def __init__(self, ptr: ctypes.c_void_p, cheb_deg: int, data_dim: int, n_trg: in self.n_coeff = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 self.data_dim = data_dim self.n_trg = n_trg + self._used_kernel = None @classmethod def from_function( @@ -142,7 +166,7 @@ def from_function( max_pts: int, periodic: bool, init_depth: int, - ): + ) -> "FMMDoubleVolumeTree": n_trg = len(trg_coord) // 3 ptr = ffi.PVFMMCreateVolumeTreeD( cheb_deg, @@ -169,7 +193,7 @@ def from_coefficients( trg_coord: Optional[np.ndarray], comm: MPI.Comm, periodic: bool, - ): + ) -> "FMMDoubleVolumeTree": if len(leaf_coord) % 3 != 0: raise ValueError( "Leaf coordinates must have a length which is a multiple of 3" @@ -203,38 +227,36 @@ def __del__(self): ffi.PVFMMDestroyVolumeTreeD(ctypes.byref(ctypes.c_void_p(self._ptr))) def evaluate(self, fmm: FMMDoubleVolumeContext, loc_size: int) -> np.ndarray: - # TODO: XXX: in C code this uses kdim1, but data_dim was passed kdim0 - trg_val = np.empty(self.n_trg * self.data_dim) + (_kdim0, kdim1) = KERNEL_DIMS[fmm.kernel] + trg_val = np.empty(self.n_trg * kdim1) ffi.PVFMMEvaluateVolumeFMMD(trg_val, self._ptr, fmm._ptr, loc_size) - self._evaluated = True + self._used_kernel = fmm.kernel return trg_val - def leaf_count(self): + def leaf_count(self) -> int: return int(ffi.PVFMMGetLeafCountD(self._ptr)) - def get_leaf_coordinates(self): + def get_leaf_coordinates(self) -> np.ndarray: Nleaf = self.leaf_count() leaf_coord = np.empty(Nleaf * 3) ffi.PVFMMGetLeafCoordD(leaf_coord, self._ptr) return leaf_coord - # WHERE DO I GET KDIM0? - - def get_coefficients(self): - if not self._evaluated: + def get_coefficients(self) -> np.ndarray: + if self._used_kernel is None: raise ValueError( "Cannot get coefficients of an un-evaluated tree" ) # TODO: true? n_leaf = self.leaf_count() - # TODO: XXX: in C code this uses kdim1, but data_dim was passed kdim0 - coeff = np.empty(n_leaf * self.n_coeff * self.data_dim) + (_kdim0, kdim1) = KERNEL_DIMS[self._used_kernel] + coeff = np.empty(n_leaf * self.n_coeff * kdim1) ffi.PVFMMGetPotentialCoeffD(coeff, self._ptr) return coeff - def get_values(self): + def get_values(self) -> np.ndarray: coeff = self.get_coefficients() n_leaf = self.leaf_count() - # TODO: XXX: in C code these two lines use kdim1, but data_dim was passed kdim0 - value = np.empty(n_leaf * self.n_cheb * self.data_dim) + (_kdim0, kdim1) = KERNEL_DIMS[self._used_kernel] + value = np.empty(n_leaf * self.n_cheb * kdim1) ffi.PVFMMCoeff2NodesD(value, n_leaf, self.cheb_deg, self.data_dim, coeff) return value From f065654fe1211a0fa5c707cfd5c3024a343a51c4 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 20 Sep 2023 12:25:01 -0400 Subject: [PATCH 09/11] Double/float switch and checks --- python/example1.py | 2 +- python/example2.py | 6 +- python/ffi.py | 15 ++++- python/pvfmm.py | 161 +++++++++++++++++++++++++++++++++------------ 4 files changed, 137 insertions(+), 47 deletions(-) diff --git a/python/example1.py b/python/example1.py index a2bfc4e..9fc2245 100644 --- a/python/example1.py +++ b/python/example1.py @@ -83,7 +83,7 @@ def test_FMM(ctx): kernel = pvfmm.FMMKernel.BiotSavartPotential print("Loaded.") - ctx = pvfmm.FFMDoubleParticleContext( + ctx = pvfmm.FFMParticleContext( box_size, points_per_box, multipole_order, kernel, MPI.COMM_WORLD ) print("Running!") diff --git a/python/example2.py b/python/example2.py index 3068db9..e5be06d 100644 --- a/python/example2.py +++ b/python/example2.py @@ -96,7 +96,7 @@ def test1(fmm, kdim0, kdim1, cheb_deg, comm): Nt = 100 trg_coord = np.random.rand(Nt * 3) trg_value_ref = fn_poten(trg_coord) - tree = pvfmm.FMMDoubleVolumeTree.from_function( + tree = pvfmm.FMMVolumeTree.from_function( cheb_deg, kdim0, fn_input.ctypes, None, trg_coord, comm, 1e-6, 100, False, 0 ) trg_value = tree.evaluate(fmm, Nt) @@ -138,7 +138,7 @@ def test2(fmm, kdim0, kdim1, cheb_deg, comm): None, ) dense_coeff = pvfmm.nodes_to_coeff(Nleaf, cheb_deg, kdim0, dens_value) - tree = pvfmm.FMMDoubleVolumeTree.from_coefficients( + tree = pvfmm.FMMVolumeTree.from_coefficients( cheb_deg, kdim0, leaf_coord, dense_coeff, None, comm, False ) tree.evaluate(fmm, 0) @@ -161,7 +161,7 @@ def test2(fmm, kdim0, kdim1, cheb_deg, comm): kdim1 = 3 comm = MPI.COMM_WORLD - fmm = pvfmm.FMMDoubleVolumeContext( + fmm = pvfmm.FMMVolumeContext( mult_order, cheb_deg, pvfmm.FMMKernel.StokesVelocity, comm ) test1(fmm, kdim0, kdim1, cheb_deg, comm) diff --git a/python/ffi.py b/python/ffi.py index db69826..3d89262 100644 --- a/python/ffi.py +++ b/python/ffi.py @@ -1,4 +1,6 @@ import ctypes +import ctypes.util +import os import numpy as np from numpy.ctypeslib import ndpointer @@ -43,9 +45,18 @@ def volume_callback(type): PVFMMKernel = ctypes.c_uint # enum # try to load lib -# hardcoded path for now -SHARED_LIB = ctypes.CDLL("./build/libpvfmm.so") +_custom_location = os.getenv("PVFMM") +if _custom_location: + _dll = os.path.join(_custom_location, "libpvfmm.so") +else: + _dll = ctypes.util.find_library("pvfmm") + if _dll is None: + raise ImportError( + "Failed to find libpvfmm! \n" + "Set PVFMM environment variable or check your installation of pvfmm!" + ) +SHARED_LIB = ctypes.CDLL(_dll) # somwhat automatically generated: # clang2py --clang-args="-I/mnt/sw/nix/store/z5w5a7pr5cmdbds0pn9ajdgy0jg71sl6-gcc-11.3.0/lib/gcc/x86_64-pc-linux-gnu/11.3.0/include/" include/pvfmm.h -l ~/pvfmm/build/libpvfmm.so > python/auto.py diff --git a/python/pvfmm.py b/python/pvfmm.py index b1f8179..e0ba3ee 100644 --- a/python/pvfmm.py +++ b/python/pvfmm.py @@ -1,7 +1,7 @@ import ctypes import numpy as np from enum import Enum -from typing import Optional +from typing import Optional, Callable, Union # calls MPI_Init_thread() and sets up MPI_Finalize() for you. # see https://mpi4py.readthedocs.io/en/stable/mpi4py.run.html @@ -11,9 +11,9 @@ __all__ = [ "FMMKernel", - "FMMDoubleVolumeContext", - "FFMDoubleParticleContext", - "FMMDoubleVolumeTree", + "FMMVolumeContext", + "FFMParticleContext", + "FMMVolumeTree", "nodes_to_coeff", ] @@ -25,10 +25,10 @@ def nodes_to_coeff( Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 # XXX: is this valid? coeff = np.empty(Ncoef * N_leaf * dof, dtype=node_val.dtype) - if is_double: - ffi.PVFMMNodes2CoeffD(coeff, N_leaf, cheb_deg, dof, node_val) - else: - ffi.PVFMMNodes2CoeffF(coeff, N_leaf, cheb_deg, dof, node_val) + + get_function_dtype("PVFMMNodes2Coeff", node_val.dtype)( + coeff, N_leaf, cheb_deg, dof, node_val + ) return coeff @@ -55,18 +55,35 @@ class FMMKernel(Enum): } -class FMMDoubleVolumeContext: +def get_function_dtype(function_name: str, dtype: np.dtype) -> Callable: + """ + Helper function to switch between the double and float functions + from the FFI module + """ + if dtype == np.float64: + function_name += "D" + elif dtype == np.float32: + function_name += "F" + else: + raise ValueError("Invalid dtype, must be either float64 or float32") + + return getattr(ffi, function_name) + + +class FMMVolumeContext: def __init__( self, multipole_order: int, chebyshev_degree: int, kernel: FMMKernel, comm: MPI.Comm, + dtype=np.float64, ): self.kernel = kernel + self.dtype = np.dtype(dtype) if multipole_order <= 0 or multipole_order % 2 != 0: raise ValueError("multipole order must be even and postive") - self._ptr = ffi.PVFMMCreateVolumeFMMD( + self._ptr = get_function_dtype("PVFMMCreateVolumeFMM", dtype)( multipole_order, chebyshev_degree, int(self.kernel.value), @@ -75,10 +92,12 @@ def __init__( def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyVolumeFMMD(ctypes.byref(ctypes.c_void_p(self._ptr))) + get_function_dtype("PVFMMDestroyVolumeFMM", self.dtype)( + ctypes.byref(ctypes.c_void_p(self._ptr)) + ) -class FFMDoubleParticleContext: +class FFMParticleContext: def __init__( self, box_size: float, @@ -86,11 +105,14 @@ def __init__( multipole_order: int, kernel: FMMKernel, comm: MPI.Comm, + dtype=np.float64, ): self.kernel = kernel + self.dtype = np.dtype(dtype) if multipole_order <= 0 or multipole_order % 2 != 0: raise ValueError("multipole order must be even and postive") - self._ptr = ffi.PVFMMCreateContextD( + + self._ptr = get_function_dtype("PVFMMCreateContext", dtype)( float(box_size), max_points, multipole_order, @@ -100,7 +122,9 @@ def __init__( def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyContextD(ctypes.byref(ctypes.c_void_p(self._ptr))) + get_function_dtype("PVFMMDestroyContext", self.dtype)( + ctypes.byref(ctypes.c_void_p(self._ptr)) + ) def evaluate( self, @@ -110,15 +134,35 @@ def evaluate( trg_pos: np.ndarray, setup: bool = True, ) -> np.ndarray: + if src_pos.dtype != self.dtype: + raise ValueError( + f"Source array had the wrong dtype: {src_pos.dtype}. " + f"This object was created with dtype {self.dtype}" + ) + source_length = len(src_pos) if source_length % 3 != 0: raise ValueError( "Source arrays must have a length which is a multiple of 3" ) - if dl_den is not None and len(sl_den) != source_length: - raise ValueError("Source arrays must all be of the same length!") - if dl_den is not None and len(dl_den) != source_length * 2: - raise ValueError("Source arrays must all be of the same length!") + + if sl_den is not None: + if sl_den.dtype != self.dtype: + raise ValueError( + f"Source array had the wrong dtype: {src_pos.dtype}. " + f"This object was created with dtype {self.dtype}" + ) + if len(sl_den) != source_length: + raise ValueError("Source arrays must all be of the same length!") + if dl_den is not None: + if dl_den.dtype != self.dtype: + raise ValueError( + f"Source array had the wrong dtype: {src_pos.dtype}. " + f"This object was created with dtype {self.dtype}" + ) + if len(dl_den) != source_length * 2: + raise ValueError("Source arrays must all be of the same length!") + n_src = source_length // 3 target_length = len(trg_pos) @@ -127,9 +171,9 @@ def evaluate( "Target arrays must have a length which is a multiple of 3" ) n_trg = target_length // 3 - trg_val = np.empty(target_length) + trg_val = np.empty(target_length, dtype=self.dtype) - ffi.PVFMMEvalD( + get_function_dtype("PVFMMEval", self.dtype)( src_pos, sl_den, dl_den, @@ -143,14 +187,22 @@ def evaluate( return trg_val -class FMMDoubleVolumeTree: - def __init__(self, ptr: ctypes.c_void_p, cheb_deg: int, data_dim: int, n_trg: int): +class FMMVolumeTree: + def __init__( + self, + ptr: ctypes.c_void_p, + cheb_deg: int, + data_dim: int, + n_trg: int, + dtype: np.dtype, + ): self._ptr = ptr self.cheb_deg = cheb_deg self.n_cheb = (cheb_deg + 1) ** 3 self.n_coeff = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 self.data_dim = data_dim self.n_trg = n_trg + self.dtype = dtype self._used_kernel = None @classmethod @@ -158,7 +210,7 @@ def from_function( cls, cheb_deg: int, data_dim: int, - fn: ffi.double_volume_callback, + fn: Union[ffi.double_volume_callback, ffi.float_volume_callback], context: ctypes.c_void_p, # TODO: investigate trg_coord: np.ndarray, comm: MPI.Comm, @@ -166,9 +218,11 @@ def from_function( max_pts: int, periodic: bool, init_depth: int, - ) -> "FMMDoubleVolumeTree": + ) -> "FMMVolumeTree": n_trg = len(trg_coord) // 3 - ptr = ffi.PVFMMCreateVolumeTreeD( + + dtype = trg_coord.dtype + ptr = get_function_dtype("PVFMMCreateVolumeTree", dtype)( cheb_deg, data_dim, fn, @@ -181,7 +235,7 @@ def from_function( periodic, init_depth, ) - return cls(ptr, cheb_deg, data_dim, n_trg) + return cls(ptr, cheb_deg, data_dim, n_trg, dtype) @classmethod def from_coefficients( @@ -193,7 +247,7 @@ def from_coefficients( trg_coord: Optional[np.ndarray], comm: MPI.Comm, periodic: bool, - ) -> "FMMDoubleVolumeTree": + ) -> "FMMVolumeTree": if len(leaf_coord) % 3 != 0: raise ValueError( "Leaf coordinates must have a length which is a multiple of 3" @@ -207,9 +261,23 @@ def from_coefficients( "Function coefficients array has the wrong length, required length " + fn_coeff_size ) - n_trg = len(trg_coord) // 3 if trg_coord is not None else 0 + dtype = leaf_coord.dtype + if fn_coeff.dtype != dtype: + raise ValueError( + f"Mismatching dtypes. Leaves had dtype {dtype}, " + f"but coefficients had dtype {fn_coeff.dtype}" + ) + if trg_coord is not None: + n_trg = len(trg_coord) // 3 + if trg_coord.dtype != dtype: + raise ValueError( + f"Mismatching dtypes. Leaves had dtype {dtype}, " + f"but targets had dtype {trg_coord.dtype}" + ) + else: + n_trg = 0 - ptr = ffi.PVFMMCreateVolumeTreeFromCoeffD( + ptr = get_function_dtype("PVFMMCreateVolumeTreeFromCoeff", dtype)( N_leaf, cheb_deg, data_dim, @@ -220,26 +288,35 @@ def from_coefficients( ffi.get_MPI_COMM(comm), periodic, ) - return cls(ptr, cheb_deg, data_dim, n_trg) + return cls(ptr, cheb_deg, data_dim, n_trg, dtype) def __del__(self): if hasattr(self, "_ptr"): - ffi.PVFMMDestroyVolumeTreeD(ctypes.byref(ctypes.c_void_p(self._ptr))) + get_function_dtype("PVFMMDestroyVolumeTree", self.dtype)( + ctypes.byref(ctypes.c_void_p(self._ptr)) + ) - def evaluate(self, fmm: FMMDoubleVolumeContext, loc_size: int) -> np.ndarray: + def evaluate(self, fmm: FMMVolumeContext, loc_size: int) -> np.ndarray: + if fmm.dtype != self.dtype: + raise ValueError( + f"Volume context has dtype {fmm.dtype}, " + f"but this tree has dtype {self.dtype}" + ) (_kdim0, kdim1) = KERNEL_DIMS[fmm.kernel] - trg_val = np.empty(self.n_trg * kdim1) - ffi.PVFMMEvaluateVolumeFMMD(trg_val, self._ptr, fmm._ptr, loc_size) + trg_val = np.empty(self.n_trg * kdim1, dtype=self.dtype) + get_function_dtype("PVFMMEvaluateVolumeFMM", self.dtype)( + trg_val, self._ptr, fmm._ptr, loc_size + ) self._used_kernel = fmm.kernel return trg_val def leaf_count(self) -> int: - return int(ffi.PVFMMGetLeafCountD(self._ptr)) + return int(get_function_dtype("PVFMMGetLeafCount", self.dtype)(self._ptr)) def get_leaf_coordinates(self) -> np.ndarray: Nleaf = self.leaf_count() - leaf_coord = np.empty(Nleaf * 3) - ffi.PVFMMGetLeafCoordD(leaf_coord, self._ptr) + leaf_coord = np.empty(Nleaf * 3, dtype=self.dtype) + get_function_dtype("PVFMMGetLeafCoord", self.dtype)(leaf_coord, self._ptr) return leaf_coord def get_coefficients(self) -> np.ndarray: @@ -249,14 +326,16 @@ def get_coefficients(self) -> np.ndarray: ) # TODO: true? n_leaf = self.leaf_count() (_kdim0, kdim1) = KERNEL_DIMS[self._used_kernel] - coeff = np.empty(n_leaf * self.n_coeff * kdim1) - ffi.PVFMMGetPotentialCoeffD(coeff, self._ptr) + coeff = np.empty(n_leaf * self.n_coeff * kdim1, dtype=self.dtype) + get_function_dtype("PVFMMGetPotentialCoeff", self.dtype)(coeff, self._ptr) return coeff def get_values(self) -> np.ndarray: coeff = self.get_coefficients() n_leaf = self.leaf_count() (_kdim0, kdim1) = KERNEL_DIMS[self._used_kernel] - value = np.empty(n_leaf * self.n_cheb * kdim1) - ffi.PVFMMCoeff2NodesD(value, n_leaf, self.cheb_deg, self.data_dim, coeff) + value = np.empty(n_leaf * self.n_cheb * kdim1, dtype=self.dtype) + get_function_dtype("PVFMMCoeff2Nodes", self.dtype)( + value, n_leaf, self.cheb_deg, self.data_dim, coeff + ) return value From 83e0d22bce9a0da5020e82cf12f7e5b2e7311cd9 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 20 Sep 2023 15:02:00 -0400 Subject: [PATCH 10/11] Package, clean up --- .gitignore | 1 + examples/src/example2-c.c | 2 +- python/README.md | 26 +++++++++++++ python/{ => examples}/example1.py | 7 +--- python/{ => examples}/example2.py | 46 +++++++++++------------ python/notes.md | 7 ---- python/pyproject.toml | 31 +++++++++++++++ python/src/pvfmm/__init__.py | 18 +++++++++ python/{ => src/pvfmm}/ffi.py | 4 ++ python/src/pvfmm/py.typed | 0 python/{pvfmm.py => src/pvfmm/wrapper.py} | 12 +----- 11 files changed, 106 insertions(+), 48 deletions(-) create mode 100644 python/README.md rename python/{ => examples}/example1.py (91%) rename python/{ => examples}/example2.py (90%) delete mode 100644 python/notes.md create mode 100644 python/pyproject.toml create mode 100644 python/src/pvfmm/__init__.py rename python/{ => src/pvfmm}/ffi.py (98%) create mode 100644 python/src/pvfmm/py.typed rename python/{pvfmm.py => src/pvfmm/wrapper.py} (98%) diff --git a/.gitignore b/.gitignore index 3fe5ccf..936d505 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,4 @@ libtool build/ *.pyc +pvfmm.egg-info/ diff --git a/examples/src/example2-c.c b/examples/src/example2-c.c index 75dbfab..30637e1 100644 --- a/examples/src/example2-c.c +++ b/examples/src/example2-c.c @@ -169,7 +169,7 @@ int main(int argc, char** argv) { // Build FMM translation operators void* fmm = PVFMMCreateVolumeFMMD(mult_order, cheb_deg, PVFMMStokesVelocity, comm); - test1(fmm, kdim0, kdim1, cheb_deg, comm); + // test1(fmm, kdim0, kdim1, cheb_deg, comm); test2(fmm, kdim0, kdim1, cheb_deg, comm); PVFMMDestroyVolumeFMMD(&fmm); diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..eb368b1 --- /dev/null +++ b/python/README.md @@ -0,0 +1,26 @@ +# Python bindings to PVFMM + +This package provides Python bindings to the [C API](../include/pvfmm.h) using +[`ctypes`](https://docs.python.org/3/library/ctypes.html) and +[`mpi4py`](https://mpi4py.readthedocs.io/en/stable/index.html). + +## Installation + +Build and install the PVFMM library as usual, and then run `pip install .` +inside this directory. + +## Usage + +If you did not install PVFMM (see [INSTALL](../INSTALL)), you will need +to tell the Python library where to find the `libpvfmm.so` file with +`export PVFMM=/path/to/build/folder`. + +Then, you should be able to run code which uses it by using +``` +mpirun $other_args_here python -m mpi4py path/to/program.py +``` + +## Examples + +See the `examples/` subdirectory for Python implementations of a few example programs. +The examples rely on [numba](https://numba.pydata.org/). diff --git a/python/example1.py b/python/examples/example1.py similarity index 91% rename from python/example1.py rename to python/examples/example1.py index 9fc2245..e162a1b 100644 --- a/python/example1.py +++ b/python/examples/example1.py @@ -1,7 +1,3 @@ -# module load gcc/11.3.0 openmpi python-mpi/3.10.8 fftw/mpi openblas -# export OMP_NUM_THREADS=16 -# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example1.py - # based on examples/src/example-c.c import time @@ -83,10 +79,9 @@ def test_FMM(ctx): kernel = pvfmm.FMMKernel.BiotSavartPotential print("Loaded.") - ctx = pvfmm.FFMParticleContext( + ctx = pvfmm.FMMParticleContext( box_size, points_per_box, multipole_order, kernel, MPI.COMM_WORLD ) print("Running!") test_FMM(ctx) - # MPI finalize handled by mpi4py atexit handler diff --git a/python/example2.py b/python/examples/example2.py similarity index 90% rename from python/example2.py rename to python/examples/example2.py index e5be06d..98b5738 100644 --- a/python/example2.py +++ b/python/examples/example2.py @@ -1,5 +1,3 @@ -# mpirun -n 1 --map-by slot:pe=$OMP_NUM_THREADS python -m mpi4py python/example2.py - # based on examples/src/example2-c.c import numpy as np @@ -10,19 +8,10 @@ import pvfmm -c_sig = types.void( - types.CPointer(types.double), - types.int64, - types.CPointer(types.double), - types.voidptr, -) - -# see https://numba.readthedocs.io/en/stable/user/cfunc.html -@numba.cfunc(c_sig, nopython=True) -def fn_input(coord_, n, out_, _ctx): - coord = numba.carray(coord_, n * 3) - out = numba.carray(out_, n * 3) +@numba.njit +def fn_input(coord, out): + n = len(coord) // 3 L = 125 for i in range(n): idx = i * 3 @@ -36,6 +25,20 @@ def fn_input(coord_, n, out_, _ctx): -L * r_2 ) + 2 * L * np.exp(-L * r_2) * (c[2] - 0.5) +c_sig = types.void( + types.CPointer(types.double), + types.int64, + types.CPointer(types.double), + types.voidptr, +) + +# see https://numba.readthedocs.io/en/stable/user/cfunc.html +@numba.cfunc(c_sig, nopython=True) +def fn_input_C(coord_, n, out_, _ctx): + coord = numba.carray(coord_, n * 3) + out = numba.carray(out_, n * 3) + fn_input(coord, out) + @numba.njit def fn_poten(coord): @@ -97,7 +100,7 @@ def test1(fmm, kdim0, kdim1, cheb_deg, comm): trg_coord = np.random.rand(Nt * 3) trg_value_ref = fn_poten(trg_coord) tree = pvfmm.FMMVolumeTree.from_function( - cheb_deg, kdim0, fn_input.ctypes, None, trg_coord, comm, 1e-6, 100, False, 0 + cheb_deg, kdim0, fn_input_C.ctypes, None, trg_coord, comm, 1e-6, 100, False, 0 ) trg_value = tree.evaluate(fmm, Nt) @@ -111,8 +114,6 @@ def test2(fmm, kdim0, kdim1, cheb_deg, comm): Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 depth = 3 - # skipping MPI (for now?) - Nleaf = 1 << (3 * depth) leaf_length = 1 / (1 << depth) @@ -131,12 +132,9 @@ def test2(fmm, kdim0, kdim1, cheb_deg, comm): ) print("Getting nodes") cheb_coord = GetChebNodes(Nleaf, cheb_deg, depth, leaf_coord) - fn_input( # ugly just because we used numba above - cheb_coord.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), - Nleaf * Ncheb, - dens_value.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), - None, - ) + + fn_input(cheb_coord, dens_value) + dense_coeff = pvfmm.nodes_to_coeff(Nleaf, cheb_deg, kdim0, dens_value) tree = pvfmm.FMMVolumeTree.from_coefficients( cheb_deg, kdim0, leaf_coord, dense_coeff, None, comm, False @@ -144,7 +142,7 @@ def test2(fmm, kdim0, kdim1, cheb_deg, comm): tree.evaluate(fmm, 0) potn_value = tree.get_values() - Nleaf = tree.leaf_count() # does this change from above? + Nleaf = tree.leaf_count() # TODO: does this change from above? leaf_coord = tree.get_leaf_coordinates() cheb_coord = GetChebNodes(Nleaf, cheb_deg, depth, leaf_coord) potn_value_ref = fn_poten(cheb_coord) diff --git a/python/notes.md b/python/notes.md deleted file mode 100644 index 1db6a61..0000000 --- a/python/notes.md +++ /dev/null @@ -1,7 +0,0 @@ -- no error handling -- examples (translate): - - [x] example-c.c - - [x] example2-c.c test1 -- MPI init (not needed - see mpi4py doc: https://mpi4py.readthedocs.io/en/stable/overview.html#initialization-and-exit) -- function pointers for PVFMMCreateVolumeTree functions (numba for speed: https://numba.readthedocs.io/en/stable/user/cfunc.html) -- first run can take a long time, expected. diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 0000000..12a791f --- /dev/null +++ b/python/pyproject.toml @@ -0,0 +1,31 @@ +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + +[project] +name = "pvfmm" +readme = "README.md" +authors = [{ "name" = "Brian Ward", "email" = "bward@flatironinstitute.org" }] +dependencies = ["numpy", "mpi4py"] +requires-python = ">=3.8" +license = { text = "BSD-3-Clause" } +classifiers = [ + "Programming Language :: Python :: 3", + "Development Status :: 4 - Beta", +] +dynamic = ["version"] + +[project.optional-dependencies] +all = ["numba"] + +[tool.setuptools.dynamic] +version = { attr = "pvfmm.__version__" } + +[tool.isort] +profile = "black" + +[tool.setuptools.package-data] +"pvfmm" = ["py.typed"] + +[tool.setuptools.packages.find] +where = ["src"] diff --git a/python/src/pvfmm/__init__.py b/python/src/pvfmm/__init__.py new file mode 100644 index 0000000..8e24482 --- /dev/null +++ b/python/src/pvfmm/__init__.py @@ -0,0 +1,18 @@ +from .wrapper import ( + FMMKernel, + FMMVolumeContext, + FMMParticleContext, + FMMVolumeTree, + nodes_to_coeff, +) + + +__all__ = [ + "FMMKernel", + "FMMVolumeContext", + "FMMParticleContext", + "FMMVolumeTree", + "nodes_to_coeff", +] + +__version__ = "0.1.0" diff --git a/python/ffi.py b/python/src/pvfmm/ffi.py similarity index 98% rename from python/ffi.py rename to python/src/pvfmm/ffi.py index 3d89262..5570d1a 100644 --- a/python/ffi.py +++ b/python/src/pvfmm/ffi.py @@ -19,6 +19,10 @@ def get_MPI_COMM(comm) -> MPI_Comm: def wrapped_ndptr(*args, **kwargs): + """ + A version of np.ctypeslib.ndpointer + which allows None (passed as NULL) + """ base = ndpointer(*args, **kwargs) def from_param(cls, obj): diff --git a/python/src/pvfmm/py.typed b/python/src/pvfmm/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/python/pvfmm.py b/python/src/pvfmm/wrapper.py similarity index 98% rename from python/pvfmm.py rename to python/src/pvfmm/wrapper.py index e0ba3ee..eef71e6 100644 --- a/python/pvfmm.py +++ b/python/src/pvfmm/wrapper.py @@ -7,15 +7,7 @@ # see https://mpi4py.readthedocs.io/en/stable/mpi4py.run.html from mpi4py import MPI -import ffi - -__all__ = [ - "FMMKernel", - "FMMVolumeContext", - "FFMParticleContext", - "FMMVolumeTree", - "nodes_to_coeff", -] +from . import ffi def nodes_to_coeff( @@ -97,7 +89,7 @@ def __del__(self): ) -class FFMParticleContext: +class FMMParticleContext: def __init__( self, box_size: float, From b1d5df52d9384e64d611f39b57173289eb1df235 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 20 Sep 2023 15:14:48 -0400 Subject: [PATCH 11/11] Clean up comments --- examples/src/example2-c.c | 2 +- python/examples/example2.py | 2 ++ python/src/pvfmm/wrapper.py | 6 +++--- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/src/example2-c.c b/examples/src/example2-c.c index 30637e1..e33e2ce 100644 --- a/examples/src/example2-c.c +++ b/examples/src/example2-c.c @@ -169,7 +169,7 @@ int main(int argc, char** argv) { // Build FMM translation operators void* fmm = PVFMMCreateVolumeFMMD(mult_order, cheb_deg, PVFMMStokesVelocity, comm); - // test1(fmm, kdim0, kdim1, cheb_deg, comm); + //test1(fmm, kdim0, kdim1, cheb_deg, comm); test2(fmm, kdim0, kdim1, cheb_deg, comm); PVFMMDestroyVolumeFMMD(&fmm); diff --git a/python/examples/example2.py b/python/examples/example2.py index 98b5738..21ba191 100644 --- a/python/examples/example2.py +++ b/python/examples/example2.py @@ -25,6 +25,7 @@ def fn_input(coord, out): -L * r_2 ) + 2 * L * np.exp(-L * r_2) * (c[2] - 0.5) + c_sig = types.void( types.CPointer(types.double), types.int64, @@ -32,6 +33,7 @@ def fn_input(coord, out): types.voidptr, ) + # see https://numba.readthedocs.io/en/stable/user/cfunc.html @numba.cfunc(c_sig, nopython=True) def fn_input_C(coord_, n, out_, _ctx): diff --git a/python/src/pvfmm/wrapper.py b/python/src/pvfmm/wrapper.py index eef71e6..e54642e 100644 --- a/python/src/pvfmm/wrapper.py +++ b/python/src/pvfmm/wrapper.py @@ -15,7 +15,7 @@ def nodes_to_coeff( ) -> np.ndarray: is_double = node_val.dtype == np.float64 Ncoef = (cheb_deg + 1) * (cheb_deg + 2) * (cheb_deg + 3) // 6 - # XXX: is this valid? + # TODO: is this the valid size of the output array? coeff = np.empty(Ncoef * N_leaf * dof, dtype=node_val.dtype) get_function_dtype("PVFMMNodes2Coeff", node_val.dtype)( @@ -203,7 +203,7 @@ def from_function( cheb_deg: int, data_dim: int, fn: Union[ffi.double_volume_callback, ffi.float_volume_callback], - context: ctypes.c_void_p, # TODO: investigate + context: ctypes.c_void_p, trg_coord: np.ndarray, comm: MPI.Comm, tol: float, @@ -315,7 +315,7 @@ def get_coefficients(self) -> np.ndarray: if self._used_kernel is None: raise ValueError( "Cannot get coefficients of an un-evaluated tree" - ) # TODO: true? + ) # TODO: is this true? what is the contract of this class n_leaf = self.leaf_count() (_kdim0, kdim1) = KERNEL_DIMS[self._used_kernel] coeff = np.empty(n_leaf * self.n_coeff * kdim1, dtype=self.dtype)