Transition to pyAMReX (#3474)

* pyAMReX: Build System

* CI Updates (Changed Options)

* Callback modernization (#4)

* refactor callbacks.py
* added binding code in `pyWarpX.cpp` to add or remove keys from the callback dictionary
* minor PR cleanups

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Added Python level reference to fetch the multifabs (#3)

* pyAMReX: Build System

* Added Python level reference to Ex_aux

* Now uses the multifab map

* Fix typo

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Add initialization and finalize routines (#5)

A basic PICMI input file will now run to completion.

* Regression Tests: WarpX_PYTHON=ON

* Update Imports to nD pyAMReX

* IPO/LTO Control

Although pybind11 relies heavily on IPO/LTO to create low-latency,
small-binary bindings, some compilers will have troubles with that.

Thus, we add a compile-time option to optionally disable it when
needed.

* Fix: Link Legacy WarpXWrappers.cpp

* Wrap WarpX instance and start multi particle container

* Fix test Python_pass_mpi_comm
* Start wrapper for multiparticle container
* Add return policy

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Update fields to use MultiFabs directly

Remove EOL white space

Removed old routines accessing MultiFabs

Update to use "node_centered"

* Fix compilation with Python

* Update fields.py to use modified MultiFab tag names

* Remove incorrect, unused code

* Add function to extract the WarpX MultiParticleContainer

* Complete class WarpXParticleContainer

* Wrap functions getNprocs / getMyProc

* restore `install___` callback API - could remove later if we want but should maintain backward compatibility for now

* add `gett_new` and `getistep` functions wrappers; fix typos in `callbacks.py`; avoid error in getting `rho` from `fields.py`

* Update callback call and `getNproc`/`getMyProc` function

* Replace function add_n_particles

* Fix setitem in fields.py for 1d and 2d

* also update `gett_new()` in `_libwarpx.py` in case we want to keep that API

* added binding for `WarpXParIter` - needed to port `libwarpx.depositChargeDensity()` which is an ongoing effort

* Wrap function num_real_comp

* added binding for `TotalNumberOfParticles` and continue progress on enabling 1d MCC test to run

* add `SyncRho()` binding and create helper function in `libwarpx.depositChargeDensity` to manage scope of the particle iter

* Clean up issues in fields.py

* update bindings for `get_particle_structs`

* Fix setitem in fields.py

* switch order of initialization for particle container and particle iterator

* switch deposit_charge loop to C++ code; bind `ApplyInverseVolumeScalingToChargeDensity`

* move `WarpXParticleContainer.cpp` and `MultiParticleContainer.cpp` to new Particles folder

* added binding for `ParticleBoundaryBuffer`

* More fixes for fields.py

* Fix: Backtraces from Python

Add the Python executable name with an absolute path, so backtraces
in AMReX work. See linked AMReX issue for details.

* Cleaning

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Fix: Backtraces from Python Part II

Do not add Python script name - it confuses the AMReX ParmParser to
build its table.

* Fix: CMake Dependencies for Wheel

This fixes a racecondition during `pip_install`: it was possible
that not all dimensions where yet build from pybind before we
start packing them in the wheel for pip install.

* MCC Test: Install Callbacks before Run

Otherwise hangs in aquiring the gil during shutdown.

* addition of `Python/pywarpx/particle_containers.py` and various associated bindings

* Fix: CMake Superbuild w/ Shared AMReX

We MUST build AMReX as a shared (so/dll/dylib) library, otherwise
all the global state in it will cause split-brain situations, where
our Python modules operate on different stacks.

* add `clear_all()` to callbacks in order to remove all callbacks at finalize

* add `-DWarpX_PYTHON=ON` to CI tests that failed to build

* add `get_comp_index` and continue to port particle data bindings

* Add AMReX Module as `libwarpx_so.amr`

Attribute to pass through the already loaded AMReX module with the
matching dimensionality to the simulation.

* Fix for fields accounting for guard cells

* Fix handling of ghost cells in fields

* Update & Test: Particle Boundary Scraping

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* CI: Python Updates

- modernize Python setups
- drop CUDA 11.0 for good and go 11.3+ as documented already
```
Error #3246: Internal Compiler Error (codegen): "there was an error in verifying the lgenfe output!"
```

* CI: Python Updates (chmod)

* Add support for cupy in fields.py

* Add MultiFab reduction routines

* CI: CUDA 11.3 is <= Ubuntu 20.04

* changed `AddNParticles` to take `amrex::Vector` arguments

* setup.py: WarpX_PYTHON=ON

* update various 2d and rz tests with new APIs

* add `-DWarpX_PYTHON_IPO=OFF` to compile string for 2d and 3d Python CI tests to speed up linking

* CI: -DpyAMReX_IPO=OFF

* CI: -DpyAMReX_IPO=OFF

actually adding `=OFF`

* CI: Intel Python

* CI: macOS Python Executable

Ensure we always use the same `python3` executable, as specified
by the `PATH` priority.

* CMake: Python Multi-Config Build

Add support for multi-config generators, especially on Windows.

* __init__.py: Windows DLL Support

Python 3.8+ on Windows: DLL search paths for dependent
shared libraries
Refs.:
- https://github.com/python/cpython/issues/80266
- https://docs.python.org/3.8/library/os.html#os.add_dll_directory

* CI: pywarpx Update

our setup.py cannot install pyamrex yet as a dependency.

* ABLASTR: `ablastr/export.H`

Add a new header to export public globals that are not covered by
`WINDOWS_EXPORT_ALL_SYMBOLS`.
https://stackoverflow.com/questions/54560832/cmake-windows-export-all-symbols-does-not-cover-global-variables/54568678#54568678

* WarpX: EXPORT Globals in `.dll` files

WarpX still uses a lot of globals:
- `static` member variables
- `extern` global variables

These globals cannot be auto-exported with CMake's
`WINDOWS_EXPORT_ALL_SYMBOLS` helper and thus we need to mark them
manually for DLL export (and import) via the new ABLASTR
`ablastr/export.H` helper macros.

This starts to export symbols in the:
- WarpX and particle container classes
- callback hook database map
- ES solver

* CI: pywarpx Clang CXXFLAGS Down

Move CXXFLAGS (`-Werror ...`) down until deps are installed.

* GNUmake: Generate `ablastr/export.H`

* CMake: More Symbol Exports for Windows

* `WarpX-tests.ini`: Simplify Python no-IPO

Also avoids subtle differences in compilation that increase
compile time.

* Update PICMI_inputs_EB_API.py for embedded_boundary_python_API CI test

* Fix Python_magnetostatic_eb_3d

* Update: Python_restart_runtime_components

New Python APIs

* Windows: no dllimport for now

* CI: Skip `PYINSTALLOPTIONS` For Now

* CMake: Dependency Bump Min-Versions

for external packages picked up by `find_package`.

* Fix F and G_fp names in fields.py

* Tests: Disable `Python_pass_mpi_comm`

* Wrappers: Cleanup

* pyWarpX: Include Cleaning

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fields.py: Fix F and G Wrappers

Correct MultiFab names (w/o components).

* Remove unused in fields.py

* Windows: New Export Headers

- ABLASTR: `ablastr/export.H`
- WarpX: `Utils/export.H`

* removed `WarpInterface.py` since that functionality is now in `particle_containers.py`; removed parts of `WarpXWrappers.cpp` that have been ported to pyamrex

* CMake: Link OBJECT Target PRIVATE

* CMake: Remove OBJECT Target

Simplify and make `app` link `lib` (default: static). Remove OBJECT
target.

* Fix in fields.py for the components index

* Update get_particle_id/cpu

As implemented in pyAMReX with
https://github.com/AMReX-Codes/pyamrex/pull/165

* WarpX: Update for Private Constructor

* Import AMReX Before pyd DLL Call

Importing AMReX will add the `add_dll_directory` to a potentially
shared amrex DLL on Windows.

* Windows CI: Set PATH to amrex_Nd.dll

* CMake: AMReX_INSTALL After Python

In superbuild, Python can modify `AMReX_BUILD_SHARED_LIBS`.

* Clang Win CI: Manually Install requirements

Sporadic error is:
```
...
Installing collected packages: pyparsing, numpy, scipy, periodictable, picmistandard
ERROR: Could not install packages due to an OSError: [WinError 32] The process cannot access the file because it is being used by another process: 'C:\\hostedtoolcache\\windows\\Python\\3.11.4\\x64\\Lib\\site-packages\\numpy\\polynomial\\__init__.py'
Consider using the `--user` option or check the permissions.
```

* Hopefully final fixes to fields.py

* Update getProbLo/getProbHi

* Set plasma length strength

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Fix fields method to remove CodeQL notice

* Update Comments & Some Finalize

* Move: set_plasma_lens_strength

to MPC

---------

Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>
Co-authored-by: David Grote <dpgrote@lbl.gov>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Dave Grote <grote1@llnl.gov>
Co-authored-by: Roelof Groenewald <regroenewald@gmail.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

13 files changed, 679 insertions, 838 deletions

diff --git a/Source/Python/CMakeLists.txt b/Source/Python/CMakeLists.txt
index e10b8486c..4b4be199c 100644
--- a/Source/Python/CMakeLists.txt
+++ b/Source/Python/CMakeLists.txt
@@ -1,8 +1,24 @@
+###############################################################################
+# These are the files equivalent to the primary C++ implementation - but here
+# we define how they will appear in our Python module (aka Python bindings).
+#
 foreach(D IN LISTS WarpX_DIMS)
     warpx_set_suffix_dims(SD ${D})
-    target_sources(WarpX_${SD}
+    target_sources(lib_${SD}
       PRIVATE
+        # callback hooks
         WarpX_py.cpp
+
+        # legacy C wrapper APIs
         WarpXWrappers.cpp
     )
+    if(WarpX_PYTHON)
+        target_sources(pyWarpX_${SD}
+          PRIVATE
+            # pybind11
+            WarpX.cpp
+        )
+    endif()
 endforeach()
+
+add_subdirectory(Particles)
diff --git a/Source/Python/Particles/CMakeLists.txt b/Source/Python/Particles/CMakeLists.txt
new file mode 100644
index 000000000..eed1bb07c
--- /dev/null
+++ b/Source/Python/Particles/CMakeLists.txt
@@ -0,0 +1,17 @@
+###############################################################################
+# These are the files equivalent to the primary C++ implementation - but here
+# we define how they will appear in our Python module (aka Python bindings).
+#
+foreach(D IN LISTS WarpX_DIMS)
+    warpx_set_suffix_dims(SD ${D})
+    if(WarpX_PYTHON)
+        target_sources(pyWarpX_${SD}
+          PRIVATE
+            # pybind11
+            ParticleBoundaryBuffer.cpp
+            MultiParticleContainer.cpp
+            PinnedMemoryParticleContainer.cpp
+            WarpXParticleContainer.cpp
+        )
+    endif()
+endforeach()
diff --git a/Source/Python/Particles/MultiParticleContainer.cpp b/Source/Python/Particles/MultiParticleContainer.cpp
new file mode 100644
index 000000000..e709f0950
--- /dev/null
+++ b/Source/Python/Particles/MultiParticleContainer.cpp
@@ -0,0 +1,46 @@
+/* Copyright 2021-2022 The WarpX Community
+ *
+ * Authors: Axel Huebl, Remi Lehe
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Python/pyWarpX.H"
+
+#include <Particles/MultiParticleContainer.H>
+
+#include <AMReX_GpuContainers.H>
+#include <AMReX_REAL.H>
+
+
+void init_MultiParticleContainer (py::module& m)
+{
+    py::class_<MultiParticleContainer>(m, "MultiParticleContainer")
+        .def("get_particle_container_from_name",
+            &MultiParticleContainer::GetParticleContainerFromName,
+            py::arg("name"),
+            py::return_value_policy::reference_internal
+        )
+
+        .def("set_plasma_lens_strength",
+             [](MultiParticleContainer& mpc, int i_lens, amrex::Real strength_E, amrex::Real strength_B) {
+                 mpc.h_repeated_plasma_lens_strengths_E.at(i_lens) = strength_E;
+                 mpc.h_repeated_plasma_lens_strengths_B.at(i_lens) = strength_B;
+                 amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                                       mpc.h_repeated_plasma_lens_strengths_E.begin(), mpc.h_repeated_plasma_lens_strengths_E.end(),
+                                       mpc.d_repeated_plasma_lens_strengths_E.begin());
+                 amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
+                                       mpc.h_repeated_plasma_lens_strengths_B.begin(), mpc.h_repeated_plasma_lens_strengths_B.end(),
+                                       mpc.d_repeated_plasma_lens_strengths_B.begin());
+                 amrex::Gpu::synchronize();
+             },
+             py::arg("i_lens"), py::arg("strength_E"), py::arg("strength_B"),
+             R"pbdoc(Set the strength of the `i_lens`-th lens
+Parameters
+----------
+i_lens: int
+  Index of the lens to be modified
+strength_E, strength_B: floats
+  The electric and magnetic focusing strength of the lens)pbdoc"
+        )
+    ;
+}
diff --git a/Source/Python/Particles/ParticleBoundaryBuffer.cpp b/Source/Python/Particles/ParticleBoundaryBuffer.cpp
new file mode 100644
index 000000000..2a35faece
--- /dev/null
+++ b/Source/Python/Particles/ParticleBoundaryBuffer.cpp
@@ -0,0 +1,59 @@
+/* Copyright 2021-2023 The WarpX Community
+ *
+ * Authors: Axel Huebl, Remi Lehe, Roelof Groenewald
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Python/pyWarpX.H"
+
+#include <Particles/ParticleBoundaryBuffer.H>
+
+namespace warpx {
+    class BoundaryBufferParIter
+        : public amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+    {
+    public:
+        using amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ParIter;
+
+        BoundaryBufferParIter(ContainerType& pc, int level) :
+            amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>(pc, level) {}
+    };
+}
+
+void init_BoundaryBufferParIter (py::module& m)
+{
+    py::class_<
+        warpx::BoundaryBufferParIter,
+        amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+    >(m, "BoundaryBufferParIter")
+        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ContainerType&, int>(),
+            py::arg("particle_container"), py::arg("level")
+        )
+    ;
+}
+
+void init_ParticleBoundaryBuffer (py::module& m)
+{
+    py::class_<ParticleBoundaryBuffer>(m, "ParticleBoundaryBuffer")
+        .def(py::init<>())
+        .def("clear_particles",
+            [](ParticleBoundaryBuffer& pbb) { pbb.clearParticles(); }
+        )
+        .def("get_particle_container",
+            [](ParticleBoundaryBuffer& pbb,
+            const std::string species_name, int boundary) {
+                return &pbb.getParticleBuffer(species_name, boundary);
+            },
+            py::arg("species_name"), py::arg("boundary"),
+            py::return_value_policy::reference_internal
+        )
+        .def("get_num_particles_in_container",
+            [](ParticleBoundaryBuffer& pbb,
+            const std::string species_name, int boundary, bool local)
+            {
+                return pbb.getNumParticlesInContainer(species_name, boundary, local);
+            },
+            py::arg("species_name"), py::arg("boundary"), py::arg("local")=true
+        )
+    ;
+}
diff --git a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp
new file mode 100644
index 000000000..d048b4d9c
--- /dev/null
+++ b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp
@@ -0,0 +1,23 @@
+/* Copyright 2021-2023 The WarpX Community
+ *
+ * Authors: Axel Huebl, Remi Lehe, Roelof Groenewald
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Python/pyWarpX.H"
+
+#include <Particles/PinnedMemoryParticleContainer.H>
+
+
+void init_PinnedMemoryParticleContainer (py::module& m)
+{
+    py::class_<
+        PinnedMemoryParticleContainer,
+        amrex::ParticleContainer<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+    > pmpc (m, "PinnedMemoryParticleContainer");
+    pmpc
+        .def("num_int_comps",
+            [](PinnedMemoryParticleContainer& pc) { return pc.NumIntComps(); }
+        )
+    ;
+}
diff --git a/Source/Python/Particles/WarpXParticleContainer.cpp b/Source/Python/Particles/WarpXParticleContainer.cpp
new file mode 100644
index 000000000..b5d3b1626
--- /dev/null
+++ b/Source/Python/Particles/WarpXParticleContainer.cpp
@@ -0,0 +1,119 @@
+/* Copyright 2021-2022 The WarpX Community
+ *
+ * Authors: Axel Huebl, Remi Lehe
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Python/pyWarpX.H"
+
+#include <Particles/WarpXParticleContainer.H>
+
+
+void init_WarpXParIter (py::module& m)
+{
+    py::class_<
+        WarpXParIter, amrex::ParIter<0,0,PIdx::nattribs>
+    >(m, "WarpXParIter")
+        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs>::ContainerType&, int>(),
+            py::arg("particle_container"), py::arg("level"))
+        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs>::ContainerType&, int, amrex::MFItInfo&>(),
+            py::arg("particle_container"), py::arg("level"),
+            py::arg("info"))
+    ;
+}
+
+void init_WarpXParticleContainer (py::module& m)
+{
+    py::class_<
+        WarpXParticleContainer,
+        amrex::ParticleContainer<0, 0, PIdx::nattribs, 0>
+    > wpc (m, "WarpXParticleContainer");
+    wpc
+        .def("add_real_comp",
+            [](WarpXParticleContainer& pc, const std::string& name, bool const comm) { pc.AddRealComp(name, comm); },
+            py::arg("name"), py::arg("comm")
+        )
+        .def("add_n_particles",
+            [](WarpXParticleContainer& pc, int lev,
+                int n, py::array_t<double> &x,
+                py::array_t<double> &y,
+                py::array_t<double> &z,
+                py::array_t<double> &ux,
+                py::array_t<double> &uy,
+                py::array_t<double> &uz,
+                const int nattr_real, py::array_t<double> &attr_real,
+                const int nattr_int, py::array_t<int> &attr_int,
+                int uniqueparticles, int id
+            ) {
+                amrex::Vector<amrex::ParticleReal> xp(x.data(), x.data() + n);
+                amrex::Vector<amrex::ParticleReal> yp(y.data(), y.data() + n);
+                amrex::Vector<amrex::ParticleReal> zp(z.data(), z.data() + n);
+                amrex::Vector<amrex::ParticleReal> uxp(ux.data(), ux.data() + n);
+                amrex::Vector<amrex::ParticleReal> uyp(uy.data(), uy.data() + n);
+                amrex::Vector<amrex::ParticleReal> uzp(uz.data(), uz.data() + n);
+
+                // create 2d arrays of real and in attributes
+                amrex::Vector<amrex::Vector<amrex::ParticleReal>> attr;
+                const double *attr_data = attr_real.data();
+                for (int ii=0; ii<nattr_real; ii++) {
+                    amrex::Vector<amrex::ParticleReal> attr_ii(n);
+                    for (int jj=0; jj<n; jj++) {
+                        attr_ii[jj] = attr_data[ii + jj*nattr_real];
+                    }
+                    attr.push_back(attr_ii);
+                }
+
+                amrex::Vector<amrex::Vector<int>> iattr;
+                const int *iattr_data = attr_int.data();
+                for (int ii=0; ii<nattr_int; ii++) {
+                    amrex::Vector<int> attr_ii(n);
+                    for (int jj=0; jj<n; jj++) {
+                        attr_ii[jj] = iattr_data[ii + jj*nattr_int];
+                    }
+                    iattr.push_back(attr_ii);
+                }
+
+                pc.AddNParticles(
+                    lev, n, xp, yp, zp, uxp, uyp, uzp, nattr_real, attr,
+                    nattr_int, iattr, uniqueparticles, id
+                );
+            },
+            py::arg("lev"), py::arg("n"),
+            py::arg("x"), py::arg("y"), py::arg("z"),
+            py::arg("ux"), py::arg("uy"), py::arg("uz"),
+            py::arg("nattr_real"), py::arg("attr_real"),
+            py::arg("nattr_int"), py::arg("attr_int"),
+            py::arg("uniqueparticles"), py::arg("id")=-1
+        )
+        .def("num_real_comps", &WarpXParticleContainer::NumRealComps)
+        .def("get_comp_index",
+            [](WarpXParticleContainer& pc, std::string comp_name)
+            {
+                auto particle_comps = pc.getParticleComps();
+                return particle_comps.at(comp_name);
+            },
+            py::arg("comp_name")
+        )
+        .def("num_local_tiles_at_level",
+            &WarpXParticleContainer::numLocalTilesAtLevel,
+            py::arg("level")
+        )
+        .def("total_number_of_particles",
+            &WarpXParticleContainer::TotalNumberOfParticles,
+            py::arg("valid_particles_only"), py::arg("local")
+        )
+        .def("deposit_charge",
+            [](WarpXParticleContainer& pc,
+            amrex::MultiFab* rho, const int lev)
+            {
+                for (WarpXParIter pti(pc, lev); pti.isValid(); ++pti)
+                {
+                    const long np = pti.numParticles();
+                    auto& wp = pti.GetAttribs(PIdx::w);
+                    pc.DepositCharge(pti, wp, nullptr, rho, 0, 0, np, 0, lev, lev);
+                }
+            },
+            py::arg("rho"), py::arg("lev")
+        )
+    ;
+}
diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp
new file mode 100644
index 000000000..9ef5281b3
--- /dev/null
+++ b/Source/Python/WarpX.cpp
@@ -0,0 +1,201 @@
+/* Copyright 2021-2022 The WarpX Community
+ *
+ * Authors: Axel Huebl
+ * License: BSD-3-Clause-LBNL
+ */
+#include "pyWarpX.H"
+
+#include <WarpX.H>
+// see WarpX.cpp - full includes for _fwd.H headers
+#include <BoundaryConditions/PML.H>
+#include <Diagnostics/MultiDiagnostics.H>
+#include <Diagnostics/ReducedDiags/MultiReducedDiags.H>
+#include <EmbeddedBoundary/WarpXFaceInfoBox.H>
+#include <FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H>
+#include <FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H>
+#include <FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H>
+#ifdef WARPX_USE_PSATD
+#   include <FieldSolver/SpectralSolver/SpectralKSpace.H>
+#   ifdef WARPX_DIM_RZ
+#       include <FieldSolver/SpectralSolver/SpectralSolverRZ.H>
+#       include <BoundaryConditions/PML_RZ.H>
+#   else
+#       include <FieldSolver/SpectralSolver/SpectralSolver.H>
+#   endif // RZ ifdef
+#endif // use PSATD ifdef
+#include <FieldSolver/WarpX_FDTD.H>
+#include <Filter/NCIGodfreyFilter.H>
+#include <Particles/MultiParticleContainer.H>
+#include <Particles/ParticleBoundaryBuffer.H>
+#include <AcceleratorLattice/AcceleratorLattice.H>
+#include <Utils/TextMsg.H>
+#include <Utils/WarpXAlgorithmSelection.H>
+#include <Utils/WarpXConst.H>
+#include <Utils/WarpXProfilerWrapper.H>
+#include <Utils/WarpXUtil.H>
+
+#include <AMReX.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_ParallelDescriptor.H>
+
+#if defined(AMREX_DEBUG) || defined(DEBUG)
+#   include <cstdio>
+#endif
+#include <string>
+
+
+//using namespace warpx;
+
+namespace warpx {
+    struct Config {};
+}
+
+void init_WarpX (py::module& m)
+{
+    // Expose the WarpX instance
+    m.def("get_instance",
+        [] () { return &WarpX::GetInstance(); },
+        "Return a reference to the WarpX object.");
+
+    m.def("finalize", &WarpX::Finalize,
+        "Close out the WarpX related data");
+
+    py::class_<WarpX> warpx(m, "WarpX");
+    warpx
+        // WarpX is a Singleton Class with a private constructor
+        //   https://github.com/ECP-WarpX/WarpX/pull/4104
+        //   https://pybind11.readthedocs.io/en/stable/advanced/classes.html?highlight=singleton#custom-constructors
+        .def(py::init([]() {
+            return &WarpX::GetInstance();
+        }))
+        .def_static("get_instance",
+            [] () { return &WarpX::GetInstance(); },
+            "Return a reference to the WarpX object."
+        )
+        .def_static("finalize", &WarpX::Finalize,
+            "Close out the WarpX related data"
+        )
+
+        .def("initialize_data", &WarpX::InitData,
+            "Initializes the WarpX simulation"
+        )
+        .def("evolve", &WarpX::Evolve,
+            "Evolve the simulation the specified number of steps"
+        )
+
+        // from AmrCore->AmrMesh
+        .def("Geom",
+            //[](WarpX const & wx, int const lev) { return wx.Geom(lev); },
+            py::overload_cast< int >(&WarpX::Geom, py::const_),
+            py::arg("lev")
+        )
+        .def("DistributionMap",
+            [](WarpX const & wx, int const lev) { return wx.DistributionMap(lev); },
+            //py::overload_cast< int >(&WarpX::DistributionMap, py::const_),
+            py::arg("lev")
+        )
+        .def("boxArray",
+            [](WarpX const & wx, int const lev) { return wx.boxArray(lev); },
+            //py::overload_cast< int >(&WarpX::boxArray, py::const_),
+            py::arg("lev")
+        )
+        .def("multifab",
+            [](WarpX const & wx, std::string const multifab_name) {
+                if (wx.multifab_map.count(multifab_name) > 0) {
+                    return wx.multifab_map.at(multifab_name);
+                } else {
+                    throw std::runtime_error("The MultiFab '" + multifab_name + "' is unknown or is not allocated!");
+                }
+            },
+            py::arg("multifab_name"),
+            py::return_value_policy::reference_internal,
+            "Return MultiFabs by name, e.g., 'Efield_aux[x][l=0]', 'Efield_cp[x][l=0]', ..."
+        )
+        .def("multi_particle_container",
+            [](WarpX& wx){ return &wx.GetPartContainer(); },
+            py::return_value_policy::reference_internal
+        )
+        .def("get_particle_boundary_buffer",
+            [](WarpX& wx){ return &wx.GetParticleBoundaryBuffer(); },
+            py::return_value_policy::reference_internal
+        )
+
+        // Expose functions used to sync the charge density multifab
+        // accross tiles and apply appropriate boundary conditions
+        .def("sync_rho",
+            [](WarpX& wx){ wx.SyncRho(); }
+        )
+#ifdef WARPX_DIM_RZ
+        .def("apply_inverse_volume_scaling_to_charge_density",
+            [](WarpX& wx, amrex::MultiFab* rho, int const lev) {
+                wx.ApplyInverseVolumeScalingToChargeDensity(rho, lev);
+            },
+            py::arg("rho"), py::arg("lev")
+        )
+#endif
+
+        // Expose functions to get the current simulation step and time
+        .def("getistep",
+            [](WarpX const & wx, int lev){ return wx.getistep(lev); },
+            py::arg("lev")
+        )
+        .def("gett_new",
+            [](WarpX const & wx, int lev){ return wx.gett_new(lev); },
+            py::arg("lev")
+        )
+
+        .def("set_potential_on_eb",
+            [](WarpX& wx, std::string potential) {
+                wx.m_poisson_boundary_handler.setPotentialEB(potential);
+            },
+            py::arg("potential")
+        )
+    ;
+
+    py::class_<warpx::Config>(m, "Config")
+//        .def_property_readonly_static(
+//            "warpx_version",
+//            [](py::object) { return Version(); },
+//            "WarpX version")
+        .def_property_readonly_static(
+            "have_mpi",
+            [](py::object){
+#ifdef AMREX_USE_MPI
+                return true;
+#else
+                return false;
+#endif
+            })
+        .def_property_readonly_static(
+            "have_gpu",
+            [](py::object){
+#ifdef AMREX_USE_GPU
+                return true;
+#else
+                return false;
+#endif
+            })
+        .def_property_readonly_static(
+            "have_omp",
+            [](py::object){
+#ifdef AMREX_USE_OMP
+                return true;
+#else
+                return false;
+#endif
+            })
+        .def_property_readonly_static(
+            "gpu_backend",
+            [](py::object){
+#ifdef AMREX_USE_CUDA
+                return "CUDA";
+#elif defined(AMREX_USE_HIP)
+                return "HIP";
+#elif defined(AMREX_USE_DPCPP)
+                return "SYCL";
+#else
+                return py::none();
+#endif
+            })
+        ;
+}
diff --git a/Source/Python/WarpXWrappers.H b/Source/Python/WarpXWrappers.H
index cd1e41cf4..38ba8b10f 100644
--- a/Source/Python/WarpXWrappers.H
+++ b/Source/Python/WarpXWrappers.H
@@ -3,6 +3,10 @@
  *
  * This file is part of WarpX.
  *
+ * This file is a legacy file and will be removed soon.
+ * Please do NOT add new bindings here! Please see the other files
+ * in this directory for the new pybind11-based bindings!
+ *
  * License: BSD-3-Clause-LBNL
  */
 #ifndef WARPX_WRAPPERS_H_
@@ -30,104 +34,20 @@ extern "C" {
 
     int warpx_galerkin_interpolation();
 
-    int warpx_nComps();
-
-    int warpx_nCompsSpecies(const char* char_species_name);
-
-    int warpx_SpaceDim();
-
-    void amrex_init (int argc, char* argv[]);
-
     void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm);
 
-    void amrex_finalize (int finalize_mpi);
-
-    void warpx_init ();
-
-    void warpx_finalize ();
-
     typedef void(*WARPX_CALLBACK_PY_FUNC_0)();
 
-    void warpx_set_callback_py (const char* char_callback_name,
-                                WARPX_CALLBACK_PY_FUNC_0 callback);
-    void warpx_clear_callback_py (const char* char_callback_name);
-
-    void warpx_evolve (int numsteps);  // -1 means the inputs parameter will be used.
-
-    void warpx_addNParticles(const char* char_species_name,
-                             int lenx,
-                             amrex::ParticleReal const * x,
-                             amrex::ParticleReal const * y,
-                             amrex::ParticleReal const * z,
-                             amrex::ParticleReal const * vx,
-                             amrex::ParticleReal const * vy,
-                             amrex::ParticleReal const * vz,
-                             const int nattr_real,
-                             amrex::ParticleReal const * attr_real,
-                             const int nattr_int,
-                             int const * attr_int,
-                             int uniqueparticles);
-
     void warpx_ConvertLabParamsToBoost();
 
     void warpx_ReadBCParams();
 
     void warpx_CheckGriddingForRZSpectral();
 
-    amrex::Real warpx_getProbLo(int dir);
-
-    amrex::Real warpx_getProbHi(int dir);
-
     amrex::Real warpx_getCellSize(int dir, int lev);
 
-    long warpx_getNumParticles(const char* char_species_name, const bool local);
-
-    void warpx_convert_id_to_long (amrex::Long* ids, const WarpXParticleContainer::ParticleType* pstructs, int size);
-
-    void warpx_convert_cpu_to_int (int* cpus, const WarpXParticleContainer::ParticleType* pstructs, int size);
-
-    amrex::ParticleReal** warpx_getParticleStructs(
-        const char* char_species_name, int lev, int* num_tiles,
-        int** particles_per_tile);
-
-    amrex::ParticleReal** warpx_getParticleArrays(
-        const char* char_species_name, const char* char_comp_name, int lev,
-        int* num_tiles, int** particles_per_tile);
-
-    int warpx_getParticleCompIndex(
-        const char* char_species_name, const char* char_comp_name);
-
-    void warpx_addRealComp(
-        const char* char_species_name, const char* char_comp_name, bool comm);
-
     amrex::Real warpx_sumParticleCharge(const char* char_species_name, const bool local);
 
-    int warpx_getParticleBoundaryBufferSize(const char* species_name, int boundary, bool local);
-
-    int** warpx_getParticleBoundaryBufferScrapedSteps(
-        const char* species_name, int boundary, int lev,
-        int* num_tiles, int** particles_per_tile);
-
-    amrex::ParticleReal** warpx_getParticleBoundaryBuffer(
-        const char* species_name, int boundary, int lev,
-        int* num_tiles, int** particles_per_tile, const char* comp_name);
-
-    amrex::ParticleReal** warpx_getParticleBoundaryBufferStructs(
-            const char* species_name, int boundary, int lev,
-            int* num_tiles, int** particles_per_tile);
-
-    void warpx_clearParticleBoundaryBuffer ();
-
-    /**
-     * \brief This function is used to deposit a given species' charge density
-     * in the rho_fp multifab which can then be accessed from python via
-     * pywarpx.fields.RhoFPWrapper()
-     *
-     * @param[in] species_name specifying the name of the species to deposit
-     * @param[in] lev mesh refinement level
-     */
-    void warpx_depositChargeDensity (const char* species_name, int lev);
-
   void warpx_ComputeDt ();
   void warpx_MoveWindow (int step, bool move_j);
 
@@ -135,7 +55,6 @@ extern "C" {
   void warpx_EvolveB (amrex::Real dt, DtType a_dt_type);
   void warpx_FillBoundaryE ();
   void warpx_FillBoundaryB ();
-  void warpx_SyncRho ();
   void warpx_SyncCurrent (
       const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp,
       const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp,
@@ -143,9 +62,7 @@ extern "C" {
   void warpx_UpdateAuxilaryData ();
   void warpx_PushParticlesandDepose (amrex::Real cur_time);
 
-  int warpx_getistep (int lev);
   void warpx_setistep (int lev, int ii);
-  amrex::Real warpx_gett_new (int lev);
   void warpx_sett_new (int lev, amrex::Real time);
   amrex::Real warpx_getdt (int lev);
 
@@ -154,118 +71,8 @@ extern "C" {
 
   int warpx_finestLevel ();
 
-  int warpx_getMyProc ();
-  int warpx_getNProcs ();
-
-  void warpx_setPotentialEB (const char * char_potential);
-
-  void warpx_setPlasmaLensStrength (const int i_lens, const amrex::Real strength_E, const amrex::Real strength_B);
-
   void mypc_Redistribute ();
 
-  amrex::Real** warpx_getEfield (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getEfieldCP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getEfieldFP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  amrex::Real** warpx_getBfield (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getBfieldCP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getBfieldFP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  amrex::Real** warpx_getCurrentDensity (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getCurrentDensityCP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getCurrentDensityFP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  amrex::Real** warpx_getCurrentDensityFP_Ampere (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  amrex::Real** warpx_getVectorPotentialFP (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  int* warpx_getEfieldLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getEfieldCPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getEfieldFPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  int* warpx_getBfieldLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getBfieldCPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getBfieldFPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  int* warpx_getCurrentDensityLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getCurrentDensityCPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getCurrentDensityFPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  int* warpx_getVectorPotentialFPLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getEdgeLengths (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getEdgeLengthsLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getFaceAreas (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getFaceAreasLoVects (int lev, int direction, int *return_size, int **ngrowvect);
-
-  int* warpx_getEx_nodal_flag ();
-  int* warpx_getEy_nodal_flag ();
-  int* warpx_getEz_nodal_flag ();
-  int* warpx_getBx_nodal_flag ();
-  int* warpx_getBy_nodal_flag ();
-  int* warpx_getBz_nodal_flag ();
-  int* warpx_getAx_nodal_flag ();
-  int* warpx_getAy_nodal_flag ();
-  int* warpx_getAz_nodal_flag ();
-  int* warpx_getJx_nodal_flag ();
-  int* warpx_getJy_nodal_flag ();
-  int* warpx_getJz_nodal_flag ();
-  int* warpx_getRho_nodal_flag ();
-  int* warpx_getPhi_nodal_flag ();
-  int* warpx_getF_nodal_flag ();
-  int* warpx_getG_nodal_flag ();
-  int* warpx_get_edge_lengths_x_nodal_flag ();
-  int* warpx_get_edge_lengths_y_nodal_flag ();
-  int* warpx_get_edge_lengths_z_nodal_flag ();
-  int* warpx_get_face_areas_x_nodal_flag ();
-  int* warpx_get_face_areas_y_nodal_flag ();
-  int* warpx_get_face_areas_z_nodal_flag ();
-
-  amrex::Real** warpx_getChargeDensityCP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getChargeDensityFP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getChargeDensityCPLoVects (int lev, int *return_size, int **ngrowvect);
-  int* warpx_getChargeDensityFPLoVects (int lev, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getPhiFP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  int* warpx_getPhiFPLoVects (int lev, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getFfieldCP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getFfieldFP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getFfieldCPLoVects (int lev, int *return_size, int **ngrowvect);
-  int* warpx_getFfieldFPLoVects (int lev, int *return_size, int **ngrowvect);
-  amrex::Real** warpx_getGfieldCP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getGfieldFP (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getGfieldCPLoVects (int lev, int *return_size, int **ngrowvect);
-  int* warpx_getGfieldFPLoVects (int lev, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getEfieldCP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getEfieldFP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getBfieldCP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getBfieldFP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getCurrentDensityCP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getCurrentDensityFP_PML (int lev, int direction, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-
-  int* warpx_getEfieldCPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getEfieldFPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getBfieldCPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getBfieldFPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getCurrentDensityCPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-  int* warpx_getCurrentDensityFPLoVects_PML (int lev, int direction, int *return_size, int **ngrowvect);
-
-  amrex::Real** warpx_getFfieldCP_PML (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getFfieldFP_PML (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getFfieldCPLoVects_PML (int lev, int *return_size, int **ngrowvect);
-  int* warpx_getFfieldFPLoVects_PML (int lev, int *return_size, int **ngrowvect);
-  amrex::Real** warpx_getGfieldCP_PML (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  amrex::Real** warpx_getGfieldFP_PML (int lev, int *return_size, int *ncomps, int **ngrowvect, int **shapes);
-  int* warpx_getGfieldCPLoVects_PML (int lev, int *return_size, int **ngrowvect);
-  int* warpx_getGfieldFPLoVects_PML (int lev, int *return_size, int **ngrowvect);
-
-  int* warpx_getF_pml_nodal_flag ();
-  int* warpx_getG_pml_nodal_flag ();
-
 #ifdef __cplusplus
 }
 #endif
diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp
index 02cee52e0..44820dd1f 100644
--- a/Source/Python/WarpXWrappers.cpp
+++ b/Source/Python/WarpXWrappers.cpp
@@ -39,71 +39,6 @@
 #include <array>
 #include <cstdlib>
 
-namespace
-{
-    amrex::Real** getMultiFabPointers (amrex::MultiFab& mf, int *num_boxes, int *ncomps, int **ngrowvect, int **shapes)
-    {
-        *ncomps = mf.nComp();
-        *num_boxes = mf.local_size();
-        int shapesize = AMREX_SPACEDIM;
-        *ngrowvect = static_cast<int*>(malloc(sizeof(int)*shapesize));
-        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
-            (*ngrowvect)[j] = mf.nGrow(j);
-        }
-        if (mf.nComp() > 1) shapesize += 1;
-        *shapes = static_cast<int*>(malloc(sizeof(int)*shapesize * (*num_boxes)));
-        auto data =
-            static_cast<amrex::Real**>(malloc((*num_boxes) * sizeof(amrex::Real*)));
-
-#ifdef AMREX_USE_OMP
-#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
-#endif
-        for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi ) {
-            const int i = mfi.LocalIndex();
-            data[i] = mf[mfi].dataPtr();
-            for (int j = 0; j < AMREX_SPACEDIM; ++j) {
-                (*shapes)[shapesize*i+j] = mf[mfi].box().length(j);
-            }
-            if (mf.nComp() > 1) (*shapes)[shapesize*i+AMREX_SPACEDIM] = mf.nComp();
-        }
-        return data;
-    }
-    int* getMultiFabLoVects (const amrex::MultiFab& mf, int *num_boxes, int **ngrowvect)
-    {
-        const int shapesize = AMREX_SPACEDIM;
-        *ngrowvect = static_cast<int*>(malloc(sizeof(int)*shapesize));
-        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
-            (*ngrowvect)[j] = mf.nGrow(j);
-        }
-        *num_boxes = mf.local_size();
-        auto loVects = static_cast<int*>(malloc((*num_boxes)*AMREX_SPACEDIM * sizeof(int)));
-
-        int i = 0;
-        for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) {
-            const int* loVect = mf[mfi].loVect();
-            for (int j = 0; j < AMREX_SPACEDIM; ++j) {
-                loVects[AMREX_SPACEDIM*i+j] = loVect[j];
-            }
-        }
-        return loVects;
-    }
-    // Copy the nodal flag data and return the copy:
-    // the nodal flag data should not be modifiable from Python.
-    int* getFieldNodalFlagData ( const amrex::MultiFab* mf )
-    {
-        if (mf == nullptr) return nullptr;
-        const amrex::IntVect nodal_flag( mf->ixType().toIntVect() );
-        auto *nodal_flag_data = static_cast<int*>(malloc(AMREX_SPACEDIM * sizeof(int)));
-
-        constexpr int NODE = amrex::IndexType::NODE;
-
-        for (int i=0 ; i < AMREX_SPACEDIM ; i++) {
-            nodal_flag_data[i] = (nodal_flag[i] == NODE ? 1 : 0);
-        }
-        return nodal_flag_data;
-    }
-}
-
     int warpx_Real_size()
     {
         return (int)sizeof(amrex::Real);
@@ -130,85 +65,9 @@ namespace
         return WarpX::galerkin_interpolation;
     }
 
-    int warpx_nComps()
-    {
-        return PIdx::nattribs;
-    }
-
-    int warpx_nCompsSpecies(const char* char_species_name)
-    {
-        auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-        return myspc.NumRealComps();
-    }
-
-    int warpx_SpaceDim()
-    {
-        return AMREX_SPACEDIM;
-    }
-
-    void amrex_init (int argc, char* argv[])
-    {
-        warpx::initialization::amrex_init(argc, argv);
-    }
-
-    void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm)
-    {
-        warpx::initialization::amrex_init(argc, argv, true, mpicomm);
-    }
-
-    void amrex_finalize (int /*finalize_mpi*/)
-    {
-        amrex::Finalize();
-    }
-
-    void warpx_init ()
-    {
-        WarpX& warpx = WarpX::GetInstance();
-        warpx.InitData();
-        ExecutePythonCallback("afterinit");
-        ExecutePythonCallback("particleloader");
-    }
-
-    void warpx_finalize ()
+    void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm /* mpicomm */)
     {
-        WarpX::Finalize();
-    }
-
-    void warpx_set_callback_py (
-        const char* char_callback_name, WARPX_CALLBACK_PY_FUNC_0 callback)
-    {
-        const std::string callback_name(char_callback_name);
-        warpx_callback_py_map[callback_name] = callback;
-    }
-
-    void warpx_clear_callback_py (const char* char_callback_name)
-    {
-        const std::string callback_name(char_callback_name);
-        warpx_callback_py_map.erase(callback_name);
-    }
-
-    void warpx_evolve (int numsteps)
-    {
-        WarpX& warpx = WarpX::GetInstance();
-        warpx.Evolve(numsteps);
-    }
-
-    void warpx_addNParticles(
-        const char* char_species_name, int lenx, amrex::ParticleReal const * x,
-        amrex::ParticleReal const * y, amrex::ParticleReal const * z,
-        amrex::ParticleReal const * vx, amrex::ParticleReal const * vy,
-        amrex::ParticleReal const * vz, const int nattr_real,
-        amrex::ParticleReal const * attr_real, const int nattr_int,
-        int const * attr_int, int uniqueparticles)
-    {
-        auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-        const int lev = 0;
-        myspc.AddNParticles(lev, lenx, x, y, z, vx, vy, vz, nattr_real, attr_real,
-                            nattr_int, attr_int, uniqueparticles);
+      warpx::initialization::amrex_init(argc, argv, true);
     }
 
     void warpx_ConvertLabParamsToBoost()
@@ -226,340 +85,11 @@ namespace
         CheckGriddingForRZSpectral();
     }
 
-    amrex::Real warpx_getProbLo(int dir)
-    {
-      WarpX& warpx = WarpX::GetInstance();
-      const amrex::Geometry& geom = warpx.Geom(0);
-      return geom.ProbLo(dir);
-    }
-
-    amrex::Real warpx_getProbHi(int dir)
-    {
-      WarpX& warpx = WarpX::GetInstance();
-      const amrex::Geometry& geom = warpx.Geom(0);
-      return geom.ProbHi(dir);
-    }
-
     amrex::Real warpx_getCellSize(int dir, int lev) {
         const std::array<amrex::Real,3>& dx = WarpX::CellSize(lev);
         return dx[dir];
     }
 
-    long warpx_getNumParticles(const char* char_species_name, const bool local) {
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-        // the first argument below is to only count valid particles
-        return myspc.TotalNumberOfParticles(true, local);
-    }
-
-#define WARPX_GET_FIELD(FIELD, GETTER) \
-    amrex::Real** FIELD(int lev, int direction, \
-                        int *return_size, int *ncomps, int **ngrowvect, int **shapes) { \
-        auto * mf = GETTER(lev, direction); \
-        if (mf != nullptr) { \
-            return getMultiFabPointers(*mf, return_size, ncomps, ngrowvect, shapes); \
-        } else { \
-            return nullptr; \
-        } \
-    }
-
-#define WARPX_GET_LOVECTS(FIELD, GETTER) \
-    int* FIELD(int lev, int direction, \
-               int *return_size, int **ngrowvect) { \
-        auto * mf = GETTER(lev, direction); \
-        if (mf != nullptr) { \
-            return getMultiFabLoVects(*mf, return_size, ngrowvect); \
-        } else { \
-            return nullptr; \
-        } \
-    }
-
-    WARPX_GET_FIELD(warpx_getEfield, WarpX::GetInstance().get_pointer_Efield_aux)
-    WARPX_GET_FIELD(warpx_getEfieldCP, WarpX::GetInstance().get_pointer_Efield_cp)
-    WARPX_GET_FIELD(warpx_getEfieldFP, WarpX::GetInstance().get_pointer_Efield_fp)
-
-    WARPX_GET_FIELD(warpx_getBfield, WarpX::GetInstance().get_pointer_Bfield_aux)
-    WARPX_GET_FIELD(warpx_getBfieldCP, WarpX::GetInstance().get_pointer_Bfield_cp)
-    WARPX_GET_FIELD(warpx_getBfieldFP, WarpX::GetInstance().get_pointer_Bfield_fp)
-
-    WARPX_GET_FIELD(warpx_getEdgeLengths, WarpX::GetInstance().get_pointer_edge_lengths)
-    WARPX_GET_FIELD(warpx_getFaceAreas, WarpX::GetInstance().get_pointer_face_areas)
-
-    WARPX_GET_FIELD(warpx_getCurrentDensity, WarpX::GetInstance().get_pointer_current_fp)
-    WARPX_GET_FIELD(warpx_getCurrentDensityCP, WarpX::GetInstance().get_pointer_current_cp)
-    WARPX_GET_FIELD(warpx_getCurrentDensityFP, WarpX::GetInstance().get_pointer_current_fp)
-
-    WARPX_GET_FIELD(warpx_getCurrentDensityFP_Ampere, WarpX::GetInstance().GetHybridPICModel().get_pointer_current_fp_ampere)
-
-    WARPX_GET_FIELD(warpx_getVectorPotentialFP, WarpX::GetInstance().get_pointer_vector_potential_fp)
-
-    WARPX_GET_LOVECTS(warpx_getEfieldLoVects, WarpX::GetInstance().get_pointer_Efield_aux)
-    WARPX_GET_LOVECTS(warpx_getEfieldCPLoVects, WarpX::GetInstance().get_pointer_Efield_cp)
-    WARPX_GET_LOVECTS(warpx_getEfieldFPLoVects, WarpX::GetInstance().get_pointer_Efield_fp)
-
-    WARPX_GET_LOVECTS(warpx_getBfieldLoVects, WarpX::GetInstance().get_pointer_Bfield_aux)
-    WARPX_GET_LOVECTS(warpx_getBfieldCPLoVects, WarpX::GetInstance().get_pointer_Bfield_cp)
-    WARPX_GET_LOVECTS(warpx_getBfieldFPLoVects, WarpX::GetInstance().get_pointer_Bfield_fp)
-
-    WARPX_GET_LOVECTS(warpx_getCurrentDensityLoVects, WarpX::GetInstance().get_pointer_current_fp)
-    WARPX_GET_LOVECTS(warpx_getCurrentDensityCPLoVects, WarpX::GetInstance().get_pointer_current_cp)
-    WARPX_GET_LOVECTS(warpx_getCurrentDensityFPLoVects, WarpX::GetInstance().get_pointer_current_fp)
-    WARPX_GET_LOVECTS(warpx_getVectorPotentialFPLoVects, WarpX::GetInstance().get_pointer_vector_potential_fp)
-
-    WARPX_GET_LOVECTS(warpx_getEdgeLengthsLoVects, WarpX::GetInstance().get_pointer_edge_lengths)
-    WARPX_GET_LOVECTS(warpx_getFaceAreasLoVects, WarpX::GetInstance().get_pointer_face_areas)
-
-    int* warpx_getEx_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Efield_aux(0,0) );}
-    int* warpx_getEy_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Efield_aux(0,1) );}
-    int* warpx_getEz_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Efield_aux(0,2) );}
-    int* warpx_getBx_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Bfield_aux(0,0) );}
-    int* warpx_getBy_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Bfield_aux(0,1) );}
-    int* warpx_getBz_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_Bfield_aux(0,2) );}
-    int* warpx_getJx_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_current_fp(0,0) );}
-    int* warpx_getJy_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_current_fp(0,1) );}
-    int* warpx_getJz_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_current_fp(0,2) );}
-    int* warpx_getAx_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_vector_potential_fp(0,0) );}
-    int* warpx_getAy_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_vector_potential_fp(0,1) );}
-    int* warpx_getAz_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_vector_potential_fp(0,2) );}
-    int* warpx_getRho_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_rho_fp(0) );}
-    int* warpx_getPhi_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_phi_fp(0) );}
-    int* warpx_getF_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_F_fp(0) );}
-    int* warpx_getG_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_G_fp(0) );}
-    int* warpx_get_edge_lengths_x_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_edge_lengths(0, 0) );}
-    int* warpx_get_edge_lengths_y_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_edge_lengths(0, 1) );}
-    int* warpx_get_edge_lengths_z_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_edge_lengths(0, 2) );}
-    int* warpx_get_face_areas_x_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_face_areas(0, 0) );}
-    int* warpx_get_face_areas_y_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_face_areas(0, 1) );}
-    int* warpx_get_face_areas_z_nodal_flag() {return getFieldNodalFlagData( WarpX::GetInstance().get_pointer_face_areas(0, 2) );}
-
-#define WARPX_GET_SCALAR(SCALAR, GETTER) \
-    amrex::Real** SCALAR(int lev, \
-                         int *return_size, int *ncomps, int **ngrowvect, int **shapes) { \
-        auto * mf = GETTER(lev); \
-        if (mf != nullptr) { \
-            return getMultiFabPointers(*mf, return_size, ncomps, ngrowvect, shapes); \
-        } else { \
-            return nullptr; \
-        } \
-    }
-
-#define WARPX_GET_LOVECTS_SCALAR(SCALAR, GETTER) \
-    int* SCALAR(int lev, \
-                int *return_size, int **ngrowvect) { \
-        auto * mf = GETTER(lev); \
-        if (mf != nullptr) { \
-            return getMultiFabLoVects(*mf, return_size, ngrowvect); \
-        } else { \
-            return nullptr; \
-        } \
-    }
-
-    WARPX_GET_SCALAR(warpx_getChargeDensityCP, WarpX::GetInstance().get_pointer_rho_cp)
-    WARPX_GET_SCALAR(warpx_getChargeDensityFP, WarpX::GetInstance().get_pointer_rho_fp)
-
-    WARPX_GET_LOVECTS_SCALAR(warpx_getChargeDensityCPLoVects, WarpX::GetInstance().get_pointer_rho_cp)
-    WARPX_GET_LOVECTS_SCALAR(warpx_getChargeDensityFPLoVects, WarpX::GetInstance().get_pointer_rho_fp)
-
-    WARPX_GET_SCALAR(warpx_getPhiFP, WarpX::GetInstance().get_pointer_phi_fp)
-
-    WARPX_GET_LOVECTS_SCALAR(warpx_getPhiFPLoVects, WarpX::GetInstance().get_pointer_phi_fp)
-
-    // F and G
-    WARPX_GET_SCALAR(warpx_getFfieldCP, WarpX::GetInstance().get_pointer_F_cp)
-    WARPX_GET_SCALAR(warpx_getFfieldFP, WarpX::GetInstance().get_pointer_F_fp)
-    WARPX_GET_LOVECTS_SCALAR(warpx_getFfieldCPLoVects, WarpX::GetInstance().get_pointer_F_cp)
-    WARPX_GET_LOVECTS_SCALAR(warpx_getFfieldFPLoVects, WarpX::GetInstance().get_pointer_F_fp)
-    WARPX_GET_SCALAR(warpx_getGfieldCP, WarpX::GetInstance().get_pointer_G_cp)
-    WARPX_GET_SCALAR(warpx_getGfieldFP, WarpX::GetInstance().get_pointer_G_fp)
-    WARPX_GET_LOVECTS_SCALAR(warpx_getGfieldCPLoVects, WarpX::GetInstance().get_pointer_G_cp)
-    WARPX_GET_LOVECTS_SCALAR(warpx_getGfieldFPLoVects, WarpX::GetInstance().get_pointer_G_fp)
-
-#define WARPX_GET_FIELD_PML(FIELD, GETTER) \
-    amrex::Real** FIELD(int lev, int direction, \
-                        int *return_size, int *ncomps, int **ngrowvect, int **shapes) { \
-        auto * pml = WarpX::GetInstance().GetPML(lev); \
-        if (!pml) return nullptr; \
-        auto * mf = (pml->GETTER()[direction]); \
-        if (!mf) return nullptr; \
-        return getMultiFabPointers(*mf, return_size, ncomps, ngrowvect, shapes); \
-    }
-
-#define WARPX_GET_LOVECTS_PML(FIELD, GETTER) \
-    int* FIELD(int lev, int direction, \
-               int *return_size, int **ngrowvect) { \
-        auto * pml = WarpX::GetInstance().GetPML(lev); \
-        if (!pml) return nullptr; \
-        auto * mf = (pml->GETTER()[direction]); \
-        if (!mf) return nullptr; \
-        return getMultiFabLoVects(*mf, return_size, ngrowvect); \
-    }
-
-    WARPX_GET_FIELD_PML(warpx_getEfieldCP_PML, GetE_cp)
-    WARPX_GET_FIELD_PML(warpx_getEfieldFP_PML, GetE_fp)
-    WARPX_GET_FIELD_PML(warpx_getBfieldCP_PML, GetB_cp)
-    WARPX_GET_FIELD_PML(warpx_getBfieldFP_PML, GetB_fp)
-    WARPX_GET_FIELD_PML(warpx_getCurrentDensityCP_PML, Getj_cp)
-    WARPX_GET_FIELD_PML(warpx_getCurrentDensityFP_PML, Getj_fp)
-    WARPX_GET_LOVECTS_PML(warpx_getEfieldCPLoVects_PML, GetE_cp)
-    WARPX_GET_LOVECTS_PML(warpx_getEfieldFPLoVects_PML, GetE_fp)
-    WARPX_GET_LOVECTS_PML(warpx_getBfieldCPLoVects_PML, GetB_cp)
-    WARPX_GET_LOVECTS_PML(warpx_getBfieldFPLoVects_PML, GetB_fp)
-    WARPX_GET_LOVECTS_PML(warpx_getCurrentDensityCPLoVects_PML, Getj_cp)
-    WARPX_GET_LOVECTS_PML(warpx_getCurrentDensityFPLoVects_PML, Getj_fp)
-
-#define WARPX_GET_SCALAR_PML(SCALAR, GETTER) \
-    amrex::Real** SCALAR(int lev, \
-                        int *return_size, int *ncomps, int **ngrowvect, int **shapes) { \
-        auto * pml = WarpX::GetInstance().GetPML(lev); \
-        if (!pml) return nullptr; \
-        auto * mf = pml->GETTER(); \
-        if (!mf) return nullptr; \
-        return getMultiFabPointers(*mf, return_size, ncomps, ngrowvect, shapes); \
-    }
-
-#define WARPX_GET_LOVECTS_PML_SCALAR(SCALAR, GETTER) \
-    int* SCALAR(int lev, \
-               int *return_size, int **ngrowvect) { \
-        auto * pml = WarpX::GetInstance().GetPML(lev); \
-        if (!pml) return nullptr; \
-        auto * mf = pml->GETTER(); \
-        if (!mf) return nullptr; \
-        return getMultiFabLoVects(*mf, return_size, ngrowvect); \
-    }
-
-    // F and G
-    WARPX_GET_SCALAR_PML(warpx_getFfieldCP_PML, GetF_cp)
-    WARPX_GET_SCALAR_PML(warpx_getFfieldFP_PML, GetF_fp)
-    WARPX_GET_LOVECTS_PML_SCALAR(warpx_getFfieldCPLoVects_PML, GetF_cp)
-    WARPX_GET_LOVECTS_PML_SCALAR(warpx_getFfieldFPLoVects_PML, GetF_fp)
-    WARPX_GET_SCALAR_PML(warpx_getGfieldCP_PML, GetG_cp)
-    WARPX_GET_SCALAR_PML(warpx_getGfieldFP_PML, GetG_fp)
-    WARPX_GET_LOVECTS_PML_SCALAR(warpx_getGfieldCPLoVects_PML, GetG_cp)
-    WARPX_GET_LOVECTS_PML_SCALAR(warpx_getGfieldFPLoVects_PML, GetG_fp)
-
-    int* warpx_getF_pml_nodal_flag ()
-    {
-        auto * pml = WarpX::GetInstance().GetPML(0);
-        if (!pml) return nullptr;
-        return getFieldNodalFlagData(pml->GetF_fp());
-    }
-
-    int* warpx_getG_pml_nodal_flag ()
-    {
-        auto * pml = WarpX::GetInstance().GetPML(0);
-        if (!pml) return nullptr;
-        return getFieldNodalFlagData(pml->GetG_fp());
-    }
-
-    amrex::ParticleReal** warpx_getParticleStructs(
-            const char* char_species_name, int lev,
-            int* num_tiles, int** particles_per_tile) {
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-
-        *num_tiles = myspc.numLocalTilesAtLevel(lev);
-        *particles_per_tile = static_cast<int*>(malloc(*num_tiles*sizeof(int)));
-        memset(*particles_per_tile, 0, *num_tiles*sizeof(int));
-
-        auto data = static_cast<amrex::ParticleReal**>(malloc(*num_tiles*sizeof(typename WarpXParticleContainer::ParticleType*)));
-        int i = 0;
-        for (WarpXParIter pti(myspc, lev); pti.isValid(); ++pti, ++i) {
-            auto& aos = pti.GetArrayOfStructs();
-            data[i] = (amrex::ParticleReal*) aos.data();
-            (*particles_per_tile)[i] = pti.numParticles();
-        }
-        return data;
-    }
-
-    amrex::ParticleReal** warpx_getParticleArrays (
-            const char* char_species_name, const char* char_comp_name,
-            int lev, int* num_tiles, int** particles_per_tile ) {
-
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-
-        const int comp = warpx_getParticleCompIndex(char_species_name, char_comp_name);
-
-        *num_tiles = myspc.numLocalTilesAtLevel(lev);
-        *particles_per_tile = static_cast<int*>(malloc(*num_tiles*sizeof(int)));
-        memset(*particles_per_tile, 0, *num_tiles*sizeof(int));
-
-        auto data = static_cast<amrex::ParticleReal**>(malloc(*num_tiles*sizeof(amrex::ParticleReal*)));
-        int i = 0;
-        for (WarpXParIter pti(myspc, lev); pti.isValid(); ++pti, ++i) {
-            auto& soa = pti.GetStructOfArrays();
-            data[i] = (amrex::ParticleReal*) soa.GetRealData(comp).dataPtr();
-            (*particles_per_tile)[i] = pti.numParticles();
-        }
-        return data;
-    }
-
-    void warpx_convert_id_to_long (amrex::Long* ids, const WarpXParticleContainer::ParticleType* pstructs, int size)
-    {
-        amrex::Long* d_ptr = nullptr;
-#ifdef AMREX_USE_GPU
-        amrex::Gpu::DeviceVector<amrex::Long> d_ids(size);
-        d_ptr = d_ids.data();
-#else
-        d_ptr = ids;
-#endif
-        amrex::ParallelFor(size, [=] AMREX_GPU_DEVICE (int i) noexcept
-        {
-            d_ptr[i] = pstructs[i].id();
-        });
-#ifdef AMREX_USE_GPU
-        amrex::Gpu::dtoh_memcpy(ids, d_ptr, size*sizeof(amrex::Long));
-#endif
-    }
-
-    void warpx_convert_cpu_to_int (int* cpus, const WarpXParticleContainer::ParticleType* pstructs, int size)
-    {
-        int* d_ptr = nullptr;
-#ifdef AMREX_USE_GPU
-        amrex::Gpu::DeviceVector<int> d_cpus(size);
-        d_ptr = d_cpus.data();
-#else
-        d_ptr = cpus;
-#endif
-        amrex::ParallelFor(size, [=] AMREX_GPU_DEVICE (int i) noexcept
-        {
-            d_ptr[i] = pstructs[i].cpu();
-        });
-#ifdef AMREX_USE_GPU
-        amrex::Gpu::dtoh_memcpy(cpus, d_ptr, size*sizeof(int));
-#endif
-    }
-
-    int warpx_getParticleCompIndex (
-         const char* char_species_name, const char* char_comp_name )
-    {
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-
-        const std::string comp_name(char_comp_name);
-        auto particle_comps = myspc.getParticleComps();
-
-        return particle_comps.at(comp_name);
-    }
-
-    void warpx_addRealComp(const char* char_species_name,
-        const char* char_comp_name, bool comm=true)
-    {
-        auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-
-        const std::string comp_name(char_comp_name);
-        myspc.AddRealComp(comp_name, comm);
-
-        mypc.defineAllParticleTiles();
-    }
-
     amrex::Real warpx_sumParticleCharge(const char* char_species_name, const bool local)
     {
         auto & mypc = WarpX::GetInstance().GetPartContainer();
@@ -568,117 +98,6 @@ namespace
         return myspc.sumParticleCharge(local);
     }
 
-    int warpx_getParticleBoundaryBufferSize(const char* species_name, int boundary, bool local)
-    {
-        const std::string name(species_name);
-        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
-        return particle_buffers.getNumParticlesInContainer(species_name, boundary, local);
-    }
-
-    int** warpx_getParticleBoundaryBufferScrapedSteps(const char* species_name, int boundary, int lev,
-                     int* num_tiles, int** particles_per_tile)
-    {
-        const std::string name(species_name);
-        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
-        auto& particle_buffer = particle_buffers.getParticleBuffer(species_name, boundary);
-
-        const int comp = particle_buffer.NumIntComps() - 1;
-
-        *num_tiles = particle_buffer.numLocalTilesAtLevel(lev);
-        *particles_per_tile = static_cast<int*>(malloc(*num_tiles*sizeof(int)));
-        memset(*particles_per_tile, 0, *num_tiles*sizeof(int));
-
-        auto data = static_cast<int**>(malloc(*num_tiles*sizeof(int*)));
-        int i = 0;
-        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {
-            auto& soa = pti.GetStructOfArrays();
-            data[i] = (int*) soa.GetIntData(comp).dataPtr();
-            (*particles_per_tile)[i] = pti.numParticles();
-        }
-
-        return data;
-    }
-
-    amrex::ParticleReal** warpx_getParticleBoundaryBuffer(const char* species_name, int boundary, int lev,
-                     int* num_tiles, int** particles_per_tile, const char* comp_name)
-    {
-        const std::string name(species_name);
-        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
-        auto& particle_buffer = particle_buffers.getParticleBuffer(species_name, boundary);
-
-        const int comp = warpx_getParticleCompIndex(species_name, comp_name);
-
-        *num_tiles = particle_buffer.numLocalTilesAtLevel(lev);
-        *particles_per_tile = static_cast<int*>(malloc(*num_tiles*sizeof(int)));
-        memset(*particles_per_tile, 0, *num_tiles*sizeof(int));
-
-        auto data = static_cast<amrex::ParticleReal**>(malloc(*num_tiles*sizeof(amrex::ParticleReal*)));
-        int i = 0;
-        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {
-            auto& soa = pti.GetStructOfArrays();
-            data[i] = (amrex::ParticleReal*) soa.GetRealData(comp).dataPtr();
-            (*particles_per_tile)[i] = pti.numParticles();
-        }
-
-        return data;
-    }
-
-    amrex::ParticleReal** warpx_getParticleBoundaryBufferStructs(const char* species_name, int boundary, int lev,
-                     int* num_tiles, int** particles_per_tile)
-    {
-        const std::string name(species_name);
-        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
-        auto& particle_buffer = particle_buffers.getParticleBuffer(species_name, boundary);
-
-        *num_tiles = particle_buffer.numLocalTilesAtLevel(lev);
-        *particles_per_tile = static_cast<int*>(malloc(*num_tiles*sizeof(int)));
-        memset(*particles_per_tile, 0, *num_tiles*sizeof(int));
-
-        auto data = static_cast<amrex::ParticleReal**>(malloc(*num_tiles*sizeof(typename WarpXParticleContainer::ParticleType*)));
-        int i = 0;
-        for (amrex::ParIter<0,0,PIdx::nattribs, 0, amrex::PinnedArenaAllocator> pti(particle_buffer, lev); pti.isValid(); ++pti, ++i) {
-            auto& aos = pti.GetArrayOfStructs();
-            data[i] = (amrex::ParticleReal*) aos.data();
-            (*particles_per_tile)[i] = pti.numParticles();
-        }
-        return data;
-    }
-
-    void warpx_clearParticleBoundaryBuffer () {
-        auto& particle_buffers = WarpX::GetInstance().GetParticleBoundaryBuffer();
-        particle_buffers.clearParticles();
-    }
-
-    void warpx_depositChargeDensity (const char* char_species_name, int lev) {
-        // this function is used to deposit a given species' charge density
-        // in the rho_fp multifab which can then be accessed from python via
-        // pywarpx.fields.RhoFPWrapper()
-        WarpX& warpx = WarpX::GetInstance();
-        const auto & mypc = warpx.GetPartContainer();
-        const std::string species_name(char_species_name);
-        auto & myspc = mypc.GetParticleContainerFromName(species_name);
-        auto * rho_fp = warpx.get_pointer_rho_fp(lev);
-
-        if (rho_fp == nullptr) {
-            ablastr::warn_manager::WMRecordWarning(
-                "WarpXWrappers", "rho_fp is not allocated",
-                ablastr::warn_manager::WarnPriority::low
-            );
-            return;
-        }
-
-        for (WarpXParIter pti(myspc, lev); pti.isValid(); ++pti)
-        {
-            const long np = pti.numParticles();
-            auto& wp = pti.GetAttribs(PIdx::w);
-            // Do this unconditionally, ignoring myspc.do_not_deposit, to support diagnostic uses
-            myspc.DepositCharge(pti, wp, nullptr, rho_fp, 0, 0, np, 0, lev, lev);
-        }
-#ifdef WARPX_DIM_RZ
-        warpx.ApplyInverseVolumeScalingToChargeDensity(rho_fp, lev);
-#endif
-    }
-
     void warpx_ComputeDt () {
         WarpX& warpx = WarpX::GetInstance();
         warpx.ComputeDt();
@@ -704,10 +123,6 @@ namespace
         WarpX& warpx = WarpX::GetInstance();
         warpx.FillBoundaryB(warpx.getngEB());
     }
-    void warpx_SyncRho () {
-        WarpX& warpx = WarpX::GetInstance();
-        warpx.SyncRho();
-    }
     void warpx_SyncCurrent (
         const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp,
         const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp,
@@ -724,18 +139,10 @@ namespace
         warpx.PushParticlesandDepose(cur_time);
     }
 
-    int warpx_getistep (int lev) {
-        const WarpX& warpx = WarpX::GetInstance();
-        return warpx.getistep(lev);
-    }
     void warpx_setistep (int lev, int ii) {
         WarpX& warpx = WarpX::GetInstance();
         warpx.setistep(lev, ii);
     }
-    amrex::Real warpx_gett_new (int lev) {
-        const WarpX& warpx = WarpX::GetInstance();
-        return warpx.gett_new(lev);
-    }
     void warpx_sett_new (int lev, amrex::Real time) {
         WarpX& warpx = WarpX::GetInstance();
         warpx.sett_new(lev, time);
@@ -759,35 +166,6 @@ namespace
         return warpx.finestLevel();
     }
 
-    int warpx_getMyProc () {
-        return amrex::ParallelDescriptor::MyProc();
-    }
-
-    int warpx_getNProcs () {
-        return amrex::ParallelDescriptor::NProcs();
-    }
-
-    void warpx_setPotentialEB (const char * char_potential) {
-        WarpX& warpx = WarpX::GetInstance();
-        const std::string potential(char_potential);
-        warpx.m_poisson_boundary_handler.setPotentialEB(potential);
-    }
-
-    void warpx_setPlasmaLensStrength (const int i_lens, const amrex::Real strength_E, const amrex::Real strength_B ) {
-        auto & mypc = WarpX::GetInstance().GetPartContainer();
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE( i_lens >= mypc.h_repeated_plasma_lens_strengths_E.size(),
-            "i_lens is larger than the number of plasma lenses.");
-        mypc.h_repeated_plasma_lens_strengths_E[i_lens] = strength_E;
-        mypc.h_repeated_plasma_lens_strengths_B[i_lens] = strength_B;
-        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
-                       mypc.h_repeated_plasma_lens_strengths_E.begin(), mypc.h_repeated_plasma_lens_strengths_E.end(),
-                       mypc.d_repeated_plasma_lens_strengths_E.begin());
-        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
-                       mypc.h_repeated_plasma_lens_strengths_B.begin(), mypc.h_repeated_plasma_lens_strengths_B.end(),
-                       mypc.d_repeated_plasma_lens_strengths_B.begin());
-        amrex::Gpu::synchronize();
-    }
-
     void mypc_Redistribute () {
         auto & mypc = WarpX::GetInstance().GetPartContainer();
         mypc.Redistribute();
diff --git a/Source/Python/WarpX_py.H b/Source/Python/WarpX_py.H
index 83fc91f21..a18f9a7f1 100644
--- a/Source/Python/WarpX_py.H
+++ b/Source/Python/WarpX_py.H
@@ -1,40 +1,48 @@
-/* Copyright 2019 David Grote, Maxence Thevenet, Weiqun Zhang
- *
+/* Copyright 2019-2022 The WarpX Community
  *
  * This file is part of WarpX.
  *
+ * Authors: David Grote, Maxence Thevenet, Weiqun Zhang, Roelof Groenewald
+ *
  * License: BSD-3-Clause-LBNL
  */
 #ifndef WARPX_PY_H_
 #define WARPX_PY_H_
 
-#include "WarpXWrappers.H"
+#include "Utils/export.H"
 #include "Utils/WarpXProfilerWrapper.H"
 
 #include <map>
 #include <string>
 
+
 /**
  * Declare global map to hold python callback functions.
  *
  * The keys of the map describe at what point in the simulation the python
- * functions will be called. Currently supported keys (callback points) are
- * afterinit, beforecollisions, aftercollisions, beforeEsolve, poissonsolver,
- * afterEsolve, beforeInitEsolve, afterInitEsolve,
- * beforedeposition, afterdeposition, particlescraper,
- * particleloader, beforestep, afterstep, afterdiagnostics, afterrestart,
- * oncheckpointsignal, onbreaksignal, particleinjection and appliedfields.
+ * functions will be called. See ``WarpX/Python/pywarpx/callbacks.py`` for a
+ * list of currently supported callback names.
 */
-extern std::map< std::string, WARPX_CALLBACK_PY_FUNC_0 > warpx_callback_py_map;
+extern WARPX_EXPORT std::map< std::string, std::function<void()> > warpx_callback_py_map;
+
+/**
+ * \brief Function to install the given name and function in warpx_callback_py_map
+ */
+void InstallPythonCallback ( std::string name, std::function<void()> callback );
 
 /**
  * \brief Function to check if the given name is a key in warpx_callback_py_map
  */
-bool IsPythonCallBackInstalled ( std::string name );
+bool IsPythonCallbackInstalled ( std::string name );
 
 /**
  * \brief Function to look for and execute Python callbacks
  */
 void ExecutePythonCallback ( std::string name );
 
+/**
+ * \brief Function to clear the given callback name from warpx_callback_py_map
+ */
+void ClearPythonCallback ( std::string name );
+
 #endif
diff --git a/Source/Python/WarpX_py.cpp b/Source/Python/WarpX_py.cpp
index 7e0e491f8..49328d0a4 100644
--- a/Source/Python/WarpX_py.cpp
+++ b/Source/Python/WarpX_py.cpp
@@ -1,15 +1,21 @@
-/* Copyright 2019 David Grote, Maxence Thevenet, Weiqun Zhang
- *
+/* Copyright 2019-2022 The WarpX Community
  *
  * This file is part of WarpX.
  *
+ * Authors: David Grote, Maxence Thevenet, Weiqun Zhang, Roelof Groenewald
+ *
  * License: BSD-3-Clause-LBNL
  */
 #include "WarpX_py.H"
 
-std::map< std::string, WARPX_CALLBACK_PY_FUNC_0 > warpx_callback_py_map;
+std::map< std::string, std::function<void()> > warpx_callback_py_map;
+
+void InstallPythonCallback ( std::string name, std::function<void()> callback )
+{
+    warpx_callback_py_map[name] = callback;
+}
 
-bool IsPythonCallBackInstalled ( std::string name )
+bool IsPythonCallbackInstalled ( std::string name )
 {
     return (warpx_callback_py_map.count(name) == 1u);
 }
@@ -17,8 +23,13 @@ bool IsPythonCallBackInstalled ( std::string name )
 // Execute Python callbacks of the type given by the input string
 void ExecutePythonCallback ( std::string name )
 {
-    if ( IsPythonCallBackInstalled(name) ) {
+    if ( IsPythonCallbackInstalled(name) ) {
         WARPX_PROFILE("warpx_py_"+name);
         warpx_callback_py_map[name]();
     }
 }
+
+void ClearPythonCallback ( std::string name )
+{
+    warpx_callback_py_map.erase(name);
+}
diff --git a/Source/Python/pyWarpX.H b/Source/Python/pyWarpX.H
new file mode 100644
index 000000000..51740f75c
--- /dev/null
+++ b/Source/Python/pyWarpX.H
@@ -0,0 +1,20 @@
+/* Copyright 2021-2022 The WarpX Community
+ *
+ * This header is used to centrally define classes that shall not violate the
+ * C++ one-definition-rule (ODR) for various Python translation units.
+ *
+ * Authors: Axel Huebl
+ * License: BSD-3-Clause-LBNL
+ */
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+
+#include <list>
+
+namespace py = pybind11;
+//using namespace warpx;
+
+// PYBIND11_MAKE_OPAQUE(std::list<...>)
diff --git a/Source/Python/pyWarpX.cpp b/Source/Python/pyWarpX.cpp
new file mode 100644
index 000000000..37d910f03
--- /dev/null
+++ b/Source/Python/pyWarpX.cpp
@@ -0,0 +1,136 @@
+/* Copyright 2021-2022 The ImpactX Community
+ *
+ * Authors: Axel Huebl
+ * License: BSD-3-Clause-LBNL
+ */
+#include "pyWarpX.H"
+#include "WarpX_py.H"
+
+#include <WarpX.H>  // todo: move this out to Python/WarpX.cpp
+#include <Utils/WarpXUtil.H>  // todo: move to its own Python/Utils.cpp
+#include <Initialization/WarpXAMReXInit.H>
+
+#define STRINGIFY(x) #x
+#define MACRO_STRINGIFY(x) STRINGIFY(x)
+#define CONCAT_NAME(PRE, SUF) PRE ## SUF
+
+// see: CMakeLists.txt, setup.py and __init__.py
+#if defined(WARPX_DIM_1D_Z)
+#  define PYWARPX_MODULE_NAME CONCAT_NAME(warpx_pybind_, 1d)
+#elif defined(WARPX_DIM_XZ)
+#  define PYWARPX_MODULE_NAME CONCAT_NAME(warpx_pybind_, 2d)
+#elif defined(WARPX_DIM_RZ)
+#  define PYWARPX_MODULE_NAME CONCAT_NAME(warpx_pybind_, rz)
+#elif defined(WARPX_DIM_3D)
+#  define PYWARPX_MODULE_NAME CONCAT_NAME(warpx_pybind_, 3d)
+#endif
+
+//using namespace warpx;
+
+
+// forward declarations of exposed classes
+void init_BoundaryBufferParIter (py::module&);
+void init_MultiParticleContainer (py::module&);
+void init_ParticleBoundaryBuffer (py::module&);
+void init_PinnedMemoryParticleContainer (py::module&);
+void init_WarpXParIter (py::module&);
+void init_WarpXParticleContainer (py::module&);
+void init_WarpX(py::module&);
+
+PYBIND11_MODULE(PYWARPX_MODULE_NAME, m) {
+    // make sure AMReX types are known
+#if defined(WARPX_DIM_3D)
+    auto amr = py::module::import("amrex.space3d");
+#elif defined(WARPX_DIM_1D_Z)
+    auto amr = py::module::import("amrex.space1d");
+#else
+    auto amr = py::module::import("amrex.space2d");
+#endif
+
+    m.doc() = R"pbdoc(
+            warpx_pybind
+            --------------
+            .. currentmodule:: warpx_pybind_(1d|2d|3d|rz)
+
+            .. autosummary::
+               :toctree: _generate
+               WarpX
+    )pbdoc";
+
+    // note: order from parent to child classes
+    init_PinnedMemoryParticleContainer(m);
+    init_WarpXParticleContainer(m);
+    init_WarpXParIter(m);
+    init_BoundaryBufferParIter(m);
+    init_ParticleBoundaryBuffer(m);
+    init_MultiParticleContainer(m);
+    init_WarpX(m);
+
+    // expose our amrex module
+    m.attr("amr") = amr;
+
+    // API runtime version
+    //   note PEP-440 syntax: x.y.zaN but x.y.z.devN
+#ifdef PYIMPACTX_VERSION_INFO
+    m.attr("__version__") = MACRO_STRINGIFY(PYIMPACTX_VERSION_INFO);
+#else
+    m.attr("__version__") = "dev";
+#endif
+
+    // authors
+    m.attr("__author__") =
+        "Jean-Luc Vay, David P. Grote, Maxence Thevenet, Remi Lehe, Andrew Myers, Weiqun Zhang, Axel Huebl, et al.";
+
+    // API runtime build-time feature variants
+    // m.attr("variants") = warpx::getVariants();
+    // TODO allow to query runtime versions of all dependencies
+
+    // license SPDX identifier
+    m.attr("__license__") = "BSD-3-Clause-LBNL";
+
+    // TODO broken numpy if not at least v1.15.0: raise warning
+    // auto numpy = py::module::import("numpy");
+    // auto npversion = numpy.attr("__version__");
+    // std::cout << "numpy version: " << py::str(npversion) << std::endl;
+
+    m.def("amrex_init",
+        [](const py::list args) {
+            amrex::Vector<std::string> cargs;
+            amrex::Vector<char*> argv;
+
+            // Populate the "command line"
+            for (const auto& v: args)
+                cargs.push_back(v.cast<std::string>());
+            for (auto& v: cargs)
+                argv.push_back(&v[0]);
+            int argc = argv.size();
+
+            // note: +1 since there is an extra char-string array element,
+            //       that ANSII C requires to be a simple NULL entry
+            //       https://stackoverflow.com/a/39096006/2719194
+            argv.push_back(NULL);
+            char** tmp = argv.data();
+
+            const bool build_parm_parse = (cargs.size() > 1);
+            // TODO: handle version with MPI
+            return warpx::initialization::amrex_init(argc, tmp, build_parm_parse);
+        }, py::return_value_policy::reference,
+        "Initialize AMReX library");
+    m.def("amrex_finalize", [] () { amrex::Finalize(); },
+        "Close out the amrex related data");
+    m.def("convert_lab_params_to_boost",  &ConvertLabParamsToBoost,
+        "Convert input parameters from the lab frame to the boosted frame");
+    m.def("read_BC_params", &ReadBCParams,
+        "Read the boundary condition parametes and check for consistency");
+    m.def("check_gridding_for_RZ_spectral", &CheckGriddingForRZSpectral,
+        "Ensure that the grid is setup appropriately with using the RZ spectral solver");
+
+    // Expose functions to get the processor number
+    m.def("getNProcs", [](){return amrex::ParallelDescriptor::NProcs();} );
+    m.def("getMyProc", [](){return amrex::ParallelDescriptor::MyProc();} );
+
+    // Expose the python callback function installation and removal functions
+    m.def("add_python_callback", &InstallPythonCallback);
+    m.def("remove_python_callback", &ClearPythonCallback);
+    m.def("execute_python_callback", &ExecutePythonCallback, py::arg("name"));
+}