Add unstructured finite volume spatial method

CHANGELOG.md

@@ -8,6 +8,7 @@
 
 ## Features
 
+- Added `FiniteVolumeUnstructured` spatial method and unstructured processed-variable support for cell-centered data on arbitrary meshes. The TPFA Laplacian carries an implicit non-orthogonal correction (`"non-orthogonal correction"` option: `"over-relaxed"` or `"minimum"`) and gradients use a least-squares reconstruction, so both are exact on linear fields and second-order on skewed triangle and tetrahedral meshes. Diffusion coefficients reach faces through the distance-weighted harmonic mean, as in `FiniteVolume`, so material interfaces carry the exact series flux. ([#5688](https://github.com/pybamm-team/PyBaMM/pull/5688))
 - Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. Hexahedra must have planar faces (warped hexes raise a `GeometryError`), and `UserSuppliedUnstructuredMesh` accepts tetrahedral, triangular, and quadrilateral cells only. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687))
 - Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields. ([#5686](https://github.com/pybamm-team/PyBaMM/pull/5686))
 - Removed the left sidebar from the documentation home page for a cleaner landing experience. ([#5699](https://github.com/pybamm-team/PyBaMM/pull/5699))

docs/source/api/expression_tree/index.rst

@@ -12,6 +12,7 @@ Expression Tree
   matrix
   vector
   state_vector
+  vector_field
   binary_operator
   unary_operator
   concatenations

docs/source/api/expression_tree/vector_field.rst

@@ -0,0 +1,5 @@
+Vector Field
+============
+
+.. autoclass:: pybamm.VectorField
+  :members:

docs/source/api/spatial_methods/finite_volume_unstructured.rst

@@ -0,0 +1,5 @@
+Unstructured Finite Volume
+==========================
+
+.. autoclass:: pybamm.FiniteVolumeUnstructured
+    :members:

docs/source/api/spatial_methods/index.rst

@@ -10,3 +10,4 @@ Discretisation and spatial methods
   scikit_finite_element
   zero_dimensional_method
   scikit_finite_element_3d
+  finite_volume_unstructured

packages/pybamm/src/pybamm/__init__.py

@@ -192,6 +192,7 @@
 from .spatial_methods.spectral_volume import SpectralVolume
 from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 from .spatial_methods.scikit_finite_element_3d import ScikitFiniteElement3D
+from .spatial_methods.finite_volume_unstructured import FiniteVolumeUnstructured
 
 # Solver classes
 from .solvers.solution import (
@@ -202,7 +203,7 @@
     make_cycle_solution,
 )
 from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral
-from .solvers.processed_variable import ProcessedVariable, ProcessedVariable2DFVM, process_variable
+from .solvers.processed_variable import ProcessedVariable, ProcessedVariable2DFVM, ProcessedVariableUnstructuredFVM, ProcessedVariableVectorFieldUnstructuredFVM, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
 from .solvers.processed_variable import ProcessedVariableUnstructured
 from .solvers.summary_variable import SummaryVariables

packages/pybamm/src/pybamm/discretisations/discretisation.py

@@ -15,6 +15,11 @@ def has_bc_of_form(symbol, side, bcs, form):
     return (symbol in bcs) and (bcs[symbol][side][1] == form)
 
 
+# legacy current-collector tab BC side names, converted to left/right for
+# 1D meshes by Discretisation.check_tab_conditions
+LEGACY_TAB_SIDES = frozenset({"negative tab", "positive tab", "no tab"})
+
+
 class Discretisation:
     """The discretisation class, with methods to process a model and replace
     Spatial Operators with Matrices and Variables with StateVectors
@@ -492,9 +497,43 @@ def boundary_gradient(left_symbol, right_symbol):
                 continue
             children = var.orphans
 
+            # Dispatch hook: a spatial method may own its own internal-BC
+            # logic (e.g. graph-traversal for arbitrary topology); a non-None
+            # return replaces the default 1D-stack pairwise routine.
+            primary_method = self.spatial_methods.get(children[0].domain[0])
+            if primary_method is not None:
+                handled = primary_method.set_internal_bcs_for_concat(
+                    self, var, children, self.bcs[var]
+                )
+                if handled is not None:
+                    # Only adopt entries for children not already user-supplied.
+                    for child, child_bcs in handled.items():
+                        if child in bc_keys:
+                            continue
+                        if not child_bcs:
+                            # adopting an empty dict would strip the child of
+                            # BCs entirely; surface it instead
+                            pybamm.logger.warning(
+                                f"No internal or external boundary conditions "
+                                f"were found for {child.name!r} in domain "
+                                f"{child.domain}; it will be discretised "
+                                "without boundary conditions."
+                            )
+                            continue
+                        internal_bcs[child] = child_bcs
+                    continue
+                # else fall through to legacy 1D-stack pairwise logic
+
             first_child = children[0]
             next_child = children[1]
 
+            if "left" not in self.bcs[var] or "right" not in self.bcs[var]:
+                raise pybamm.DiscretisationError(
+                    f"Boundary conditions for the concatenated variable "
+                    f"{var.name!r} must include 'left' and 'right' entries "
+                    f"(got {sorted(self.bcs[var])}); other sides are not "
+                    "supported by the 1D-stack internal-BC routine."
+                )
             lbc = self.bcs[var]["left"]
             rbc = (boundary_gradient(first_child, next_child), "Neumann")
 
@@ -581,8 +620,8 @@ def process_boundary_conditions(self, model):
                         f"Neumann condition for {self.mesh[subdomain].coord_sys} coordinates"
                     )
 
-            # Handle any boundary conditions applied on the tabs
-            if any("tab" in side for side in list(bcs.keys())):
+            # Handle legacy tab boundary conditions ("negative tab", etc.)
+            if LEGACY_TAB_SIDES & set(bcs.keys()):
                 bcs = self.check_tab_conditions(key, bcs)
 
             # Process boundary conditions
@@ -944,7 +983,7 @@ def _process_symbol(self, symbol):
             # If boundary conditions are provided, need to check for BCs on tabs
             if self.bcs:
                 key_id = next(iter(self.bcs.keys()))
-                if any("tab" in side for side in list(self.bcs[key_id].keys())):
+                if LEGACY_TAB_SIDES & set(self.bcs[key_id].keys()):
                     self.bcs[key_id] = self.check_tab_conditions(
                         symbol, self.bcs[key_id]
                     )
@@ -964,6 +1003,16 @@ def _process_symbol(self, symbol):
                     isinstance(left, (pybamm.VectorField, pybamm.Gradient))
                 ):
                     right = pybamm.VectorField(right, right)
+            elif isinstance(spatial_method, pybamm.FiniteVolumeUnstructured):
+                dim = self.mesh[symbol.domain[0]].dimension
+                if isinstance(left, pybamm.Scalar) and isinstance(
+                    right, pybamm.VectorField | pybamm.Gradient
+                ):
+                    left = pybamm.VectorField(*[left] * dim)
+                elif isinstance(right, pybamm.Scalar) and isinstance(
+                    left, pybamm.VectorField | pybamm.Gradient
+                ):
+                    right = pybamm.VectorField(*[right] * dim)
             disc_left = self.process_symbol(left)
             disc_right = self.process_symbol(right)
             if symbol.domain == []:
@@ -1000,6 +1049,33 @@ def _process_symbol(self, symbol):
         elif isinstance(symbol, pybamm.UnaryOperator):
             child = symbol.child
 
+            # Intercept div(grad(u)) and div(D*grad(u)) before processing
+            # children, to avoid the expensive Green-Gauss gradient assembly.
+            if isinstance(symbol, pybamm.Divergence) and child.domain != []:
+                child_spatial_method = self.spatial_methods[child.domain[0]]
+                if isinstance(child_spatial_method, pybamm.FiniteVolumeUnstructured):
+                    grad_sym = None
+                    coeff_sym = None
+                    if isinstance(child, pybamm.Gradient):
+                        grad_sym = child
+                        coeff_sym = pybamm.Scalar(1)
+                    elif isinstance(child, pybamm.Multiplication):
+                        left_c, right_c = child.children
+                        if isinstance(right_c, pybamm.Gradient):
+                            grad_sym, coeff_sym = right_c, left_c
+                        elif isinstance(left_c, pybamm.Gradient):
+                            grad_sym, coeff_sym = left_c, right_c
+                    if grad_sym is not None:
+                        disc_coeff = self.process_symbol(coeff_sym)
+                        disc_u = self.process_symbol(grad_sym.child)
+                        return child_spatial_method.div_D_grad(
+                            symbol,
+                            grad_sym.child,
+                            disc_coeff,
+                            disc_u,
+                            self.bcs,
+                        )
+
             disc_child = self.process_symbol(child)
             if child.domain != []:
                 child_spatial_method = self.spatial_methods[child.domain[0]]
@@ -1119,6 +1195,11 @@ def _process_symbol(self, symbol):
                     raise pybamm.DiscretisationError(
                         "Component can only be applied to a VectorField"
                     )
+                if symbol.index >= disc_child.n_components:
+                    raise pybamm.DiscretisationError(
+                        f"Component index {symbol.index} is out of range for a "
+                        f"VectorField with {disc_child.n_components} components"
+                    )
                 return disc_child.components[symbol.index]
             elif isinstance(symbol, pybamm.Norm):
                 if not isinstance(disc_child, pybamm.VectorField):

packages/pybamm/src/pybamm/meshes/unstructured_submesh.py

@@ -551,6 +551,30 @@ def optimize_ordering(self):
                 if mirror.get("other_mesh") is self:
                     mirror["right_cells"] = permuted
 
+    def contains_points(self, query_pts):
+        """Boolean mask of ``query_pts`` lying inside the domain.
+
+        In 2D the boundary loops are cached on the mesh and containment uses
+        the even-odd rule, so holes are excluded, disconnected components
+        kept, and on-boundary points count as inside.  In 3D this is
+        :meth:`contains_points_3d`.  Returns ``None`` for a 2D mesh without
+        boundary edges.
+        """
+        query_pts = np.asarray(query_pts, dtype=np.float64)
+        if self.dimension == 3:
+            return self.contains_points_3d(query_pts)
+        if not hasattr(self, "_cached_boundary_loops"):
+            self._cached_boundary_loops = self.boundary_loops()
+        loops = self._cached_boundary_loops
+        if loops is None or len(loops) == 0:
+            return None
+        radius = 1e-9 * max(np.ptp(self.vertices, axis=0).max(), np.finfo(float).tiny)
+        containment_count = sum(
+            path.contains_points(query_pts[:, :2], radius=radius).astype(int)
+            for path in loops
+        )
+        return (containment_count % 2) == 1
+
     def boundary_loops(self):
         """Return boundary loops as a list of ``matplotlib.path.Path`` (2D only).
 
@@ -1159,8 +1183,8 @@ def compute_interface_data(left_mesh, right_mesh, left_name=None, right_name=Non
     Returns
     -------
     dict
-        Keys: ``"left_cells"``, ``"right_cells"``, ``"face_areas"``,
-        ``"cell_distances"``.
+        Keys: ``"left_cells"``, ``"right_cells"``, ``"left_faces"``,
+        ``"right_faces"``, ``"face_areas"``, ``"cell_distances"``.
     """
     left_bnd = left_mesh.boundary_faces.get("right", np.array([], dtype=int))
     right_bnd = right_mesh.boundary_faces.get("left", np.array([], dtype=int))
@@ -1211,9 +1235,12 @@ def compute_interface_data(left_mesh, right_mesh, left_name=None, right_name=Non
     right_cell_centroids = right_mesh.cell_centroids[right_cells]
     cell_distances = np.linalg.norm(right_cell_centroids - left_cell_centroids, axis=1)
 
+    right_faces = right_bnd[right_indices]
     result = {
         "left_cells": left_cells,
         "right_cells": right_cells,
+        "left_faces": left_bnd,
+        "right_faces": right_faces,
         "face_areas": face_areas,
         "cell_distances": cell_distances,
         "other_mesh": right_mesh,
@@ -1225,6 +1252,8 @@ def compute_interface_data(left_mesh, right_mesh, left_name=None, right_name=Non
         right_mesh.interface_data[left_name] = {
             "left_cells": right_cells,
             "right_cells": left_cells,
+            "left_faces": right_faces,
+            "right_faces": left_bnd,
             "face_areas": face_areas,
             "cell_distances": cell_distances,
             "other_mesh": left_mesh,

packages/pybamm/src/pybamm/parameters/parameter_substitutor.py

@@ -580,47 +580,26 @@ def process_boundary_conditions(
 
         Boundary conditions are dictionaries {"left": left bc, "right": right bc}
         in general, but may be imposed on the tabs (or *not* on the tab) for a
-        small number of variables.
+        small number of variables, or on arbitrary named boundary regions of an
+        unstructured mesh (e.g. Gmsh physical groups).
+
+        Every side present in the model is processed: a fixed list of side
+        names would silently discard boundary conditions on any other tag.
         """
         new_boundary_conditions: dict[
             pybamm.Symbol, dict[str, tuple[pybamm.Symbol, str]]
         ] = {}
-        sides = [
-            "left",
-            "right",
-            "negative tab",
-            "positive tab",
-            "no tab",
-            "top",
-            "bottom",
-            "x_min",
-            "x_max",
-            "y_min",
-            "y_max",
-            "z_min",
-            "z_max",
-            "r_min",
-            "r_max",
-        ]
         for variable, bcs in model.boundary_conditions.items():
             processed_variable = self.process_symbol(variable)
             new_boundary_conditions[processed_variable] = {}
-            for side in sides:
-                try:
-                    bc, typ = bcs[side]
-                    pybamm.logger.verbose(
-                        f"Processing parameters for {variable!r} ({side} bc)"
-                    )
-                    processed_bc = (self.process_symbol(bc), typ)
-                    new_boundary_conditions[processed_variable][side] = processed_bc
-                except KeyError as err:
-                    # don't raise error if the key error comes from the side not being
-                    # found
-                    if err.args[0] in side:
-                        pass
-                    # do raise error otherwise (e.g. can't process symbol)
-                    else:
-                        raise
+            for side, (bc, typ) in bcs.items():
+                pybamm.logger.verbose(
+                    f"Processing parameters for {variable!r} ({side} bc)"
+                )
+                new_boundary_conditions[processed_variable][side] = (
+                    self.process_symbol(bc),
+                    typ,
+                )
 
         return new_boundary_conditions
 

packages/pybamm/src/pybamm/plotting/quick_plot.py

@@ -317,6 +317,16 @@ def set_output_variables(self, output_variables, solutions):
             # just use the first solution to check this
             first_solution = variables[0]
             first_variable = first_solution[0]
+            if isinstance(
+                first_variable,
+                pybamm.ProcessedVariableUnstructuredFVM
+                | pybamm.ProcessedVariableVectorFieldUnstructuredFVM,
+            ):
+                raise NotImplementedError(
+                    f"QuickPlot cannot plot '{variable_tuple[0]}': variables on "
+                    "unstructured meshes have no plotting support yet. Query the "
+                    "variable at points with solution[name](t, x=..., z=...) instead."
+                )
             domain = first_variable.domain
             # check all other solutions against the first solution
             for idx, variable in enumerate(first_solution):