Unstructured finite volume

.github/workflows/_nox.yml

@@ -81,7 +81,7 @@ jobs:
         uses: awalsh128/cache-apt-pkgs-action@553a35bb8ebd9fcabcb1c9451aa4c98e1b4ca8a9 # v1.6.3
         if: startsWith(matrix.leg.os, 'ubuntu')
         with:
-          packages: gfortran gcc graphviz pandoc
+          packages: gfortran gcc graphviz pandoc libosmesa6 libegl1
           execute_install_scripts: true
 
       # dot -c is for registering graphviz fonts and plugins
@@ -92,6 +92,11 @@ jobs:
           sudo dot -c
           sudo apt-get install libopenblas-dev
 
+      # VTK off-screen GIF export needs a software OpenGL backend on headless runners.
+      - name: Prefer OSMesa for VTK on Linux
+        if: startsWith(matrix.leg.os, 'ubuntu')
+        run: echo "VTK_DEFAULT_OPENGL_WINDOW=vtkOSOpenGLRenderWindow" >> "$GITHUB_ENV"
+
       # Kept separate and opt-out: texlive-latex-extra is large and uncached.
       - name: Install TeXLive for Linux
         if: ${{ startsWith(matrix.leg.os, 'ubuntu') && inputs.texlive }}
@@ -115,6 +120,15 @@ jobs:
         if: startsWith(matrix.leg.os, 'windows')
         run: winget install --id Graphviz.Graphviz --exact --accept-source-agreements --accept-package-agreements
 
+      # VTK save_gif needs OSMesa on headless Windows runners (osmesa.dll on PATH).
+      - name: Setup headless OpenGL on Windows
+        if: startsWith(matrix.leg.os, 'windows')
+        uses: pyvista/setup-headless-display-action@5bc8de3bc71fcda7a96439571287a554901541a0 # v4
+        with:
+          pyvista: "false"
+          mesa3d-release: "24.3.0"
+          install-mesa3d-offscreen: "true"
+
       - name: Set up Python ${{ matrix.leg.python }}
         uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0
         with:

CHANGELOG.md

@@ -1,7 +1,14 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- Added unstructured finite volume support: `pybamm.UnstructuredSubMesh` (cell-centred meshes of triangles, quadrilaterals, tetrahedra, or hexahedra), the `pybamm.FiniteVolumeUnstructured` spatial method, and the `pybamm.lithium_ion.BasicDFN2DUnstructured`/`BasicDFN3DUnstructured` models. Meshes can be read from gmsh files via `pybamm.UserSuppliedUnstructuredMesh` or `pybamm.TaggedSubMeshGenerator`, and interfaces between adjacent submeshes are discovered automatically for arbitrary topologies rather than assuming a 1D stack. ([#5397](https://github.com/pybamm-team/PyBaMM/pull/5397))
+- Added `pybamm.VTKQuickPlot`, a VTK-based interactive alternative to `QuickPlot` for 2D and 3D unstructured mesh solutions. Requires the new `vtk` extra (`pip install pybamm[vtk]`). ([#5397](https://github.com/pybamm-team/PyBaMM/pull/5397))
+
 ## Bug fixes
 
+- `pybamm.max` and `pybamm.min` now clear their child's domains, reflecting that a reduction over a spatial field is a scalar. ([#5397](https://github.com/pybamm-team/PyBaMM/pull/5397))
+
 - `BatchStudy.solve` no longer ignores its `solver` argument: previously the loop over study inputs shadowed it, so a caller-supplied solver was silently dropped. A solver from `BatchStudy(solvers=...)` still takes precedence. ([#5677](https://github.com/pybamm-team/PyBaMM/pull/5677))
 - `pybamm.citations.register` now names the citation the caller passed in when a BibTeX string fails to parse, instead of whichever entry the parser had reached. ([#5677](https://github.com/pybamm-team/PyBaMM/pull/5677))
 - Deserialising a parameter set whose interpolant specification is invalid now logs a warning naming the offending parameter, instead of printing the bare exception to stdout with no indication of which parameter fell back to zero. ([#5679](https://github.com/pybamm-team/PyBaMM/pull/5679))

docs/source/api/meshes/index.rst

@@ -8,3 +8,4 @@ Meshes
   one_dimensional_submeshes
   two_dimensional_submeshes
   three_dimensional_submeshes
+  unstructured_submeshes

docs/source/api/meshes/unstructured_submeshes.rst

@@ -0,0 +1,16 @@
+Unstructured Sub Meshes
+=======================
+
+.. autoclass:: pybamm.UnstructuredSubMesh
+    :members:
+
+.. autoclass:: pybamm.UnstructuredMeshGenerator
+    :members:
+
+.. autoclass:: pybamm.UserSuppliedUnstructuredMesh
+    :members:
+
+.. autoclass:: pybamm.TaggedSubMeshGenerator
+    :members:
+
+.. autofunction:: pybamm.compute_interface_data

docs/source/api/models/lithium_ion/dfn.rst

@@ -13,4 +13,10 @@ Doyle-Fuller-Newman (DFN)
 .. autoclass:: pybamm.lithium_ion.BasicDFNHalfCell
     :members:
 
+.. autoclass:: pybamm.lithium_ion.BasicDFN2DUnstructured
+    :members:
+
+.. autoclass:: pybamm.lithium_ion.BasicDFN3DUnstructured
+    :members:
+
 .. footbibliography::

docs/source/api/plotting/index.rst

@@ -10,3 +10,4 @@ Plotting
   plot_summary_variables
   plot_3d_cross_section
   plot_3d_heatmap
+  plot_vtk

docs/source/api/plotting/plot_vtk.rst

@@ -0,0 +1,5 @@
+VTK Quick Plot
+==============
+
+.. autoclass:: pybamm.VTKQuickPlot
+    :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/pyproject.toml

@@ -66,11 +66,13 @@ bpx = ["bpx>=1.1.1,<1.2.0"]
 # Low-overhead progress bars
 tqdm = ["tqdm"]
 jax = ["jax>=0.7.0, <0.9.0; python_version >= '3.11' and (sys_platform != 'darwin' or platform_machine != 'x86_64')"]
+# VTK-based interactive visualization for unstructured meshes
+vtk = ["vtk>=9.0.0"]
 # Contains all optional dependencies, except for jax, and dev dependencies
 all = [
     "scikit-fem>=8.1.0",
     "meshio>=5.3.0",
-    "pybamm[examples,plot,cite,bpx,tqdm]",
+    "pybamm[examples,plot,cite,bpx,tqdm,vtk]",
 ]
 
 [dependency-groups]

packages/pybamm/src/pybamm/__init__.py

@@ -168,6 +168,14 @@
     UserSuppliedSubmesh3D,
 )
 
+from .meshes.unstructured_submesh import (
+    UnstructuredSubMesh,
+    UnstructuredMeshGenerator,
+    UserSuppliedUnstructuredMesh,
+    TaggedSubMeshGenerator,
+    compute_interface_data,
+)
+
 # Serialisation
 from .models.base_model import load_model
 
@@ -183,6 +191,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 (
@@ -193,7 +202,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
@@ -227,6 +236,7 @@
 from .plotting.dynamic_plot import dynamic_plot
 from .plotting.plot_3d_cross_section import plot_3d_cross_section
 from .plotting.plot_3d_heatmap import plot_3d_heatmap
+from .plotting.plot_vtk import VTKQuickPlot
 from .plotting.nyquist_plot import nyquist_plot
 
 # Simulation

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

@@ -492,6 +492,26 @@ 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).  If the
+            # spatial method on the first child's domain provides
+            # ``set_internal_bcs_for_concat``, defer to it and skip the
+            # default 1D-stack pairwise routine.
+            primary_method = self.spatial_methods.get(children[0].domain[0])
+            if primary_method is not None and hasattr(
+                primary_method, "set_internal_bcs_for_concat"
+            ):
+                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 not in bc_keys:
+                            internal_bcs[child] = child_bcs
+                    continue
+                # else fall through to legacy 1D-stack pairwise logic
+
             first_child = children[0]
             next_child = children[1]
 
@@ -581,8 +601,9 @@ 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.)
+            legacy_tab_sides = {"negative tab", "positive tab", "no tab"}
+            if legacy_tab_sides & set(bcs.keys()):
                 bcs = self.check_tab_conditions(key, bcs)
 
             # Process boundary conditions
@@ -882,7 +903,11 @@ def process_symbol(self, symbol):
 
         # Assign mesh as an attribute to the processed variable
         if symbol.domain != []:
-            discretised_symbol.mesh = self.mesh[symbol.domain]
+            mesh_for_symbol = self.mesh[symbol.domain]
+            discretised_symbol.mesh = mesh_for_symbol
+            if isinstance(discretised_symbol, pybamm.VectorField):
+                for comp in discretised_symbol._components:
+                    comp.mesh = mesh_for_symbol
         else:
             discretised_symbol.mesh = None
 
@@ -928,29 +953,47 @@ 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 == []:
                 if isinstance(disc_left, pybamm.VectorField) or isinstance(
                     disc_right, pybamm.VectorField
                 ):
+                    if isinstance(disc_left, pybamm.VectorField):
+                        n = disc_left.n_components
+                    else:
+                        n = disc_right.n_components
                     if not isinstance(disc_right, pybamm.VectorField):
-                        disc_right = pybamm.VectorField(disc_right, disc_right)
+                        disc_right = pybamm.VectorField(*[disc_right] * n)
                     if not isinstance(disc_left, pybamm.VectorField):
-                        disc_left = pybamm.VectorField(disc_left, disc_left)
-                    else:  # both are vector fields already
-                        pass
-                    disc_lr = pybamm.simplify_if_constant(
-                        symbol.create_copy(
-                            new_children=[disc_left.lr_field, disc_right.lr_field]
-                        )
-                    )
-                    disc_tb = pybamm.simplify_if_constant(
-                        symbol.create_copy(
-                            new_children=[disc_left.tb_field, disc_right.tb_field]
+                        disc_left = pybamm.VectorField(*[disc_left] * n)
+                    new_comps = [
+                        pybamm.simplify_if_constant(
+                            symbol.create_copy(
+                                new_children=[
+                                    disc_left._components[k],
+                                    disc_right._components[k],
+                                ]
+                            )
                         )
-                    )
-                    return pybamm.VectorField(disc_lr, disc_tb)
+                        for k in range(n)
+                    ]
+                    result = pybamm.VectorField(*new_comps)
+                    for src in (disc_left, disc_right):
+                        if hasattr(src, "_disc_state_vector"):
+                            result._disc_state_vector = src._disc_state_vector
+                            break
+                    return result
 
                 return pybamm.simplify_if_constant(
                     symbol.create_copy(new_children=[disc_left, disc_right])
@@ -978,6 +1021,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]]
@@ -1092,6 +1162,18 @@ def _process_symbol(self, symbol):
             elif isinstance(symbol, pybamm.NotConstant):
                 # After discretisation, we can make the symbol constant
                 return disc_child
+            elif isinstance(symbol, pybamm.Component):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Component can only be applied to a VectorField")
+                return disc_child._components[symbol.index]
+            elif isinstance(symbol, pybamm.Norm):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Norm can only be applied to a VectorField")
+                result = None
+                for comp in disc_child._components:
+                    sq = comp**2
+                    result = sq if result is None else result + sq
+                return result**0.5
             elif isinstance(symbol, pybamm.Magnitude):
                 if not isinstance(disc_child, pybamm.VectorField):
                     raise ValueError("Magnitude can only be applied to a vector field")
@@ -1104,10 +1186,14 @@ def _process_symbol(self, symbol):
                     raise ValueError("Invalid direction")
             else:
                 if isinstance(disc_child, pybamm.VectorField):
-                    return pybamm.VectorField(
-                        symbol.create_copy(new_children=[disc_child.lr_field]),
-                        symbol.create_copy(new_children=[disc_child.tb_field]),
-                    )
+                    new_comps = [
+                        symbol.create_copy(new_children=[c])
+                        for c in disc_child._components
+                    ]
+                    result = pybamm.VectorField(*new_comps)
+                    if hasattr(disc_child, "_disc_state_vector"):
+                        result._disc_state_vector = disc_child._disc_state_vector
+                    return result
                 else:
                     return symbol.create_copy(new_children=[disc_child])
 
@@ -1181,10 +1267,8 @@ def _process_symbol(self, symbol):
             )
 
         elif isinstance(symbol, pybamm.VectorField):
-            # VectorField is a subclass of TensorField, handle it first for specificity
-            left_symbol = self.process_symbol(symbol.lr_field)
-            right_symbol = self.process_symbol(symbol.tb_field)
-            return symbol.create_copy(new_children=[left_symbol, right_symbol])
+            processed = [self.process_symbol(c) for c in symbol._components]
+            return symbol.create_copy(new_children=processed)
 
         elif isinstance(symbol, pybamm.TensorField):
             # General TensorField handling (rank-2 tensors)