Merge pull request #568 from dpgrote/pml_python

Adds Python wrapper of PML fields

6 files changed, 576 insertions, 128 deletions

diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py
index 6c3e54619..0444140e4 100755
--- a/Python/pywarpx/_libwarpx.py
+++ b/Python/pywarpx/_libwarpx.py
@@ -130,18 +130,30 @@ libwarpx.warpx_getEfieldCP.restype = _LP_LP_c_real
 libwarpx.warpx_getEfieldCPLoVects.restype = _LP_c_int
 libwarpx.warpx_getEfieldFP.restype = _LP_LP_c_real
 libwarpx.warpx_getEfieldFPLoVects.restype = _LP_c_int
+libwarpx.warpx_getEfieldCP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getEfieldCPLoVects_PML.restype = _LP_c_int
+libwarpx.warpx_getEfieldFP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getEfieldFPLoVects_PML.restype = _LP_c_int
 libwarpx.warpx_getBfield.restype = _LP_LP_c_real
 libwarpx.warpx_getBfieldLoVects.restype = _LP_c_int
 libwarpx.warpx_getBfieldCP.restype = _LP_LP_c_real
 libwarpx.warpx_getBfieldCPLoVects.restype = _LP_c_int
 libwarpx.warpx_getBfieldFP.restype = _LP_LP_c_real
 libwarpx.warpx_getBfieldFPLoVects.restype = _LP_c_int
+libwarpx.warpx_getBfieldCP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getBfieldCPLoVects_PML.restype = _LP_c_int
+libwarpx.warpx_getBfieldFP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getBfieldFPLoVects_PML.restype = _LP_c_int
 libwarpx.warpx_getCurrentDensity.restype = _LP_LP_c_real
 libwarpx.warpx_getCurrentDensityLoVects.restype = _LP_c_int
 libwarpx.warpx_getCurrentDensityCP.restype = _LP_LP_c_real
 libwarpx.warpx_getCurrentDensityCPLoVects.restype = _LP_c_int
 libwarpx.warpx_getCurrentDensityFP.restype = _LP_LP_c_real
 libwarpx.warpx_getCurrentDensityFPLoVects.restype = _LP_c_int
+libwarpx.warpx_getCurrentDensityCP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getCurrentDensityCPLoVects_PML.restype = _LP_c_int
+libwarpx.warpx_getCurrentDensityFP_PML.restype = _LP_LP_c_real
+libwarpx.warpx_getCurrentDensityFPLoVects_PML.restype = _LP_c_int
 
 #libwarpx.warpx_getPMLSigma.restype = _LP_c_real
 #libwarpx.warpx_getPMLSigmaStar.restype = _LP_c_real
@@ -747,6 +759,66 @@ def get_mesh_electric_field_fp(level, direction, include_ghosts=True):
     return _get_mesh_field_list(libwarpx.warpx_getEfieldFP, level, direction, include_ghosts)
 
 
+def get_mesh_electric_field_cp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh electric field
+    data on each grid for this process. This version returns the field on
+    the coarse patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getEfieldCP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
+def get_mesh_electric_field_fp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh electric field
+    data on each grid for this process. This version returns the field on
+    the fine patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getEfieldFP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
 def get_mesh_magnetic_field(level, direction, include_ghosts=True):
     '''
 
@@ -829,6 +901,66 @@ def get_mesh_magnetic_field_fp(level, direction, include_ghosts=True):
     return _get_mesh_field_list(libwarpx.warpx_getBfieldFP, level, direction, include_ghosts)
 
 
+def get_mesh_magnetic_field_cp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh magnetic field
+    data on each grid for this process. This version returns the field on
+    the coarse patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getBfieldCP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
+def get_mesh_magnetic_field_fp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh magnetic field
+    data on each grid for this process. This version returns the field on
+    the fine patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getBfieldFP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
 def get_mesh_current_density(level, direction, include_ghosts=True):
     '''
 
@@ -909,6 +1041,66 @@ def get_mesh_current_density_fp(level, direction, include_ghosts=True):
     return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityFP, level, direction, include_ghosts)
 
 
+def get_mesh_current_density_cp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh current density
+    data on each grid for this process. This version returns the density for
+    the coarse patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityCP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
+def get_mesh_current_density_fp_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of numpy arrays containing the mesh current density
+    data on each grid for this process. This version returns the density on
+    the fine patch for the PML for the given level.
+
+    The data for the numpy arrays are not copied, but share the underlying
+    memory buffer with WarpX. The numpy arrays are fully writeable.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A List of numpy arrays.
+
+    '''
+
+    try:
+        return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityFP_PML, level, direction, include_ghosts)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
 def _get_mesh_array_lovects(level, direction, include_ghosts=True, getarrayfunc=None):
     assert(0 <= level and level <= libwarpx.warpx_finestLevel())
 
@@ -998,6 +1190,56 @@ def get_mesh_electric_field_fp_lovects(level, direction, include_ghosts=True):
     return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldFPLoVects)
 
 
+def get_mesh_electric_field_cp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldCPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
+def get_mesh_electric_field_fp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldFPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
 def get_mesh_magnetic_field_lovects(level, direction, include_ghosts=True):
     '''
 
@@ -1066,6 +1308,56 @@ def get_mesh_magnetic_field_fp_lovects(level, direction, include_ghosts=True):
     return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldFPLoVects)
 
 
+def get_mesh_magnetic_field_cp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldCPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
+def get_mesh_magnetic_field_fp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldFPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+
 def get_mesh_current_density_lovects(level, direction, include_ghosts=True):
     '''
 
@@ -1129,3 +1421,52 @@ def get_mesh_current_density_fp_lovects(level, direction, include_ghosts=True):
 
     '''
     return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityFPLoVects)
+
+
+def get_mesh_current_density_cp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityCPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
+
+def get_mesh_current_density_fp_lovects_pml(level, direction, include_ghosts=True):
+    '''
+
+    This returns a list of the lo vectors of the arrays containing the mesh electric field
+    data on each PML grid for this process.
+
+    Parameters
+    ----------
+
+        level          : the AMR level to get the data for
+        direction      : the component of the data you want
+        include_ghosts : whether to include ghost zones or not
+
+    Returns
+    -------
+
+        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+
+    '''
+    try:
+        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityFPLoVects_PML)
+    except ValueError:
+        raise Exception('PML not initialized')
diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py
index fa247db07..920c109d6 100644
--- a/Python/pywarpx/fields.py
+++ b/Python/pywarpx/fields.py
@@ -87,14 +87,28 @@ class _MultiFABWrapper(object):
     def _getitem3d(self, index):
         """Returns slices of a 3D decomposed array,
         """
-        ix = index[0]
-        iy = index[1]
-        iz = index[2]
 
         lovects = self._getlovects()
         hivects = self._gethivects()
         fields = self._getfields()
 
+        ix = index[0]
+        iy = index[1]
+        iz = index[2]
+
+        if len(fields[0].shape) > self.dim:
+            ncomps = fields[0].shape[-1]
+        else:
+            ncomps = 1
+
+        if len(index) > self.dim:
+            if ncomps > 1:
+                ic = index[-1]
+            else:
+                raise Exception('Too many indices given')
+        else:
+            ic = None
+
         nx = hivects[0,:].max() - self.nghosts
         ny = hivects[1,:].max() - self.nghosts
         nz = hivects[2,:].max() - self.nghosts
@@ -124,9 +138,12 @@ class _MultiFABWrapper(object):
             izstop = iz + 1
 
         # --- Setup the size of the array to be returned and create it.
+        # --- Space is added for multiple components if needed.
         sss = (max(0, ixstop - ixstart),
                max(0, iystop - iystart),
                max(0, izstop - izstart))
+        if ncomps > 1 and ic is None:
+            sss = tuple(list(sss) + [ncomps])
         resultglobal = np.zeros(sss, dtype=_libwarpx._numpy_real_dtype)
 
         datalist = []
@@ -145,6 +162,8 @@ class _MultiFABWrapper(object):
                 sss = (slice(ix1 - lovects[0,i], ix2 - lovects[0,i]),
                        slice(iy1 - lovects[1,i], iy2 - lovects[1,i]),
                        slice(iz1 - lovects[2,i], iz2 - lovects[2,i]))
+                if ic is not None:
+                    sss = tuple(list(sss) + [ic])
 
                 vslice = (slice(ix1 - ixstart, ix2 - ixstart),
                           slice(iy1 - iystart, iy2 - iystart),
@@ -295,6 +314,14 @@ class _MultiFABWrapper(object):
         hivects = self._gethivects()
         fields = self._getfields()
 
+        if len(index) > self.dim:
+            if ncomps > 1:
+                ic = index[-1]
+            else:
+                raise Exception('Too many indices given')
+        else:
+            ic = None
+
         nx = hivects[0,:].max() - self.nghosts
         ny = hivects[1,:].max() - self.nghosts
         nz = hivects[2,:].max() - self.nghosts
@@ -343,6 +370,8 @@ class _MultiFABWrapper(object):
                 sss = (slice(ix1 - lovects[0,i], ix2 - lovects[0,i]),
                        slice(iy1 - lovects[1,i], iy2 - lovects[1,i]),
                        slice(iz1 - lovects[2,i], iz2 - lovects[2,i]))
+                if ic is not None:
+                    sss = tuple(list(sss) + [ic])
 
                 if isinstance(value, np.ndarray):
                     vslice = (slice(ix1 - ixstart, ix2 - ixstart),
@@ -593,3 +622,119 @@ def JzFPWrapper(level=0, include_ghosts=False):
                             get_lovects=_libwarpx.get_mesh_current_density_fp_lovects,
                             get_fabs=_libwarpx.get_mesh_current_density_fp,
                             level=level, include_ghosts=include_ghosts)
+def ExCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=0, overlaps=[0,1,1],
+                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def EyCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=1, overlaps=[1,0,1],
+                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def EzCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=2, overlaps=[1,1,0],
+                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def BxCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=0, overlaps=[1,0,0],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def ByCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=1, overlaps=[0,1,0],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def BzCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=2, overlaps=[0,0,1],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JxCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=0, overlaps=[0,1,1],
+                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JyCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=1, overlaps=[1,0,1],
+                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JzCPPMLWrapper(level=1, include_ghosts=False):
+    assert level>0, Exception('Coarse patch only available on levels > 0')
+    return _MultiFABWrapper(direction=2, overlaps=[1,1,0],
+                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def ExFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=0, overlaps=[0,1,1],
+                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def EyFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=1, overlaps=[1,0,1],
+                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def EzFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=2, overlaps=[1,1,0],
+                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def BxFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=0, overlaps=[1,0,0],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def ByFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=1, overlaps=[0,1,0],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def BzFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=2, overlaps=[0,0,1],
+                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JxFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=0, overlaps=[0,1,1],
+                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JyFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=1, overlaps=[1,0,1],
+                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
+
+def JzFPPMLWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(direction=2, overlaps=[1,1,0],
+                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
+                            level=level, include_ghosts=include_ghosts)
diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp
index 3074b1990..54dfe1955 100644
--- a/Source/Python/WarpXWrappers.cpp
+++ b/Source/Python/WarpXWrappers.cpp
@@ -208,113 +208,80 @@ extern "C"
         return myspc.TotalNumberOfParticles();
     }
 
-    amrex::Real** warpx_getEfield(int lev, int direction,
-                             int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getEfield(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getEfieldLoVects(int lev, int direction,
-                                int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getEfield(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getEfieldCP(int lev, int direction,
-                               int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getEfieldCPLoVects(int lev, int direction,
-                                  int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getEfieldFP(int lev, int direction,
-                               int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getEfieldFPLoVects(int lev, int direction,
-                                  int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getBfield(int lev, int direction,
-                             int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getBfield(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getBfieldLoVects(int lev, int direction,
-                                int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getBfield(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getBfieldCP(int lev, int direction,
-                               int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getBfieldCPLoVects(int lev, int direction,
-                                  int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getBfieldFP(int lev, int direction,
-                               int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getBfieldFPLoVects(int lev, int direction,
-                                  int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getCurrentDensity(int lev, int direction,
-                                     int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getcurrent(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getCurrentDensityLoVects(int lev, int direction,
-                                        int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getcurrent(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getCurrentDensityCP(int lev, int direction,
-                                       int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getCurrentDensityCPLoVects(int lev, int direction,
-                                          int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
-
-    amrex::Real** warpx_getCurrentDensityFP(int lev, int direction,
-                                       int *return_size, int *ncomps, int *ngrow, int **shapes) {
-        auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction);
-        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes);
-    }
-
-    int* warpx_getCurrentDensityFPLoVects(int lev, int direction,
-                                          int *return_size, int *ngrow) {
-        auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction);
-        return getMultiFabLoVects(mf, return_size, ngrow);
-    }
+#define WARPX_GET_FIELD(FIELD, GETTER) \
+    amrex::Real** FIELD(int lev, int direction, \
+                        int *return_size, int *ncomps, int *ngrow, int **shapes) { \
+        auto & mf = GETTER(lev, direction); \
+        return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); \
+    }
+
+#define WARPX_GET_LOVECTS(FIELD, GETTER) \
+    int* FIELD(int lev, int direction, \
+               int *return_size, int *ngrow) { \
+        auto & mf = GETTER(lev, direction); \
+        return getMultiFabLoVects(mf, return_size, ngrow); \
+    }
+
+    WARPX_GET_FIELD(warpx_getEfield, WarpX::GetInstance().getEfield);
+    WARPX_GET_FIELD(warpx_getEfieldCP, WarpX::GetInstance().getEfield_cp);
+    WARPX_GET_FIELD(warpx_getEfieldFP, WarpX::GetInstance().getEfield_fp);
+
+    WARPX_GET_FIELD(warpx_getBfield, WarpX::GetInstance().getBfield);
+    WARPX_GET_FIELD(warpx_getBfieldCP, WarpX::GetInstance().getBfield_cp);
+    WARPX_GET_FIELD(warpx_getBfieldFP, WarpX::GetInstance().getBfield_fp);
+
+    WARPX_GET_FIELD(warpx_getCurrentDensity, WarpX::GetInstance().getcurrent);
+    WARPX_GET_FIELD(warpx_getCurrentDensityCP, WarpX::GetInstance().getcurrent_cp);
+    WARPX_GET_FIELD(warpx_getCurrentDensityFP, WarpX::GetInstance().getcurrent_fp);
+
+    WARPX_GET_LOVECTS(warpx_getEfieldLoVects, WarpX::GetInstance().getEfield);
+    WARPX_GET_LOVECTS(warpx_getEfieldCPLoVects, WarpX::GetInstance().getEfield_cp);
+    WARPX_GET_LOVECTS(warpx_getEfieldFPLoVects, WarpX::GetInstance().getEfield_fp);
+
+    WARPX_GET_LOVECTS(warpx_getBfieldLoVects, WarpX::GetInstance().getBfield);
+    WARPX_GET_LOVECTS(warpx_getBfieldCPLoVects, WarpX::GetInstance().getBfield_cp);
+    WARPX_GET_LOVECTS(warpx_getBfieldFPLoVects, WarpX::GetInstance().getBfield_fp);
+
+    WARPX_GET_LOVECTS(warpx_getCurrentDensityLoVects, WarpX::GetInstance().getcurrent);
+    WARPX_GET_LOVECTS(warpx_getCurrentDensityCPLoVects, WarpX::GetInstance().getcurrent_cp);
+    WARPX_GET_LOVECTS(warpx_getCurrentDensityFPLoVects, WarpX::GetInstance().getcurrent_fp);
+
+#define WARPX_GET_FIELD_PML(FIELD, GETTER) \
+    amrex::Real** FIELD(int lev, int direction, \
+                        int *return_size, int *ncomps, int *ngrow, int **shapes) { \
+        auto * pml = WarpX::GetInstance().GetPML(lev); \
+        if (pml) { \
+            auto & mf = *(pml->GETTER()[direction]); \
+            return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); \
+        } else { \
+            return nullptr; \
+        } \
+    }
+
+#define WARPX_GET_LOVECTS_PML(FIELD, GETTER) \
+    int* FIELD(int lev, int direction, \
+               int *return_size, int *ngrow) { \
+        auto * pml = WarpX::GetInstance().GetPML(lev); \
+        if (pml) { \
+            auto & mf = *(pml->GETTER()[direction]); \
+            return getMultiFabLoVects(mf, return_size, ngrow); \
+        } else { \
+            return nullptr; \
+        } \
+    }
+
+    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);
 
     amrex::ParticleReal** warpx_getParticleStructs(int speciesnumber, int lev,
                                       int* num_tiles, int** particles_per_tile) {
diff --git a/Source/Python/WarpXWrappers.h b/Source/Python/WarpXWrappers.h
index a272b9250..4de885b88 100644
--- a/Source/Python/WarpXWrappers.h
+++ b/Source/Python/WarpXWrappers.h
@@ -67,24 +67,6 @@ extern "C" {
 
     long warpx_getNumParticles(int speciesnumber);
 
-    amrex::Real** warpx_getEfield(int lev, int direction,
-                                  int *return_size, int* ncomps, int* ngrow, int **shapes);
-
-    int* warpx_getEfieldLoVects(int lev, int direction,
-                                int *return_size, int* ngrow);
-
-    amrex::Real** warpx_getBfield(int lev, int direction,
-                                  int *return_size, int* ncomps, int* ngrow, int **shapes);
-
-    int* warpx_getBfieldLoVects(int lev, int direction,
-                                int *return_size, int* ngrow);
-
-    amrex::Real** warpx_getCurrentDensity(int lev, int direction,
-                                          int *return_size, int* ncomps, int* ngrow, int **shapes);
-
-    int* warpx_getCurrentDensityLoVects(int lev, int direction,
-                                        int *return_size, int* ngrow);
-
     amrex::ParticleReal** warpx_getParticleStructs(int speciesnumber, int lev,
                                                    int* num_tiles, int** particles_per_tile);
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 7fa93c156..fe8b2840d 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -207,6 +207,8 @@ public:
 
     void CopyJPML ();
 
+    PML* GetPML (int lev);
+
     void PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full);
     void PushParticlesandDepose (         amrex::Real cur_time);
 
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 53728b54f..377d103d1 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -1251,6 +1251,17 @@ WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp,
     }
 }
 
+PML*
+WarpX::GetPML (int lev)
+{
+    if (do_pml) {
+        // This should check if pml was initialized
+        return pml[lev].get();
+    } else {
+        return nullptr;
+    }
+}
+
 void
 WarpX::BuildBufferMasks ()
 {