Merge branch 'develop' into GSoC

Author: medha-14

CHANGELOG.md

@@ -1,6 +1,13 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
+
+- Generalises `set_initial_soc` to `set_initial_state` and adds Thevenin initial state setting. ([#5129](https://github.com/pybamm-team/PyBaMM/pull/5129))
+
+# [v25.8.0](https://github.com/pybamm-team/PyBaMM/tree/v25.8.0) - 2025-08-04
+
+## Features
+
 - Added `plot_3d_cross_section` & `plot_3d_heatmap` functions to support plotting for 3D thermal simulations. ([#5130](https://github.com/pybamm-team/PyBaMM/pull/5130))
 - Added a `Basic3DThermalSPM` with two way coupling. ([#5112](https://github.com/pybamm-team/PyBaMM/pull/5112))
 - Enables the passing of `inputs` throughout `set_initial_soc`. ([#5122](https://github.com/pybamm-team/PyBaMM/pull/5122))
@@ -18,6 +25,7 @@
 
 ## Bug fixes
 
+- Fixed non-deterministic plotting CI issues ([#5150](https://github.com/pybamm-team/PyBaMM/pull/5150))
 - Fix non-deterministic ShapeError in 3D FEM gradient method ([#5143](https://github.com/pybamm-team/PyBaMM/pull/5143))
 - Fixes negative electrode boundary values for half-cell voltage contributions. ([#5139](https://github.com/pybamm-team/PyBaMM/pull/5139))
 - Makes `A_cc` L_z * L_y * number of layers ([#5138](https://github.com/pybamm-team/PyBaMM/pull/5138))

CITATION.cff

@@ -24,6 +24,6 @@ keywords:
   - "expression tree"
   - "python"
   - "symbolic differentiation"
-version: "25.6.0"
+version: "25.8.0"
 repository-code: "https://github.com/pybamm-team/PyBaMM"
 title: "Python Battery Mathematical Modelling (PyBaMM)"

src/pybamm/__init__.py

@@ -39,6 +39,7 @@
 from .expression_tree.independent_variable import *
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
+from .expression_tree.vector_field import VectorField
 from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot
 
 from .expression_tree.exceptions import *
@@ -135,7 +136,10 @@
     SpectralVolume1DSubMesh,
     SymbolicUniform1DSubMesh,
 )
-
+from .meshes.two_dimensional_submeshes import (
+    SubMesh2D,
+    Uniform2DSubMesh,
+)
 from .meshes.scikit_fem_submeshes import (
     ScikitSubMesh2D,
     ScikitUniform2DSubMesh,
@@ -156,14 +160,15 @@
 from .spatial_methods.spatial_method import SpatialMethod
 from .spatial_methods.zero_dimensional_method import ZeroDimensionalSpatialMethod
 from .spatial_methods.finite_volume import FiniteVolume
+from .spatial_methods.finite_volume_2d import FiniteVolume2D
 from .spatial_methods.spectral_volume import SpectralVolume
 from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 from .spatial_methods.scikit_finite_element_3d import ScikitFiniteElement3D
 
 # Solver classes
 from .solvers.solution import Solution, EmptySolution, make_cycle_solution
 from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral
-from .solvers.processed_variable import ProcessedVariable, process_variable
+from .solvers.processed_variable import ProcessedVariable, ProcessedVariable2DFVM, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
 from .solvers.processed_variable import ProcessedVariableUnstructured
 from .solvers.summary_variable import SummaryVariables

src/pybamm/discretisations/discretisation.py

@@ -273,26 +273,59 @@ def set_variable_slices(self, variables):
         end = 0
         lower_bounds = []
         upper_bounds = []
+
         # Iterate through unpacked variables, adding appropriate slices to y_slices
         for variable in variables:
+            if variable in y_slices:
+                continue
             # Add up the size of all the domains in variable.domain
             if isinstance(variable, pybamm.ConcatenationVariable):
-                start_ = start
                 spatial_method = self.spatial_methods[variable.domain[0]]
+                dimension = spatial_method.mesh[variable.domain[0]].dimension
+                start_ = start
                 children = variable.children
                 meshes = OrderedDict()
+                lr_points = OrderedDict()
+                tb_points = OrderedDict()
                 for child in children:
                     meshes[child] = [spatial_method.mesh[dom] for dom in child.domain]
+                    if dimension == 2:
+                        lr_points[child] = sum(
+                            spatial_method.mesh[dom].npts_lr for dom in child.domain
+                        )
+                        tb_points[child] = sum(
+                            spatial_method.mesh[dom].npts_tb for dom in child.domain
+                        )
                 sec_points = spatial_method._get_auxiliary_domain_repeats(
                     variable.domains
                 )
                 for _ in range(sec_points):
+                    start_this_child = start_
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             end += domain_mesh.npts_for_broadcast_to_nodes
                         # Add to slices
-                        y_slices[child].append(slice(start_, end))
-                        y_slices_explicit[child].append(slice(start_, end))
+                        if dimension == 2:
+                            other_children = set(meshes.keys()) - {child}
+                            num_pts_to_skip = sum(
+                                lr_points[other_child] for other_child in other_children
+                            )
+                            for row in range(tb_points[child]):
+                                start_this_row = (
+                                    start_this_child
+                                    + (lr_points[child] + num_pts_to_skip) * row
+                                )
+                                end_this_row = start_this_row + lr_points[child]
+                                y_slices[child].append(
+                                    slice(start_this_row, end_this_row)
+                                )
+                                y_slices_explicit[child].append(
+                                    slice(start_this_row, end_this_row)
+                                )
+                            start_this_child += lr_points[child]
+                        else:
+                            y_slices[child].append(slice(start_, end))
+                            y_slices_explicit[child].append(slice(start_, end))
                         # Increment start_
                         start_ = end
             else:
@@ -737,6 +770,7 @@ def process_dict(self, var_eqn_dict, ics=False):
                     eqn = pybamm.FullBroadcast(eqn, broadcast_domains=eqn_key.domains)
 
             pybamm.logger.debug(f"Discretise {eqn_key!r}")
+
             processed_eqn = self.process_symbol(eqn)
             # Calculate scale if the key has a scale
             scale = getattr(eqn_key, "scale", 1)
@@ -811,9 +845,46 @@ def _process_symbol(self, symbol):
         if isinstance(symbol, pybamm.BinaryOperator):
             # Pre-process children
             left, right = symbol.children
+            # Catch case where diffusion is a scalar and turn it into an identity matrix vector field.
+            if len(symbol.domain) != 0:
+                spatial_method = self.spatial_methods[symbol.domain[0]]
+            else:
+                spatial_method = None
+            if isinstance(spatial_method, pybamm.FiniteVolume2D):
+                if isinstance(left, pybamm.Scalar) and (
+                    isinstance(right, pybamm.VectorField)
+                    or isinstance(right, pybamm.Gradient)
+                ):
+                    left = pybamm.VectorField(left, left)
+                elif isinstance(right, pybamm.Scalar) and (
+                    isinstance(left, pybamm.VectorField)
+                    or isinstance(left, pybamm.Gradient)
+                ):
+                    right = pybamm.VectorField(right, right)
             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 not isinstance(disc_right, pybamm.VectorField):
+                        disc_right = pybamm.VectorField(disc_right, disc_right)
+                    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]
+                        )
+                    )
+                    return pybamm.VectorField(disc_lr, disc_tb)
+
                 return pybamm.simplify_if_constant(
                     symbol.create_copy(new_children=[disc_left, disc_right])
                 )
@@ -878,7 +949,10 @@ def _process_symbol(self, symbol):
                     symbol.integration_variable[0].domain[0]
                 ]
                 out = integral_spatial_method.integral(
-                    child, disc_child, symbol._integration_dimension
+                    child,
+                    disc_child,
+                    symbol._integration_dimension,
+                    symbol.integration_variable,
                 )
                 out.copy_domains(symbol)
                 return out
@@ -888,6 +962,17 @@ def _process_symbol(self, symbol):
                     child, vector_type=symbol.vector_type
                 )
 
+            elif isinstance(symbol, pybamm.OneDimensionalIntegral):
+                child_spatial_method = self.spatial_methods[
+                    symbol.integration_domain[0]
+                ]
+                return child_spatial_method.one_dimensional_integral(
+                    symbol,
+                    child,
+                    disc_child,
+                    symbol.integration_domain,
+                    symbol.direction,
+                )
             elif isinstance(symbol, pybamm.BoundaryIntegral):
                 return child_spatial_method.boundary_integral(
                     child, disc_child, symbol.region
@@ -918,6 +1003,14 @@ def _process_symbol(self, symbol):
                 return child_spatial_method.evaluate_at(
                     symbol, disc_child, symbol.position
                 )
+            elif isinstance(symbol, pybamm.UpwindDownwind2D):
+                return spatial_method.upwind_or_downwind(
+                    child,
+                    disc_child,
+                    self.bcs,
+                    symbol.lr_direction,
+                    symbol.tb_direction,
+                )
             elif isinstance(symbol, pybamm.UpwindDownwind):
                 direction = symbol.name  # upwind or downwind
                 return spatial_method.upwind_or_downwind(
@@ -926,8 +1019,24 @@ 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.Magnitude):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Magnitude can only be applied to a vector field")
+                direction = symbol.direction
+                if direction == "lr":
+                    return disc_child.lr_field
+                elif direction == "tb":
+                    return disc_child.tb_field
+                else:
+                    raise ValueError("Invalid direction")
             else:
-                return symbol.create_copy(new_children=[disc_child])
+                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]),
+                    )
+                else:
+                    return symbol.create_copy(new_children=[disc_child])
 
         elif isinstance(symbol, pybamm.Function):
             disc_children = [self.process_symbol(child) for child in symbol.children]
@@ -996,6 +1105,12 @@ def _process_symbol(self, symbol):
         elif isinstance(symbol, pybamm.CoupledVariable):
             new_symbol = self.process_symbol(symbol.children[0])
             return new_symbol
+
+        elif isinstance(symbol, pybamm.VectorField):
+            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])
+
         elif isinstance(symbol, pybamm.Constant):
             # after discretisation we just care about the value, not the name
             return self.process_symbol(pybamm.Scalar(symbol.value))

src/pybamm/expression_tree/__init__.py

@@ -2,4 +2,4 @@
            'concatenations', 'exceptions', 'functions', 'independent_variable',
            'input_parameter', 'interpolant', 'matrix', 'operations',
            'parameter', 'printing', 'scalar', 'state_vector', 'symbol',
-           'unary_operators', 'variable', 'vector', 'discrete_time_sum']
+           'unary_operators', 'variable', 'vector', 'discrete_time_sum', 'vector_field', 'magnitude' ]

src/pybamm/expression_tree/concatenations.py

@@ -454,15 +454,15 @@ class SparseStack(Concatenation):
         The equations to concatenate
     """
 
-    def __init__(self, *children):
+    def __init__(self, *children, name="sparse_stack"):
         children = list(children)
         if not any(issparse(child.evaluate_for_shape()) for child in children):
             concatenation_function = np.vstack
         else:
             concatenation_function = vstack
         super().__init__(
             *children,
-            name="sparse_stack",
+            name=name,
             check_domain=False,
             concat_fun=concatenation_function,
         )

src/pybamm/expression_tree/independent_variable.py

@@ -144,12 +144,14 @@ def __init__(
         auxiliary_domains: AuxiliaryDomainType = None,
         domains: DomainsType = None,
         coord_sys=None,
+        direction=None,
     ) -> None:
         self.coord_sys = coord_sys
         super().__init__(
             name, domain=domain, auxiliary_domains=auxiliary_domains, domains=domains
         )
         domain = self.domain
+        self.direction = direction
 
         if domain == []:
             raise ValueError("domain must be provided")

src/pybamm/expression_tree/symbol.py

@@ -911,11 +911,13 @@ def evaluates_to_number(self):
     def evaluates_to_constant_number(self):
         return self.evaluates_to_number() and self.is_constant()
 
-    def evaluates_on_edges(self, dimension: str) -> bool:
+    def evaluates_on_edges(self, dimension: str) -> bool | str:
         """
         Returns True if a symbol evaluates on an edge, i.e. symbol contains a gradient
         operator, but not a divergence operator, and is not an IndefiniteIntegral.
-        Caches the solution for faster results.
+        Caches the solution for faster results. If the symbol is a component of a
+        vector field, returns a string indicating the direction that it evaluates on
+        edges.
 
         Parameters
         ----------