From 98cc140cadcbfb5b48c6052ac31de7bd2d24544f Mon Sep 17 00:00:00 2001 From: David Zwicker Date: Sat, 23 Nov 2024 14:55:11 +0100 Subject: [PATCH] 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 | 45 +- docs/source/manual/contributing.rst | 4 +- examples/boundary_conditions.py | 8 +- examples/heterogeneous_bcs.py | 5 +- examples/jupyter/Discretized Fields.ipynb | 2 +- examples/jupyter/Solve PDEs.ipynb | 2 +- examples/laplace_eq_2d.py | 4 +- examples/poisson_eq_1d.py | 2 +- .../Tutorial 1 - Grids and fields.ipynb | 8 +- ...fined partial differential equations.ipynb | 6 +- pde/fields/datafield_base.py | 14 +- pde/grids/_mesh.py | 24 +- pde/grids/base.py | 20 +- pde/grids/boundaries/__init__.py | 94 ++-- pde/grids/boundaries/axes.py | 520 +++++++++++++----- pde/grids/boundaries/axis.py | 60 +- pde/grids/boundaries/local.py | 14 +- pde/grids/cartesian.py | 21 +- pde/grids/cylindrical.py | 1 - pde/grids/operators/cartesian.py | 23 +- pde/grids/operators/cylindrical_sym.py | 10 +- pde/grids/operators/polar_sym.py | 10 +- pde/grids/operators/spherical_sym.py | 10 +- pde/grids/spherical.py | 6 +- pde/pdes/allen_cahn.py | 8 +- pde/pdes/cahn_hilliard.py | 14 +- pde/pdes/diffusion.py | 7 +- pde/pdes/kpz_interface.py | 7 +- pde/pdes/kuramoto_sivashinsky.py | 16 +- pde/pdes/laplace.py | 3 +- pde/pdes/pde.py | 7 +- pde/pdes/swift_hohenberg.py | 9 +- pde/pdes/wave.py | 7 +- pde/tools/docstrings.py | 24 +- scripts/format_code.sh | 6 +- scripts/performance_laplace.py | 6 +- tests/fields/test_generic_fields.py | 2 +- .../grids/boundaries/test_axes_boundaries.py | 108 +++- .../boundaries/test_axes_boundaries_legacy.py | 81 +++ .../operators/test_cartesian_operators.py | 2 +- tests/grids/test_cartesian_grids.py | 4 +- tests/test_integration.py | 22 + tests/tools/test_expressions.py | 4 +- 43 files changed, 873 insertions(+), 377 deletions(-) create mode 100644 tests/grids/boundaries/test_axes_boundaries_legacy.py diff --git a/docs/source/manual/advanced_usage.rst b/docs/source/manual/advanced_usage.rst index 8c7a9f5d..7f12d2cb 100644 --- a/docs/source/manual/advanced_usage.rst +++ b/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: diff --git a/docs/source/manual/contributing.rst b/docs/source/manual/contributing.rst index 3e97f96a..60173c90 100644 --- a/docs/source/manual/contributing.rst +++ b/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 diff --git a/examples/boundary_conditions.py b/examples/boundary_conditions.py index c5722ee0..1b54926b 100644 --- a/examples/boundary_conditions.py +++ b/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() diff --git a/examples/heterogeneous_bcs.py b/examples/heterogeneous_bcs.py index f31cc79e..822b3b05 100644 --- a/examples/heterogeneous_bcs.py +++ b/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() diff --git a/examples/jupyter/Discretized Fields.ipynb b/examples/jupyter/Discretized Fields.ipynb index 74c08e6e..caf243e6 100644 --- a/examples/jupyter/Discretized Fields.ipynb +++ b/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);" ] }, diff --git a/examples/jupyter/Solve PDEs.ipynb b/examples/jupyter/Solve PDEs.ipynb index 0afe6097..e4193a67 100644 --- a/examples/jupyter/Solve PDEs.ipynb +++ b/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", diff --git a/examples/laplace_eq_2d.py b/examples/laplace_eq_2d.py index 82990ef5..c75ff0eb 100644 --- a/examples/laplace_eq_2d.py +++ b/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() diff --git a/examples/poisson_eq_1d.py b/examples/poisson_eq_1d.py index 1a6674a1..db9e9649 100644 --- a/examples/poisson_eq_1d.py +++ b/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() diff --git a/examples/tutorial/Tutorial 1 - Grids and fields.ipynb b/examples/tutorial/Tutorial 1 - Grids and fields.ipynb index 635ba43a..ecf45f1f 100644 --- a/examples/tutorial/Tutorial 1 - Grids and fields.ipynb +++ b/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" ] }, diff --git a/examples/tutorial/Tutorial 2 - Solving pre-defined partial differential equations.ipynb b/examples/tutorial/Tutorial 2 - Solving pre-defined partial differential equations.ipynb index 61ad7cb6..0615689e 100644 --- a/examples/tutorial/Tutorial 2 - Solving pre-defined partial differential equations.ipynb +++ b/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", diff --git a/pde/fields/datafield_base.py b/pde/fields/datafield_base.py index 9762189f..7e42246b 100644 --- a/pde/fields/datafield_base.py +++ b/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 diff --git a/pde/grids/_mesh.py b/pde/grids/_mesh.py index e722a2a2..928d2b0a 100644 --- a/pde/grids/_mesh.py +++ b/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 diff --git a/pde/grids/base.py b/pde/grids/base.py index e2e23e36..38837e48 100644 --- a/pde/grids/base.py +++ b/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,8 +1060,8 @@ 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: @@ -1069,17 +1069,17 @@ def get_boundary_conditions( 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 diff --git a/pde/grids/boundaries/__init__.py b/pde/grids/boundaries/__init__.py index fe9a5c0e..0b893ad7 100644 --- a/pde/grids/boundaries/__init__.py +++ b/pde/grids/boundaries/__init__.py @@ -8,30 +8,38 @@ The mathematical details of boundary conditions for partial differential equations are treated in more detail in the :download:`documentation document `. -Since the :mod:`pde` package only supports orthogonal grids, boundary conditions need to -be applied at the end of each axis. -Consequently, methods expecting boundary conditions typically receive a list of +Since the :mod:`pde` package only supports orthogonal grids, boundary conditions +generally need to be applied at both ends of each axis. +Consequently, methods expecting boundary conditions typically receive a dictionary of conditions for each axes: .. code-block:: python field = ScalarField(UnitGrid([16, 16], periodic=[True, False])) - field.laplace(bc=[bc_x, bc_y]) + field.laplace(bc={"x": bc_x, "y-": bc_y_lower, "y+": bc_y_upper}) + +If both sides of an axis have the same boundary condition, they can be specified +together, e.g., if `bc_y_lower == bc_y_upper`, one could have used +:code:`{"x": bc_x, "y": bc_y_lower}` instead of the example above. Moreover, it is +possible to specify boundary conditions for all sides that have not a specific condition +specified using :code:`{"*": default_bc}`. Similarly, boundary conditions for entire +axes can be overwritten by conditions specified on one side. Finally, the boundary sides +often have aliases defined by the grid, so one can use `left` instead of `x-` and so on. If an axis is periodic (like the first one in the example above), the only valid boundary conditions are 'periodic' and its cousin 'anti-periodic', which imposes -opposite signs on both sides. For non-periodic axes (e.g., the second axis), -different boundary conditions can be specified for the lower and upper end of the axis, -which is done using a tuple of two conditions. Typical choices for individual -conditions are Dirichlet conditions that enforce a value NUM (specified by -`{'value': NUM}`) and Neumann conditions that enforce the value DERIV for the -derivative in the normal direction (specified by `{'derivative': DERIV}`). The specific -choices for the example above could be +opposite signs on both sides. For non-periodic axes (e.g., the second axis), different +boundary conditions can be specified for the lower and upper end of the axis, as in the +example above. Typical choices for individual conditions are Dirichlet conditions that +enforce a value NUM (specified by `{'value': NUM}`) and Neumann conditions that enforce +the value DERIV for the derivative in the normal direction (specified by +`{'derivative': DERIV}`). The specific choices for the example above could be .. code-block:: python bc_x = "periodic" - bc_y = ({"value": 2}, {"derivative": -1}) + bc_y_lower = {"value": 2} + bc_y_upper = {"derivative": -1} which enforces a value of `2` at the lower side of the y-axis and a derivative (in outward normal direction) of `-1` on the upper side. Instead of plain @@ -42,12 +50,12 @@ .. code-block:: python - bc_y = ({"value": "y**2"}, {"derivative": "-sin(x)"}) - + bc_y_lower = {"value": "y**2"} + bc_y_upper = {"derivative": "-sin(x)"} Warning: - To interpret arbitrary expressions, the package uses :func:`exec`. It - should therefore not be used in a context where malicious input could occur. + 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 needs to be compatible with the boundary, i.e., it needs to have the same shape as the @@ -58,18 +66,19 @@ .. code-block:: python - bc_y = ({"type": "mixed", "value": 2, "const": 7}, {"curvature": 2}) + bc_y_lower = {"type": "mixed", "value": 2, "const": 7} + bc_y_upper = {"curvature": 2} which enforces the condition :math:`\partial_n c + 2 c = 7` and -:math:`\partial^2_n c = 2` onto the field :math:`c` on the lower and upper side -of the axis, respectively. +:math:`\partial^2_n c = 2` onto the field :math:`c` on the lower and upper side of the +axis, respectively. -Beside the full specification of the boundary conditions, various short-hand -notations are supported. If both sides of an axis have the same boundary -condition, only one needs to be specified (instead of the tuple). For instance, -:code:`bc_y = {'value': 2}` imposes a value of `2` on both sides of the y-axis. -Similarly, if all axes have the same boundary conditions, only one axis needs to -be specified (instead of the list). For instance, the following example +Beside the full specification of boundary conditions, various short-hand notations +are supported. If both sides of an axis have the same boundary condition, only one needs +to be specified. For instance, :code:`{"x": {"value": 2}}` is equivalent to +:code:`{"x-": {"value": 2}, "x+": {"value": 2}}` and imposes a value of `2` on both +sides of the x-axis. In the special case where all sides have the same boundary +conditions, only this condition can be specified instead of the full dictionary, e.g. .. code-block:: python @@ -77,8 +86,9 @@ field.laplace(bc={"value": 2}) imposes a value of `2` on all sides of the grid. Finally, the special values -'auto_periodic_neumann' and 'auto_periodic_dirichlet' impose periodic boundary -conditions for periodic axis and a vanishing derivative or value otherwise. For example, +:code:`"auto_periodic_neumann"` and :code:`"auto_periodic_dirichlet"` impose periodic +boundary conditions for periodic axis and a vanishing derivative or value otherwise. +For example, .. code-block:: python @@ -92,6 +102,24 @@ Derivatives are given relative to the outward normal vector, such that positive derivatives correspond to a function that increases across the boundary. +If more complex boundary conditions are required, a custom function that directly sets +the boundary conditions can also be supplied. This special approach comes with few +checks, so only use it in exceptional circumstances. The following example shows a +setter function, which sets specific boundary conditions in the x-direction and periodic +conditions in the y-direction of a grid with two axes. + +.. code-block:: python + + def setter(data, args=None): + data[0, :] = data[1, :] # Vanishing derivative at left side + data[-1, :] = 2 - data[-2, :] # Fixed value `1` at right side + data[:, 0] = data[:, -2] # Periodic BC at top + data[:, -1] = data[:, 1] # Periodic BC at bottom + + + field = ScalarField(UnitGrid([16, 16], periodic=[False, True])) + field.laplace(bc=setter) + Boundaries overview ^^^^^^^^^^^^^^^^^^^ @@ -136,15 +164,19 @@ * :class:`~pde.grids.boundaries.axis.BoundaryPeriodic`: Indicates that an axis has periodic boundary conditions -**Boundaries for all axes of a grid:** +**BoundariesList for all axes of a grid:** -* :class:`~pde.grids.boundaries.axes.Boundaries`: +* :class:`~pde.grids.boundaries.axes.BoundariesList`: Collection of boundaries to describe conditions for all axes +* :class:`~pde.grids.boundaries.axes.BoundariesSetter`: + Describes custom function setting virtual points to impose boundary conditions **Inheritance structure of the classes:** -.. inheritance-diagram:: pde.grids.boundaries.axes pde.grids.boundaries.axis +.. inheritance-diagram:: + pde.grids.boundaries.axes + pde.grids.boundaries.axis pde.grids.boundaries.local :parts: 2 @@ -154,7 +186,7 @@ """ from ..base import DomainError, PeriodicityError -from .axes import Boundaries +from .axes import BoundariesBase, BoundariesList, set_default_bc from .local import ( registered_boundary_condition_classes, registered_boundary_condition_names, diff --git a/pde/grids/boundaries/axes.py b/pde/grids/boundaries/axes.py index a323d0e3..8ea7b567 100644 --- a/pde/grids/boundaries/axes.py +++ b/pde/grids/boundaries/axes.py @@ -1,36 +1,117 @@ -r""" -.. codeauthor:: David Zwicker +r"""This module handles the boundaries of all axes of a grid. + +.. autosummary:: + :nosignatures: + + ~BoundariesBase + ~BoundariesList + ~BoundariesSetter + ~set_default_bc -This module handles the boundaries of all axes of a grid. It only defines -:class:`Boundaries`, which acts as a list of -:class:`~pde.grids.boundaries.axis.BoundaryAxisBase`. +.. codeauthor:: David Zwicker """ from __future__ import annotations import itertools import logging +import warnings from collections.abc import Iterator, Sequence -from typing import Union +from typing import Any, Callable, Union import numpy as np from numba.extending import register_jitable +from ...tools.numba import jit from ...tools.typing import GhostCellSetter from ..base import GridBase, PeriodicityError -from .axis import BoundaryPair, BoundaryPairData, get_boundary_axis -from .local import BCBase, BCDataError +from .axis import BoundaryAxisBase, BoundaryPair, BoundaryPairData, get_boundary_axis +from .local import BCBase, BCDataError, BoundaryData + +BoundariesData = Union[ + BoundaryPairData, Sequence[BoundaryPairData], Callable, "BoundariesBase" +] + +BC_LOCAL_KEYS = ["type", "value"] + list(BCBase._conditions.keys()) + + +def _is_local_bc_data(data: dict[str, Any]) -> bool: + """Tries to identify whether data specifies a local boundary condition.""" + return any(key in data for key in BC_LOCAL_KEYS) + + +class BoundariesBase: + """Base class keeping information about how to set conditions on all boundaries.""" + + def set_ghost_cells(self, data_full: np.ndarray, *, args=None) -> None: + """Set the ghost cells for all boundaries. + + Args: + data_full (:class:`~numpy.ndarray`): + The full field data including ghost points + set_corners (bool): + Determines whether the corner cells are set using interpolation + args: + Additional arguments that might be supported by special boundary + conditions. + """ + raise NotImplementedError -BoundariesData = Union[BoundaryPairData, Sequence[BoundaryPairData]] + def make_ghost_cell_setter(self) -> GhostCellSetter: + """Return function that sets the ghost cells on a full array. + Returns: + Callable with signature :code:`(data_full: np.ndarray, args=None)`, which + sets the ghost cells of the full data, potentially using additional + information in `args` (e.g., the time `t` during solving a PDE) + """ + raise NotImplementedError -class Boundaries(list): - """Class that bundles all boundary conditions for all axes.""" + @classmethod + def from_data(cls, data, **kwargs) -> BoundariesBase: + r"""Creates all boundaries from given data. - grid: GridBase - """:class:`~pde.grids.base.GridBase`: grid for which boundaries are defined.""" + Args: + data (str or dict or callable): + Data that describes the boundaries. If this is a callable, we create + :class:`~pde.grids.boundaries.axes.BoundariesSetter`. In all other, + cases :class:`~pde.grids.boundaries.axes.BoundariesList` is created and + `data` can either be string denoting a specific boundary condition + applied to all sides or a dictionary with detailed information. + **kwargs: + In some cases additional data can be specified or is even required. For + instance, :class:`~pde.grids.boundaries.axes.BoundariesList` expects + a `grid` (:class:`~pde.grids.base.GridBase`): to which the boundary + condition are associated, and it can use a `rank` (int), which sets + the tensorial rank of the field for this boundary condition. + """ + # check whether this is already of the correct type + if isinstance(data, BoundariesBase): + return data.__class__.from_data(data, **kwargs) + + # best guess based on the data: + if callable(data): + return BoundariesSetter.from_data(data) + else: + return BoundariesList.from_data(data, **kwargs) - def __init__(self, boundaries): + @classmethod + def get_help(cls) -> str: + """Return information on how boundary conditions can be set.""" + return ( + 'Boundary conditions for each axis are set using a dictionary: {"x": bc_x, ' + '"y-": bc_y_lower, "y+": bc_y_upper}. If the associated axis is periodic, ' + 'the boundary condition needs to be set to "periodic". Otherwise, ' + + BCBase.get_help() + ) + + +class BoundariesList(BoundariesBase): + """Defines boundary conditions for all axes individually.""" + + _logger = logging.getLogger(__qualname__) + + def __init__(self, boundaries: list[BoundaryAxisBase]): """Initialize with a list of boundaries.""" if len(boundaries) == 0: raise BCDataError("List of boundaries must not be empty") @@ -45,10 +126,10 @@ def __init__(self, boundaries): # check consistency for axis, boundary in enumerate(boundaries): if boundary.grid != self.grid: - raise BCDataError("Boundaries are not defined on the same grid") + raise BCDataError("BoundariesList are not defined on the same grid") if boundary.axis != axis: raise BCDataError( - "Boundaries need to be ordered like the respective axes" + "BoundariesList need to be ordered like the respective axes" ) if boundary.periodic != self.grid.periodic[axis]: raise PeriodicityError( @@ -58,123 +139,195 @@ def __init__(self, boundaries): ) # create the list of boundaries - super().__init__(boundaries) - - def __str__(self): - items = ", ".join(str(item) for item in self) - return f"[{items}]" + self._axes: list[BoundaryAxisBase] = boundaries @classmethod - def from_data(cls, grid: GridBase, boundaries, rank: int = 0) -> Boundaries: + def from_data( # type: ignore + cls, data, *, grid: GridBase, rank: int = 0, **kwargs + ) -> BoundariesList: """Creates all boundaries from given data. Args: grid (:class:`~pde.grids.base.GridBase`): The grid with which the boundary condition is associated - boundaries (str or list or tuple or dict): - Data that describes the boundaries. This can either be a list of - specifications for each dimension or a single one, which is then - applied to all dimensions. The boundary for a dimensions can be - specified by one of the following formats: - - * string specifying a single type for all boundaries - * dictionary specifying the type and values for all boundaries - * tuple pair specifying the low and high boundary individually + data (str or dict): + Data that describes the boundaries. This should either be a string + naming a boundary condition or a dictionary with detailed information. rank (int): The tensorial rank of the field for this boundary condition """ - # check whether this is already the correct class - if isinstance(boundaries, Boundaries): + # distinguish different possible data formats based on their type + if isinstance(data, BoundariesList): # boundaries are already in the correct format - if boundaries.grid._mesh is not None: - # we need to exclude this case since otherwise we get into a rabit hole + if data.grid._mesh is not None: + # we need to exclude this case since otherwise we get into a rabbit hole # where it is not clear what grid boundary conditions belong to. The # idea is that users only create boundary conditions for the full grid # and that the splitting onto subgrids is only done once, automatically, # and without involving calls to `from_data` raise ValueError("Cannot create MPI subgrid BC from data") - if boundaries.grid != grid: + if data.grid != grid: raise ValueError( "The grid of the supplied boundary condition is incompatible with " - f"the current grid ({boundaries.grid!r} != {grid!r})" + f"the current grid ({data.grid!r} != {grid!r})" ) - boundaries.check_value_rank(rank) - return boundaries - - # convert natural boundary conditions if present - if ( - boundaries == "natural" - or boundaries == "auto_periodic_neumann" - or boundaries == "auto_periodic_derivative" - ): - # set the respective natural conditions for all axes - boundaries = [] - for periodic in grid.periodic: - if periodic: - boundaries.append("periodic") - else: - boundaries.append("derivative") - - elif ( - boundaries == "auto_periodic_dirichlet" - or boundaries == "auto_periodic_value" - ): - # set the respective natural conditions (with vanishing values) for all axes - boundaries = [] - for periodic in grid.periodic: - if periodic: - boundaries.append("periodic") - else: - boundaries.append("value") - - elif boundaries == "auto_periodic_curvature": - # set the respective natural conditions (with vanishing curvature) for all axes - boundaries = [] - for periodic in grid.periodic: - if periodic: - boundaries.append("periodic") - else: - boundaries.append("curvature") - - # create the list of BoundaryAxis objects - bcs = None - if isinstance(boundaries, (str, dict)): - # one specification for all axes - bcs = [ - get_boundary_axis(grid, i, boundaries, rank=rank) - for i in range(grid.num_axes) - ] - - elif hasattr(boundaries, "__len__"): - # handle cases that look like sequences - if len(boundaries) == grid.num_axes: + data.check_value_rank(rank) + return data + + elif isinstance(data, BoundariesBase): + # data seems to be given as another base class, which indicates problems + raise TypeError( + "Can only use type `BoundariesList`. Use `BoundariesBase.from_data` " + "for more general data." + ) + + elif isinstance(data, str): + # a string implies the same boundary condition for all axes + + if data.startswith("auto_periodic_"): + # initialize boundary condition that could be periodic + bc = data[len("auto_periodic_") :] + bcs = [ + get_boundary_axis(grid, i, "periodic" if per else bc, rank=rank) + for i, per in enumerate(grid.periodic) + ] + + else: + # assume the same boundary condition for all axes + bcs = [ + get_boundary_axis(grid, i, data, rank=rank) + for i in range(grid.num_axes) + ] + + elif isinstance(data, dict): + # dictionaries can either specify boundary conditions for separate sides or + # they can specify a single boundary condition that is used on all sides + + if "low" in data or "high" in data: + # specifying conditions using this have been deprecated on 2024-11-23 + warnings.warn( + "Deprecated format for boundary conditions." + cls.get_help(), + DeprecationWarning, + ) + bcs = [ + get_boundary_axis(grid, i, data, rank=rank) + for i in range(grid.num_axes) + ] + + if _is_local_bc_data(data): + # detected identifier signifying that a single condition was specified + bcs = [ + get_boundary_axis(grid, i, data, rank=rank) + for i in range(grid.num_axes) + ] + + else: + # assume that boundary conditions are given for separate axes + + # initialize boundary data with wildcard default + data = data.copy() # we want to modify this dictionary + bc_all = data.pop("*", None) + bc_data = [[bc_all, bc_all] for _ in range(grid.num_axes)] + bc_seen = [[False, False] for _ in range(grid.num_axes)] + + # check specific boundary conditions for all axes + for ax, ax_name in enumerate(grid.axes): + # overwrite boundaries whose axes are given + if bc_axes := data.pop(ax_name, None): + bc_data[ax] = [bc_axes, bc_axes] + + # overwrite specific conditions for one side + if bc_lower := data.pop(ax_name + "-", None): + bc_data[ax][0] = bc_lower + bc_seen[ax][0] = True + if bc_upper := data.pop(ax_name + "+", None): + bc_data[ax][1] = bc_upper + bc_seen[ax][1] = True + + # overwrite conditions for named boundaries + for name, (ax, upper) in grid.boundary_names.items(): + if bc := data.pop(name, None): + i = int(upper) + if bc_seen[ax][i]: + cls._logger.warning( + "Duplicate BC data for axis %s%s", ax, "-+"[i] + ) + bc_data[ax][i] = bc + bc_seen[ax][i] = True + + # warn if some keys were left over + if data: + cls._logger.warning("Didn't use BC data from %s", list(data.keys())) + # find boundary conditions that have not been specified + bcs_unspecified = [] + for ax, bc_ax in enumerate(bc_data): + for i, bc_side in enumerate(bc_ax): + if bc_side is None: + bcs_unspecified.append(grid.axes[ax] + "-+"[int(i)]) + if bcs_unspecified: + cls._logger.warning("Didn't specified BCs for %s", bcs_unspecified) + + # create the actual boundary conditions + cls._logger.debug("Parsed BCs as %s", bc_data) + bcs = [ + get_boundary_axis(grid, i, tuple(boundary), rank=rank) + for i, boundary in enumerate(bc_data) + ] + + elif hasattr(data, "__len__"): + # sequences have been deprecated on 2024-11-23 + warnings.warn( + "Deprecated format for boundary conditions." + cls.get_help(), + DeprecationWarning, + ) + if len(data) == grid.num_axes: # assume that data is given for each boundary bcs = [ get_boundary_axis(grid, i, boundary, rank=rank) - for i, boundary in enumerate(boundaries) + for i, boundary in enumerate(data) ] - elif grid.num_axes == 1 and len(boundaries) == 2: + elif grid.num_axes == 1 and len(data) == 2: # special case where the two sides can be specified directly - bcs = [get_boundary_axis(grid, 0, boundaries, rank=rank)] + bcs = [get_boundary_axis(grid, 0, data, rank=rank)] + else: + raise BCDataError( + f"Got {len(data)} boundary conditions, but grid has " + f"{grid.num_axes} axes." + cls.get_help() + ) - if bcs is None: - # none of the logic worked + else: + # unknown format raise BCDataError( - f"Unsupported boundary format: `{boundaries}`. " + cls.get_help() + f"Unsupported boundary format: `{data}`. " + cls.get_help() ) - return cls(bcs) + return BoundariesList(bcs) + + def __str__(self): + items = ", ".join(str(item) for item in self) + return f"[{items}]" + + def __len__(self): + return len(self._axes) + + def __iter__(self) -> Iterator[BoundaryAxisBase]: + yield from self._axes def __eq__(self, other): - if not isinstance(other, Boundaries): + if not isinstance(other, BoundariesList): + return NotImplemented + return self.grid == other.grid and self._axes == other._axes + + def __ne__(self, other): + if not isinstance(other, BoundariesList): return NotImplemented - return super().__eq__(other) and self.grid == other.grid + return self.grid != other.grid or self._axes != other._axes @property def boundaries(self) -> Iterator[BCBase]: """Iterator over all non-periodic boundaries.""" - for boundary_axis in self: # iterate all axes + for boundary_axis in self._axes: # iterate all axes if not boundary_axis.periodic: # skip periodic axes yield from boundary_axis @@ -188,21 +341,12 @@ def check_value_rank(self, rank: int) -> None: Throws: RuntimeError: if any value does not have rank `rank` """ - for b in self: + for b in self._axes: b.check_value_rank(rank) - @classmethod - def get_help(cls) -> str: - """Return information on how boundary conditions can be set.""" - return ( - "Boundary conditions for each axis are set using a list: [bc_x, bc_y, " - "bc_z]. If the associated axis is periodic, the boundary condition needs " - f"to be set to 'periodic'. Otherwise, {BoundaryPair.get_help()}" - ) - - def copy(self) -> Boundaries: + def copy(self) -> BoundariesList: """Create a copy of the current boundaries.""" - return self.__class__([bc.copy() for bc in self]) + return self.__class__([bc.copy() for bc in self._axes]) @property def periodic(self) -> list[bool]: @@ -222,15 +366,25 @@ def __getitem__(self, index): if isinstance(index, str): # assume that the index is a known identifier if index in self.grid.boundary_names: + # found a known boundary axis, upper = self.grid.boundary_names[index] - return self[axis][upper] - elif index in self.grid.axes: - return self[self.grid.axes.index(index)] - else: - raise KeyError(index) + return self._axes[axis][upper] + + # check all axes + for ax, ax_name in enumerate(self.grid.axes): + if index == ax_name: + return self._axes[ax] + if index == ax_name + "-": + return self._axes[ax][False] + if index == ax_name + "+": + return self._axes[ax][True] + + # found nothing + raise KeyError(index) + else: # handle all other cases, in particular integer indices - return super().__getitem__(index) + return self._axes[index] def __setitem__(self, index, data) -> None: """Set specific boundary conditions. @@ -245,31 +399,41 @@ def __setitem__(self, index, data) -> None: """ if isinstance(index, str): # assume that the index is a known identifier + if index in self.grid.boundary_names: # set a specific boundary side - axis, upper = self.grid.boundary_names[index] - self[axis][upper] = data - - elif index in self.grid.axes: - # set both sides of the boundary condition - axis = self.grid.axes.index(index) - self[axis] = get_boundary_axis( - grid=self.grid, - axis=axis, - data=data, - rank=self[axis].rank, - ) + ax, upper = self.grid.boundary_names[index] + self._axes[ax][upper] = data + else: - raise KeyError(index) + # check all axes + for ax, ax_name in enumerate(self.grid.axes): + if index == ax_name: + # found just the axis -> set both sides + self._axes[ax] = get_boundary_axis( + grid=self.grid, axis=ax, data=data, rank=self[ax].rank + ) + break + if index == ax_name + "-": + # found lower part of the axis + self._axes[ax][False] = data + break + if index == ax_name + "+": + # found upper part of the axis + self._axes[ax][True] = data + break + + else: + raise KeyError(index) else: # handle all other cases, in particular integer indices - super().__setitem__(index, data) + self._axes[index] = data def get_mathematical_representation(self, field_name: str = "C") -> str: """Return mathematical representation of the boundary condition.""" - result = [] - for b in self: + result: list[str] = [] + for b in self._axes: try: result.extend(b.get_mathematical_representation(field_name)) except NotImplementedError: @@ -364,3 +528,97 @@ def wrap(data_full: np.ndarray, args=None) -> None: return wrap # type: ignore return chain(ghost_cell_setters) + + +class BoundariesSetter(BoundariesBase): + """Represents a function that sets ghost cells to determine boundary conditions. + + The function must have accept a :class:`~numpy.ndarray`, which contains the full + field data including the ghost points, and a second, optional argument, which is a + dictionary containing additional parameters, like the current time point `t` in case + of a simulation. + + Example: + Here is an example for a simple boundary setter, which sets specific boundary + conditions in the x-direction and periodic conditions in the y-direction of a + grid with two axes. Note that this boundary condition will not work for grids + with other number of axes and no additional checks are performed. + + .. code-block:: python + + def setter(data, args=None): + data[0, :] = data[1, :] # Vanishing derivative at left side + data[-1, :] = 2 - data[-2, :] # Fixed value `1` at right side + data[:, 0] = data[:, -2] # Periodic BC at top + data[:, -1] = data[:, 1] # Periodic BC at bottom + """ + + def __init__(self, setter: GhostCellSetter): + self._setter = setter + + @classmethod + def from_data(cls, data, **kwargs) -> BoundariesSetter: + """Creates all boundaries from given data. + + Args: + data (callable): + Function that sets the ghost cells + """ + # check whether this is already the correct class + if isinstance(data, BoundariesSetter): + # boundaries are already in the correct format + return data + + elif isinstance(data, BoundariesBase): + raise TypeError( + "Can only use type `BoundariesSetter`. Use `BoundariesBase.from_data` " + "for more general data." + ) + + return BoundariesSetter(data) + + def set_ghost_cells(self, data_full: np.ndarray, *, args=None) -> None: + """Set the ghost cells for all boundaries. + + Args: + data_full (:class:`~numpy.ndarray`): + The full field data including ghost points + set_corners (bool): + Determines whether the corner cells are set using interpolation + args: + Additional arguments that might be supported by special boundary + conditions. + """ + self._setter(data_full, args=args) + + def make_ghost_cell_setter(self) -> GhostCellSetter: + """Return function that sets the ghost cells on a full array. + + Returns: + Callable with signature :code:`(data_full: np.ndarray, args=None)`, which + sets the ghost cells of the full data, potentially using additional + information in `args` (e.g., the time `t` during solving a PDE) + """ + return jit(self._setter) # type: ignore + + +def set_default_bc( + bc_data: BoundariesData | None, default: BoundaryData +) -> BoundariesData: + """Set a default boundary condition. + + Args: + bc_data (str or list or tuple or dict or callable): + User-supplied data specifying boundary conditions + default: + Default condition that should be imposed where user conditions are not given + + Returns: + Modified `bc_data` with added defaults + """ + if bc_data is None: + bc_data = default + elif isinstance(bc_data, dict) and not _is_local_bc_data(bc_data): + # set default when boundary conditions for axes are specified + bc_data.setdefault("*", default) + return bc_data diff --git a/pde/grids/boundaries/axis.py b/pde/grids/boundaries/axis.py index 61805c45..4a65e578 100644 --- a/pde/grids/boundaries/axis.py +++ b/pde/grids/boundaries/axis.py @@ -1,12 +1,19 @@ -r""" -.. codeauthor:: David Zwicker - -This module handles the boundaries of a single axis of a grid. There are -generally only two options, depending on whether the axis of the underlying -grid is defined as periodic or not. If it is periodic, the class -:class:`~pde.grids.boundaries.axis.BoundaryPeriodic` should be used, while -non-periodic axes have more option, which are represented by +r"""This module handles the boundaries of a single axis of a grid. There are generally +only two options, depending on whether the axis of the underlying grid is defined as +periodic or not. If it is periodic, the class +:class:`~pde.grids.boundaries.axis.BoundaryPeriodic` should be used, while non-periodic +axes have more option, which are represented by :class:`~pde.grids.boundaries.axis.BoundaryPair`. + +.. autosummary:: + :nosignatures: + + ~BoundaryAxisBase + ~BoundaryPair + ~BoundaryPeriodic + + +.. codeauthor:: David Zwicker """ from __future__ import annotations @@ -89,6 +96,15 @@ def copy(self) -> BoundaryAxisBase: """Return a copy of itself, but with a reference to the same grid.""" return self.__class__(self.low.copy(), self.high.copy()) + def check_value_rank(self, rank: int) -> None: + """Check whether the values at the boundaries have the correct rank. + + Args: + rank (int): + The tensorial rank of the field for this boundary condition + """ + raise NotImplementedError + def __eq__(self, other): if not isinstance(other, self.__class__): return NotImplemented @@ -286,11 +302,16 @@ def from_data(cls, grid: GridBase, axis: int, data, rank: int = 0) -> BoundaryPa except TypeError: # if len is not supported, the format must be wrong raise BCDataError( - f"Unsupported boundary format: `{data}`. " + cls.get_help() + f"Unsupported boundary format: `{data}`." + cls.get_help() ) from None else: if data_len == 2: # assume that data is given for each boundary + if data[0] == "periodic" or data[1] == "periodic": + raise BCDataError( + f"Only one side of {grid.axes[axis]} axis was set to have " + "periodic boundary conditions." + ) low = BCBase.from_data( grid, axis, upper=False, data=data[0], rank=rank ) @@ -387,24 +408,9 @@ def get_boundary_axis( :class:`~pde.grids.boundaries.axis.BoundaryAxisBase`: Appropriate boundary condition for the axis """ - # handle special constructs that describe boundary conditions - if ( - data == "natural" - or data == "auto_periodic_neumann" - or data == "auto_periodic_derivative" - ): - # automatic choice between periodic and Neumann condition - if grid.periodic[axis]: - data = "periodic" - else: - data = "derivative" - - elif data == "auto_periodic_dirichlet" or data == "auto_periodic_value": - # automatic choice between periodic and Dirichlet condition - if grid.periodic[axis]: - data = "periodic" - else: - data = "value" + # handle special case describing potentially periodic boundary conditions + if isinstance(data, str) and data.startswith("auto_periodic_"): + data = "periodic" if grid.periodic[axis] else data[len("auto_periodic_") :] # handle different types of data that specify boundary conditions if isinstance(data, BoundaryAxisBase): diff --git a/pde/grids/boundaries/local.py b/pde/grids/boundaries/local.py index 20a1c42d..eeb9d479 100644 --- a/pde/grids/boundaries/local.py +++ b/pde/grids/boundaries/local.py @@ -1,7 +1,4 @@ -r""" -.. codeauthor:: David Zwicker - -This module contains classes for handling a single boundary of a non-periodic axis. +r"""This module contains classes for handling a single boundary of a non-periodic axis. Since an axis has two boundary, we simply distinguish them by a boolean flag `upper`, which is `True` for the side of the axis with the larger coordinate. @@ -55,6 +52,8 @@ .. inheritance-diagram:: pde.grids.boundaries.local :parts: 1 + +.. codeauthor:: David Zwicker """ from __future__ import annotations @@ -531,8 +530,13 @@ def from_data( # create a specific condition given by a string bc = cls.from_str(grid, axis, upper=upper, condition=data, rank=rank) + elif data is None: + raise BCDataError( + f"Unspecified condition for boundary {grid.axes[axis]}{'-+'[int(upper)]}" + ) + else: - raise BCDataError(f"Unsupported format: `{data}`. " + cls.get_help()) + raise BCDataError(f"Unsupported BC format: `{data}`. " + cls.get_help()) # check consistency if bc.periodic != grid.periodic[axis]: diff --git a/pde/grids/cartesian.py b/pde/grids/cartesian.py index 187561de..9431bf82 100644 --- a/pde/grids/cartesian.py +++ b/pde/grids/cartesian.py @@ -7,7 +7,7 @@ import itertools from collections.abc import Generator, Sequence -from typing import TYPE_CHECKING, Any +from typing import Any import numpy as np @@ -22,9 +22,6 @@ ) from .coordinates import CartesianCoordinates -if TYPE_CHECKING: - from .boundaries.axes import Boundaries, BoundariesData - class CartesianGrid(GridBase): r"""D-dimensional Cartesian grid with uniform discretization for each axis. @@ -54,15 +51,6 @@ class CartesianGrid(GridBase): cuboid: Cuboid - boundary_names = { # name all the boundaries - "left": (0, False), - "right": (0, True), - "bottom": (1, False), - "top": (1, True), - "back": (2, False), - "front": (2, True), - } - def __init__( self, bounds: Sequence[tuple[float, float]], @@ -138,6 +126,13 @@ def __init__( self._axes_coords = tuple(axes_coords) self._axes_bounds = tuple(self.cuboid.bounds) + # name all the boundaries + self.boundary_names = {"left": (0, False), "right": (0, True)} + if self.num_axes > 1: + self.boundary_names.update({"bottom": (1, False), "top": (1, True)}) + if self.num_axes > 2: + self.boundary_names.update({"back": (2, False), "front": (2, True)}) + @property def state(self) -> dict[str, Any]: """dict: the state of the grid""" diff --git a/pde/grids/cylindrical.py b/pde/grids/cylindrical.py index a684467c..7e2b85f4 100644 --- a/pde/grids/cylindrical.py +++ b/pde/grids/cylindrical.py @@ -22,7 +22,6 @@ from .coordinates import CylindricalCoordinates if TYPE_CHECKING: - from .boundaries.axes import Boundaries, BoundariesData from .spherical import PolarSymGrid diff --git a/pde/grids/operators/cartesian.py b/pde/grids/operators/cartesian.py index b8fbda87..e0abf8dc 100644 --- a/pde/grids/operators/cartesian.py +++ b/pde/grids/operators/cartesian.py @@ -26,7 +26,8 @@ from ...tools.misc import module_available from ...tools.numba import jit from ...tools.typing import OperatorType -from ..boundaries import Boundaries +from ..boundaries.axes import BoundariesBase, BoundariesList +from ..boundaries.axis import BoundaryAxisBase from ..cartesian import CartesianGrid from .common import make_derivative as _make_derivative from .common import make_derivative2 as _make_derivative2 @@ -36,11 +37,11 @@ # deprecated since 2023-12-06 -def _get_laplace_matrix_1d(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix_1d(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for Laplace operator on a 1d Cartesian grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -78,11 +79,11 @@ def _get_laplace_matrix_1d(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: return matrix, vector -def _get_laplace_matrix_2d(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix_2d(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for Laplace operator on a 2d Cartesian grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -147,11 +148,11 @@ def i(x, y): return matrix, vector -def _get_laplace_matrix_3d(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix_3d(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for Laplace operator on a 3d Cartesian grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -234,11 +235,11 @@ def i(x, y, z): return matrix, vector -def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for Laplace operator on a Cartesian grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -1286,12 +1287,12 @@ def make_tensor_divergence( @CartesianGrid.register_operator("poisson_solver", rank_in=0, rank_out=0) def make_poisson_solver( - bcs: Boundaries, *, method: Literal["auto", "scipy"] = "auto" + bcs: BoundariesList, *, method: Literal["auto", "scipy"] = "auto" ) -> OperatorType: """Make a operator that solves Poisson's equation. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} method (str): Method used for calculating the tensor divergence operator. diff --git a/pde/grids/operators/cylindrical_sym.py b/pde/grids/operators/cylindrical_sym.py index 0594fe80..5b9418f7 100644 --- a/pde/grids/operators/cylindrical_sym.py +++ b/pde/grids/operators/cylindrical_sym.py @@ -25,16 +25,16 @@ from ...tools.docstrings import fill_in_docstring from ...tools.numba import jit from ...tools.typing import OperatorType -from ..boundaries import Boundaries +from ..boundaries.axes import BoundariesBase, BoundariesList from ..cylindrical import CylindricalSymGrid from .common import make_general_poisson_solver -def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for Laplace operator on a cylindrical grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -432,14 +432,14 @@ def tensor_divergence(arr: np.ndarray, out: np.ndarray) -> None: @CylindricalSymGrid.register_operator("poisson_solver", rank_in=0, rank_out=0) @fill_in_docstring def make_poisson_solver( - bcs: Boundaries, *, method: Literal["auto", "scipy"] = "auto" + bcs: BoundariesList, *, method: Literal["auto", "scipy"] = "auto" ) -> OperatorType: """Make a operator that solves Poisson's equation. {DESCR_CYLINDRICAL_GRID} Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} method (str): The chosen method for implementing the operator diff --git a/pde/grids/operators/polar_sym.py b/pde/grids/operators/polar_sym.py index 078c2120..60c5207a 100644 --- a/pde/grids/operators/polar_sym.py +++ b/pde/grids/operators/polar_sym.py @@ -22,7 +22,7 @@ from ...tools.docstrings import fill_in_docstring from ...tools.numba import jit from ...tools.typing import OperatorType -from ..boundaries import Boundaries +from ..boundaries.axes import BoundariesBase, BoundariesList from ..spherical import PolarSymGrid from .common import make_general_poisson_solver @@ -266,11 +266,11 @@ def tensor_divergence(arr: np.ndarray, out: np.ndarray) -> None: @fill_in_docstring -def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for laplace operator on a polar grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -325,14 +325,14 @@ def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: @PolarSymGrid.register_operator("poisson_solver", rank_in=0, rank_out=0) @fill_in_docstring def make_poisson_solver( - bcs: Boundaries, *, method: Literal["auto", "scipy"] = "auto" + bcs: BoundariesList, *, method: Literal["auto", "scipy"] = "auto" ) -> OperatorType: """Make a operator that solves Poisson's equation. {DESCR_POLAR_GRID} Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} method (str): The chosen method for implementing the operator diff --git a/pde/grids/operators/spherical_sym.py b/pde/grids/operators/spherical_sym.py index f5972b47..db256683 100644 --- a/pde/grids/operators/spherical_sym.py +++ b/pde/grids/operators/spherical_sym.py @@ -22,7 +22,7 @@ from ...tools.docstrings import fill_in_docstring from ...tools.numba import jit from ...tools.typing import OperatorType -from ..boundaries import Boundaries +from ..boundaries.axes import BoundariesBase, BoundariesList from ..spherical import SphericalSymGrid from .common import make_general_poisson_solver @@ -578,11 +578,11 @@ def tensor_double_divergence(arr: np.ndarray, out: np.ndarray) -> None: @fill_in_docstring -def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: +def _get_laplace_matrix(bcs: BoundariesList) -> tuple[np.ndarray, np.ndarray]: """Get sparse matrix for laplace operator on a spherical grid. Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} Returns: @@ -642,14 +642,14 @@ def _get_laplace_matrix(bcs: Boundaries) -> tuple[np.ndarray, np.ndarray]: @SphericalSymGrid.register_operator("poisson_solver", rank_in=0, rank_out=0) @fill_in_docstring def make_poisson_solver( - bcs: Boundaries, *, method: Literal["auto", "scipy"] = "auto" + bcs: BoundariesList, *, method: Literal["auto", "scipy"] = "auto" ) -> OperatorType: """Make a operator that solves Poisson's equation. {DESCR_POLAR_GRID} Args: - bcs (:class:`~pde.grids.boundaries.axes.Boundaries`): + bcs (:class:`~pde.grids.boundaries.axes.BoundariesList`): {ARG_BOUNDARIES_INSTANCE} method (str): The chosen method for implementing the operator diff --git a/pde/grids/spherical.py b/pde/grids/spherical.py index 3b9d2732..8b94986b 100644 --- a/pde/grids/spherical.py +++ b/pde/grids/spherical.py @@ -10,7 +10,7 @@ from __future__ import annotations from abc import ABCMeta -from typing import TYPE_CHECKING, Any, Literal, TypeVar +from typing import Any, Literal, TypeVar import numpy as np @@ -20,10 +20,6 @@ from .cartesian import CartesianGrid from .coordinates import PolarCoordinates, SphericalCoordinates -if TYPE_CHECKING: - from .boundaries.axes import Boundaries - - TNumArr = TypeVar("TNumArr", float, np.ndarray) diff --git a/pde/pdes/allen_cahn.py b/pde/pdes/allen_cahn.py index 2b4d0106..f1bb6d61 100644 --- a/pde/pdes/allen_cahn.py +++ b/pde/pdes/allen_cahn.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -30,6 +31,8 @@ class AllenCahnPDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" interface_width: float @@ -38,7 +41,8 @@ def __init__( self, interface_width: float = 1, mobility: float = 1, - bc: BoundariesData = "auto_periodic_neumann", + *, + bc: BoundariesData | None = None, ): """ Args: @@ -54,7 +58,7 @@ def __init__( self.interface_width = interface_width self.mobility = mobility - self.bc = bc + self.bc = set_default_bc(bc, self.default_bc) @property def expression(self) -> str: diff --git a/pde/pdes/cahn_hilliard.py b/pde/pdes/cahn_hilliard.py index 8fd870a9..4d31fec4 100644 --- a/pde/pdes/cahn_hilliard.py +++ b/pde/pdes/cahn_hilliard.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -30,13 +31,18 @@ class CahnHilliardPDE(PDEBase): """ explicit_time_dependence = False + default_bc_c = "auto_periodic_neumann" + """Default boundary condition for order parameter.""" + default_bc_mu = "auto_periodic_neumann" + """Default boundary condition for chemical potential.""" @fill_in_docstring def __init__( self, interface_width: float = 1, - bc_c: BoundariesData = "auto_periodic_neumann", - bc_mu: BoundariesData = "auto_periodic_neumann", + *, + bc_c: BoundariesData | None = None, + bc_mu: BoundariesData | None = None, ): """ Args: @@ -54,8 +60,8 @@ def __init__( super().__init__() self.interface_width = interface_width - self.bc_c = bc_c - self.bc_mu = bc_mu + self.bc_c = set_default_bc(bc_c, self.default_bc_c) + self.bc_mu = set_default_bc(bc_mu, self.default_bc_mu) @property def expression(self) -> str: diff --git a/pde/pdes/diffusion.py b/pde/pdes/diffusion.py index a3032448..71473129 100644 --- a/pde/pdes/diffusion.py +++ b/pde/pdes/diffusion.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -29,13 +30,15 @@ class DiffusionPDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" @fill_in_docstring def __init__( self, diffusivity: float = 1, *, - bc: BoundariesData = "auto_periodic_neumann", + bc: BoundariesData | None = None, noise: float = 0, rng: np.random.Generator | None = None, ): @@ -59,7 +62,7 @@ def __init__( super().__init__(noise=noise, rng=rng) self.diffusivity = diffusivity - self.bc = bc + self.bc = set_default_bc(bc, self.default_bc) @property def expression(self) -> str: diff --git a/pde/pdes/kpz_interface.py b/pde/pdes/kpz_interface.py index 6d1197c2..c6c8e51b 100644 --- a/pde/pdes/kpz_interface.py +++ b/pde/pdes/kpz_interface.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -33,6 +34,8 @@ class KPZInterfacePDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" @fill_in_docstring def __init__( @@ -40,7 +43,7 @@ def __init__( nu: float = 0.5, lmbda: float = 1, *, - bc: BoundariesData = "auto_periodic_neumann", + bc: BoundariesData | None = None, noise: float = 0, rng: np.random.Generator | None = None, ): @@ -67,7 +70,7 @@ def __init__( self.nu = nu self.lmbda = lmbda - self.bc = bc + self.bc = set_default_bc(bc, self.default_bc) @property def expression(self) -> str: diff --git a/pde/pdes/kuramoto_sivashinsky.py b/pde/pdes/kuramoto_sivashinsky.py index 700411c6..d669c749 100644 --- a/pde/pdes/kuramoto_sivashinsky.py +++ b/pde/pdes/kuramoto_sivashinsky.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -32,13 +33,15 @@ class KuramotoSivashinskyPDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" @fill_in_docstring def __init__( self, nu: float = 1, *, - bc: BoundariesData = "auto_periodic_neumann", + bc: BoundariesData | None = None, bc_lap: BoundariesData | None = None, noise: float = 0, rng: np.random.Generator | None = None, @@ -51,10 +54,9 @@ def __init__( The boundary conditions applied to the field. {ARG_BOUNDARIES} bc_lap: - The boundary conditions applied to the second derivative of the - scalar field :math:`c`. If `None`, the same boundary condition - as `bc` is chosen. Otherwise, this supports the same options as - `bc`. + The boundary conditions applied to the second derivative of the scalar + field :math:`c`. If `None`, the same boundary condition as `bc` is + chosen. Otherwise, this supports the same options as `bc`. noise (float): Variance of the (additive) noise term rng (:class:`~numpy.random.Generator`): @@ -68,8 +70,8 @@ def __init__( super().__init__(noise=noise, rng=rng) self.nu = nu - self.bc = bc - self.bc_lap = bc if bc_lap is None else bc_lap + self.bc = set_default_bc(bc, self.default_bc) + self.bc_lap = self.bc if bc_lap is None else bc_lap @property def expression(self) -> str: diff --git a/pde/pdes/laplace.py b/pde/pdes/laplace.py index da6e17ae..abc9581d 100644 --- a/pde/pdes/laplace.py +++ b/pde/pdes/laplace.py @@ -15,6 +15,7 @@ def solve_poisson_equation( rhs: ScalarField, bc: BoundariesData, + *, label: str = "Solution to Poisson's equation", **kwargs, ) -> ScalarField: @@ -81,7 +82,7 @@ def solve_poisson_equation( @fill_in_docstring def solve_laplace_equation( - grid: GridBase, bc: BoundariesData, label: str = "Solution to Laplace's equation" + grid: GridBase, bc: BoundariesData, *, label: str = "Solution to Laplace's equation" ) -> ScalarField: """Solve Laplace's equation on a given grid. diff --git a/pde/pdes/pde.py b/pde/pdes/pde.py index 5c52a788..a7913763 100644 --- a/pde/pdes/pde.py +++ b/pde/pdes/pde.py @@ -19,6 +19,7 @@ from ..fields import FieldCollection, VectorField from ..fields.base import FieldBase from ..fields.datafield_base import DataFieldBase +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..grids.boundaries.local import BCDataError from ..pdes.base import PDEBase, TState @@ -59,12 +60,15 @@ class PDE(PDEBase): appear in the `state`. """ + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" + @fill_in_docstring def __init__( self, rhs: dict[str, str], *, - bc: BoundariesData = "auto_periodic_neumann", + bc: BoundariesData | None = None, bc_ops: dict[str, BoundariesData] | None = None, user_funcs: dict[str, Callable] | None = None, consts: dict[str, NumberOrArray] | None = None, @@ -196,6 +200,7 @@ def __init__( self.complex_valued = complex_valued # setup boundary conditions + bc = set_default_bc(bc, self.default_bc) if bc_ops is None: bcs = {"*:*": bc} elif isinstance(bc_ops, dict): diff --git a/pde/pdes/swift_hohenberg.py b/pde/pdes/swift_hohenberg.py index c14fb1b6..3e37ec2a 100644 --- a/pde/pdes/swift_hohenberg.py +++ b/pde/pdes/swift_hohenberg.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -32,6 +33,8 @@ class SwiftHohenbergPDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" @fill_in_docstring def __init__( @@ -40,7 +43,7 @@ def __init__( kc2: float = 1.0, delta: float = 1.0, *, - bc: BoundariesData = "auto_periodic_neumann", + bc: BoundariesData | None = None, bc_lap: BoundariesData | None = None, ): r""" @@ -64,8 +67,8 @@ def __init__( self.rate = rate self.kc2 = kc2 self.delta = delta - self.bc = bc - self.bc_lap = bc if bc_lap is None else bc_lap + self.bc = set_default_bc(bc, self.default_bc) + self.bc_lap = self.bc if bc_lap is None else bc_lap @property def expression(self) -> str: diff --git a/pde/pdes/wave.py b/pde/pdes/wave.py index 2f986bc6..2b731770 100644 --- a/pde/pdes/wave.py +++ b/pde/pdes/wave.py @@ -11,6 +11,7 @@ import numpy as np from ..fields import FieldCollection, ScalarField +from ..grids.boundaries import set_default_bc from ..grids.boundaries.axes import BoundariesData from ..tools.docstrings import fill_in_docstring from ..tools.numba import jit @@ -34,9 +35,11 @@ class WavePDE(PDEBase): """ explicit_time_dependence = False + default_bc = "auto_periodic_neumann" + """Default boundary condition used when no specific conditions are chosen.""" @fill_in_docstring - def __init__(self, speed: float = 1, bc: BoundariesData = "auto_periodic_neumann"): + def __init__(self, speed: float = 1, *, bc: BoundariesData | None = None): """ Args: speed (float): @@ -48,7 +51,7 @@ def __init__(self, speed: float = 1, bc: BoundariesData = "auto_periodic_neumann super().__init__() self.speed = speed - self.bc = bc + self.bc = set_default_bc(bc, self.default_bc) def get_initial_condition(self, u: ScalarField, v: ScalarField | None = None): """Create a suitable initial condition. diff --git a/pde/tools/docstrings.py b/pde/tools/docstrings.py index 95114017..9e919df0 100644 --- a/pde/tools/docstrings.py +++ b/pde/tools/docstrings.py @@ -21,21 +21,21 @@ # description of function arguments "ARG_BOUNDARIES_INSTANCE": """ Specifies the boundary conditions applied to the field. This must be an - instance of :class:`~pde.grids.boundaries.axes.Boundaries`, which can be + instance of :class:`~pde.grids.boundaries.axes.BoundariesList`, which can be created from various data formats using the class method - :func:`~pde.grids.boundaries.axes.Boundaries.from_data`. + :meth:`~pde.grids.boundaries.axes.BoundariesBase.from_data`. """, "ARG_BOUNDARIES": """ - Boundary conditions are generally given as a list with one condition for each - axis. For periodic axes, only periodic boundary conditions are allowed + Boundary conditions are generally given as a dictionary with one condition for + each axis side. For periodic axes, only periodic boundary conditions are allowed (indicated by 'periodic' and 'anti-periodic'). For non-periodic axes, different - boundary conditions can be specified for the lower and upper end (using a tuple - of two conditions). For instance, Dirichlet conditions enforcing a value NUM - (specified by `{'value': NUM}`) and Neumann conditions enforcing the value DERIV - for the derivative in the normal direction (specified by `{'derivative': - DERIV}`) are supported. Note that the special value 'natural' imposes periodic - boundary conditions for periodic axis and a vanishing derivative otherwise. - More information can be found in the + boundary conditions can be specified for the lower and upper end (using specific + identifiers, like `x-` and `y+`). For instance, Dirichlet conditions enforcing a + value NUM (specified by `{'value': NUM}`) and Neumann conditions enforcing the + value DERIV for the derivative in the normal direction (specified by + `{'derivative': DERIV}`) are supported. Note that the special value + 'auto_periodic_neumann' imposes periodic boundary conditions for periodic axis + and a vanishing derivative otherwise. More information can be found in the :ref:`boundaries documentation `. """, "ARG_TRACKER_INTERRUPT": """ @@ -43,7 +43,7 @@ determines an interval (measured in the simulation time unit) of regular interruption. A string is interpreted as a duration in real time assuming the format 'hh:mm:ss'. A list of numbers is taken as explicit simulation time points. - More fine-grained contol is possible by passing an instance of classes defined + More fine-grained control is possible by passing an instance of classes defined in :mod:`~pde.trackers.interrupts`. """, "ARG_PLOT_QUANTITIES": """ diff --git a/scripts/format_code.sh b/scripts/format_code.sh index 8608a544..07d72676 100755 --- a/scripts/format_code.sh +++ b/scripts/format_code.sh @@ -6,11 +6,11 @@ pushd .. > /dev/null find . -name '*.py' -exec pyupgrade --py39-plus {} + popd > /dev/null -echo "Formating import statements..." +echo "Formatting import statements..." ruff check --fix --config=../pyproject.toml .. -echo "Formating docstrings..." +echo "Formatting docstrings..." docformatter --in-place --black --recursive .. -echo "Formating source code..." +echo "Formatting source code..." ruff format --config=../pyproject.toml .. \ No newline at end of file diff --git a/scripts/performance_laplace.py b/scripts/performance_laplace.py index 68201176..9896f78f 100755 --- a/scripts/performance_laplace.py +++ b/scripts/performance_laplace.py @@ -12,7 +12,7 @@ import numpy as np from pde import CylindricalSymGrid, ScalarField, SphericalSymGrid, UnitGrid, config -from pde.grids.boundaries import Boundaries +from pde.grids.boundaries import BoundariesBase from pde.tools.misc import estimate_computation_speed from pde.tools.numba import jit @@ -193,7 +193,7 @@ def main(): grid = CylindricalSymGrid(shape[0], [0, shape[1]], shape) print(f"Cylindrical grid, shape={shape}") field = ScalarField.random_normal(grid) - bcs = Boundaries.from_data(grid, "derivative") + bcs = BoundariesBase.from_data(grid, "derivative") expected = field.laplace(bcs) for method in ["CUSTOM", "numba"]: @@ -214,7 +214,7 @@ def main(): grid = SphericalSymGrid(shape, shape) print(grid) field = ScalarField.random_normal(grid) - bcs = Boundaries.from_data(grid, "derivative") + bcs = BoundariesBase.from_data(grid, "derivative") for conservative in [True, False]: laplace = grid.make_operator("laplace", bcs, conservative=conservative) diff --git a/tests/fields/test_generic_fields.py b/tests/fields/test_generic_fields.py index 9c395f7a..b1442eab 100644 --- a/tests/fields/test_generic_fields.py +++ b/tests/fields/test_generic_fields.py @@ -42,7 +42,7 @@ def test_set_label(field_class): @pytest.mark.parametrize("grid", iter_grids()) @pytest.mark.parametrize("field_class", [ScalarField, Tensor2Field]) def test_interpolation_natural(grid, field_class, rng): - """Test some interpolation for natural boundary conditions.""" + """Test some interpolation for auto_periodic_neumann boundary conditions.""" msg = f"grid={grid}, field={field_class}" f = field_class.random_uniform(grid, rng=rng) diff --git a/tests/grids/boundaries/test_axes_boundaries.py b/tests/grids/boundaries/test_axes_boundaries.py index 2a889488..0eb4cf86 100644 --- a/tests/grids/boundaries/test_axes_boundaries.py +++ b/tests/grids/boundaries/test_axes_boundaries.py @@ -9,7 +9,12 @@ from pde import ScalarField, UnitGrid from pde.grids.base import PeriodicityError -from pde.grids.boundaries.axes import BCDataError, Boundaries +from pde.grids.boundaries.axes import ( + BCDataError, + BoundariesBase, + BoundariesList, + BoundariesSetter, +) from pde.grids.boundaries.axis import BoundaryPair, BoundaryPeriodic, get_boundary_axis from pde.grids.boundaries.local import NeumannBC @@ -21,7 +26,7 @@ def test_boundaries(): periodic = [b == "periodic" for b in (bx, by)] g = UnitGrid([2, 2], periodic=periodic) - bcs = Boundaries.from_data(g, [bx, by]) + bcs = BoundariesBase.from_data({"x": bx, "y": by}, grid=g) bc_x = get_boundary_axis(g, 0, bx) bc_y = get_boundary_axis(g, 1, by) @@ -33,16 +38,18 @@ def test_boundaries(): assert isinstance(str(bcs), str) assert isinstance(repr(bcs), str) - assert bcs == Boundaries.from_data(g, [bc_x, bc_y]) + assert bcs == BoundariesBase.from_data([bc_x, bc_y], grid=g) if bx == by: - assert bcs == Boundaries.from_data(g, bx) + assert bcs == BoundariesBase.from_data(bx, grid=g) bc2 = bcs.copy() assert bcs == bc2 assert bcs is not bc2 - b1 = Boundaries.from_data(UnitGrid([2, 2]), "auto_periodic_neumann") - b2 = Boundaries.from_data(UnitGrid([3, 3]), "auto_periodic_neumann") + b1 = BoundariesBase.from_data("auto_periodic_neumann", grid=UnitGrid([2, 2])) + b2 = BoundariesBase.from_data("auto_periodic_neumann", grid=UnitGrid([3, 3])) + assert isinstance(b1, BoundariesList) + assert isinstance(b2, BoundariesList) assert b1 != b2 @@ -51,32 +58,34 @@ def test_boundaries_edge_cases(): grid = UnitGrid([3, 3]) bcs = grid.get_boundary_conditions("auto_periodic_neumann") with pytest.raises(BCDataError): - Boundaries([]) + BoundariesList([]) with pytest.raises(BCDataError): - Boundaries([bcs[0]]) + BoundariesList([bcs[0]]) with pytest.raises(BCDataError): - Boundaries([bcs[0], bcs[0]]) + BoundariesList([bcs[0], bcs[0]]) - assert bcs == Boundaries([bcs[0], bcs[1]]) + assert bcs == BoundariesList([bcs[0], bcs[1]]) bc0 = get_boundary_axis(grid.copy(), 0, "auto_periodic_neumann") - assert bcs == Boundaries([bc0, bcs[1]]) + assert bcs == BoundariesList([bc0, bcs[1]]) bc0 = get_boundary_axis(UnitGrid([4, 3]), 0, "auto_periodic_neumann") with pytest.raises(BCDataError): - Boundaries([bc0, bcs[1]]) + BoundariesList([bc0, bcs[1]]) bc0 = get_boundary_axis(UnitGrid([3, 3], periodic=True), 0, "auto_periodic_neumann") with pytest.raises(BCDataError): - Boundaries([bc0, bcs[1]]) + BoundariesList([bc0, bcs[1]]) def test_boundary_specifications(): """Test different ways of specifying boundary conditions.""" g = UnitGrid([2]) - bc1 = Boundaries.from_data( - g, [{"type": "derivative", "value": 0}, {"type": "value", "value": 0}] - ) - assert bc1 == Boundaries.from_data(g, [{"type": "derivative"}, {"type": "value"}]) - assert bc1 == Boundaries.from_data(g, [{"derivative": 0}, {"value": 0}]) - assert bc1 == Boundaries.from_data(g, ["neumann", "dirichlet"]) + ref = BoundariesBase.from_data({"x-": "neumann", "x+": "dirichlet"}, grid=g) + BC_DATA = [ + {"x-": {"type": "derivative", "value": 0}, "x+": {"type": "value", "value": 0}}, + {"left": {"type": "derivative"}, "right": {"type": "value"}}, + {"x-": {"derivative": 0}, "x+": {"value": 0}}, + ] + for data in BC_DATA: + assert ref == BoundariesBase.from_data(data, grid=g) def test_mixed_boundary_condition(rng): @@ -99,8 +108,8 @@ def test_natural_boundary_conditions(cond, is_value): """Test special automatic boundary conditions.""" g = UnitGrid([2, 2], periodic=[True, False]) for bc in [ - Boundaries.from_data(g, cond), - Boundaries.from_data(g, ["periodic", cond]), + BoundariesBase.from_data(cond, grid=g), + BoundariesBase.from_data(["periodic", cond], grid=g), ]: assert isinstance(bc[0], BoundaryPeriodic) if is_value: @@ -120,7 +129,7 @@ def test_special_cases(): def test_bc_values(): """Test setting the values of boundary conditions.""" g = UnitGrid([5]) - bc = g.get_boundary_conditions([{"value": 2}, {"derivative": 3}]) + bc = g.get_boundary_conditions({"x-": {"value": 2}, "x+": {"derivative": 3}}) assert bc[0].low.value == 2 assert bc[0].high.value == 3 @@ -159,9 +168,13 @@ def test_setting_specific_bcs(): assert str(bcs["x"]) == '"value"' bcs["left"] = "derivative" assert str(bcs["left"]) == '"derivative"' + assert str(bcs["x-"]) == '"derivative"' assert str(bcs["right"]) == '"value"' + assert str(bcs["x+"]) == '"value"' bcs["right"] = "derivative" assert str(bcs["x"]) == '"derivative"' + bcs["x-"] = bcs["x+"] = "value" + assert str(bcs["x"]) == '"value"' with pytest.raises(PeriodicityError): bcs["x"] = "periodic" @@ -174,6 +187,8 @@ def test_setting_specific_bcs(): bcs["y"] = "value" with pytest.raises(PeriodicityError): bcs["top"] = "value" + with pytest.raises(PeriodicityError): + bcs["y+"] = "value" # test wrong input with pytest.raises(KeyError): @@ -187,13 +202,56 @@ def test_setting_specific_bcs(): def test_boundaries_property(): """Test boundaries property.""" g = UnitGrid([2, 2]) - bc = Boundaries.from_data(g, ["neumann", "dirichlet"]) + bc = BoundariesBase.from_data({"x": "neumann", "y": "dirichlet"}, grid=g) assert len(list(bc.boundaries)) == 4 - bc = Boundaries.from_data(g, "neumann") + bc = BoundariesBase.from_data("neumann", grid=g) for b in bc.boundaries: assert isinstance(b, NeumannBC) g = UnitGrid([2, 2], periodic=[True, False]) - bc = Boundaries.from_data(g, "auto_periodic_neumann") + bc = BoundariesBase.from_data("auto_periodic_neumann", grid=g) assert len(list(bc.boundaries)) == 2 + + +@pytest.mark.parametrize("periodic", [True, False]) +def test_boundaries_setter_1d(periodic, rng): + """Test BoundariesSetter class for 1d grids.""" + + def setter(data, args=None): + if periodic: + data[0] = data[-2] + data[-1] = data[1] + else: + data[0] = data[1] # Neumann + data[-1] = -data[-2] # Dirichlet + + f1 = ScalarField.random_normal(UnitGrid([4], periodic=periodic)) + f2 = f1.copy() + + f1.set_ghost_cells(bc=BoundariesSetter(setter)) + if periodic: + f2.set_ghost_cells(bc="periodic") + else: + f2.set_ghost_cells(bc={"x-": "neumann", "x+": "dirichlet"}) + np.testing.assert_allclose(f1._data_full, f2._data_full) + + +def test_boundaries_setter_2d(rng): + """Test BoundariesSetter class for 2d grids.""" + + def setter(data, args=None): + data[0, :] = data[1, :] # Neumann + data[-1, :] = -data[-2, :] # Dirichlet + data[:, 0] = data[:, -2] # periodic + data[:, -1] = data[:, 1] # periodic + + f1 = ScalarField.random_normal(UnitGrid([4, 4], periodic=[False, True])) + f2 = f1.copy() + + f1.set_ghost_cells(bc=BoundariesSetter(setter)) + f2.set_ghost_cells(bc={"x-": "neumann", "x+": "dirichlet", "y": "periodic"}) + # compare full fields without corner points + mask = np.ones((6, 6), dtype=bool) + mask[0, 0] = mask[-1, 0] = mask[0, -1] = mask[-1, -1] = False + np.testing.assert_allclose(f1._data_full[mask], f2._data_full[mask]) diff --git a/tests/grids/boundaries/test_axes_boundaries_legacy.py b/tests/grids/boundaries/test_axes_boundaries_legacy.py new file mode 100644 index 00000000..d10a0565 --- /dev/null +++ b/tests/grids/boundaries/test_axes_boundaries_legacy.py @@ -0,0 +1,81 @@ +""" +.. codeauthor:: David Zwicker +""" + +import itertools + +import numpy as np +import pytest + +from pde import ScalarField, UnitGrid +from pde.grids.base import PeriodicityError +from pde.grids.boundaries.axes import ( + BCDataError, + BoundariesBase, + BoundariesList, + BoundariesSetter, +) +from pde.grids.boundaries.axis import BoundaryPair, BoundaryPeriodic, get_boundary_axis +from pde.grids.boundaries.local import NeumannBC + + +def test_boundaries_legacy(): + """Test setting boundaries for multiple systems.""" + b = ["periodic", "value", {"type": "derivative", "value": 1}] + for bx, by in itertools.product(b, b): + periodic = [b == "periodic" for b in (bx, by)] + g = UnitGrid([2, 2], periodic=periodic) + + bcs = BoundariesBase.from_data([bx, by], grid=g) + bc_x = get_boundary_axis(g, 0, bx) + bc_y = get_boundary_axis(g, 1, by) + + assert bcs.grid.num_axes == 2 + assert bcs.periodic == periodic + assert bcs[0] == bc_x + assert bcs[1] == bc_y + assert "field" in bcs.get_mathematical_representation("field") + assert isinstance(str(bcs), str) + assert isinstance(repr(bcs), str) + + assert bcs == BoundariesBase.from_data([bc_x, bc_y], grid=g) + if bx == by: + assert bcs == BoundariesBase.from_data(bx, grid=g) + + bc2 = bcs.copy() + assert bcs == bc2 + assert bcs is not bc2 + + b1 = BoundariesBase.from_data("auto_periodic_neumann", grid=UnitGrid([2, 2])) + b2 = BoundariesBase.from_data("auto_periodic_neumann", grid=UnitGrid([3, 3])) + assert isinstance(b1, BoundariesList) + assert isinstance(b2, BoundariesList) + assert b1 != b2 + + +def test_boundary_specifications_legacy(): + """Test different ways of specifying boundary conditions.""" + g = UnitGrid([2]) + bc1 = BoundariesBase.from_data( + [{"type": "derivative", "value": 0}, {"type": "value", "value": 0}], grid=g + ) + assert bc1 == BoundariesBase.from_data( + [{"type": "derivative"}, {"type": "value"}], grid=g + ) + assert bc1 == BoundariesBase.from_data([{"derivative": 0}, {"value": 0}], grid=g) + assert bc1 == BoundariesBase.from_data(["neumann", "dirichlet"], grid=g) + + +def test_boundaries_property_legacy(): + """Test boundaries property.""" + g = UnitGrid([2, 2]) + bc = BoundariesBase.from_data(["neumann", "dirichlet"], grid=g) + assert len(list(bc.boundaries)) == 4 + + bc = BoundariesBase.from_data("neumann", grid=g) + for b in bc.boundaries: + assert isinstance(b, NeumannBC) + + g = UnitGrid([2, 2], periodic=[True, False]) + bc = BoundariesBase.from_data("auto_periodic_neumann", grid=g) + assert len(list(bc.boundaries)) == 2 diff --git a/tests/grids/operators/test_cartesian_operators.py b/tests/grids/operators/test_cartesian_operators.py index f1315f29..f4775e0b 100644 --- a/tests/grids/operators/test_cartesian_operators.py +++ b/tests/grids/operators/test_cartesian_operators.py @@ -24,7 +24,7 @@ def _get_random_grid_bcs(ndim: int, dx="random", periodic="random", rank=0): - """Create a random Cartesian grid with natural bcs.""" + """Create a random Cartesian grid with auto_periodic_neumann bcs.""" rng = np.random.default_rng(0) shape = tuple(rng.integers(2, 5, ndim)) diff --git a/tests/grids/test_cartesian_grids.py b/tests/grids/test_cartesian_grids.py index 74e5e0c0..65f25635 100644 --- a/tests/grids/test_cartesian_grids.py +++ b/tests/grids/test_cartesian_grids.py @@ -9,7 +9,7 @@ import pytest from pde import CartesianGrid, UnitGrid -from pde.grids.boundaries import Boundaries, PeriodicityError +from pde.grids.boundaries import BoundariesBase, PeriodicityError from pde.grids.operators.common import make_derivative @@ -252,7 +252,7 @@ def test_setting_boundary_conditions(): grid.get_boundary_conditions("auto_periodic_neumann"), grid.get_boundary_conditions(["auto_periodic_neumann", "derivative"]), ]: - assert isinstance(bc, Boundaries) + assert isinstance(bc, BoundariesBase) for bc in ["periodic", "value"]: with pytest.raises(PeriodicityError): diff --git a/tests/test_integration.py b/tests/test_integration.py index f538d9ce..d357df45 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -10,6 +10,7 @@ import pytest from pde import CartesianGrid, DiffusionPDE, FileStorage, PDEBase, ScalarField, UnitGrid +from pde.grids.boundaries.axes import BoundariesSetter from pde.tools import misc, mpi, numba from pde.tools.misc import module_available @@ -317,3 +318,24 @@ def test_modelrunner_storage_many(tmp_path): # delete temporary files for path in tmp_path.iterdir(): path.unlink() + + +@pytest.mark.parametrize("backend", ["numpy", "numba"]) +def test_pde_with_bc_setter(backend): + """Test a PDE with boundary conditions controlled by an explicit setter.""" + + def setter(data, args=None): + data[0, :] = data[1, :] # Neumann + data[-1, :] = 2 * args["t"] - data[-2, :] # Dirichlet ~ t + data[:, 0] = data[:, -2] # periodic + data[:, -1] = data[:, 1] # periodic + + field = ScalarField.random_normal(UnitGrid([4, 4], periodic=[False, True])) + + eq1 = DiffusionPDE(bc=setter) + res1 = eq1.solve(field, t_range=1.01, dt=0.1, backend=backend) + + eq2 = DiffusionPDE(bc=[["neumann", {"value_expression": "t"}], "periodic"]) + res2 = eq2.solve(field, t_range=1.01, dt=0.1, backend=backend) + + np.testing.assert_allclose(res1.data, res2.data) diff --git a/tests/tools/test_expressions.py b/tests/tools/test_expressions.py index 462d6699..a81997e6 100644 --- a/tests/tools/test_expressions.py +++ b/tests/tools/test_expressions.py @@ -349,7 +349,7 @@ def test_evaluate_func_scalar(): field = ScalarField.from_expression(grid, "x") res1 = evaluate("laplace(a)", {"a": field}) - res2 = field.laplace("natural") + res2 = field.laplace("auto_periodic_neumann") np.testing.assert_almost_equal(res1.data, res2.data) res = evaluate("a - x", {"a": field}) @@ -402,7 +402,7 @@ def test_evaluate_func_invalid(): evaluate("a + b", {"a": field, "b": field2}) with pytest.raises(BCDataError): - evaluate("laplace(a)", {"a": field}, bc=RuntimeError) + evaluate("laplace(a)", {"a": field}, bc=Ellipsis) with pytest.raises(RuntimeError): evaluate("a + b", {"a": field})