diff --git a/pde/fields/base.py b/pde/fields/base.py
index f8c3f3b6..f2a61642 100644
--- a/pde/fields/base.py
+++ b/pde/fields/base.py
@@ -210,6 +210,9 @@ def from_state(
                 The attributes that describe the current instance
             data (:class:`~numpy.ndarray`, optional):
                 Data values at the support points of the grid defining the field
+
+        Returns:
+            :class:`FieldBase`: The field created from the state
         """
         # base class was chosen => select correct class from attributes
         class_name = attributes.pop("class")
@@ -237,7 +240,8 @@ def from_file(cls, filename: str) -> FieldBase:
                 field_copy = pde.FieldBase.from_file("test.hdf5")
 
         Args:
-            filename (str): Path to the file being read
+            filename (str):
+                Path to the file being read
 
         Returns:
             :class:`FieldBase`: The field with the appropriate sub-class
@@ -280,10 +284,10 @@ def _from_hdf_dataset(cls, dataset) -> FieldBase:
 
     @property
     def grid(self) -> GridBase:
-        """GridBase: The grid on which the field is defined"""
+        """:class:`~pde.grids.base,GridBase`: The grid on which the field is defined"""
         return self._grid
 
-    def to_file(self, filename: str, **kwargs):
+    def to_file(self, filename: str, **kwargs) -> None:
         r"""store field in a file
 
         The extension of the filename determines what format is being used. If it ends
@@ -324,7 +328,7 @@ def to_file(self, filename: str, **kwargs):
         else:
             raise ValueError(f"Do not know how to save data to `*{extension}`")
 
-    def _write_hdf_dataset(self, hdf_path, key: str = "data"):
+    def _write_hdf_dataset(self, hdf_path, key: str = "data") -> None:
         """write data to a given hdf5 path `hdf_path`"""
         # write the data
         dataset = hdf_path.create_dataset(key, data=self.data)
@@ -347,7 +351,9 @@ def copy(
     ) -> TField:
         ...
 
-    def assert_field_compatible(self, other: FieldBase, accept_scalar: bool = False):
+    def assert_field_compatible(
+        self, other: FieldBase, accept_scalar: bool = False
+    ) -> None:
         """checks whether `other` is compatible with the current field
 
         Args:
@@ -437,7 +443,7 @@ def _unary_operation(self: TField, op: Callable) -> TField:
                 A function calculating the result
 
         Returns:
-            FieldBase: An field that contains the result of the operation.
+            :class:`FieldBase`: An field that contains the result of the operation.
         """
         return self.__class__(grid=self.grid, data=op(self.data), label=self.label)
 
@@ -452,11 +458,18 @@ def imag(self: TField) -> TField:
         return self._unary_operation(np.imag)
 
     def conjugate(self: TField) -> TField:
-        """returns complex conjugate of the field"""
+        """returns complex conjugate of the field
+
+        Returns:
+            :class:`FieldBase`: the complex conjugated field
+        """
         return self._unary_operation(np.conjugate)
 
     def __neg__(self):
-        """return the negative of the current field"""
+        """return the negative of the current field
+
+        :class:`FieldBase`: The negative of the current field
+        """
         return self._unary_operation(np.negative)
 
     def _binary_operation(
@@ -473,7 +486,7 @@ def _binary_operation(
                 Flag determining whether the second operator must be a scalar
 
         Returns:
-            FieldBase: An field that contains the result of the operation. If
+            :class:`FieldBase`: An field that contains the result of the operation. If
             `scalar_second == True`, the type of FieldBase is the same as `self`
         """
         # determine the dtype of the output
@@ -526,7 +539,7 @@ def _binary_operation_inplace(
                 Flag determining whether the second operator must be a scalar.
 
         Returns:
-            FieldBase: The field `self` with updated data
+            :class:`FieldBase`: The field `self` with updated data
         """
         if isinstance(other, FieldBase):
             # right operator is a field
@@ -638,7 +651,7 @@ def apply(
                 Only used when `func` is a string.
 
         Returns:
-            Field with new data. This is identical to `out` if it was given.
+            :class:`FieldBase`: Field with new data. Identical to `out` if given.
         """
         if isinstance(func, str):
             # function is given as an expression that will be evaluated
@@ -694,7 +707,9 @@ def plot(self, *args, **kwargs):
     def _get_napari_data(self, **kwargs) -> Dict[str, Dict[str, Any]]:
         ...
 
-    def plot_interactive(self, viewer_args: Optional[Dict[str, Any]] = None, **kwargs):
+    def plot_interactive(
+        self, viewer_args: Optional[Dict[str, Any]] = None, **kwargs
+    ) -> None:
         """create an interactive plot of the field using :mod:`napari`
 
         For a detailed description of the launched program, see the
@@ -878,7 +893,7 @@ def random_uniform(
             vmax (float):
                 Largest random value
             label (str, optional):
-                Name of the field
+                Name of the returned field
             dtype (numpy dtype):
                 The data type of the field. If omitted, it defaults to `double` if both
                 `vmin` and `vmax` are real, otherwise it is `complex`.
@@ -937,7 +952,7 @@ def random_normal(
                 scaled by the inverse volume of the grid cell; this is for instance
                 useful for concentration fields, which vary less in larger volumes).
             label (str, optional):
-                Name of the field
+                Name of the returned field
             dtype (numpy dtype):
                 The data type of the field. If omitted, it defaults to `double` if both
                 `mean` and `std` are real, otherwise it is `complex`.
@@ -1016,7 +1031,7 @@ def random_harmonic(
                 resulting in products and sums of the values along axes,
                 respectively.
             label (str, optional):
-                Name of the field
+                Name of the returned field
             dtype (numpy dtype):
                 The data type of the field. If omitted, it defaults to `double`.
             rng (:class:`~numpy.random.Generator`):
@@ -1077,7 +1092,7 @@ def random_colored(
             scale (float):
                 Scaling factor :math:`\Gamma` determining noise strength
             label (str, optional):
-                Name of the field
+                Name of the returned field
             dtype (numpy dtype):
                 The data type of the field. If omitted, it defaults to `double`.
             rng (:class:`~numpy.random.Generator`):
@@ -1105,6 +1120,9 @@ def get_class_by_rank(cls, rank: int) -> Type[DataFieldBase]:
 
         Args:
             rank (int): The rank of the tensor field
+
+        Returns:
+            The DataField class that corresponds to the rank
         """
         for field_cls in cls._subclasses.values():
             if (
@@ -1128,6 +1146,9 @@ def from_state(
                 The attributes that describe the current instance
             data (:class:`~numpy.ndarray`, optional):
                 Data values at the support points of the grid defining the field
+
+        Returns:
+            :class:`DataFieldBase`: The instance created from the stored state
         """
         if "class" in attributes:
             class_name = attributes.pop("class")
@@ -1142,7 +1163,7 @@ def copy(
         label: Optional[str] = None,
         dtype: Optional[DTypeLike] = None,
     ) -> TDataField:
-        """return a copy of the data, but not of the grid
+        """return a new field with the data (but not the grid) copied
 
         Args:
             label (str, optional):
@@ -1150,6 +1171,9 @@ def copy(
             dtype (numpy dtype):
                 The data type of the field. If omitted, it will be determined from
                 `data` automatically or the dtype of the current field is used.
+
+        Returns:
+            :class:`DataFieldBase`: A copy of the current field
         """
         if label is None:
             label = self.label
@@ -1188,7 +1212,7 @@ def unserialize_attributes(cls, attributes: Dict[str, str]) -> Dict[str, Any]:
                 results[key] = json.loads(value)
         return results
 
-    def _write_to_image(self, filename: str, **kwargs):
+    def _write_to_image(self, filename: str, **kwargs) -> None:
         r"""write data to image
 
         Args:
@@ -1519,7 +1543,7 @@ def to_scalar(
 
     @property
     def average(self) -> NumberOrArray:
-        """determine the average of data
+        """float or :class:`~numpy.ndarray`: the average of data
 
         This is calculated by integrating each component of the field over space
         and dividing by the grid volume
@@ -1528,7 +1552,7 @@ def average(self) -> NumberOrArray:
 
     @property
     def fluctuations(self) -> NumberOrArray:
-        """:class:`~numpy.ndarray`: fluctuations over the entire space.
+        """float or :class:`~numpy.ndarray`: quantification of the average fluctuations
 
         The fluctuations are defined as the standard deviation of the data scaled by the
         cell volume. This definition makes the fluctuations independent of the
@@ -1547,7 +1571,7 @@ def fluctuations(self) -> NumberOrArray:
 
     @property
     def magnitude(self) -> float:
-        """float: determine the magnitude of the field.
+        """float: determine the (scalar) magnitude of the field
 
         This is calculated by getting a scalar field using the default arguments of the
         :func:`to_scalar` method, averaging the result over the whole grid, and taking
@@ -1584,12 +1608,14 @@ def apply_operator(
                 Optional field to which the  result is written.
             label (str, optional):
                 Name of the returned field
+            args (dict):
+                Additional arguments for the boundary conditions
             **kwargs:
                 Additional arguments affecting how the operator behaves.
 
         Returns:
-            Field data after applying the operator. This field is identical to `out` if
-            this argument was specified.
+            :class:`DataFieldBase`: Field data after applying the operator. This field
+            is identical to `out` if this argument was specified.
         """
         # get information about the operator
         operator_info = self.grid._get_operator_info(operator)
@@ -1637,6 +1663,8 @@ def _apply_operator(
                 Optional field to which the  result is written.
             label (str, optional):
                 Name of the returned field
+            args (dict):
+                Additional arguments for the boundary conditions
             **kwargs:
                 Additional arguments affecting how the operator behaves.
 
@@ -1829,17 +1857,18 @@ def smooth(
     ) -> TDataField:
         """applies Gaussian smoothing with the given standard deviation
 
-        This function respects periodic boundary conditions of the underlying
-        grid, using reflection when no periodicity is specified.
-
-        sigma (float):
-            Gives the standard deviation of the smoothing in real length units
-            (default: 1)
-        out (FieldBase, optional):
-            Optional field into which the smoothed data is stored. Setting this
-            to the input field enables in-place smoothing.
-        label (str, optional):
-            Name of the returned field
+        This function respects periodic boundary conditions of the underlying grid,
+        using reflection when no periodicity is specified.
+
+        Args:
+            sigma (float):
+                Gives the standard deviation of the smoothing in real length units
+                (default: 1)
+            out (FieldBase, optional):
+                Optional field into which the smoothed data is stored. Setting this
+                to the input field enables in-place smoothing.
+            label (str, optional):
+                Name of the returned field
 
         Returns:
             Field with smoothed data. This is stored at `out` if given.
diff --git a/pde/fields/collection.py b/pde/fields/collection.py
index eb4635bc..eb5eda36 100644
--- a/pde/fields/collection.py
+++ b/pde/fields/collection.py
@@ -23,6 +23,7 @@
 )
 
 import numpy as np
+from numpy.typing import DTypeLike
 
 from ..grids.base import GridBase
 from ..tools.docstrings import fill_in_docstring
@@ -49,7 +50,7 @@ def __init__(
         copy_fields: bool = False,
         label: Optional[str] = None,
         labels: Union[List[Optional[str]], _FieldLabels, None] = None,
-        dtype=None,
+        dtype: DTypeLike = None,
     ):
         """
         Args:
@@ -263,7 +264,7 @@ def from_dict(
         *,
         copy_fields: bool = False,
         label: Optional[str] = None,
-        dtype=None,
+        dtype: DTypeLike = None,
     ) -> FieldCollection:
         """create a field collection from a dictionary of fields
 
@@ -326,7 +327,7 @@ def from_data(
         with_ghost_cells: bool = True,
         label: Optional[str] = None,
         labels: Union[List[Optional[str]], _FieldLabels, None] = None,
-        dtype=None,
+        dtype: DTypeLike = None,
     ):
         """create a field collection from classes and data
 
@@ -423,7 +424,7 @@ def from_scalar_expressions(
         consts: Optional[Dict[str, NumberOrArray]] = None,
         label: Optional[str] = None,
         labels: Optional[Sequence[str]] = None,
-        dtype=None,
+        dtype: DTypeLike = None,
     ) -> FieldCollection:
         """create a field collection on a grid from given expressions
 
@@ -570,7 +571,7 @@ def copy(
         self: FieldCollection,
         *,
         label: Optional[str] = None,
-        dtype=None,
+        dtype: DTypeLike = None,
     ) -> FieldCollection:
         """return a copy of the data, but not of the grid
 
diff --git a/pde/fields/vectorial.py b/pde/fields/vectorial.py
index c34481b1..abfaa883 100644
--- a/pde/fields/vectorial.py
+++ b/pde/fields/vectorial.py
@@ -173,7 +173,7 @@ def dot(
             other (VectorField or Tensor2Field):
                 the second field
             out (ScalarField or VectorField, optional):
-                Optional field to which the  result is written.
+                Optional field to which the result is written.
             conjugate (bool):
                 Whether to use the complex conjugate for the second operand
             label (str, optional):
diff --git a/pde/grids/_mesh.py b/pde/grids/_mesh.py
index 11e453f1..3bb31d1d 100644
--- a/pde/grids/_mesh.py
+++ b/pde/grids/_mesh.py
@@ -50,8 +50,10 @@ def _subdivide(num: int, chunks: int) -> np.ndarray:
     r"""subdivide `num` intervals in `chunk` chunks
 
     Args:
-        num (int): Number of intervals
-        chunks (int): Number of chunks
+        num (int):
+            Number of intervals
+        chunks (int):
+            Number of chunks
 
     Returns:
         list: The number of intervals per chunk
@@ -65,8 +67,10 @@ def _subdivide_along_axis(grid: GridBase, axis: int, chunks: int) -> List[GridBa
     """subdivide the grid along a given axis
 
     Args:
-        axis (int): The axis along which the subdivision will happen
-        chunks (int): The number of chunks along this axis
+        axis (int):
+            The axis along which the subdivision will happen
+        chunks (int):
+            The number of chunks along this axis
 
     Returns:
         list: A list of subgrids
@@ -152,7 +156,7 @@ def from_grid(
                 nodes.
 
         Returns:
-            :class:`GridMesh`
+            :class:`GridMesh`: The grid mesh created from the grid
         """
         # parse `decomposition`
         try:
@@ -678,7 +682,8 @@ def broadcast(self, data: TData) -> TData:
         """distribute a value from the main node to all nodes
 
         Args:
-            data: The data that will be broadcasted from the main node
+            data:
+                The data that will be broadcasted from the main node
 
         Returns:
             The same data, but on all nodes
@@ -691,7 +696,8 @@ def gather(self, data: TData) -> Optional[List[TData]]:
         """gather a value from all nodes
 
         Args:
-            data: The data that will be sent to the main node
+            data:
+                The data that will be sent to the main node
 
         Returns:
             None on all nodes, except the main node, which receives an ordered list with
@@ -705,7 +711,8 @@ def allgather(self, data: TData) -> List[TData]:
         """gather a value from reach node and sends them to all nodes
 
         Args:
-            data: The data that will be sent to the main node
+            data:
+                The data that will be sent to the main node
 
         Returns:
             list: data from all nodes.
@@ -715,17 +722,18 @@ def allgather(self, data: TData) -> List[TData]:
         return COMM_WORLD.allgather(data)
 
     @plot_on_axes()
-    def plot(self, ax, **kwargs):
+    def plot(self, ax, **kwargs) -> None:
         r"""visualize the grid mesh
 
         Args:
             {PLOT_ARGS}
-            \**kwargs: Extra arguments are passed on the to the matplotlib
-                plotting routines, e.g., to set the color of the lines
+            \**kwargs:
+                Extra arguments are passed on the to the matplotlib plotting routines,
+                e.g., to set the color of the lines
         """
         if self.num_axes not in {1, 2}:
             raise NotImplementedError(
-                f"Plotting is not implemented for grids of dimension {self.dim}"
+                f"Cannot plot data of dimension {self.basegrid.dim}"
             )
 
         kwargs.setdefault("color", "k")
diff --git a/pde/grids/base.py b/pde/grids/base.py
index bdee8b88..dc03e934 100644
--- a/pde/grids/base.py
+++ b/pde/grids/base.py
@@ -2,7 +2,6 @@
 Bases classes
 
 .. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
- 
 """
 
 from __future__ import annotations
@@ -23,6 +22,7 @@
     Generator,
     Iterator,
     List,
+    Literal,
     NamedTuple,
     Optional,
     Sequence,
@@ -189,6 +189,9 @@ def from_state(cls, state: Union[str, Dict[str, Any]]) -> GridBase:
             state (`str` or `dict`):
                 The state from which the grid is reconstructed. If `state` is a
                 string, it is decoded as JSON, which should yield a `dict`.
+
+        Returns:
+            :class:`GridBase`: Grid re-created from the state
         """
         # decode the json data
         if isinstance(state, str):
@@ -229,8 +232,10 @@ def get_axis_index(self, key: Union[int, str], allow_symmetric: bool = True) ->
         """return the index belonging to an axis
 
         Args:
-            key (int or str): The index or name of an axis
-            allow_symmetric (bool): Whether axes with assumed symmetry are included
+            key (int or str):
+                The index or name of an axis
+            allow_symmetric (bool):
+                Whether axes with assumed symmetry are included
 
         Returns:
             int: The index of the axis
@@ -264,7 +269,7 @@ def _get_boundary_index(
                 e.g., :code:`(1, True)`.
 
         Returns:
-            Tuple: axis index perpendicular to the boundary and a boolean specifying
+            tuple: axis index perpendicular to the boundary and a boolean specifying
                 whether the boundary is at the upper side of the axis or not.
         """
         if isinstance(index, str):
@@ -300,7 +305,12 @@ def _idx_valid(self) -> Tuple[slice, ...]:
         return tuple(slice(1, s + 1) for s in self.shape)
 
     def _make_get_valid(self) -> Callable[[np.ndarray], np.ndarray]:
-        """callable: function to extract the valid part of a full data array"""
+        """create a function to extract the valid part of a full data array
+
+        Returns:
+            callable: Mapping a numpy array containing the full data of the grid to a
+                numpy array of only the valid data
+        """
         num_axes = self.num_axes
 
         @jit
@@ -318,7 +328,12 @@ def get_valid(arr: np.ndarray) -> np.ndarray:
         return get_valid  # type: ignore
 
     def _make_set_valid(self) -> Callable[[np.ndarray, np.ndarray], None]:
-        """callable: function to extract the valid part of a full data array"""
+        """create a function to set the valid part of a full data array
+
+        Returns:
+            callable: Takes two numpy arrays, setting the valid data in the first one,
+                using the second array
+        """
         num_axes = self.num_axes
 
         @jit
@@ -456,8 +471,10 @@ def difference_vector_real(self, p1: np.ndarray, p2: np.ndarray) -> np.ndarray:
             for curvilinear coordinates.
 
         Args:
-            p1 (:class:`~numpy.ndarray`): First point(s)
-            p2 (:class:`~numpy.ndarray`): Second point(s)
+            p1 (:class:`~numpy.ndarray`):
+                First point(s)
+            p2 (:class:`~numpy.ndarray`):
+                Second point(s)
 
         Returns:
             :class:`~numpy.ndarray`: The difference vectors between the points with
@@ -482,8 +499,10 @@ def distance_real(self, p1: np.ndarray, p2: np.ndarray) -> float:
             results for coordinates given in curvilinear coordinate systems.
 
         Args:
-            p1 (:class:`~numpy.ndarray`): First position
-            p2 (:class:`~numpy.ndarray`): Second position
+            p1 (:class:`~numpy.ndarray`):
+                First position
+            p2 (:class:`~numpy.ndarray`):
+                Second position
 
         Returns:
             float: Distance between the two positions
@@ -640,7 +659,10 @@ def _grid_to_cell(
             return cells  # type: ignore
 
     def transform(
-        self, coordinates: np.ndarray, source: str, target: str
+        self,
+        coordinates: np.ndarray,
+        source: Literal["cartesian", "cell", "grid"],
+        target: Literal["cartesian", "cell", "grid"],
     ) -> np.ndarray:
         """converts coordinates from one coordinate system to another
 
@@ -659,9 +681,12 @@ def transform(
             originating at the origin.
 
         Args:
-            coordinates (:class:`~numpy.ndarray`): The coordinates to convert
-            source (str): The source coordinate system
-            target (str): The target coordinate system
+            coordinates (:class:`~numpy.ndarray`):
+                The coordinates to convert
+            source (str):
+                The source coordinate system
+            target (str):
+                The target coordinate system
 
         Returns:
             :class:`~numpy.ndarray`: The transformed coordinates
@@ -725,7 +750,10 @@ def polar_coordinates_real(
         ...
 
     def contains_point(
-        self, points: np.ndarray, *, coords: str = "cartesian", wrap: bool = True
+        self,
+        points: np.ndarray,
+        *,
+        coords: Literal["cartesian", "cell", "grid"] = "cartesian",
     ) -> np.ndarray:
         """check whether the point is contained in the grid
 
@@ -1105,12 +1133,20 @@ def apply_op_compiled(
             raise NotImplementedError(f"Undefined backend '{backend}'")
 
     def slice(self, indices: Sequence[int]) -> GridBase:
-        """return a subgrid of only the specified axes"""
+        """return a subgrid of only the specified axes
+
+        Args:
+            indices (list):
+                Indices indicating the axes that are retained in the subgrid
+
+        Returns:
+            :class:`GridBase`: The subgrid
+        """
         raise NotImplementedError(
             f"Slicing is not implemented for class {self.__class__.__name__}"
         )
 
-    def plot(self):
+    def plot(self) -> None:
         """visualize the grid"""
         raise NotImplementedError(
             f"Plotting is not implemented for class {self.__class__.__name__}"
@@ -1286,9 +1322,9 @@ def _make_interpolation_axis_data(
                 point coordinates.
 
         Returns:
-            A function that is called with a coordinate value for the axis. The function
-            returns the indices of the neighboring support points as well as the
-            associated weights
+            callable: A function that is called with a coordinate value for the axis.
+            The function returns the indices of the neighboring support points as well
+            as the associated weights.
         """
         # obtain information on how this axis is discretized
         size = self.shape[axis]
@@ -1377,10 +1413,10 @@ def _make_interpolator_compiled(
                 point coordinates.
 
         Returns:
-            A function which returns interpolated values when called with arbitrary
-            positions within the space of the grid. The signature of this function is
-            (data, point), where `data` is the numpy array containing the field data and
-            position denotes the position in grid coordinates.
+            callable: A function which returns interpolated values when called with
+            arbitrary positions within the space of the grid. The signature of this
+            function is (data, point), where `data` is the numpy array containing the
+            field data and position denotes the position in grid coordinates.
         """
         args = {"with_ghost_cells": with_ghost_cells, "cell_coords": cell_coords}
 
@@ -1515,9 +1551,10 @@ def make_inserter_compiled(
                 boundaries are not checked and the coordinates are used as is.
 
         Returns:
-            A function with signature (data, position, amount), where `data` is the numpy
-            array containing the field data, position is denotes the position in grid
-            coordinates, and `amount` is the  that is to be added to the field.
+            callable: A function with signature (data, position, amount), where `data`
+            is the numpy array containing the field data, position is denotes the
+            position in grid coordinates, and `amount` is the  that is to be added to
+            the field.
         """
         cell_volume = self.make_cell_volume_compiled()
 
@@ -1663,7 +1700,7 @@ def insert(
         return insert  # type: ignore
 
     def make_integrator(self) -> Callable[[np.ndarray], NumberOrArray]:
-        """Return function that can be used to integrates discretized data over the grid
+        """return function that can be used to integrates discretized data over the grid
 
         If this function is used in a multiprocessing run (using MPI), the integrals are
         performed on all subgrids and then accumulated. Each process then receives the
@@ -1749,7 +1786,6 @@ def integrate_global(arr: np.ndarray) -> NumberOrArray:
 def registered_operators() -> Dict[str, List[str]]:
     """returns all operators that are currently defined
 
-
     Returns:
         dict: a dictionary with the names of the operators defined for each grid class
     """
diff --git a/pde/grids/cartesian.py b/pde/grids/cartesian.py
index d01c72a7..4ccbb1f2 100644
--- a/pde/grids/cartesian.py
+++ b/pde/grids/cartesian.py
@@ -2,7 +2,6 @@
 Cartesian grids of arbitrary dimension.
 
 .. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
- 
 """
 
 from __future__ import annotations
@@ -169,6 +168,9 @@ def from_state(cls, state: Dict[str, Any]) -> CartesianGrid:  # type: ignore
         Args:
             state (dict):
                 The state from which the grid is reconstructed.
+
+        Returns:
+            :class:`CartesianGrid`: the grid re-created from the state data
         """
         state_copy = state.copy()
         obj = cls(
@@ -202,7 +204,7 @@ def from_bounds(
                 or a single boolean value specifying the same periodicity for all axes.
 
         Returns:
-            :class:`CartesianGrid` representing the region chosen by bounds
+            :class:`CartesianGrid`: representing the region chosen by bounds
         """
         return CartesianGrid(bounds, shape, periodic)
 
@@ -222,9 +224,12 @@ def iter_mirror_points(
         """generates all mirror points corresponding to `point`
 
         Args:
-            point (:class:`~numpy.ndarray`): the point within the grid
-            with_self (bool): whether to include the point itself
-            only_periodic (bool): whether to only mirror along periodic axes
+            point (:class:`~numpy.ndarray`):
+                The point within the grid
+            with_self (bool):
+                Whether to include the point itself
+            only_periodic (bool):
+                Whether to only mirror along periodic axes
 
         Returns:
             A generator yielding the coordinates that correspond to mirrors
@@ -305,8 +310,8 @@ def get_line_data(self, data: np.ndarray, extract: str = "auto") -> Dict[str, An
                 a letter denoting the axis.
 
         Returns:
-            A dictionary with information about the line cut, which is
-            convenient for plotting.
+            dict: A dictionary with information about the line cut, which is convenient
+            for plotting.
         """
         if data.shape[-self.dim :] != self.shape:
             raise ValueError(
@@ -365,10 +370,11 @@ def get_image_data(self, data: np.ndarray) -> Dict[str, Any]:
         """return a 2d-image of the data
 
         Args:
-            data (:class:`~numpy.ndarray`): The values at the grid points
+            data (:class:`~numpy.ndarray`):
+                The values at the grid points
 
         Returns:
-            A dictionary with information about the image, which is  convenient
+            dict: A dictionary with information about the image, which is  convenient
             for plotting.
         """
         if data.shape[-self.dim :] != self.shape:
@@ -405,8 +411,10 @@ def point_to_cartesian(
         """convert coordinates of a point to Cartesian coordinates
 
         Args:
-            points (:class:`~numpy.ndarray`): Points given in grid coordinates
-            full (bool): Compatibility option not used in this method
+            points (:class:`~numpy.ndarray`):
+                Points given in grid coordinates
+            full (bool):
+                Compatibility option not used in this method
 
         Returns:
             :class:`~numpy.ndarray`: The Cartesian coordinates of the point
@@ -418,7 +426,8 @@ def point_from_cartesian(self, coords: np.ndarray) -> np.ndarray:
         """convert points given in Cartesian coordinates to this grid
 
         Args:
-            coords (:class:`~numpy.ndarray`): Points in Cartesian coordinates.
+            coords (:class:`~numpy.ndarray`):
+                Points in Cartesian coordinates.
 
         Returns:
             :class:`~numpy.ndarray`: Points given in the coordinates of the grid
@@ -443,6 +452,10 @@ def polar_coordinates_real(
                 point and the origin, so that angles can either be 1 or -1. For 2d
                 systems and 3d systems, polar coordinates and spherical coordinates are
                 used, respectively.
+
+        Returns:
+            :class:`~numpy.ndarray` or tuple of :class:`~numpy.ndarray`:
+                Coordinates values in polar coordinates
         """
         origin = np.array(origin, dtype=np.double, ndmin=1)
         if len(origin) != self.dim:
@@ -622,7 +635,11 @@ def from_state(cls, state: Dict[str, Any]) -> UnitGrid:  # type: ignore
         return obj
 
     def to_cartesian(self) -> CartesianGrid:
-        """convert unit grid to :class:`CartesianGrid`"""
+        """convert unit grid to :class:`CartesianGrid`
+
+        Returns:
+            :class:`CartesianGrid`: The equivalent cartesian grid
+        """
         return CartesianGrid(
             self.cuboid.bounds, shape=self.shape, periodic=self.periodic
         )
diff --git a/pde/grids/cylindrical.py b/pde/grids/cylindrical.py
index 3cd3bf60..43f80411 100644
--- a/pde/grids/cylindrical.py
+++ b/pde/grids/cylindrical.py
@@ -2,7 +2,6 @@
 Cylindrical grids with azimuthal symmetry
 
 .. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
- 
 """
 
 from __future__ import annotations
@@ -21,7 +20,7 @@
 
 
 class CylindricalSymGrid(GridBase):
-    r""" 3-dimensional cylindrical grid assuming polar symmetry 
+    r"""3-dimensional cylindrical grid assuming polar symmetry 
 
     The polar symmetry implies that states only depend on the radial and axial
     coordinates :math:`r` and :math:`z`, respectively. These are discretized uniformly as
@@ -138,7 +137,7 @@ def state(self) -> Dict[str, Any]:
         }
 
     @classmethod
-    def from_state(cls, state: Dict[str, Any]) -> "CylindricalSymGrid":  # type: ignore
+    def from_state(cls, state: Dict[str, Any]) -> CylindricalSymGrid:  # type: ignore
         """create a field from a stored `state`.
 
         Args:
@@ -177,7 +176,7 @@ def from_bounds(
                 entry is ignored.
 
         Returns:
-            CylindricalGrid representing the region chosen by bounds
+            :class:`CylindricalSymGrid`: grid representing the region chosen by bounds
         """
         radii, bounds_z = bounds
         if radii[0] != 0:
@@ -186,7 +185,7 @@ def from_bounds(
 
     @property
     def has_hole(self) -> bool:
-        """returns whether the inner radius is larger than zero"""
+        """bool: whether the inner radius is larger than zero"""
         return self.axes_bounds[0][0] > 0
 
     @property
@@ -283,8 +282,8 @@ def get_line_data(self, data: np.ndarray, extract: str = "auto") -> Dict[str, An
                 * `project_r` or `project_radial`: average values for each
                   radial position (axial average)
         Returns:
-            A dictionary with information about the line cut, which is
-            convenient for plotting.
+            dict: A dictionary with information about the line cut, which is convenient
+            for plotting.
         """
         if extract == "auto":
             extract = "cut_axial"
@@ -322,10 +321,11 @@ def get_image_data(self, data: np.ndarray) -> Dict[str, Any]:
         """return a 2d-image of the data
 
         Args:
-            data (:class:`~numpy.ndarray`): The values at the grid points
+            data (:class:`~numpy.ndarray`):
+                The values at the grid points
 
         Returns:
-            A dictionary with information about the image, which is  convenient
+            dict: A dictionary with information about the image, which is  convenient
             for plotting.
         """
         bounds_r, bounds_z = self.axes_bounds
@@ -352,9 +352,12 @@ def iter_mirror_points(
         """generates all mirror points corresponding to `point`
 
         Args:
-            point (:class:`~numpy.ndarray`): the point within the grid
-            with_self (bool): whether to include the point itself
-            only_periodic (bool): whether to only mirror along periodic axes
+            point (:class:`~numpy.ndarray`):
+                The point within the grid
+            with_self (bool):
+                Whether to include the point itself
+            only_periodic (bool):
+                Whether to only mirror along periodic axes
 
         Returns:
             A generator yielding the coordinates that correspond to mirrors
@@ -383,8 +386,10 @@ def point_to_cartesian(
         """convert coordinates of a point to Cartesian coordinates
 
         Args:
-            points (:class:`~numpy.ndarray`): The grid coordinates of the points
-            full (bool): Flag indicating whether angular coordinates are specified
+            points (:class:`~numpy.ndarray`):
+                The grid coordinates of the points
+            full (bool):
+                Flag indicating whether angular coordinates are specified
 
         Returns:
             :class:`~numpy.ndarray`: The Cartesian coordinates of the point
@@ -430,13 +435,15 @@ def polar_coordinates_real(
         """return spherical coordinates associated with the grid
 
         Args:
-            origin (:class:`~numpy.ndarray`): Coordinates of the origin at which the polar
-                coordinate system is anchored. Note that this must be of the
-                form `[0, 0, z_val]`, where only `z_val` can be chosen freely.
-            ret_angle (bool): Determines whether the azimuthal angle is returned
-                alongside the distance. If `False` only the distance to the
-                origin is  returned for each support point of the grid.
-                If `True`, the distance and angles are returned.
+            origin (:class:`~numpy.ndarray`):
+                Coordinates of the origin at which the polar coordinate system is
+                anchored. Note that this must be of the form `[0, 0, z_val]`, where only
+                `z_val` can be chosen freely.
+            ret_angle (bool):
+                Determines whether the azimuthal angle is returned alongside the
+                distance. If `False` only the distance to the origin is  returned for
+                each support point of the grid. If `True`, the distance and angles are
+                returned.
         """
         origin = np.array(origin, dtype=np.double, ndmin=1)
         if len(origin) != self.dim:
@@ -459,10 +466,10 @@ def get_cartesian_grid(self, mode: str = "valid") -> CartesianGrid:
 
         Args:
             mode (str):
-                Determines how the grid is determined. Setting it to 'valid'
-                only returns points that are fully resolved in the cylindrical
-                grid, e.g., the cylinder is circumscribed. Conversely, 'full'
-                returns all data, so the cylinder is inscribed.
+                Determines how the grid is determined. Setting it to 'valid' only
+                returns points that are fully resolved in the cylindrical grid, e.g.,
+                the cylinder is circumscribed. Conversely, 'full' returns all data, so
+                the cylinder is inscribed.
 
         Returns:
             :class:`pde.grids.cartesian.CartesianGrid`: The requested grid
diff --git a/pde/grids/spherical.py b/pde/grids/spherical.py
index b13f2308..9dba45c0 100644
--- a/pde/grids/spherical.py
+++ b/pde/grids/spherical.py
@@ -7,7 +7,6 @@
 vanishes.
 
 .. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
- 
 """
 
 from __future__ import annotations
@@ -38,8 +37,10 @@ def volume_from_radius(radius: TNumArr, dim: int) -> TNumArr:
     """Return the volume of a sphere with a given radius
 
     Args:
-        radius (float or :class:`~numpy.ndarray`): Radius of the sphere
-        dim (int): Dimension of the space
+        radius (float or :class:`~numpy.ndarray`):
+            Radius of the sphere
+        dim (int):
+            Dimension of the space
 
     Returns:
         float or :class:`~numpy.ndarray`: Volume of the sphere
@@ -83,11 +84,11 @@ def __init__(
         r"""
         Args:
             radius (float or tuple of floats):
-                radius :math:`R_\mathrm{outer}` in case a simple float is given.
-                If a tuple is supplied it is interpreted as the inner and outer
-                radius, :math:`(R_\mathrm{inner}, R_\mathrm{outer})`.
-            shape (tuple or int): A single number setting the number :math:`N`
-                of support points along the radial coordinate
+                Radius :math:`R_\mathrm{outer}` in case a simple float is given. If a
+                tuple is supplied it is interpreted as the inner and outer radius,
+                :math:`(R_\mathrm{inner}, R_\mathrm{outer})`.
+            shape (tuple or int):
+                The number :math:`N` of support points along the radial coordinate.
         """
         super().__init__()
         shape_list = _check_shape(shape)
@@ -140,12 +141,13 @@ def from_bounds(  # type: ignore
     ) -> SphericalSymGridBase:
         """
         Args:
-            bounds (tuple): Give the coordinate range for the radial axis.
-            shape (tuple): The number of support points for the radial axis
-            periodic (bool or list): Not used
+            bounds (tuple):
+                Give the coordinate range for the radial axis.
+            shape (tuple):
+                The number of support points for the radial axis
 
         Returns:
-            SphericalGridBase representing the region chosen by bounds
+            :class:`SphericalGridBase`: represents the region chosen by bounds
         """
         if len(bounds) != 1:
             raise ValueError(
@@ -258,9 +260,8 @@ def get_line_data(self, data: np.ndarray, extract: str = "auto") -> Dict[str, An
             data (:class:`~numpy.ndarray`):
                 The values at the grid points
             extract (str):
-                Determines which cut is done through the grid. This parameter is
-                mainly supplied for a consistent interface and has no effect for
-                polar grids.
+                Determines which cut is done through the grid. This parameter is mainly
+                supplied for a consistent interface and has no effect for polar grids.
 
         Returns:
             A dictionary with information about the line cut, which is
@@ -289,17 +290,17 @@ def get_image_data(
             data (:class:`~numpy.ndarray`):
                 The values at the grid points
             performance_goal (str):
-                Determines the method chosen for interpolation. Possible options
-                are `speed` and `quality`.
+                Determines the method chosen for interpolation. Possible options are
+                `speed` and `quality`.
             fill_value (float):
-                The value assigned to invalid positions (those inside the hole
-                or outside the region).
+                The value assigned to invalid positions (those inside the hole or
+                outside the region).
             masked (bool):
-                Whether a :class:`numpy.ma.MaskedArray` is returned for the data
-                instead of the normal :class:`~numpy.ndarray`.
+                Whether a :class:`numpy.ma.MaskedArray` is returned for the data instead
+                of the normal :class:`~numpy.ndarray`.
 
         Returns:
-            A dictionary with information about the image, which is  convenient
+            :dict: A dictionary with information about the image, which is  convenient
             for plotting.
         """
         from scipy import interpolate
@@ -356,9 +357,12 @@ def iter_mirror_points(
         """generates all mirror points corresponding to `point`
 
         Args:
-            point (:class:`~numpy.ndarray`): the point within the grid
-            with_self (bool): whether to include the point itself
-            only_periodic (bool): whether to only mirror along periodic axes
+            point (:class:`~numpy.ndarray`):
+                The point within the grid
+            with_self (bool):
+                Whether to include the point itself
+            only_periodic (bool):
+                Whether to only mirror along periodic axes
 
         Returns:
             A generator yielding the coordinates that correspond to mirrors
@@ -449,8 +453,9 @@ def plot(self, ax, **kwargs):
 
         Args:
             {PLOT_ARGS}
-            \**kwargs: Extra arguments are passed on the to the matplotlib
-                plotting routines, e.g., to set the color of the lines
+            \**kwargs:
+                Extra arguments are passed on the to the matplotlib plotting routines,
+                e.g., to set the color of the lines
         """
         from matplotlib import collections, patches
 
@@ -508,8 +513,10 @@ def point_to_cartesian(
         coordinates will be zero.
 
         Args:
-            points (:class:`~numpy.ndarray`): The grid coordinates of the points
-            full (bool): Flag indicating whether angular coordinates are specified
+            points (:class:`~numpy.ndarray`):
+                The grid coordinates of the points
+            full (bool):
+                Flag indicating whether angular coordinates are specified
 
         Returns:
             :class:`~numpy.ndarray`: The Cartesian coordinates of the point
@@ -573,8 +580,10 @@ def point_to_cartesian(
         coordinates will be zero.
 
         Args:
-            points (:class:`~numpy.ndarray`): The grid coordinates of the points
-            full (bool): Flag indicating whether angular coordinates are specified
+            points (:class:`~numpy.ndarray`):
+                The grid coordinates of the points
+            full (bool):
+                Flag indicating whether angular coordinates are specified
 
         Returns:
             :class:`~numpy.ndarray`: The Cartesian coordinates of the point
diff --git a/pde/tools/cuboid.py b/pde/tools/cuboid.py
index f3ed8e4d..7428202e 100644
--- a/pde/tools/cuboid.py
+++ b/pde/tools/cuboid.py
@@ -11,6 +11,7 @@
 from typing import List, Tuple
 
 import numpy as np
+from numpy.typing import DTypeLike
 
 from .typing import FloatNumerical
 
@@ -225,7 +226,7 @@ def contains_point(self, points: np.ndarray) -> np.ndarray:
         return np.all(c1 <= points, axis=-1) & np.all(points <= c2, axis=-1)  # type: ignore
 
 
-def asanyarray_flags(data: np.ndarray, dtype=None, writeable: bool = True):
+def asanyarray_flags(data: np.ndarray, dtype: DTypeLike = None, writeable: bool = True):
     """turns data into an array and sets the respective flags.
 
     A copy is only made if necessary
diff --git a/pde/tools/misc.py b/pde/tools/misc.py
index cf8afc43..3d2c5ca0 100644
--- a/pde/tools/misc.py
+++ b/pde/tools/misc.py
@@ -31,6 +31,7 @@
 from typing import Any, Callable, Dict, Optional, Sequence, TypeVar, Union
 
 import numpy as np
+from numpy.typing import DTypeLike
 
 from .typing import ArrayLike, Number
 
@@ -352,7 +353,9 @@ def get_common_dtype(*args):
     return np.double
 
 
-def number_array(data: ArrayLike, dtype=None, copy: bool = True) -> np.ndarray:
+def number_array(
+    data: ArrayLike, dtype: DTypeLike = None, copy: bool = True
+) -> np.ndarray:
     """convert an array with arbitrary dtype either to np.double or np.cdouble
 
     Args:
diff --git a/scripts/tests_parallel.sh b/scripts/tests_parallel.sh
index 6e9db162..0914f21d 100755
--- a/scripts/tests_parallel.sh
+++ b/scripts/tests_parallel.sh
@@ -1,12 +1,12 @@
 #!/bin/bash
 
-if [ ! -z $1 ] 
-then 
-	# test pattern was specified 
-	echo 'Run unittests with pattern '$1':'
+if [ ! -z $1 ]
+then
+    # test pattern was specified 
+    echo 'Run unittests with pattern '$1':'
     ./run_tests.py --unit --runslow --num_cores auto --pattern "$1"
 else
-	# test pattern was not specified
-	echo 'Run all unittests:'
+    # test pattern was not specified
+    echo 'Run all unittests:'
     ./run_tests.py --unit --num_cores auto
 fi