Refactored boundary conditions (#619)

This PR changes the way in which boundary conditions should be normally set. In particular, it deprecates the normal nested list approach, where conditions were set for each side of all axes in a more flexible, dictionary-based approach. In the new approach, conditions are directly specified for entire axes (`x`) or for one side (`x+`).

Moreover, we refactored the `Boundaries` class, introducing a new base class and added the `BoundariesSetter` class, so ghost points can be set directly using a callback function.

We also cleaned up some code and moved around some classes, so this might break old code!

docs/source/manual/advanced_usage.rst

@@ -13,20 +13,16 @@ field to be `4` at the lower side and its derivative (in the outward direction)
 
 .. code-block:: python
 
-    bc_lower = {"value": 4}
-    bc_upper = {"derivative": 2}
-    bc = [bc_lower, bc_upper]
+    bc = {"x-": {"value": 4}, "x+": {"derivative": 2}}
 
     grid = pde.UnitGrid([16])
     field = pde.ScalarField(grid)
     field.laplace(bc)
 
 Here, the Laplace operator applied to the field in the last line will respect
-the boundary conditions.
-Note that it suffices to give one condition if both sides of the axis require the same
-condition.
-For instance, to enforce a value of `3` on both side, one could simply use
-:code:`bc = {'value': 3}`.
+the boundary conditions. Note that both sides of the axis can be specified together if
+their conditions are the same. For instance, to enforce a value of `3` on both side, one
+could simply use :code:`bc = {"x": {"value": 3}}`.
 Vectorial boundary conditions, e.g., to calculate the vector gradient or tensor
 divergence, can have vectorial values for the boundary condition.
 Generally, only the normal components at a boundary need to be specified if an operator
@@ -39,33 +35,42 @@ which may depend on the coordinates of all axes:
 .. code-block:: python
 
     # two different conditions for lower and upper end of x-axis
-    bc_x = [{"derivative": 0.1}, {"value": "sin(y / 2)"}] 
+    bc_left = {"derivative": 0.1}
+    bc_right = {"value": "sin(y / 2)"}
     # the same condition on the lower and upper end of the y-axis
     bc_y = {"value": "sqrt(1 + cos(x))"}
 
     grid = UnitGrid([32, 32])
     field = pde.ScalarField(grid)
-    field.laplace(bc=[bc_x, bc_y])
+    field.laplace(bc={"x-": bc_left, "x+": bc_right, "y": bc_y})
 
 .. warning::
     To interpret arbitrary expressions, the package uses :func:`exec`. It
     should therefore not be used in a context where malicious input could occur.
 
-Inhomogeneous values can also be specified by directly supplying an array, whose shape
+Heterogeneous values can also be specified by directly supplying an array, whose shape
 needs to be compatible with the boundary, i.e., it needs to have the same shape as the
 grid but with the dimension of the axis along which the boundary is specified removed.
 
-There exist also special boundary conditions that impose a more complex value of the
-field (:code:`bc='value_expression'`) or its derivative
-(:code:`bc='derivative_expression'`). Beyond the spatial coordinates that are already
+There also exist special boundary conditions that impose a more complex value of the
+field (:code:`bc="value_expression"`) or its derivative
+(:code:`bc="derivative_expression"`). Beyond the spatial coordinates that are already
 supported for the constant conditions above, the expressions of these boundary
 conditions can depend on the time variable :code:`t`. Moreover, these boundary
-conditions also except python functions (with signature `adjacent_value, dx, *coords, t`),
+conditions also except python functions with signature `(adjacent_value, dx, *coords, t)`,
 thus greatly enlarging the flexibility with which boundary conditions can be expressed.
 Note that PDEs need to supply the current time `t` when setting the boundary conditions,
 e.g., when applying the differential operators. The pre-defined PDEs and the general
 class :class:`~pde.pdes.pde.PDE` already support time-dependent boundary conditions.
 
+To specify the same boundary conditions for many sides, the wildcard specifier can be
+used: For example, :code:`bc = {"*": {"value": 1}, "y+": {"derivative": 0}}` specifies
+Dirichlet conditions for all axes, except the upper y-axis, where a Neumann condition is
+imposed instead. If all axes have the same condition, the outer dictionary can be
+skipped, so that :code:`bc = {"*": {"value": 1}}` imposes the same conditions as
+:code:`bc = {"value": 1}`. Moreover, many boundaries have convenient names, so that for
+instance `x-` can be replaced by `left`, and `y+` can be replaced by `top`.
+
 One important aspect about boundary conditions is that they need to respect the
 periodicity of the underlying grid. For instance, in a 2d grid with one periodic axis,
 the following boundary condition can be used:
@@ -74,8 +79,7 @@ the following boundary condition can be used:
 
     grid = pde.UnitGrid([16, 16], periodic=[True, False])
     field = pde.ScalarField(grid)
-    bc = ["periodic", {"derivative": 0}]
-    field.laplace(bc)
+    field.laplace({"x": "periodic", "y": {"derivative": 0})
 
 For convenience, this typical situation can be described with the special boundary
 condition `auto_periodic_neumann`, e.g., calling the Laplace operator using 
@@ -146,6 +150,13 @@ Here, :math:`\partial_n` denotes a derivative in outward normal direction, :math
 denotes an arbitrary function given by an expression (see next section), :math:`x`
 denotes coordinates along the boundary, :math:`t` denotes time.
 
+Finally, we support the advanced technique of setting the virtual points at the boundary
+manually. This can be achieved by passing a python function that takes as 
+its first argument a :class:`~numpy.ndarray`, which contains the full field data
+including the virtual points, and a second, optional argument, which is a dictionary
+containing additional parameters, like the current time point `t` in case of a
+simulation; see :class:`~pde.grids.boundaries.axes.BoundariesSetter` for more details.
+
 
 .. _documentation-expressions:
 

docs/source/manual/contributing.rst

@@ -32,10 +32,10 @@ New operators can be associated with grids by registering them using
 :meth:`~pde.grids.base.GridBase.register_operator`.
 For instance, to create a new operator for the cylindrical grid one needs to 
 define a factory function that creates the operator. This factory function takes
-an instance of :class:`~pde.grids.boundaries.axes.Boundaries` as an argument and
+an instance of :class:`~pde.grids.boundaries.axes.BoundariesList` as an argument and
 returns a function that takes as an argument the actual data array for the grid.
 Note that the grid itself is an attribute of
-:class:`~pde.grids.boundaries.axes.Boundaries`.
+:class:`~pde.grids.boundaries.axes.BoundariesList`.
 This operator would be registered with the grid by calling
 :code:`CylindricalSymGrid.register_operator("operator", make_operator)`, where the
 first argument is the name of the operator and the second argument is the

examples/boundary_conditions.py

@@ -11,11 +11,9 @@
 state = ScalarField.random_uniform(grid, 0.2, 0.3)  # generate initial condition
 
 # set boundary conditions `bc` for all axes
-bc_x_left = {"derivative": 0.1}
-bc_x_right = {"value": "sin(y / 2)"}
-bc_x = [bc_x_left, bc_x_right]
-bc_y = "periodic"
-eq = DiffusionPDE(bc=[bc_x, bc_y])
+eq = DiffusionPDE(
+    bc={"x-": {"derivative": 0.1}, "x+": {"value": "sin(y / 2)"}, "y": "periodic"}
+)
 
 result = eq.solve(state, t_range=10, dt=0.005)
 result.plot()

examples/heterogeneous_bcs.py

@@ -32,10 +32,7 @@ def bc_value(adjacent_value, dx, x, y, t):
     return np.sign(x)
 
 
-bc_x = "derivative"
-bc_y = ["derivative", {"value_expression": bc_value}]
-
 # define and solve a simple diffusion equation
-eq = DiffusionPDE(bc=[bc_x, bc_y])
+eq = DiffusionPDE(bc={"*": {"derivative": 0}, "y+": {"value_expression": bc_value}})
 res = eq.solve(field, t_range=10, dt=0.01, backend="numpy")
 res.plot()

examples/jupyter/Discretized Fields.ipynb

@@ -64,7 +64,7 @@
    "source": [
     "# apply operators to the field\n",
     "smoothed = scalar.smooth(1)\n",
-    "laplace = smoothed.laplace(bc=\"natural\")\n",
+    "laplace = smoothed.laplace(bc=\"auto_periodic_neumann\")\n",
     "laplace.plot(colorbar=True);"
    ]
   },

examples/jupyter/Solve PDEs.ipynb

@@ -86,7 +86,7 @@
     "    Taken from https://en.wikipedia.org/wiki/Chafee-Infante_equation\n",
     "    \"\"\"\n",
     "\n",
-    "    def __init__(self, lmbd=1, bc=\"natural\"):\n",
+    "    def __init__(self, lmbd=1, bc=\"auto_periodic_neumann\"):\n",
     "        super().__init__()\n",
     "        self.bc = bc\n",
     "        self.lmbd = lmbd\n",

examples/laplace_eq_2d.py

@@ -10,8 +10,8 @@
 
 from pde import CartesianGrid, solve_laplace_equation
 
-grid = CartesianGrid([[0, 2 * np.pi]] * 2, 64)
-bcs = [{"value": "sin(y)"}, {"value": "sin(x)"}]
+grid = CartesianGrid([[0, 2 * np.pi], [0, 2 * np.pi]], 64)
+bcs = {"x": {"value": "sin(y)"}, "y": {"value": "sin(x)"}}
 
 res = solve_laplace_equation(grid, bcs)
 res.plot()

examples/poisson_eq_1d.py

@@ -9,6 +9,6 @@
 
 grid = CartesianGrid([[0, 1]], 32, periodic=False)
 field = ScalarField(grid, 1)
-result = solve_poisson_equation(field, bc=[{"value": 0}, {"derivative": 1}])
+result = solve_poisson_equation(field, bc={"x-": {"value": 0}, "x+": {"derivative": 1}})
 
 result.plot()

examples/tutorial/Tutorial 1 - Grids and fields.ipynb

@@ -321,7 +321,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Periodic and non-periodic axes can also be mixed as in the example below. In this case, the $x$-axis is periodic and thus requires periodic boundary conditions. Conversely, the $y$-axis is non-periodic and any other boundary condition can be specified. The most generic one is a Neumann condition of vanishing derivative. For convenience, we also define `natural` boundary conditions, which indicate periodic conditions for periodic axes and Neumann conditions otherwise."
+    "Periodic and non-periodic axes can also be mixed as in the example below. In this case, the $x$-axis is periodic and thus requires periodic boundary conditions. Conversely, the $y$-axis is non-periodic and any other boundary condition can be specified. The most generic one is a Neumann condition of vanishing derivative. For convenience, we also define `auto_periodic_neumann` boundary conditions, which indicate periodic conditions for periodic axes and Neumann conditions otherwise."
    ]
   },
   {
@@ -336,7 +336,7 @@
     "field_mixed = pde.ScalarField.from_expression(grid_mixed, \"sin(x) * cos(y)\")\n",
     "\n",
     "laplace_mixed = field_mixed.laplace([\"periodic\", {\"derivative\": 0}])\n",
-    "# laplace_mixed = field_mixed.laplace('natural')\n",
+    "# laplace_mixed = field_mixed.laplace('auto_periodic_neumann')\n",
     "laplace_mixed.plot(title=\"Laplacian of mixed field\", colorbar=True);"
    ]
   },
@@ -373,7 +373,7 @@
    "outputs": [],
    "source": [
     "field_per = pde.ScalarField.from_expression(grid_per, \"sin(x) * cos(y)\")\n",
-    "field_grad = field_per.gradient(\"natural\")\n",
+    "field_grad = field_per.gradient(\"auto_periodic_neumann\")\n",
     "field_grad"
    ]
   },
@@ -463,7 +463,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "field_hess = field_grad.gradient(\"natural\", label=\"Hessian of field\")\n",
+    "field_hess = field_grad.gradient(\"auto_periodic_neumann\", label=\"Hessian of field\")\n",
     "field_hess.attributes"
    ]
   },

examples/tutorial/Tutorial 2 - Solving pre-defined partial differential equations.ipynb

@@ -375,9 +375,9 @@
     "\n",
     "    def evolution_rate(self, state, t=0):\n",
     "        \"\"\"implement the python version of the evolution equation\"\"\"\n",
-    "        state_lap = state.laplace(bc=\"natural\")\n",
-    "        state_lap2 = state_lap.laplace(bc=\"natural\")\n",
-    "        state_grad = state.gradient(bc=\"natural\")\n",
+    "        state_lap = state.laplace(bc=\"auto_periodic_neumann\")\n",
+    "        state_lap2 = state_lap.laplace(bc=\"auto_periodic_neumann\")\n",
+    "        state_grad = state.gradient(bc=\"auto_periodic_neumann\")\n",
     "        return -state_grad.to_scalar(\"squared_sum\") / 2 - state_lap - state_lap2\n",
     "\n",
     "\n",

pde/fields/datafield_base.py

@@ -738,23 +738,23 @@ def get_boundary_values(
         return (self._data_full[i_wall] + self._data_full[i_ghost]) / 2  # type: ignore
 
     @fill_in_docstring
-    def set_ghost_cells(
-        self, bc: BoundariesData, *, set_corners: bool = False, args=None
-    ) -> None:
-        """Set the boundary values on virtual points for all boundaries.
+    def set_ghost_cells(self, bc: BoundariesData, *, args=None, **kwargs) -> None:
+        r"""Set the boundary values on virtual points for all boundaries.
 
         Args:
             bc (str or list or tuple or dict):
                 The boundary conditions applied to the field.
                 {ARG_BOUNDARIES}
-            set_corners (bool):
-                Determines whether the corner cells are set using interpolation
             args:
                 Additional arguments that might be supported by special boundary
                 conditions.
+            **kwargs:
+                Some boundary conditions might provide additional control via keyword
+                arguments, e.g., `set_corners` to control whether corner cells are set
+                using interpolation.
         """
         bcs = self.grid.get_boundary_conditions(bc, rank=self.rank)
-        bcs.set_ghost_cells(self._data_full, args=args, set_corners=set_corners)
+        bcs.set_ghost_cells(self._data_full, args=args, **kwargs)
 
     @property
     @abstractmethod

pde/grids/_mesh.py

@@ -19,8 +19,8 @@
 from ..tools.cache import cached_method
 from ..tools.plotting import plot_on_axes
 from .base import GridBase
-from .boundaries.axes import Boundaries
-from .boundaries.axis import BoundaryPair
+from .boundaries.axes import BoundariesBase, BoundariesList
+from .boundaries.axis import BoundaryAxisBase, BoundaryPair
 from .boundaries.local import _MPIBC
 
 
@@ -514,18 +514,26 @@ def extract_subfield(
         else:
             raise TypeError(f"Field type {field.__class__.__name__} unsupported")
 
-    def extract_boundary_conditions(self, bcs_base: Boundaries) -> Boundaries:
+    def extract_boundary_conditions(self, bcs_base: BoundariesBase) -> BoundariesList:
         """Extract boundary conditions for current subgrid from global conditions.
 
         Args:
-            bcs_base (:class:`~pde.grids.boundaries.axes.Boundaries`):
+            bcs_base (:class:`~pde.grids.boundaries.axes.BoundariesBase`):
                 The boundary conditions that will be split
 
         Returns:
-            :class:`~pde.grids.boundaries.axes.Boundaries`: Boundary conditions of the
-            sub grid
+            :class:`~pde.grids.boundaries.axes.BoundariesList`: Boundary conditions of
+            the subgrid
         """
-        bcs = []
+        if not isinstance(bcs_base, BoundariesList):
+            raise TypeError(
+                "Simulations parallelized with MPI only work with boundary conditions "
+                "based on the `BoundariesList` class, i.e., conditions need to be "
+                "specified for each axes separately and cannot be set using a "
+                "function globally."
+            )
+
+        bcs: list[BoundaryAxisBase] = []
         for axis in range(self.num_axes):
             bcs_axis = []
             for upper in [False, True]:
@@ -539,7 +547,7 @@ def extract_boundary_conditions(self, bcs_base: Boundaries) -> Boundaries:
                 bcs_axis.append(bc)
             bcs.append(BoundaryPair(*bcs_axis))
 
-        return Boundaries(bcs)
+        return BoundariesList(bcs)
 
     def split_field_data_mpi(
         self, field_data: np.ndarray, *, with_ghost_cells: bool = False

pde/grids/base.py

@@ -37,7 +37,7 @@
 
 if TYPE_CHECKING:
     from ._mesh import GridMesh
-    from .boundaries.axes import Boundaries, BoundariesData
+    from .boundaries.axes import BoundariesBase, BoundariesData
 
 
 PI_4 = 4 * np.pi
@@ -338,14 +338,14 @@ def _make_set_valid(self) -> Callable[[np.ndarray, np.ndarray], None]: ...
 
     @overload
     def _make_set_valid(
-        self, bcs: Boundaries
+        self, bcs: BoundariesBase
     ) -> Callable[[np.ndarray, np.ndarray, dict], None]: ...
 
-    def _make_set_valid(self, bcs: Boundaries | None = None) -> Callable:
+    def _make_set_valid(self, bcs: BoundariesBase | None = None) -> Callable:
         """Create a function to set the valid part of a full data array.
 
         Args:
-            bcs (:class:`~pde.grids.boundaries.axes.Boundaries`, optional):
+            bcs (:class:`~pde.grids.boundaries.axes.BoundariesBase`, optional):
                 If supplied, the returned function also enforces boundary conditions by
                 setting the ghost cells to the correct values
 
@@ -1049,7 +1049,7 @@ def iter_mirror_points(
     @fill_in_docstring
     def get_boundary_conditions(
         self, bc: BoundariesData = "auto_periodic_neumann", rank: int = 0
-    ) -> Boundaries:
+    ) -> BoundariesBase:
         """Constructs boundary conditions from a flexible data format.
 
         Args:
@@ -1060,26 +1060,26 @@ def get_boundary_conditions(
                 The tensorial rank of the value associated with the boundary conditions.
 
         Returns:
-            :class:`~pde.grids.boundaries.axes.Boundaries`: The boundary conditions for
-            all axes.
+            :class:`~pde.grids.boundaries.axes.BoundariesBase`: The boundary conditions
+            for all axes.
 
         Raises:
             ValueError:
                 If the data given in `bc` cannot be read
             PeriodicityError:
                 If the boundaries are not compatible with the periodic axes of the grid.
         """
-        from .boundaries import Boundaries
+        from .boundaries import BoundariesBase
 
         if self._mesh is None:
             # get boundary conditions for a simple grid that is not part of a mesh
-            bcs = Boundaries.from_data(self, bc, rank=rank)
+            bcs = BoundariesBase.from_data(bc, grid=self, rank=rank)
 
         else:
             # this grid is part of a mesh and we thus need to set special conditions to
             # support parallelism via MPI. We here assume that bc is given for the full
             # system and not
-            bcs_base = Boundaries.from_data(self._mesh.basegrid, bc, rank=rank)
+            bcs_base = BoundariesBase.from_data(bc, grid=self._mesh.basegrid, rank=rank)
             bcs = self._mesh.extract_boundary_conditions(bcs_base)
 
         return bcs