diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a158f7d3f..4a6154e5f4 100644 --- a/CHANGELOG.md +++ b/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)) diff --git a/docs/source/api/expression_tree/index.rst b/docs/source/api/expression_tree/index.rst index 0a6f3d757c..1a32e7a0f9 100644 --- a/docs/source/api/expression_tree/index.rst +++ b/docs/source/api/expression_tree/index.rst @@ -12,6 +12,7 @@ Expression Tree matrix vector state_vector + vector_field binary_operator unary_operator concatenations diff --git a/docs/source/api/expression_tree/vector_field.rst b/docs/source/api/expression_tree/vector_field.rst new file mode 100644 index 0000000000..2516c9c4a2 --- /dev/null +++ b/docs/source/api/expression_tree/vector_field.rst @@ -0,0 +1,5 @@ +Vector Field +============ + +.. autoclass:: pybamm.VectorField + :members: diff --git a/docs/source/api/spatial_methods/finite_volume_unstructured.rst b/docs/source/api/spatial_methods/finite_volume_unstructured.rst new file mode 100644 index 0000000000..eeb01b8a43 --- /dev/null +++ b/docs/source/api/spatial_methods/finite_volume_unstructured.rst @@ -0,0 +1,5 @@ +Unstructured Finite Volume +========================== + +.. autoclass:: pybamm.FiniteVolumeUnstructured + :members: diff --git a/docs/source/api/spatial_methods/index.rst b/docs/source/api/spatial_methods/index.rst index f9ccacd0d4..207c2a485d 100644 --- a/docs/source/api/spatial_methods/index.rst +++ b/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 diff --git a/packages/pybamm/src/pybamm/__init__.py b/packages/pybamm/src/pybamm/__init__.py index 8370ce6f58..f9fbc95875 100644 --- a/packages/pybamm/src/pybamm/__init__.py +++ b/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 diff --git a/packages/pybamm/src/pybamm/discretisations/discretisation.py b/packages/pybamm/src/pybamm/discretisations/discretisation.py index cdb404505d..781d4b62d6 100644 --- a/packages/pybamm/src/pybamm/discretisations/discretisation.py +++ b/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): diff --git a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py index f2d344fd44..816e8c4dfb 100644 --- a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py +++ b/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, diff --git a/packages/pybamm/src/pybamm/parameters/parameter_substitutor.py b/packages/pybamm/src/pybamm/parameters/parameter_substitutor.py index ab88e133f7..000b893828 100644 --- a/packages/pybamm/src/pybamm/parameters/parameter_substitutor.py +++ b/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 diff --git a/packages/pybamm/src/pybamm/plotting/quick_plot.py b/packages/pybamm/src/pybamm/plotting/quick_plot.py index 6b3fb3b4aa..dea94bf328 100644 --- a/packages/pybamm/src/pybamm/plotting/quick_plot.py +++ b/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): diff --git a/packages/pybamm/src/pybamm/solvers/processed_variable.py b/packages/pybamm/src/pybamm/solvers/processed_variable.py index 263becd5ed..345d9f2dd8 100644 --- a/packages/pybamm/src/pybamm/solvers/processed_variable.py +++ b/packages/pybamm/src/pybamm/solvers/processed_variable.py @@ -966,6 +966,265 @@ def _shape(self, t): return [self.first_dim_size, self.second_dim_size, len(t)] +class ProcessedVariableUnstructuredFVM(ProcessedVariable): + """ + Processed variable for cell-centered data on an unstructured mesh + (triangles, quads in 2D; tetrahedra in 3D). + + Spatial interpolation uses ``scipy.interpolate.LinearNDInterpolator`` + on cell centroids; query it at arbitrary points with ``pv(t, x=, y=, z=)``. + """ + + def __init__( + self, + name: str, + base_variables, + base_variables_casadi, + solution, + time_integral: pybamm.ProcessedVariableTimeIntegral | None = None, + ): + mesh = base_variables[0].mesh + if base_variables[0].size != mesh.npts: + raise NotImplementedError( + "Post-processing of unstructured-mesh variables with " + f"auxiliary domains is not yet supported: variable {name!r} " + f"has {base_variables[0].size} entries but the mesh has " + f"{mesh.npts} cells." + ) + if time_integral is not None: + # silently returning the integrand would be wrong; the postfix + # sum assumes time on axis 0, which only holds for 0D variables + raise NotImplementedError( + "Time integrals of unstructured-mesh variables are not yet supported." + ) + self.dimensions = 3 if mesh.dimension == 3 else 2 + super().__init__( + name, + base_variables, + base_variables_casadi, + solution, + time_integral=time_integral, + ) + self._time_interpolator = None + self.internal_boundaries = [] + + def _shape(self, t): + return [self.mesh.npts, len(t)] + + def initialise(self): + if self.entries_raw_initialized: + return + self._entries_raw = self.observe_raw() + + from scipy.interpolate import interp1d + + self._time_interpolator = interp1d( + self.t_pts, + self._entries_raw, + kind="linear", + axis=1, + bounds_error=False, + fill_value="extrapolate", + ) + + def _augmented_points(self): + """Return the interpolation point cloud (cell centroids + boundary + face centroids) and the index array for mapping cell values to + boundary face values. Cached after first call.""" + if not hasattr(self, "_aug_pts"): + mesh = self.mesh + bnd_start = mesh._boundary_face_start + bnd_centroids = mesh.face_centroids[bnd_start:] + self._aug_pts = np.concatenate([mesh.cell_centroids, bnd_centroids], axis=0) + self._aug_bnd_owners = mesh.face_owner[bnd_start:] + return self._aug_pts, self._aug_bnd_owners + + def _get_triangulation(self): + """Return a cached Delaunay triangulation of the augmented point cloud.""" + if not hasattr(self, "_cached_tri"): + from scipy.spatial import Delaunay + + pts, _ = self._augmented_points() + self._cached_tri = Delaunay(pts) + return self._cached_tri + + def _get_boundary_mask(self, query_pts): + """Boolean mask of query points outside the domain, or ``None`` when + the mesh cannot decide (2D mesh without boundary edges).""" + inside = self.mesh.contains_points(query_pts) + return None if inside is None else ~inside + + def _interpolate_spatial(self, values, query_pts, fill_value=np.nan): + """Interpolate cell-centered data to query points. + + ``values`` has one entry per cell and, optionally, one column per + time step; all columns are interpolated in a single pass, so the + nearest-neighbour tree and boundary mask are built once rather + than per time step. Points outside the domain are set to + ``fill_value``. NaNs in ``values`` propagate to the output. + """ + from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator + + pts, bnd_owners = self._augmented_points() + vals = np.concatenate([values, values[bnd_owners]]) + + tri = self._get_triangulation() + linear = LinearNDInterpolator(tri, vals) + result = linear(query_pts) + + # A query point outside the convex hull is NaN in every column + # (it is a location property), so whole rows are nearest-filled; + # requiring all columns avoids clobbering valid columns when the + # input itself contains NaNs. Points outside the domain get + # fill_value instead, so they are excluded before the (costly) + # nearest-neighbour fill rather than filled and then overwritten. + outside = self._get_boundary_mask(query_pts) + mask = np.isnan(result) + if result.ndim == 2: + mask = mask.all(axis=1) + if outside is not None: + mask &= ~outside + if np.any(mask): + nearest = NearestNDInterpolator(pts, vals) + result[mask] = nearest(query_pts[mask]) + if outside is not None: + result[outside] = fill_value + + return result + + def _data_at_time(self, t): + """Return cell-centered data at time t.""" + self.initialise() + t_observe, observe_raw = self._check_observe_raw(t) + if observe_raw: + return self._entries_raw + return self._time_interpolator(t_observe) + + def __call__( + self, t=None, x=None, r=None, y=None, z=None, R=None, fill_value=np.nan + ): + if r is not None or R is not None: + raise ValueError( + f"Variable {self._name!r} is on an unstructured mesh, which " + "has no r or R coordinates." + ) + if y is not None and self.mesh.dimension == 2: + raise ValueError( + f"Variable {self._name!r} is on a 2D unstructured mesh, which " + "has no y coordinate; its in-plane coordinates are x and z." + ) + data_at_t = self._data_at_time(t) + scalar_t = t is not None and np.ndim(t) == 0 + + spatial_provided = any(c is not None for c in [x, y, z]) + if not spatial_provided: + return data_at_t + + nodes = self.mesh.vertices + + def coord(values, axis): + if values is not None: + return np.asarray(values).ravel() + # a missing coordinate defaults to the domain midplane + return np.array([0.5 * (nodes[:, axis].min() + nodes[:, axis].max())]) + + if self.mesh.dimension == 2: + axes = [coord(x, 0), coord(z, 1)] + else: + axes = [coord(x, 0), coord(y, 1), coord(z, 2)] + grid = np.meshgrid(*axes, indexing="ij") + query = np.column_stack([g.ravel() for g in grid]) + out_shape = grid[0].shape + + if data_at_t.ndim == 1: + data_at_t = data_at_t[:, np.newaxis] + n_t = data_at_t.shape[1] + result = self._interpolate_spatial( + data_at_t, query, fill_value=fill_value + ).reshape(*out_shape, n_t) + + # scalar t drops the time axis; array-valued t (any length) keeps it + if scalar_t: + result = result[..., 0] + + return result + + +class ProcessedVariableVectorFieldUnstructuredFVM: + """ + Processed variable for a VectorField on an unstructured mesh. + + Wraps N scalar ``ProcessedVariableUnstructuredFVM`` instances (one per + component) and provides a unified interface for querying vector-valued + data. + """ + + def __init__( + self, + name: str, + base_variables, + base_variables_casadi, + solution, + time_integral=None, + ): + vf = base_variables[0] + + self.name = name + self.mesh = vf.mesh + self.domain = vf.domain + self.is_vector_field = True + self.n_components = vf.n_components + self.internal_boundaries = [] + + self._component_vars = [] + for k in range(vf.n_components): + comp_base_k = [bv.components[k] for bv in base_variables] + comp_casadi_k = [] + for bvc in base_variables_casadi: + if isinstance(bvc, list): + comp_casadi_k.append(bvc[k]) + elif isinstance(bvc, pybamm.VectorField): + comp_casadi_k.append(bvc.components[k]) + else: + comp_casadi_k.append(bvc) + pv = ProcessedVariableUnstructuredFVM( + f"{name}[{k}]", + comp_base_k, + comp_casadi_k, + solution, + time_integral=time_integral, + ) + self._component_vars.append(pv) + + self.dimensions = self._component_vars[0].dimensions + + @property + def entries(self): + """Tuple of per-component entry arrays.""" + return tuple(pv.entries for pv in self._component_vars) + + @property + def data(self): + """Tuple of per-component data arrays.""" + return tuple(pv.data for pv in self._component_vars) + + def update(self, other, new_sol): + raise NotImplementedError( + f"Variable {self.name!r}: vector-valued output_variables cannot " + "yet be merged across solution segments (multi-step experiments " + "or solution addition). Post-process the components separately." + ) + + def __call__( + self, t=None, x=None, r=None, y=None, z=None, R=None, fill_value=np.nan + ): + """Return a tuple of arrays, one per component.""" + return tuple( + pv(t=t, x=x, r=r, y=y, z=z, R=R, fill_value=fill_value) + for pv in self._component_vars + ) + + class ProcessedVariableRawFVM(ProcessedVariable): def _shape(self, t): return [self.base_variables[0].size, len(t)] @@ -1395,6 +1654,20 @@ def process_variable(name: str, base_variables, *args, **kwargs): ) return ProcessedVariable2DFVM(name, base_variables, *args, **kwargs) + if isinstance(mesh, pybamm.UnstructuredSubMesh): + if isinstance(base_variables[0], pybamm.VectorField): + return ProcessedVariableVectorFieldUnstructuredFVM( + name, base_variables, *args, **kwargs + ) + # Scalar reductions (e.g. Max/Min) keep the spatial domain but + # evaluate to a single value, so they are 0D in space. A one-cell + # mesh is also size 1, hence the cell-count check takes precedence. + if base_eval_size != mesh.npts and ( + len(base_eval_shape) == 0 or base_eval_shape[0] == 1 + ): + return ProcessedVariable0D(name, base_variables, *args, **kwargs) + return ProcessedVariableUnstructuredFVM(name, base_variables, *args, **kwargs) + # check variable shape if len(base_eval_shape) == 0 or base_eval_shape[0] == 1: return ProcessedVariable0D(name, base_variables, *args, **kwargs) diff --git a/packages/pybamm/src/pybamm/spatial_methods/__init__.py b/packages/pybamm/src/pybamm/spatial_methods/__init__.py index af2b1b4a23..bc57ff310f 100644 --- a/packages/pybamm/src/pybamm/spatial_methods/__init__.py +++ b/packages/pybamm/src/pybamm/spatial_methods/__init__.py @@ -1,2 +1,3 @@ __all__ = ['finite_volume', 'scikit_finite_element', 'spatial_method', - 'spectral_volume', 'zero_dimensional_method', 'scikit_finite_element_3d', 'finite_volume_2d'] + 'spectral_volume', 'zero_dimensional_method', 'scikit_finite_element_3d', 'finite_volume_2d', + 'finite_volume_unstructured'] diff --git a/packages/pybamm/src/pybamm/spatial_methods/finite_volume_unstructured.py b/packages/pybamm/src/pybamm/spatial_methods/finite_volume_unstructured.py new file mode 100644 index 0000000000..2c7f5179a7 --- /dev/null +++ b/packages/pybamm/src/pybamm/spatial_methods/finite_volume_unstructured.py @@ -0,0 +1,1833 @@ +""" +Finite Volume spatial method for unstructured simplex meshes (2D triangles / 3D tets). + +Dimension-agnostic: the same code path handles both 2D and 3D, with +dimension inferred from the mesh. All operators are assembled from +face-cell connectivity as sparse matrices. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +from scipy.sparse import coo_matrix, csr_matrix, diags, eye, kron +from scipy.spatial import cKDTree + +import pybamm + + +class FiniteVolumeUnstructured(pybamm.SpatialMethod): + """ + Cell-centered finite volume method on unstructured meshes. + + Supports triangles and quadrilaterals (2D), tetrahedra and hexahedra + (3D). Operators: + + * **Laplacian** – Two-Point Flux Approximation (TPFA) with an implicit + non-orthogonal correction: the face normal is split as + :math:`\\hat n = \\alpha \\hat e + \\mathbf{k}` along the unit + centroid-to-centroid direction :math:`\\hat e`, so the normal + derivative is :math:`\\alpha (u_j - u_i)/d + \\mathbf{k}\\cdot\\nabla + u_f` with the cross term taken from the Green-Gauss gradient. On + orthogonal meshes :math:`\\mathbf{k} = 0` and this is plain TPFA. + * **Gradient** – Green-Gauss cell-centroid reconstruction + * **Divergence** – face-flux summation (adjoint of gradient) + * **Boundary conditions** – ghost-cell (Dirichlet) / direct injection (Neumann) + + Neumann boundary values on the named axis sides (``"left"``/``"right"``, + ``"front"``/``"back"``, ``"bottom"``/``"top"``) are coordinate-direction + derivatives (:math:`\\partial u/\\partial x`, etc.), matching + :class:`pybamm.FiniteVolume`; e.g. ``u = x`` takes value ``+1`` on both + ``"left"`` and ``"right"``. Values on any other face tag (Gmsh region + names, ``"iface_*"``) are outward-normal derivatives + :math:`\\partial u/\\partial n`. + + Parameters + ---------- + options : dict, optional + Passed through to :class:`pybamm.SpatialMethod`. Additionally + ``"non-orthogonal correction"`` selects the decomposition of the + face normal: ``"over-relaxed"`` (default, :math:`\\alpha = 1/\\cos + \\theta`, favouring diagonal dominance) or ``"minimum"`` + (:math:`\\alpha = \\cos\\theta`, the smallest cross term). Both + are exact on linear fields. + """ + + _CORRECTIONS = ("over-relaxed", "minimum") + # Floor on cos(theta) in the over-relaxed weight (as in OpenFOAM): it + # bounds alpha, and k is built from the same alpha so consistency holds. + _COS_THETA_FLOOR = 0.05 + # Common CFD mesh-quality limit; beyond it the scheme stays consistent + # but conditioning degrades. + _NON_ORTHOGONALITY_WARNING_DEG = 70.0 + # Faces with |k| (about the angle in radians) below this are orthogonal: + # centroid rounding on high-aspect-ratio cells reaches 1e-11 and must not + # switch on the wide cross-term stencil. + _ORTHOGONALITY_TOL = 1e-8 + + def __init__(self, options=None): + super().__init__(options) + self.options.setdefault("non-orthogonal correction", "over-relaxed") + correction = self.options["non-orthogonal correction"] + if correction not in self._CORRECTIONS: + raise pybamm.OptionError( + "'non-orthogonal correction' must be one of " + f"{self._CORRECTIONS}, not {correction!r}" + ) + + # ------------------------------------------------------------------ + # build + # ------------------------------------------------------------------ + + def build(self, mesh): + """See :meth:`pybamm.SpatialMethod.build`.""" + from pybamm.meshes.unstructured_submesh import UnstructuredSubMesh + + super().build(mesh) + for dom in mesh: + mesh[dom].npts_for_broadcast_to_nodes = mesh[dom].npts + sm = mesh[dom] + if not isinstance(sm, UnstructuredSubMesh): + continue + name = dom[0] if isinstance(dom, tuple) else dom + max_angle = self._face_geometry(sm)["max_angle_deg"] + if max_angle > self._NON_ORTHOGONALITY_WARNING_DEG: + pybamm.logger.warning( + f"Unstructured submesh for domain {name!r} has faces with " + f"{max_angle:.1f} degrees of non-orthogonality (angle " + "between the face normal and the centroid line). The " + "discretisation remains consistent but the linear systems " + "become poorly conditioned; consider improving the mesh." + ) + # Tags come from the generator, not the constructor: a hand-built + # mesh with none gets no BCs and is invisible to interface + # discovery, so surface that before it fails downstream. + if not sm.boundary_faces and len(sm.face_owner) > sm._boundary_face_start: + pybamm.logger.warning( + f"Unstructured submesh for domain {name!r} has exterior " + "faces but no boundary tags: boundary conditions cannot " + "be applied and interface auto-discovery will not pair " + "it with neighboring domains. Tag it (e.g. " + "detect_box_boundaries() for axis-aligned boxes) or use " + "a mesh generator that supplies tags." + ) + # Discover interfaces between all unstructured submesh pairs so + # internal BCs work for arbitrary topology, not just 1D stacks. + self._auto_compute_all_interfaces(mesh) + + # ------------------------------------------------------------------ + # interface auto-discovery (graph topology support) + # ------------------------------------------------------------------ + + @staticmethod + def _interface_face_match(a_mesh, b_mesh, tol_factor=1e-3): + """Return matched boundary-face index pairs between two submeshes. + + Boundary faces whose centroids coincide within + :func:`pybamm.meshes.unstructured_submesh._geometric_tolerance` + (``tol_factor`` of the smallest element edge) are paired. Returns + ``(a_idx, b_idx, matched)`` where ``matched`` is True iff at least + one pair was found. + """ + a_idx = ( + np.concatenate(list(a_mesh.boundary_faces.values())) + if a_mesh.boundary_faces + else np.array([], dtype=int) + ) + b_idx = ( + np.concatenate(list(b_mesh.boundary_faces.values())) + if b_mesh.boundary_faces + else np.array([], dtype=int) + ) + if ( + len(a_idx) == 0 + or len(b_idx) == 0 + # meshes of different spatial dimension can never share an interface + or a_mesh.face_centroids.shape[1] != b_mesh.face_centroids.shape[1] + ): + return np.array([], dtype=int), np.array([], dtype=int), False + from pybamm.meshes.unstructured_submesh import _geometric_tolerance + + a_c = a_mesh.face_centroids[a_idx] + b_c = b_mesh.face_centroids[b_idx] + # main's mesh module owns the geometric tolerance definition + tol = _geometric_tolerance([a_mesh, b_mesh], rel=tol_factor) + tree = cKDTree(b_c) + d, j = tree.query(a_c, distance_upper_bound=tol) + keep = np.isfinite(d) + matched_b = j[keep] + if len(np.unique(matched_b)) != len(matched_b): + raise pybamm.GeometryError( + f"Interface between meshes is not one-to-one: multiple faces " + f"matched the same neighbor face within tolerance {tol:.2e}. " + "The meshes are non-conforming at the interface." + ) + return a_idx[keep], b_idx[matched_b], bool(keep.any()) + + def _compute_pair_interface(self, a_mesh, b_mesh, a_name, b_name): + """Populate ``interface_data`` and ``iface_`` face buckets for + a pair of submeshes that share a non-empty conformal interface. + + If either mesh already has an interface entry for the other (e.g. set + up by 1D-stack auto-pairing in :class:`pybamm.Mesh` or by a manual + ``compute_interface_data`` call), this method is a no-op so existing + models keep their original face-tag scheme. + """ + if b_name in a_mesh.interface_data or a_name in b_mesh.interface_data: + return False + a_match, b_match, ok = self._interface_face_match(a_mesh, b_mesh) + if not ok: + return False + + a_cells = a_mesh.face_owner[a_match] + b_cells = b_mesh.face_owner[b_match] + face_areas = a_mesh.face_areas[a_match] + cell_distances = np.linalg.norm( + b_mesh.cell_centroids[b_cells] - a_mesh.cell_centroids[a_cells], + axis=1, + ) + + a_mesh.interface_data[b_name] = { + "left_cells": a_cells, + "right_cells": b_cells, + "left_faces": a_match, + "right_faces": b_match, + "face_areas": face_areas, + "cell_distances": cell_distances, + "other_mesh": b_mesh, + } + b_mesh.interface_data[a_name] = { + "left_cells": b_cells, + "right_cells": a_cells, + "left_faces": b_match, + "right_faces": a_match, + "face_areas": face_areas, + "cell_distances": cell_distances, + "other_mesh": a_mesh, + } + + # Add new face-tag buckets for these interfaces. Order matches + # across both meshes so per-face BCs line up element-wise. + a_iface_tag = f"iface_{b_name}" + b_iface_tag = f"iface_{a_name}" + a_mesh.boundary_faces[a_iface_tag] = a_match + b_mesh.boundary_faces[b_iface_tag] = b_match + + # Remove interface faces from the axis-aligned buckets so external + # BCs don't double-count them. + a_match_set = {int(i) for i in a_match} + b_match_set = {int(i) for i in b_match} + for tag in list(a_mesh.boundary_faces.keys()): + if tag.startswith("iface_"): + continue + keep = np.array( + [int(i) not in a_match_set for i in a_mesh.boundary_faces[tag]], + dtype=bool, + ) + if keep.any(): + a_mesh.boundary_faces[tag] = a_mesh.boundary_faces[tag][keep] + else: + del a_mesh.boundary_faces[tag] + for tag in list(b_mesh.boundary_faces.keys()): + if tag.startswith("iface_"): + continue + keep = np.array( + [int(i) not in b_match_set for i in b_mesh.boundary_faces[tag]], + dtype=bool, + ) + if keep.any(): + b_mesh.boundary_faces[tag] = b_mesh.boundary_faces[tag][keep] + else: + del b_mesh.boundary_faces[tag] + return True + + def _auto_compute_all_interfaces(self, mesh): + """Walk every pair of unstructured submeshes; pair faces where they + coincide. Replaces the 1D-stack adjacency assumption with arbitrary + topology (star, tree, graph).""" + from pybamm.meshes.unstructured_submesh import UnstructuredSubMesh + + domains = [] + for raw in mesh: + name = raw[0] if isinstance(raw, tuple) else raw + sm = mesh[raw] + if isinstance(sm, UnstructuredSubMesh): + domains.append((name, sm)) + seen = set() + for (a, ma), (b, mb) in itertools.combinations(domains, 2): + if ma is mb or (a, b) in seen or (b, a) in seen: + continue + seen.add((a, b)) + # returns False when the pair shares no conformal interface + self._compute_pair_interface(ma, mb, a, b) + + # ------------------------------------------------------------------ + # internal BC assembly for arbitrary-topology Concatenation + # ------------------------------------------------------------------ + + def set_internal_bcs_for_concat(self, disc, var, children, outer_bcs): + """Build internal BC dict for each ``Concatenation`` child by walking + its mesh's ``interface_data`` graph instead of assuming consecutive + 1D-stack pairs. + + Returns ``None`` when no ``iface_`` face buckets exist in any + child mesh — that means there's no graph-discovered topology, so + the caller should fall through to the legacy 1D-stack pairwise + routine. + + For each child ``T_a`` on submesh ``mesh_a`` (graph case): + - Pass through any user-supplied ``outer_bcs`` whose tag matches an + external boundary tag present in ``mesh_a.boundary_faces``. + - For each interface ``mesh_a ↔ mesh_b`` (one entry per neighbor + in ``mesh_a.interface_data``), set ``iface_`` to the + discretised internal Neumann gradient between ``T_a`` and the + matching child ``T_b``. + """ + from pybamm.meshes.unstructured_submesh import UnstructuredSubMesh + + # Skip if no graph-discovered interfaces — caller falls back to the + # default 1D-stack pairwise logic. + has_iface = False + for c in children: + primary = c.domain[0] + sm = self.mesh[primary] + if isinstance(sm, UnstructuredSubMesh) and any( + k.startswith("iface_") for k in sm.boundary_faces + ): + has_iface = True + break + if not has_iface: + return None + + bcs_out = {} + name_to_child = {c.domain[0]: c for c in children} + for child in children: + primary = child.domain[0] + child_mesh = self.mesh[primary] + if not isinstance(child_mesh, UnstructuredSubMesh): + continue # leave default handling for non-unstructured children + bcs = {} + for tag, bc_value in outer_bcs.items(): + if tag in child_mesh.boundary_faces: + bcs[tag] = bc_value + for neighbor_name in child_mesh.interface_data: + neighbor_child = name_to_child.get(neighbor_name) + if neighbor_child is None: + continue + if f"iface_{neighbor_name}" not in child_mesh.boundary_faces: + pybamm.logger.warning( + f"Domain {primary!r} has interface data for " + f"{neighbor_name!r} but no 'iface_{neighbor_name}' " + "face bucket; skipping the internal BC, so these " + "domains will not be coupled." + ) + continue + left_disc = disc.process_symbol(child) + right_disc = disc.process_symbol(neighbor_child) + neighbor_mesh = self.mesh[neighbor_name] + # External conditions feed the interface gradient's cross + # term; the interface faces themselves enter as cross rows. + grad = self.internal_neumann_condition( + left_disc, + right_disc, + child_mesh, + neighbor_mesh, + left_bcs=self._external_bcs(child_mesh, outer_bcs), + right_bcs=self._external_bcs(neighbor_mesh, outer_bcs), + ) + bcs[f"iface_{neighbor_name}"] = (grad, "Neumann") + bcs_out[child] = bcs + return bcs_out + + @staticmethod + def _external_bcs(submesh, outer_bcs): + """The entries of ``outer_bcs`` on this mesh's exterior face tags.""" + return { + tag: bc + for tag, bc in outer_bcs.items() + if tag in submesh.boundary_faces and not tag.startswith("iface_") + } + + @staticmethod + def _is_scalar_value(symbol): + """Whether ``symbol`` is a single value to broadcast over faces. + + Time- or input-dependent scalars evaluate for shape to ``()`` or + ``(1,)`` rather than ``(1, 1)``, so every scalar shape counts. + """ + if isinstance(symbol, pybamm.Scalar): + return True + shape = getattr(symbol, "shape_for_testing", None) + return shape is not None and int(np.prod(shape)) == 1 + + @staticmethod + def _bc_contribution(n, n_bnd, owners, coeffs, bc_value, repeats=1): + """Build a symbolic BC contribution vector of size ``n * repeats``. + + For scalar ``bc_value``: returns ``Vector(accumulated_coeffs) * bc_value``. + For vector ``bc_value``: returns ``Matrix @ bc_value``, where the value + has one entry per boundary face (shared across auxiliary-domain + repeats) or ``n_bnd * repeats`` entries (one per face per repeat). + """ + is_scalar = FiniteVolumeUnstructured._is_scalar_value(bc_value) + if is_scalar: + row = np.zeros(n) + np.add.at(row, owners, coeffs) + if repeats > 1: + row = np.tile(row, repeats) + return pybamm.Vector(row) * bc_value + else: + M = csr_matrix((coeffs, (owners, np.arange(n_bnd))), shape=(n, n_bnd)) + if repeats > 1: + bc_shape = getattr(bc_value, "shape_for_testing", None) + if bc_shape == (n_bnd * repeats, 1): + M = csr_matrix(kron(eye(repeats, dtype=np.float64), M)) + else: + M = csr_matrix(kron(np.ones((repeats, 1)), M)) + return pybamm.Matrix(M) @ bc_value + + @staticmethod + def _tile_bc_value(bc_value, n_bnd, repeats): + """Lift a BC value to ``n_bnd * repeats`` entries. + + Scalars and already-full values (``n_bnd * repeats`` entries) pass + through; a per-face value (``n_bnd`` entries, shared across + auxiliary-domain repeats) is tiled, matching :meth:`_bc_contribution`. + """ + if repeats == 1: + return bc_value + is_scalar = FiniteVolumeUnstructured._is_scalar_value(bc_value) + if is_scalar or getattr(bc_value, "shape_for_testing", None) == ( + n_bnd * repeats, + 1, + ): + return bc_value + tile = csr_matrix(kron(np.ones((repeats, 1)), eye(n_bnd, dtype=np.float64))) + return pybamm.Matrix(tile) @ bc_value + + # ------------------------------------------------------------------ + # spatial_variable + # ------------------------------------------------------------------ + + def spatial_variable(self, symbol): + """Return a vector of cell-centroid coordinates for ``symbol``'s + direction (or its leading name token, e.g. ``x_n`` -> ``x``), tiled + over auxiliary domains. Raises :class:`pybamm.DomainError` rather + than guessing when neither identifies a coordinate.""" + symbol_mesh = self.mesh[symbol.domain] + repeats = self._get_auxiliary_domain_repeats(symbol.domains) + dim = symbol_mesh.dimension + + direction = getattr(symbol, "direction", None) + if direction is not None: + direction_cols = {"lr": 0, "fb": 1, "tb": dim - 1} + if direction not in direction_cols or (direction == "fb" and dim == 2): + valid = "'lr', 'tb'" if dim == 2 else "'lr', 'fb', 'tb'" + raise pybamm.DomainError( + f"Unknown direction {direction!r} for spatial variable " + f"{symbol.name!r} on a {dim}D unstructured mesh; valid " + f"directions are {valid}." + ) + col = direction_cols[direction] + else: + axis_name = symbol.name.split("_")[0] + name_cols = {"x": 0, "y": 1, "z": dim - 1} + if axis_name not in name_cols or (axis_name == "y" and dim == 2): + valid = "'x'/'z'" if dim == 2 else "'x'/'y'/'z'" + raise pybamm.DomainError( + f"Cannot infer a coordinate for spatial variable " + f"{symbol.name!r} on a {dim}D unstructured mesh; name it " + f"with a leading {valid} token (e.g. 'x_n') or set its " + "direction." + ) + col = name_cols[axis_name] + + entries = np.tile(symbol_mesh.cell_centroids[:, col], repeats) + return pybamm.Vector(entries, domains=symbol.domains) + + # ------------------------------------------------------------------ + # broadcast + # ------------------------------------------------------------------ + + def broadcast(self, symbol, domains, broadcast_type): + """See :meth:`pybamm.SpatialMethod.broadcast`.""" + domain = domains["primary"] + primary_pts = self.mesh[domain].npts + aux_repeats = self._get_auxiliary_domain_repeats(domains) + full_size = primary_pts * aux_repeats + + if broadcast_type.startswith("primary"): + sub_vector = np.ones((primary_pts, 1)) + if symbol.shape_for_testing == (): + out = symbol * pybamm.Vector(sub_vector) + else: + matrix = csr_matrix(kron(eye(symbol.shape_for_testing[0]), sub_vector)) + out = pybamm.Matrix(matrix) @ symbol + elif broadcast_type.startswith("full"): + out = symbol * pybamm.Vector(np.ones(full_size), domains=domains) + else: + from scipy.sparse import vstack + + # secondary/tertiary broadcast tiles the child by the size of the + # new (slower-varying) dimension, matching SpatialMethod.broadcast + if broadcast_type.startswith("secondary"): + reps = self._get_auxiliary_domain_repeats( + {"secondary": domains.get("secondary", [])} + ) + else: + reps = self._get_auxiliary_domain_repeats( + {"tertiary": domains.get("tertiary", [])} + ) + identity = eye(symbol.shape[0]) + matrix = vstack([identity for _ in range(reps)]) + out = pybamm.Matrix(matrix) @ symbol + + if out is symbol: + # simplification can hand back the child itself (e.g. ones-vector + # multiply); copy before stamping domains on a possibly shared node + out = symbol.create_copy(perform_simplifications=False) + out.domains = domains.copy() + return out + + # ================================================================== + # Core operators + # ================================================================== + + # ------------------------------------------------------------------ + # Laplacian (TPFA) + # ------------------------------------------------------------------ + + def laplacian(self, symbol, discretised_symbol, boundary_conditions): + """Laplacian ``Matrix @ discretised_symbol + bc_rhs``: the two-point + flux plus the non-orthogonal cross term, which is built from the + BC-aware Green-Gauss gradient so it stays fully implicit.""" + domain = symbol.domain + submesh = self.mesh[domain] + n = submesh.npts + d = submesh.dimension + repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + L = self._tpfa_matrix(submesh) + K = self._cross_term_matrices(submesh) + bcs = boundary_conditions.get(symbol, {}) + + # The gradient is only assembled if some face actually needs a cross + # term; on orthogonal meshes and boundaries the callable never fires. + gradient_cache = [] + + def gradient(): + if not gradient_cache: + gradient_cache.append( + self._least_squares_gradient(submesh, bcs, repeats) + ) + return gradient_cache[0] + + bc_rhs = pybamm.Vector(np.zeros(n * repeats)) + if bcs: + L, bc_rhs = self._apply_bcs_to_laplacian( + submesh, L, bc_rhs, bcs, repeats=repeats, gradient=gradient + ) + + if K is not None: + G_components, grad_bc_vecs, _ = gradient() + for k in range(d): + L = L + K[k] @ G_components[k] + if bcs: + K_full = csr_matrix(kron(eye(repeats, dtype=np.float64), K[k])) + bc_rhs = bc_rhs + pybamm.Matrix(K_full) @ grad_bc_vecs[k] + L = csr_matrix(L) + + L_full = csr_matrix(kron(eye(repeats, dtype=np.float64), L)) + result = pybamm.Matrix(L_full) @ discretised_symbol + bc_rhs + + return result + + @staticmethod + def _operator_cache(submesh): + """Per-submesh cache for assembled operator matrices. + + Keyed on the face-owner connectivity so a cell reordering (e.g. + ``optimize_ordering``) invalidates it. Cached matrices must never be + mutated in place. + """ + fingerprint = hash(submesh.face_owner.tobytes()) + cache = getattr(submesh, "_fv_operator_cache", None) + if cache is None or cache.get("fingerprint") != fingerprint: + cache = submesh._fv_operator_cache = {"fingerprint": fingerprint} + return cache + + def _face_geometry(self, submesh): + """Cached per-internal-face geometry shared by the TPFA operators. + + Returns a dict with the owner-to-neighbor centroid distance ``dist`` + and unit direction ``e_ij``, the signed ``cos_theta = n · e_ij``, + the distance-weighted owner interpolation weight ``w_owner`` for + face values, and the largest non-orthogonality angle in degrees. + + Raises + ------ + pybamm.GeometryError + If a face normal points away from the neighbor centroid: the + two-point flux is undefined on such (inverted or non-star-shaped) + cells. + """ + cache = self._operator_cache(submesh) + if "face_geometry" in cache: + return cache["face_geometry"] + n_int = submesh.n_internal_faces + owner = submesh.face_owner[:n_int] + neighbor = submesh.face_neighbor[:n_int] + centroids = submesh.cell_centroids + face_centroids = submesh.face_centroids[:n_int] + + delta = centroids[neighbor] - centroids[owner] + dist = np.linalg.norm(delta, axis=1) + e_ij = delta / dist[:, np.newaxis] + cos_theta = np.sum(submesh.face_normals[:n_int] * e_ij, axis=1) + if np.any(cos_theta <= 0): + raise pybamm.GeometryError( + f"{int(np.count_nonzero(cos_theta <= 0))} internal face(s) " + "have a normal pointing away from the neighbor centroid " + "(inverted or non-star-shaped cells), so the two-point flux " + "is undefined there. Fix the mesh." + ) + + d_owner = np.linalg.norm(face_centroids - centroids[owner], axis=1) + d_neighbor = np.linalg.norm(face_centroids - centroids[neighbor], axis=1) + w_owner = d_neighbor / (d_owner + d_neighbor) + + max_angle = np.degrees(np.arccos(np.min(cos_theta))) if n_int else 0.0 + cache["face_geometry"] = { + "dist": dist, + "e_ij": e_ij, + "cos_theta": cos_theta, + "w_owner": w_owner, + "max_angle_deg": float(max_angle), + } + return cache["face_geometry"] + + def _alpha(self, cos_theta): + """Implicit weight of the two-point difference in ``n = alpha e + k``. + + Any ``alpha`` is consistent because ``k`` is built from the same + value; the choice only sets how much flux the compact stencil + carries versus the reconstructed-gradient cross term. + """ + if self.options["non-orthogonal correction"] == "minimum": + return cos_theta + return 1.0 / np.maximum(cos_theta, self._COS_THETA_FLOOR) + + def _decomposition(self, submesh): + """``(alpha, k)`` per internal face for ``n = alpha e_ij + k``.""" + geometry = self._face_geometry(submesh) + alpha = self._alpha(geometry["cos_theta"]) + n_int = submesh.n_internal_faces + k = submesh.face_normals[:n_int] - alpha[:, np.newaxis] * geometry["e_ij"] + return alpha, self._drop_orthogonal(k) + + @classmethod + def _drop_orthogonal(cls, k): + """Zero ``k`` on faces that are orthogonal to within rounding, so + they neither enter nor widen the cross-term stencil.""" + k = k.copy() + k[np.linalg.norm(k, axis=1) < cls._ORTHOGONALITY_TOL] = 0.0 + return k + + def _boundary_decomposition(self, submesh, faces): + """``(dist, alpha, k)`` for boundary ``faces``, splitting the outward + normal along the unit vector from the owner centroid to the face + centroid: ``n = alpha e_b + k``. ``dist * cos(theta)`` is the + perpendicular distance, so ``alpha / dist`` is ``1 / (delta · n)`` + for the over-relaxed choice. + """ + delta = ( + submesh.face_centroids[faces] + - submesh.cell_centroids[submesh.face_owner[faces]] + ) + dist = np.linalg.norm(delta, axis=1) + e_b = delta / dist[:, np.newaxis] + normals = submesh.face_normals[faces] + alpha = self._alpha(np.sum(normals * e_b, axis=1)) + return dist, alpha, self._drop_orthogonal(normals - alpha[:, np.newaxis] * e_b) + + def _cross_term_matrices(self, submesh): + """Assemble (or fetch the cached) matrices ``K_k`` mapping the cell + gradient components to the cell divergence of the internal-face + cross fluxes ``A_f k_f · grad(u)_f``, where ``grad(u)_f`` is the + distance-weighted interpolation of the two cell gradients. + + Returns ``None`` when every internal face is orthogonal (``k = 0``), + so orthogonal meshes pay nothing for the correction. + """ + cache = self._operator_cache(submesh) + key = ("cross", self.options["non-orthogonal correction"]) + if key in cache: + return cache[key] + n = submesh.npts + n_int = submesh.n_internal_faces + d = submesh.dimension + _, k = self._decomposition(submesh) + if n_int == 0 or not k.any(): + cache[key] = None + return None + + owner = submesh.face_owner[:n_int] + neighbor = submesh.face_neighbor[:n_int] + areas = submesh.face_areas[:n_int] + vol = submesh.cell_volumes + w_owner = self._face_geometry(submesh)["w_owner"] + + face_rows = np.tile(np.arange(n_int), 2) + both = np.concatenate([owner, neighbor]) + # P: cell gradient -> face gradient; S: face flux -> cell divergence + # (+owner, -neighbor, so the cross flux is conservative by construction) + P = csr_matrix( + (np.concatenate([w_owner, 1.0 - w_owner]), (face_rows, both)), + shape=(n_int, n), + ) + S = csr_matrix( + ( + np.concatenate([1.0 / vol[owner], -1.0 / vol[neighbor]]), + (both, face_rows), + ), + shape=(n, n_int), + ) + matrices = [] + for kk in range(d): + matrix = csr_matrix(S @ diags(areas * k[:, kk]) @ P) + matrix.eliminate_zeros() # orthogonal faces must not widen the stencil + matrices.append(matrix) + cache[key] = matrices + return cache[key] + + def _tpfa_matrix(self, submesh): + """Assemble (or fetch the cached) two-point part of the Laplacian for + internal faces only: the ``alpha (u_j - u_i) / d`` term of the + decomposition ``n = alpha e_ij + k``. :meth:`_cross_term_matrices` + supplies the ``k · grad(u)_f`` remainder; on orthogonal meshes + ``alpha = 1`` and this is the whole operator. + """ + cache = self._operator_cache(submesh) + key = ("tpfa", self.options["non-orthogonal correction"]) + if key in cache: + return cache[key] + n = submesh.npts + n_int = submesh.n_internal_faces + + owner = submesh.face_owner[:n_int] + neighbor = submesh.face_neighbor[:n_int] + alpha, _ = self._decomposition(submesh) + coeff = ( + submesh.face_areas[:n_int] * alpha / self._face_geometry(submesh)["dist"] + ) + + vol = submesh.cell_volumes + + rows = np.concatenate([owner, neighbor, owner, neighbor]) + cols = np.concatenate([neighbor, owner, owner, neighbor]) + data = np.concatenate( + [ + coeff / vol[owner], + coeff / vol[neighbor], + -coeff / vol[owner], + -coeff / vol[neighbor], + ] + ) + + cache[key] = csr_matrix(coo_matrix((data, (rows, cols)), shape=(n, n))) + return cache[key] + + def _div_D_grad_matrices(self, submesh): + """Assemble (or fetch the cached) matrices for :meth:`div_D_grad`: + ``G`` (two-point difference per internal face), ``W`` (linear + interpolation to faces), ``S`` (face flux to cell divergence), the + geometric factor ``geo`` per internal face, the cross-term matrices + ``C`` (``None`` on orthogonal meshes) and ``W_harmonic`` (resistance + weights for the harmonic mean of ``D``). + """ + cache = self._operator_cache(submesh) + key = ("div_D_grad", self.options["non-orthogonal correction"]) + if key in cache: + return cache[key] + + n = submesh.npts + n_int = submesh.n_internal_faces + vol = submesh.cell_volumes + owner = submesh.face_owner[:n_int] + neighbor = submesh.face_neighbor[:n_int] + + alpha, _ = self._decomposition(submesh) + geo = submesh.face_areas[:n_int] * alpha / self._face_geometry(submesh)["dist"] + + # G (n_int x n): u_neighbor - u_owner per face + G = csr_matrix( + ( + np.concatenate([-np.ones(n_int), np.ones(n_int)]), + (np.tile(np.arange(n_int), 2), np.concatenate([owner, neighbor])), + ), + shape=(n_int, n), + ) + + # W (n_int x n): linear (distance-weighted) interpolation to faces, + # used for the face gradient of the cross term + w_owner = self._face_geometry(submesh)["w_owner"] + face_rows = np.tile(np.arange(n_int), 2) + both = np.concatenate([owner, neighbor]) + W = csr_matrix( + (np.concatenate([w_owner, 1.0 - w_owner]), (face_rows, both)), + shape=(n_int, n), + ) + # W_h (n_int x n): resistance weights for the harmonic mean of D, + # D_f = 1 / (W_h @ (1/D)); each cell weighs by its own centroid-to-face + # distance, so a face between two materials carries the series flux + W_harmonic = csr_matrix( + (np.concatenate([1.0 - w_owner, w_owner]), (face_rows, both)), + shape=(n_int, n), + ) + + # S (n x n_int): face flux -> cell divergence (+owner, -neighbor, /V) + S = csr_matrix( + ( + np.concatenate([1.0 / vol[owner], -1.0 / vol[neighbor]]), + (np.concatenate([owner, neighbor]), np.tile(np.arange(n_int), 2)), + ), + shape=(n, n_int), + ) + + # C[k] (n_int x n): cell gradient component -> face cross flux + # A_f k_f,k grad_k(u)_f, interpolated like D; None when orthogonal + _, k_vec = self._decomposition(submesh) + if not k_vec.any(): + C = None + else: + areas = submesh.face_areas[:n_int] + C = [] + for kk in range(submesh.dimension): + matrix = csr_matrix(diags(areas * k_vec[:, kk]) @ W) + matrix.eliminate_zeros() + C.append(matrix) + + cache[key] = (G, W, S, geo, C, W_harmonic) + return cache[key] + + def div_D_grad(self, div_symbol, grad_child, disc_D, disc_u, boundary_conditions): + """Discretise ``div(D * grad(u))`` as a single TPFA operation. + + Fully symbolic — works for both constant and state-dependent scalar + ``D``. Internal-face fluxes use the distance-weighted harmonic mean of + ``D`` (resistances in series, as :class:`pybamm.FiniteVolume` does for + coefficients of a gradient, so material interfaces carry the exact + two-cell flux) and the two-point normal derivative plus its + non-orthogonal cross term (see :meth:`_tpfa_matrix`). ``D`` must be + strictly positive. + + This method is only reached when the expression is written as + ``div(D * grad(u))`` (a single product, matched syntactically during + discretisation); other flux forms go through the generic + :meth:`gradient`/:meth:`divergence` operators, which cannot apply + boundary conditions conservatively and raise instead. + """ + if isinstance(disc_D, pybamm.VectorField): + raise pybamm.DiscretisationError( + "Anisotropic (vector-valued) diffusion coefficients are not " + "supported by the TPFA discretisation of div(D * grad(u))." + ) + domain = div_symbol.domain + submesh = self.mesh[domain] + n = submesh.npts + repeats = self._get_auxiliary_domain_repeats(div_symbol.domains) + vol = submesh.cell_volumes + + G, _, S, geo, C, W_harmonic = self._div_D_grad_matrices(submesh) + bcs = boundary_conditions.get(grad_child, {}) + + def lift(matrix): + if repeats == 1: + return matrix + return csr_matrix(kron(eye(repeats, dtype=np.float64), matrix)) + + def tile(values): + return np.tile(values, repeats) if repeats > 1 else values + + # Cell gradient components, assembled lazily: only non-orthogonal + # faces (internal or Dirichlet) need them for their cross term. + gradient_cache = [] + + def gradient(): + if not gradient_cache: + G_grad, grad_bc, _ = self._least_squares_gradient(submesh, bcs, repeats) + gradient_cache.append( + [ + pybamm.Matrix(lift(G_grad[k])) @ disc_u + grad_bc[k] + for k in range(submesh.dimension) + ] + ) + return gradient_cache[0] + + normal_grad = pybamm.Matrix(lift(G)) @ disc_u * pybamm.Vector(tile(geo)) + if C is not None: + for k, grad_k in enumerate(gradient()): + normal_grad = normal_grad + pybamm.Matrix(lift(C[k])) @ grad_k + is_scalar_D = self._is_scalar_value(disc_D) + if is_scalar_D: + D_face = disc_D + else: + if isinstance(disc_D, pybamm.Vector) and np.any(disc_D.entries <= 0): + raise pybamm.DiscretisationError( + "div(D * grad(u)) needs a strictly positive coefficient D: " + "faces take its harmonic mean." + ) + D_face = 1 / (pybamm.Matrix(lift(W_harmonic)) @ (1 / disc_D)) + result = pybamm.Matrix(lift(S)) @ (D_face * normal_grad) + + # Boundary conditions + bc_rhs = pybamm.Vector(np.zeros(n * repeats)) + if bcs: + for side, (bc_value, bc_type) in bcs.items(): + self._check_bc_type(bc_type) + fi_arr = self._boundary_faces_for_side(submesh, side) + n_bnd = len(fi_arr) + bnd_own = submesh.face_owner[fi_arr] + + E = csr_matrix( + (np.ones(n_bnd), (np.arange(n_bnd), bnd_own)), + shape=(n_bnd, n), + ) + P = csr_matrix( + (np.ones(n_bnd), (bnd_own, np.arange(n_bnd))), + shape=(n, n_bnd), + ) + E_f, P_f = lift(E), lift(P) + D_bnd = disc_D if is_scalar_D else pybamm.Matrix(E_f) @ disc_D + bc_value = self._tile_bc_value(bc_value, n_bnd, repeats) + a_over_v = submesh.face_areas[fi_arr] / vol[bnd_own] + + if bc_type == "Dirichlet": + dist, alpha, k_vec = self._boundary_decomposition(submesh, fi_arr) + u_bnd = pybamm.Matrix(E_f) @ disc_u + normal_grad_bnd = (bc_value - u_bnd) * pybamm.Vector( + tile(a_over_v * alpha / dist) + ) + if k_vec.any(): + for k, grad_k in enumerate(gradient()): + normal_grad_bnd = normal_grad_bnd + ( + pybamm.Matrix(E_f) @ grad_k + ) * pybamm.Vector(tile(a_over_v * k_vec[:, k])) + bc_rhs = bc_rhs + pybamm.Matrix(P_f) @ (D_bnd * normal_grad_bnd) + + elif bc_type == "Neumann" and bc_value != pybamm.Scalar(0): + bc_rhs = bc_rhs + pybamm.Matrix(P_f) @ ( + D_bnd + * bc_value + * pybamm.Vector(tile(self._neumann_sign(side) * a_over_v)) + ) + + return result + bc_rhs + + def _apply_bcs_to_laplacian( + self, submesh, L, bc_rhs, bcs, repeats=1, gradient=None + ): + """Return the Laplacian matrix and RHS modified for boundary + conditions. + + ``bc_rhs`` is a pybamm expression (symbolic vector of size + ``npts * repeats``). ``L`` is not mutated (it may be cached). + ``gradient`` is a zero-argument callable returning the + ``(matrices, bc_vecs)`` of the cell gradient; it is only called for + Dirichlet faces whose centroid direction is not normal to the face, + which need the cross term ``A k · grad(u)``. Without it those + faces get the two-point term only. + """ + n = submesh.npts + d = submesh.dimension + diag_correction = np.zeros(n) + cross_diag = np.zeros((d, n)) + + for side, (bc_value, bc_type) in bcs.items(): + self._check_bc_type(bc_type) + face_indices = self._boundary_faces_for_side(submesh, side) + n_bnd = len(face_indices) + owners = submesh.face_owner[face_indices] + a_over_v = submesh.face_areas[face_indices] / submesh.cell_volumes[owners] + + if bc_type == "Dirichlet": + dist, alpha, k_vec = self._boundary_decomposition(submesh, face_indices) + coeffs = a_over_v * alpha / dist + np.add.at(diag_correction, owners, -coeffs) + bc_rhs = bc_rhs + self._bc_contribution( + n, n_bnd, owners, coeffs, bc_value, repeats=repeats + ) + if gradient is not None: + for k in range(d): + np.add.at(cross_diag[k], owners, a_over_v * k_vec[:, k]) + + elif bc_type == "Neumann": + coeffs = self._neumann_sign(side) * a_over_v + bc_rhs = bc_rhs + self._bc_contribution( + n, n_bnd, owners, coeffs, bc_value, repeats=repeats + ) + + if np.any(diag_correction): + L = csr_matrix(L + diags(diag_correction)) + if cross_diag.any(): + G_components, grad_bc_vecs, _ = gradient() + for k in range(d): + scale = diags(cross_diag[k]) + L = L + scale @ G_components[k] + scale_full = csr_matrix(kron(eye(repeats, dtype=np.float64), scale)) + bc_rhs = bc_rhs + pybamm.Matrix(scale_full) @ grad_bc_vecs[k] + L = csr_matrix(L) + return L, bc_rhs + + @staticmethod + def _boundary_faces_for_side(submesh, side): + """Boundary-face indices for a BC side, raising when the tag is unknown. + + BC sides map directly onto ``submesh.boundary_faces`` keys; a missing + key means the BC cannot be applied, so failing loudly here is what + stops typos and interface-consumed sides from silently dropping BCs. + """ + if side not in submesh.boundary_faces: + raise pybamm.DiscretisationError( + f"No boundary faces tagged {side!r} on this mesh (available " + f"tags: {sorted(submesh.boundary_faces)}). The side may be " + "misspelled, or its faces were absorbed into an internal " + "interface by interface discovery." + ) + return submesh.boundary_faces[side] + + @staticmethod + def _check_bc_type(bc_type): + if bc_type not in ("Dirichlet", "Neumann"): + raise pybamm.DiscretisationError( + f"boundary condition must be Dirichlet or Neumann, not {bc_type!r}" + ) + + # Named sides whose outward normal points along the negative coordinate + # axis (see UnstructuredSubMesh._identify_boundary_faces). + _NEGATIVE_NORMAL_SIDES = frozenset({"left", "front", "bottom"}) + + @classmethod + def _neumann_sign(cls, side): + """Sign converting a Neumann boundary value to an outward-normal + derivative. + + Named axis sides carry PyBaMM coordinate-direction values, so sides + with a negative outward normal flip sign; any other face tag (Gmsh + region names, ``iface_*``) is already outward-normal. + """ + return -1.0 if side in cls._NEGATIVE_NORMAL_SIDES else 1.0 + + # ------------------------------------------------------------------ + # Gradient (Green-Gauss) + # ------------------------------------------------------------------ + + def gradient(self, symbol, discretised_symbol, boundary_conditions): + """Least-squares cell-centroid gradient, returned as a + :class:`pybamm.VectorField` with one component per dimension. + + Exact on linear fields for any cell shape: every face contributes + one directional-derivative equation — towards the neighbour + centroid (internal faces), towards the face centroid holding the + prescribed value (Dirichlet), or the normal derivative itself + (Neumann). Boundary faces without a condition contribute nothing. + """ + domain = symbol.domain + submesh = self.mesh[domain] + d = submesh.dimension + repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + bcs = boundary_conditions.get(symbol, {}) + if bcs: + missing = [tag for tag in submesh.boundary_faces if tag not in bcs] + if missing: + pybamm.logger.warning( + f"Gradient of {symbol.name!r}: boundary face buckets " + f"{missing} have no boundary condition and are fitted as " + "zero normal derivative (the operators' zero-flux " + "treatment), which is wrong if the field varies normal " + "to those boundaries." + ) + G_components, bc_vecs, _ = self._least_squares_gradient(submesh, bcs, repeats) + + components = [] + for k in range(d): + Gk = csr_matrix(kron(eye(repeats, dtype=np.float64), G_components[k])) + comp = pybamm.Matrix(Gk) @ discretised_symbol + bc_vecs[k] + components.append(comp) + + return pybamm.VectorField(*components) + + # Row kinds of the least-squares gradient fit + _ROW_INTERNAL, _ROW_DIRICHLET, _ROW_NEUMANN, _ROW_NO_BC, _ROW_INTERFACE = range(5) + + def _face_bc_kinds(self, submesh, bcs, interface_faces=None): + """Per-face row kind (see the ``_ROW_*`` constants).""" + kinds = np.full(len(submesh.face_owner), self._ROW_NO_BC, dtype=int) + kinds[: submesh.n_internal_faces] = self._ROW_INTERNAL + for side, (_, bc_type) in bcs.items(): + self._check_bc_type(bc_type) + faces = self._boundary_faces_for_side(submesh, side) + kinds[faces] = ( + self._ROW_DIRICHLET if bc_type == "Dirichlet" else self._ROW_NEUMANN + ) + if interface_faces is not None: + kinds[interface_faces] = self._ROW_INTERFACE + return kinds + + def _least_squares_matrices(self, submesh, bcs, interface=None): + """Cached matrix part of the least-squares gradient for one BC layout. + + Each cell has the same number of faces ``m``, so the per-cell normal + equations are solved in a batch. Rows are unit-direction equations + ``e · grad(u) = b``: ``e`` towards the neighbour centroid with + ``b = (u_j - u_i) / dist`` (internal faces, and interface faces + towards the other mesh's cell), towards the face centroid with + ``b = (u_b - u_i) / dist`` (Dirichlet), or the outward normal with + ``b`` the prescribed derivative (Neumann). Boundary faces without a + condition take ``b = 0``, matching the operators' zero-flux treatment + of such faces. Cells whose directions do not span the space get the + minimum-norm fit via the pseudo-inverse. + + Parameters + ---------- + submesh : UnstructuredSubMesh + bcs : dict + ``{side: (value, type)}`` boundary conditions. + interface : dict, optional + Cross-mesh rows: ``faces`` (this mesh's interface faces), + ``other_cells`` and ``other_centroids`` (the paired cells of the + other mesh), ``n_other`` and a hashable ``key`` for caching. + + Returns + ------- + tuple + ``(G, coeff, slot, length, G_cross)``: ``G[k]`` maps cell values + to gradient component ``k`` and ``G_cross[k]`` (``None`` without + an interface) maps the other mesh's values; ``coeff`` of shape + ``(n, d, m)`` holds ``grad_k(cell) = sum_m coeff[cell, k, m] + b_m``; ``slot[f]`` and ``length[f]`` are the row position within + the owner cell and the row distance of face ``f``, used to place + boundary values. + """ + cache = self._operator_cache(submesh) + signature = tuple(sorted((side, bc_type) for side, (_, bc_type) in bcs.items())) + interface_key = None if interface is None else interface["key"] + key = ("least_squares", signature, interface_key) + if key in cache: + return cache[key] + + n = submesh.npts + d = submesh.dimension + n_int = submesh.n_internal_faces + n_faces = len(submesh.face_owner) + centroids = submesh.cell_centroids + interface_faces = None if interface is None else interface["faces"] + kinds = self._face_bc_kinds(submesh, bcs, interface_faces) + + # Half-face rows: the owner side of every face, then the neighbour + # side of internal faces, so row f (< n_faces) belongs to face f. + row_face = np.concatenate([np.arange(n_faces), np.arange(n_int)]) + cell = np.concatenate([submesh.face_owner, submesh.face_neighbor[:n_int]]) + other = np.concatenate( + [ + submesh.face_neighbor[:n_int], + np.full(n_faces - n_int, -1), + submesh.face_owner[:n_int], + ] + ) + row_kind = kinds[row_face] + toward = np.where( + (row_kind == self._ROW_INTERNAL)[:, np.newaxis], + centroids[np.maximum(other, 0)], + submesh.face_centroids[row_face], + ) + if interface is not None: + toward[interface_faces] = interface["other_centroids"] + other[interface_faces] = interface["other_cells"] + delta = toward - centroids[cell] + length = np.linalg.norm(delta, axis=1) + direction = delta / length[:, np.newaxis] + normal_rows = (row_kind == self._ROW_NEUMANN) | (row_kind == self._ROW_NO_BC) + direction[normal_rows] = submesh.face_normals[row_face[normal_rows]] + length[normal_rows] = 1.0 + + counts = np.bincount(cell, minlength=n) + if np.any(counts != counts[0]): + raise pybamm.DiscretisationError( + "Least-squares gradient needs every cell to have the same " + "number of faces; the mesh connectivity is inconsistent." + ) + m = int(counts[0]) + order = np.argsort(cell, kind="stable") + dirs = direction[order].reshape(n, m, d) + normal = np.einsum("nmi,nmj->nij", dirs, dirs) + coeff = np.einsum("nij,nmj->nim", np.linalg.pinv(normal), dirs) + slot = np.empty(len(cell), dtype=int) + slot[order] = np.arange(len(cell)) % m + + def row_coefficients(rows, k): + return coeff[cell[rows], k, slot[rows]] / length[rows] + + internal = np.nonzero(row_kind == self._ROW_INTERNAL)[0] + dirichlet = np.nonzero(row_kind == self._ROW_DIRICHLET)[0] + across = np.nonzero(row_kind == self._ROW_INTERFACE)[0] + G = [] + G_cross = None if interface is None else [] + for k in range(d): + c_int = row_coefficients(internal, k) + c_dir = row_coefficients(dirichlet, k) + c_across = row_coefficients(across, k) + rows = np.concatenate( + [cell[internal], cell[internal], cell[dirichlet], cell[across]] + ) + cols = np.concatenate( + [other[internal], cell[internal], cell[dirichlet], cell[across]] + ) + data = np.concatenate([c_int, -c_int, -c_dir, -c_across]) + G.append(csr_matrix(coo_matrix((data, (rows, cols)), shape=(n, n)))) + if interface is not None: + G_cross.append( + csr_matrix( + coo_matrix( + (c_across, (cell[across], other[across])), + shape=(n, interface["n_other"]), + ) + ) + ) + + cache[key] = (G, coeff, slot[:n_faces], length[:n_faces], G_cross) + return cache[key] + + def _least_squares_gradient(self, submesh, bcs, repeats=1, interface=None): + """``(matrices, bc_vecs, cross_matrices)`` of the least-squares + gradient: component ``k`` is ``matrices[k] @ u + bc_vecs[k]``, plus + ``cross_matrices[k] @ u_other`` when ``interface`` rows are given + (sizes lifted by ``repeats`` for auxiliary domains).""" + n = submesh.npts + d = submesh.dimension + G, coeff, slot, length, G_cross = self._least_squares_matrices( + submesh, bcs, interface + ) + bc_vecs = [pybamm.Vector(np.zeros(n * repeats)) for _ in range(d)] + for side, (bc_value, bc_type) in bcs.items(): + faces = self._boundary_faces_for_side(submesh, side) + owners = submesh.face_owner[faces] + for k in range(d): + coeffs = coeff[owners, k, slot[faces]] + if bc_type == "Dirichlet": + coeffs = coeffs / length[faces] + else: + coeffs = coeffs * self._neumann_sign(side) + bc_vecs[k] = bc_vecs[k] + self._bc_contribution( + n, len(faces), owners, coeffs, bc_value, repeats=repeats + ) + return G, bc_vecs, G_cross + + def _green_gauss_matrices(self, submesh): + """ + Build (or fetch the cached) Green-Gauss gradient matrices G_k for + k = 0..d-1. + + For each cell i, the gradient component k is: + (grad u)_k,i = (1/V_i) * sum_f [u_f * n_k,f * A_f] + + where u_f is interpolated from owner/neighbor (distance-weighted + for internal faces) or just the owner value (boundary faces). + """ + cache = self._operator_cache(submesh) + if "green_gauss" in cache: + return cache["green_gauss"] + n = submesh.npts + d = submesh.dimension + n_int = submesh.n_internal_faces + + owner = submesh.face_owner + neighbor = submesh.face_neighbor + normals = submesh.face_normals + areas = submesh.face_areas + vol = submesh.cell_volumes + centroids = submesh.cell_centroids + face_centroids = submesh.face_centroids + + G = [csr_matrix((n, n)) for _ in range(d)] + + # --- internal faces: distance-weighted interpolation --- + int_owner = owner[:n_int] + int_neighbor = neighbor[:n_int] + + d_owner = np.linalg.norm(face_centroids[:n_int] - centroids[int_owner], axis=1) + d_neighbor = np.linalg.norm( + face_centroids[:n_int] - centroids[int_neighbor], axis=1 + ) + d_total = d_owner + d_neighbor + w_owner = d_neighbor / d_total # weight for owner value + w_neighbor = d_owner / d_total # weight for neighbor value + + for k in range(d): + nk_A = normals[:n_int, k] * areas[:n_int] + + # distance-weighted face value scatters to owner (+) and + # neighbor (-), each divided by that cell's volume + rows = np.concatenate([int_owner, int_owner, int_neighbor, int_neighbor]) + cols = np.concatenate([int_owner, int_neighbor, int_owner, int_neighbor]) + data = np.concatenate( + [ + w_owner * nk_A / vol[int_owner], + w_neighbor * nk_A / vol[int_owner], + -w_owner * nk_A / vol[int_neighbor], + -w_neighbor * nk_A / vol[int_neighbor], + ] + ) + + G[k] = G[k] + csr_matrix(coo_matrix((data, (rows, cols)), shape=(n, n))) + + # --- boundary faces: u_f = u_owner (zeroth-order extrapolation) --- + n_total = len(owner) + bnd_indices = np.arange(n_int, n_total) + if len(bnd_indices) > 0: + bnd_owner = owner[bnd_indices] + for k in range(d): + nk_A = normals[bnd_indices, k] * areas[bnd_indices] + rows = bnd_owner + cols = bnd_owner + data = nk_A / vol[bnd_owner] + G[k] = G[k] + csr_matrix(coo_matrix((data, (rows, cols)), shape=(n, n))) + + cache["green_gauss"] = G + return G + + # ------------------------------------------------------------------ + # Divergence + # ------------------------------------------------------------------ + + def divergence(self, symbol, discretised_symbol, boundary_conditions): + """Face-flux divergence of a cell-centred vector field. + + Boundary faces use zeroth-order flux extrapolation, so fluxes built + from BC-carrying gradients are rejected (see the raise below). + """ + domain = symbol.domain + submesh = self.mesh[domain] + n = submesh.npts + d = submesh.dimension + repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + # BC-bearing fluxes must use the div(D*grad(u)) TPFA intercept; + # here the prescribed boundary flux would be silently ignored. + bc_gradient_parents = [ + node.child + for node in symbol.pre_order() + if isinstance(node, pybamm.Gradient) and node.child in boundary_conditions + ] + if bc_gradient_parents: + names = sorted({parent.name for parent in bc_gradient_parents}) + raise pybamm.DiscretisationError( + f"Cannot discretise div of a general flux containing grad of " + f"{names} on an unstructured mesh: the boundary conditions " + "would be ignored and the result would not be conservative. " + "Write the equation as div(D * grad(u)) (a single product) so " + "the TPFA discretisation applies the boundary conditions." + ) + + if isinstance(discretised_symbol, pybamm.VectorField): + comps = discretised_symbol.components + elif isinstance(discretised_symbol, (list, tuple)): + comps = list(discretised_symbol) + else: + raise pybamm.DiscretisationError( + "FiniteVolumeUnstructured.divergence expects a VectorField or " + f"list of {d} component arrays, got {type(discretised_symbol)}" + ) + + D_components = self._divergence_matrices(submesh) + + result = pybamm.Vector(np.zeros(n * repeats)) + for k in range(d): + Dk = csr_matrix(kron(eye(repeats, dtype=np.float64), D_components[k])) + result = result + pybamm.Matrix(Dk) @ comps[k] + + return result + + def _divergence_matrices(self, submesh): + """Divergence matrices ``D_k``: ``(div F)_i = (1/V_i) sum_f F_k,f n_k,f A_f``. + + The face-value interpolation is identical to the Green-Gauss + gradient's, so the two operators share one assembly. + """ + return self._green_gauss_matrices(submesh) + + # ------------------------------------------------------------------ + # gradient_squared |grad u|^2 + # ------------------------------------------------------------------ + + def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): + """Pointwise ``|grad u|^2`` via :meth:`gradient`.""" + grad = self.gradient(symbol, discretised_symbol, boundary_conditions) + result = None + for comp in grad.components: + sq = comp**2 + result = sq if result is None else result + sq + return result + + # ------------------------------------------------------------------ + # Binary operator handling (scalar * VectorField, etc.) + # ------------------------------------------------------------------ + + def process_binary_operators(self, bin_op, left, right, disc_left, disc_right): + """Apply a binary operator componentwise when either operand is a + :class:`pybamm.VectorField`, lifting scalars to N components.""" + if isinstance(disc_left, pybamm.VectorField) or isinstance( + disc_right, pybamm.VectorField + ): + if isinstance(disc_left, pybamm.VectorField) and isinstance( + disc_right, pybamm.VectorField + ): + n = disc_left.n_components + elif isinstance(disc_left, pybamm.VectorField): + n = disc_left.n_components + disc_right = pybamm.VectorField(*[disc_right] * n) + else: + n = disc_right.n_components + disc_left = pybamm.VectorField(*[disc_left] * n) + + new_comps = [ + pybamm.simplify_if_constant( + bin_op.create_copy( + [disc_left.components[k], disc_right.components[k]] + ) + ) + for k in range(n) + ] + return pybamm.VectorField(*new_comps) + + return bin_op._binary_new_copy(disc_left, disc_right) + + # ------------------------------------------------------------------ + # Integral operators + # ------------------------------------------------------------------ + + def integral( + self, child, discretised_child, integration_dimension, integration_variable=None + ): + """Volume integral over the primary domain (cell-volume weights).""" + int_mat = self.definite_integral_matrix( + child, integration_dimension=integration_dimension + ) + return int_mat @ discretised_child + + def definite_integral_matrix( + self, child, vector_type="row", integration_dimension="primary" + ): + """Cell-volume weights of the primary domain as a + :class:`pybamm.Matrix`, one block per auxiliary-domain repeat. + + Parameters + ---------- + child : pybamm.Symbol + The symbol being integrated. + vector_type : str, optional + ``"row"`` (default) or ``"column"``. + integration_dimension : str, optional + Only ``"primary"`` is supported: cells have no secondary + structure to integrate over. + + Raises + ------ + NotImplementedError + For a non-primary ``integration_dimension``. + """ + if integration_dimension != "primary": + raise NotImplementedError( + f"Integral in the {integration_dimension!r} dimension is not " + "implemented on unstructured meshes; only the primary (cell) " + "dimension can be integrated." + ) + if vector_type not in ("row", "column"): + raise pybamm.DiscretisationError( + f"vector_type must be 'row' or 'column', not {vector_type!r}" + ) + submesh = self.mesh[child.domain] + repeats = self._get_auxiliary_domain_repeats(child.domains) + shape = (1, -1) if vector_type == "row" else (-1, 1) + block = csr_matrix(submesh.cell_volumes.reshape(shape)) + return pybamm.Matrix(csr_matrix(kron(eye(repeats, dtype=np.float64), block))) + + def boundary_integral(self, child, discretised_child, region): + """Integral of the owner-cell values of ``child`` over the boundary + faces of ``region`` (``"entire"`` = all exterior faces).""" + submesh = self.mesh[child.domain] + repeats = self._get_auxiliary_domain_repeats(child.domains) + + if region == "entire": + # every exterior boundary face; iface_* buckets are internal + # interfaces, not part of the domain boundary + iface = [ + indices + for tag, indices in submesh.boundary_faces.items() + if tag.startswith("iface_") + ] + face_indices = np.setdiff1d( + np.arange(submesh._boundary_face_start, len(submesh.face_owner)), + np.concatenate(iface) if iface else np.array([], dtype=int), + ) + else: + face_indices = self._boundary_faces_for_side(submesh, region) + n = submesh.npts + + owners = submesh.face_owner[face_indices] + face_areas = submesh.face_areas[face_indices] + + row = np.zeros(n) + np.add.at(row, owners, face_areas) + mat = csr_matrix(row.reshape(1, -1)) + mat = csr_matrix(kron(eye(repeats, dtype=np.float64), mat)) + + return pybamm.Matrix(mat) @ discretised_child + + # ------------------------------------------------------------------ + # boundary_value_or_flux + # ------------------------------------------------------------------ + + _CORNER_SIDES = { + "top-right": ("top", "right"), + "top-left": ("top", "left"), + "bottom-right": ("bottom", "right"), + "bottom-left": ("bottom", "left"), + } + + def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): + """Owner-cell values on a boundary side (zeroth-order boundary + value); corner sides return the single closest boundary cell.""" + if isinstance(symbol, pybamm.BoundaryGradient): + raise NotImplementedError( + "BoundaryGradient is not implemented for unstructured meshes; " + "returning the boundary value instead would be silently wrong." + ) + submesh = self.mesh[discretised_child.domain] + n = submesh.npts + repeats = self._get_auxiliary_domain_repeats(discretised_child.domains) + + side = symbol.side + + if side in self._CORNER_SIDES: + return self._corner_boundary_value( + submesh, n, repeats, side, discretised_child + ) + + face_indices = self._boundary_faces_for_side(submesh, side) + n_bnd = len(face_indices) + owners = submesh.face_owner[face_indices] + + sub_matrix = csr_matrix( + (np.ones(n_bnd), (np.arange(n_bnd), owners)), + shape=(n_bnd, n), + ) + + mat = csr_matrix(kron(eye(repeats, dtype=np.float64), sub_matrix)) + bv_vector = pybamm.Matrix(mat) + + out = bv_vector @ discretised_child + out.clear_domains() + return out + + def _corner_boundary_value(self, submesh, n, repeats, side, discretised_child): + """Extract the value from the boundary cell closest to a corner of + the (x, z) bounding box. + + Zeroth-order (cell value) regardless of the ``extrapolation`` option. + Candidates are restricted to cells owning boundary faces on the two + named sides, so interior cells of non-convex domains are never + picked; in 3D, ties across y pick the lowest cell index. + """ + tb_side, lr_side = self._CORNER_SIDES[side] + centroids = submesh.cell_centroids + + # cells owning boundary faces on either named side (fall back to all + # boundary-owner cells when a side bucket is missing) + candidates = np.unique( + np.concatenate( + [ + submesh.face_owner[submesh.boundary_faces[tag]] + for tag in (tb_side, lr_side) + if tag in submesh.boundary_faces + ] + or [submesh.face_owner[submesh._boundary_face_start :]] + ) + ) + + x_coords = centroids[candidates, 0] + z_coords = centroids[candidates, -1] + target_x = x_coords.max() if lr_side == "right" else x_coords.min() + target_z = z_coords.max() if tb_side == "top" else z_coords.min() + + dists = (x_coords - target_x) ** 2 + (z_coords - target_z) ** 2 + cell_idx = int(candidates[np.argmin(dists)]) + + sub_matrix = csr_matrix( + (np.ones(1), (np.zeros(1, dtype=int), [cell_idx])), + shape=(1, n), + ) + mat = csr_matrix(kron(eye(repeats, dtype=np.float64), sub_matrix)) + out = pybamm.Matrix(mat) @ discretised_child + out.clear_domains() + return out + + # ------------------------------------------------------------------ + # internal_neumann_condition + # ------------------------------------------------------------------ + + def internal_neumann_condition( + self, + left_symbol_disc, + right_symbol_disc, + left_mesh, + right_mesh, + left_bcs=None, + right_bcs=None, + ): + """Normal gradient across the interface between two submeshes, one + value per interface face (outward from ``left_mesh``). + + On unstructured meshes this is the two-point difference plus the + non-orthogonal cross term (see :meth:`_tpfa_matrix`), whose face + gradient is fitted on both sides with the interface faces as + cross-mesh rows. ``left_bcs``/``right_bcs`` are the external + boundary conditions of each side, used by that fit. + """ + from pybamm.meshes.unstructured_submesh import UnstructuredSubMesh + + repeats = self._get_auxiliary_domain_repeats(left_symbol_disc.domains) + + if repeats != self._get_auxiliary_domain_repeats(right_symbol_disc.domains): + raise pybamm.DomainError( + "Number of secondary points in subdomains do not match" + ) + + if isinstance(left_mesh, UnstructuredSubMesh): + return self._internal_neumann_unstructured( + left_symbol_disc, + right_symbol_disc, + left_mesh, + right_mesh, + repeats, + left_bcs or {}, + right_bcs or {}, + ) + else: + return self._internal_neumann_structured( + left_symbol_disc, + right_symbol_disc, + left_mesh, + right_mesh, + repeats, + ) + + def _internal_neumann_unstructured( + self, + left_symbol_disc, + right_symbol_disc, + left_mesh, + right_mesh, + repeats, + left_bcs=None, + right_bcs=None, + ): + left_bcs = left_bcs or {} + right_bcs = right_bcs or {} + # Find the interface_data entry that pairs left_mesh with right_mesh. + # Each entry stores ``other_mesh`` so multi-neighbor topologies pick + # the correct partner instead of grabbing the first dict value. + interface = next( + ( + data + for data in left_mesh.interface_data.values() + if data.get("other_mesh") is right_mesh + ), + None, + ) + + if interface is None: + rev = next( + ( + data + for data in right_mesh.interface_data.values() + if data.get("other_mesh") is left_mesh + ), + None, + ) + if rev is not None: + interface = { + "left_cells": rev["right_cells"], + "right_cells": rev["left_cells"], + "left_faces": rev["right_faces"], + "right_faces": rev["left_faces"], + "face_areas": rev["face_areas"], + "cell_distances": rev["cell_distances"], + } + + if interface is None: + raise pybamm.DiscretisationError( + "No interface data pairs these two unstructured meshes, so " + "the internal gradient between them cannot be formed and the " + "domains would be silently decoupled. Check that both meshes " + "carry boundary tags (e.g. detect_box_boundaries() for " + "axis-aligned boxes) so interface discovery can pair their " + "faces." + ) + + left_cells = interface["left_cells"] + right_cells = interface["right_cells"] + left_faces = interface["left_faces"] + n_faces = len(left_cells) + d = left_mesh.dimension + + def lift(matrix): + return pybamm.Matrix( + csr_matrix(kron(eye(repeats, dtype=np.float64), matrix)) + ) + + def without_domains(expr): + expr.clear_domains() + return expr + + def tile(values): + return pybamm.Vector(np.tile(values, repeats)) + + left_sub = csr_matrix( + (np.ones(n_faces), (np.arange(n_faces), left_cells)), + shape=(n_faces, left_mesh.npts), + ) + right_sub = csr_matrix( + (np.ones(n_faces), (np.arange(n_faces), right_cells)), + shape=(n_faces, right_mesh.npts), + ) + + # n = alpha e + k with e the unit left-to-right centroid direction and + # n the interface normal, outward from the left cells. + normals = left_mesh.face_normals[left_faces] + c_left = left_mesh.cell_centroids[left_cells] + c_right = right_mesh.cell_centroids[right_cells] + delta = c_right - c_left + dist = np.linalg.norm(delta, axis=1) + e_ij = delta / dist[:, np.newaxis] + alpha = self._alpha(np.sum(normals * e_ij, axis=1)) + k_vec = normals - alpha[:, np.newaxis] * e_ij + + two_point = diags(alpha / dist) + value = without_domains( + lift(two_point @ right_sub) @ right_symbol_disc + ) - without_domains(lift(two_point @ left_sub) @ left_symbol_disc) + + k_vec = self._drop_orthogonal(k_vec) + if not k_vec.any(): + return value + + face_centroids = left_mesh.face_centroids[left_faces] + d_left = np.linalg.norm(face_centroids - c_left, axis=1) + d_right = np.linalg.norm(face_centroids - c_right, axis=1) + w_left = d_right / (d_left + d_right) + + def side_gradient(mesh, bcs, faces, other, other_cells, own_disc, other_disc): + interface_rows = { + "key": (id(other), hash(faces.tobytes())), + "faces": faces, + "other_cells": other_cells, + "other_centroids": other.cell_centroids[other_cells], + "n_other": other.npts, + } + G, bc_vecs, G_cross = self._least_squares_gradient( + mesh, bcs, repeats, interface=interface_rows + ) + return [ + without_domains(lift(G[k]) @ own_disc) + + without_domains(lift(G_cross[k]) @ other_disc) + + without_domains(bc_vecs[k]) + for k in range(d) + ] + + grad_left = side_gradient( + left_mesh, + left_bcs, + left_faces, + right_mesh, + right_cells, + left_symbol_disc, + right_symbol_disc, + ) + grad_right = side_gradient( + right_mesh, + right_bcs, + interface["right_faces"], + left_mesh, + left_cells, + right_symbol_disc, + left_symbol_disc, + ) + for k in range(d): + face_gradient = ( + lift(diags(w_left) @ left_sub) @ grad_left[k] + + lift(diags(1.0 - w_left) @ right_sub) @ grad_right[k] + ) + value = value + face_gradient * tile(k_vec[:, k]) + return value + + def _internal_neumann_structured( + self, + left_symbol_disc, + right_symbol_disc, + left_mesh, + right_mesh, + repeats, + ): + """Fallback for structured meshes (same logic as FiniteVolume).""" + left_npts = left_mesh.npts + right_npts = right_mesh.npts + + left_sub_matrix = np.zeros((1, left_npts)) + left_sub_matrix[0][left_npts - 1] = 1 + left_matrix = pybamm.Matrix( + csr_matrix(kron(eye(repeats, dtype=np.float64), left_sub_matrix)) + ) + + right_sub_matrix = np.zeros((1, right_npts)) + right_sub_matrix[0][0] = 1 + right_matrix = pybamm.Matrix( + csr_matrix(kron(eye(repeats, dtype=np.float64), right_sub_matrix)) + ) + + # structured fallback: 1D submeshes expose ``nodes``, not ``vertices`` + right_mesh_x = right_mesh.nodes[0] + left_mesh_x = left_mesh.nodes[-1] + dx = right_mesh_x - left_mesh_x + + dy_r = (right_matrix / dx) @ right_symbol_disc + dy_r.clear_domains() + dy_l = (left_matrix / dx) @ left_symbol_disc + dy_l.clear_domains() + + return dy_r - dy_l + + # ------------------------------------------------------------------ + # concatenation + # ------------------------------------------------------------------ + + def concatenation(self, disc_children): + """See :meth:`pybamm.SpatialMethod.concatenation`.""" + return pybamm.domain_concatenation(disc_children, self.mesh) + + # ------------------------------------------------------------------ + # Not implemented + # ------------------------------------------------------------------ + + def indefinite_integral(self, child, discretised_child, direction): + raise NotImplementedError( + "Indefinite integral is not supported on unstructured meshes. " + "Use the direct PDE form instead." + ) + + def delta_function(self, symbol, discretised_symbol): + raise NotImplementedError( + "Delta function is not supported on unstructured meshes." + ) diff --git a/packages/pybamm/src/pybamm/spatial_methods/spatial_method.py b/packages/pybamm/src/pybamm/spatial_methods/spatial_method.py index c6000e86c6..8dda168d6f 100644 --- a/packages/pybamm/src/pybamm/spatial_methods/spatial_method.py +++ b/packages/pybamm/src/pybamm/spatial_methods/spatial_method.py @@ -330,6 +330,32 @@ def internal_neumann_condition( raise NotImplementedError + def set_internal_bcs_for_concat(self, disc, var, children, outer_bcs): + """ + Hook for spatial methods that own their internal-BC logic for + concatenated variables (e.g. graph topologies on unstructured + meshes). + + Parameters + ---------- + disc : :class:`pybamm.Discretisation` + The discretisation, for processing child symbols + var : :class:`pybamm.Concatenation` + The concatenated variable whose boundary conditions are being set + children : list of :class:`pybamm.Symbol` + The orphaned children of ``var`` + outer_bcs : dict + The user-supplied boundary conditions for ``var``, + ``{side: (value, type)}`` + + Returns + ------- + dict or None + ``{child: {side: (value, type)}}`` to replace the default + 1D-stack pairwise routine, or ``None`` to use it. + """ + return + def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): """ Returns the boundary value or flux using the appropriate expression for the diff --git a/packages/pybamm/tests/unit/test_discretisations/test_discretisation.py b/packages/pybamm/tests/unit/test_discretisations/test_discretisation.py index 077296e4af..2641660460 100644 --- a/packages/pybamm/tests/unit/test_discretisations/test_discretisation.py +++ b/packages/pybamm/tests/unit/test_discretisations/test_discretisation.py @@ -66,6 +66,43 @@ def test_add_internal_boundary_conditions(self): for child in c_e.children: assert child in disc.bcs + def test_internal_boundary_conditions_require_left_right(self): + model = pybamm.BaseModel() + c_e_n = pybamm.Variable("c_e_n", ["negative electrode"]) + c_e_s = pybamm.Variable("c_e_s", ["separator"]) + c_e_p = pybamm.Variable("c_e_p", ["positive electrode"]) + c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p) + bc = (pybamm.Scalar(0), "Neumann") + model.boundary_conditions = {c_e: {"top": bc, "bottom": bc}} + + mesh = get_mesh_for_testing() + spatial_methods = {"macroscale": SpatialMethodForTesting()} + disc = pybamm.Discretisation(mesh, spatial_methods) + disc.set_variable_slices([c_e_n, c_e_s, c_e_p]) + disc.bcs = model.boundary_conditions + # the base hook declines, so the legacy 1D-stack routine runs (and + # raises here because it requires left/right BCs) + assert ( + spatial_methods["macroscale"].set_internal_bcs_for_concat( + disc, c_e, c_e.orphans, disc.bcs[c_e] + ) + is None + ) + with pytest.raises(pybamm.DiscretisationError, match="'left' and 'right'"): + disc.set_internal_boundary_conditions(model) + + def test_custom_side_containing_tab_substring(self): + # arbitrary mesh region tags containing "tab" (e.g. "tab_weld") must + # not be routed into the legacy tab-condition check, which raises + # ModelError outside the current-collector domain + mesh = get_mesh_for_testing() + spatial_methods = {"macroscale": SpatialMethodForTesting()} + disc = pybamm.Discretisation(mesh, spatial_methods) + var = pybamm.Variable("var", domain=["negative electrode"]) + disc.set_variable_slices([var]) + disc.bcs = {var: {"tab_region": (pybamm.Scalar(0), "Neumann")}} + disc.process_symbol(var) + def test_add_internal_boundary_conditions_symbolic(self): submesh_types = { "left domain": pybamm.SymbolicUniform1DSubMesh, diff --git a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py index f37772781d..91faff2387 100644 --- a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py +++ b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py @@ -1515,3 +1515,30 @@ def test_optimize_ordering_preserves_interface_pairing(self): np.testing.assert_allclose( other_pre, mesh_r.cell_centroids[mirror["left_cells"]] ) + + +class TestContainsPoints: + def test_2d_even_odd_rule(self): + nodes, elements = _unit_square_two_triangles() + mesh = UnstructuredSubMesh(nodes, elements) + points = np.array([[0.5, 0.5], [2.0, 2.0], [0.5, 0.0], [-1e-3, 0.5]]) + np.testing.assert_array_equal( + mesh.contains_points(points), [True, False, True, False] + ) + # loops are cached on the mesh after the first query + assert mesh._cached_boundary_loops is not None + + def test_3d_delegates_to_winding_number(self): + nodes, elements = _unit_cube_five_tets() + mesh = UnstructuredSubMesh(nodes, elements) + points = np.array([[0.5, 0.5, 0.5], [2.0, 2.0, 2.0]]) + np.testing.assert_array_equal(mesh.contains_points(points), [True, False]) + np.testing.assert_array_equal( + mesh.contains_points(points), mesh.contains_points_3d(points) + ) + + def test_2d_without_boundary_loops_returns_none(self): + nodes, elements = _unit_square_two_triangles() + mesh = UnstructuredSubMesh(nodes, elements) + mesh._cached_boundary_loops = None + assert mesh.contains_points(np.array([[0.5, 0.5]])) is None diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py index ccee5ab4c2..a8b2e7395e 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py @@ -1,6 +1,9 @@ # # Tests for the basic lithium-ion models # +import numpy as np +import pytest + import pybamm @@ -25,3 +28,22 @@ def test_dfn_composite_well_posed(self): def test_dfn_2d(self): model = pybamm.lithium_ion.BasicDFN2D() model.check_well_posedness() + + @pytest.mark.filterwarnings("ignore:Could not determine how to combine submeshes") + def test_dfn_2d_vector_field_variable(self): + # A VectorField variable on a structured 2D mesh cannot be read + # directly, but must fail with guidance rather than an opaque error, + # and extracting a component must work. + model = pybamm.lithium_ion.BasicDFN2D() + model.variables["Electrolyte current density x [A.m-2]"] = pybamm.Component( + model.variables["Electrolyte current density [A.m-2]"], 0 + ) + var_pts = {k: 5 for k in model.default_var_pts} + sim = pybamm.Simulation(model, var_pts=var_pts) + solution = sim.solve([0, 10]) + + with pytest.raises(NotImplementedError, match=r"pybamm\.Component"): + solution["Electrolyte current density [A.m-2]"] + + component = solution["Electrolyte current density x [A.m-2]"] + assert np.all(np.isfinite(component(t=5))) diff --git a/packages/pybamm/tests/unit/test_solvers/test_processed_variable.py b/packages/pybamm/tests/unit/test_solvers/test_processed_variable.py index 82c359cb9d..94d17788ba 100644 --- a/packages/pybamm/tests/unit/test_solvers/test_processed_variable.py +++ b/packages/pybamm/tests/unit/test_solvers/test_processed_variable.py @@ -2124,3 +2124,412 @@ def test_process_variable_unstructured_detection(self): ) assert isinstance(processed_var, pybamm.ProcessedVariableUnstructured) + + +class TestProcessedVariableUnstructuredFVM: + @staticmethod + def _make_setup(dim=2, n=6): + from pybamm.meshes.unstructured_submesh import UnstructuredMeshGenerator + + domain = "negative electrode" + x = pybamm.SpatialVariable("x_n", domain=[domain], coord_sys="cartesian") + if dim == 2: + z = pybamm.SpatialVariable( + "z_2d", domain=[domain], coord_sys="cartesian", direction="tb" + ) + geometry = { + domain: {x: {"min": 0.0, "max": 1.0}, z: {"min": 0.0, "max": 1.0}} + } + var_pts = {x: n, z: n} + else: + y = pybamm.SpatialVariable("y", domain=[domain], coord_sys="cartesian") + z = pybamm.SpatialVariable("z", domain=[domain], coord_sys="cartesian") + geometry = { + domain: { + x: {"min": 0.0, "max": 1.0}, + y: {"min": 0.0, "max": 1.0}, + z: {"min": 0.0, "max": 1.0}, + } + } + var_pts = {x: n, y: n, z: n} + + mesh = pybamm.Mesh(geometry, {domain: UnstructuredMeshGenerator()}, var_pts) + disc = pybamm.Discretisation(mesh, {domain: pybamm.FiniteVolumeUnstructured()}) + var = pybamm.Variable("u", domain=[domain]) + disc.set_variable_slices([var]) + var_disc = disc.process_symbol(var) + return geometry, mesh[domain], disc, var, var_disc + + def _make_pv(self, var_disc, geometry, t_sol, y_sol): + var_casadi = to_casadi(var_disc, y_sol) + model = pybamm.BaseModel() + model._geometry = geometry + solution = pybamm.Solution(t_sol, y_sol, model, {}) + return pybamm.process_variable("u", [var_disc], [var_casadi], solution) + + def test_2d_dispatch_and_interpolation(self): + geometry, submesh, _, _, var_disc = self._make_setup(dim=2) + centroid_x = submesh.cell_centroids[:, 0] + t_sol = np.linspace(0, 1, 5) + y_sol = centroid_x[:, np.newaxis] * (1 + t_sol)[np.newaxis, :] + + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + assert isinstance(pv, pybamm.ProcessedVariableUnstructuredFVM) + assert pv.dimensions == 2 + + # At solver times with no spatial coords: raw cell data + np.testing.assert_allclose(pv(t_sol), y_sol, rtol=1e-12) + + # Linear time interpolation between solver times + np.testing.assert_allclose(pv(0.5).ravel(), centroid_x * 1.5, rtol=1e-10) + + # Spatial interpolation reproduces the linear field in the interior + x_q = np.linspace(0.4, 0.6, 3) + z_q = np.linspace(0.4, 0.6, 3) + result = pv(0.5, x=x_q, z=z_q) + assert result.shape == (3, 3) + expected = 1.5 * x_q[:, np.newaxis] * np.ones((1, 3)) + np.testing.assert_allclose(result, expected, rtol=1e-8) + + # Vector time: interpolation per time slice, time on the last axis + result_t = pv(t_sol[:2], x=x_q, z=z_q) + assert result_t.shape == (3, 3, 2) + np.testing.assert_allclose( + result_t[..., 0], x_q[:, np.newaxis] * np.ones((1, 3)), rtol=1e-8 + ) + + def test_2d_outside_domain_is_nan(self): + geometry, submesh, _, _, var_disc = self._make_setup(dim=2, n=3) + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts, 2)) + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + outside = pv(0.5, x=np.array([-0.5]), z=np.array([0.5])) + assert np.isnan(outside).all() + # Second call goes through the cached boundary mask + outside_again = pv(0.5, x=np.array([-0.5]), z=np.array([0.5])) + assert np.isnan(outside_again).all() + inside = pv(0.5, x=np.array([0.5]), z=np.array([0.5])) + np.testing.assert_allclose(inside, 1.0, rtol=1e-10) + + def test_call_coordinate_handling(self): + geometry, submesh, _, _, var_disc = self._make_setup(dim=2) + centroid_x = submesh.cell_centroids[:, 0] + t_sol = np.linspace(0, 1, 5) + y_sol = centroid_x[:, np.newaxis] * np.ones_like(t_sol)[np.newaxis, :] + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + # z only: x defaults to the domain midplane (0.5), not 0.0 + z_q = np.linspace(0.4, 0.6, 3) + result = pv(0.5, z=z_q) + assert result.shape == (1, 3) + np.testing.assert_allclose(result, 0.5, rtol=1e-8) + + # length-1 time arrays keep the time axis; scalars drop it + x_q = np.array([0.5]) + assert pv(np.array([0.5]), x=x_q, z=z_q).shape == (1, 3, 1) + assert pv(0.5, x=x_q, z=z_q).shape == (1, 3) + + # fill_value replaces NaN outside the domain + outside = pv(0.5, x=np.array([-0.5]), z=np.array([0.5]), fill_value=-7.0) + np.testing.assert_allclose(outside, -7.0) + + # r/R are not unstructured coordinates + with pytest.raises(ValueError, match="no r or R"): + pv(0.5, r=np.array([0.5])) + + # y is not a 2D-mesh coordinate; silently ignoring it would return + # midplane values for a query the user thinks is at y + with pytest.raises(ValueError, match="no y coordinate"): + pv(0.5, x=x_q, y=np.array([0.5])) + + def test_nan_time_slice_does_not_corrupt_others(self): + # One all-NaN time step must propagate as NaN without triggering + # nearest-neighbour refill of the valid time steps (rows are only + # refilled when NaN in every column, the outside-hull signature). + geometry, submesh, _, _, var_disc = self._make_setup(dim=2) + centroid_x = submesh.cell_centroids[:, 0] + t_sol = np.array([0.0, 1.0]) + y_sol = np.column_stack([centroid_x, np.full_like(centroid_x, np.nan)]) + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + x_q = np.linspace(0.4, 0.6, 3) + z_q = np.linspace(0.4, 0.6, 3) + result = pv(t_sol, x=x_q, z=z_q) + np.testing.assert_allclose( + result[..., 0], x_q[:, np.newaxis] * np.ones((1, 3)), rtol=1e-8 + ) + assert np.isnan(result[..., 1]).all() + + def test_time_integral_raises(self): + geometry, submesh, _, _, var_disc = self._make_setup(dim=2, n=3) + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts, 2)) + var_casadi = to_casadi(var_disc, y_sol) + model = pybamm.BaseModel() + model._geometry = geometry + solution = pybamm.Solution(t_sol, y_sol, model, {}) + time_integral = object() # any non-None marker + with pytest.raises(NotImplementedError, match="Time integrals"): + pybamm.process_variable( + "u", [var_disc], [var_casadi], solution, time_integral=time_integral + ) + + def test_vector_field_pv_interface(self): + geometry, submesh, disc, _, _ = self._make_setup(dim=2, n=3) + var = pybamm.Variable("u", domain=["negative electrode"]) + disc.set_variable_slices([var]) + grad_disc = disc.process_symbol(pybamm.grad(var)) + grad_disc.mesh = submesh + for comp in grad_disc.components: + comp.mesh = submesh + + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts, 2)) + comp_casadi = [to_casadi(comp, y_sol) for comp in grad_disc.components] + model = pybamm.BaseModel() + model._geometry = geometry + solution = pybamm.Solution(t_sol, y_sol, model, {}) + pv = pybamm.process_variable("grad u", [grad_disc], [comp_casadi], solution) + + assert isinstance(pv, pybamm.ProcessedVariableVectorFieldUnstructuredFVM) + # entries and data return one array per component, not component 0 + assert isinstance(pv.entries, tuple) + assert len(pv.entries) == 2 + assert isinstance(pv.data, tuple) + # merging across solution segments is not supported: clear error + with pytest.raises(NotImplementedError, match="merged across"): + pv.update(pv, solution) + + def test_scalar_reduction_routes_to_0d(self): + # Max/Min of a spatial variable keep the domain (and hence the + # unstructured mesh) but evaluate to one value: 0D in space + geometry, submesh, _, _, var_disc = self._make_setup(dim=2, n=3) + max_disc = pybamm.Max(var_disc) + max_disc.mesh = submesh + t_sol = np.array([0.0, 1.0]) + y_sol = np.arange(submesh.npts)[:, np.newaxis] * (1 + t_sol)[np.newaxis, :] + var_casadi = to_casadi(max_disc, y_sol) + model = pybamm.BaseModel() + model._geometry = geometry + solution = pybamm.Solution(t_sol, y_sol, model, {}) + pv = pybamm.process_variable("max u", [max_disc], [var_casadi], solution) + from pybamm.solvers.processed_variable import ProcessedVariable0D + + assert isinstance(pv, ProcessedVariable0D) + np.testing.assert_allclose( + pv(t_sol), (submesh.npts - 1) * (1 + t_sol), rtol=1e-12 + ) + + def test_single_cell_mesh_variable_stays_spatial(self): + # a one-cell mesh also evaluates to size 1, but it is a genuine + # spatial variable and must not be routed to the 0D PV + from pybamm.meshes.unstructured_submesh import UnstructuredSubMesh + + submesh = UnstructuredSubMesh( + np.array( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + ), + np.array([[0, 1, 2, 3]]), + ) + submesh.detect_box_boundaries() + var_disc = pybamm.StateVector(slice(0, 1)) + var_disc.mesh = submesh + t_sol = np.array([0.0, 1.0]) + y_sol = np.array([[1.0, 2.0]]) + var_casadi = to_casadi(var_disc, y_sol) + model = pybamm.BaseModel() + model._geometry = {"mesh": {}} + solution = pybamm.Solution(t_sol, y_sol, model, {}) + pv = pybamm.process_variable("u", [var_disc], [var_casadi], solution) + assert isinstance(pv, pybamm.ProcessedVariableUnstructuredFVM) + + def test_disconnected_component_is_not_masked(self): + from pybamm.meshes.unstructured_submesh import ( + UnstructuredSubMesh, + _make_quad_grid, + ) + + # two disjoint unit squares with a gap between x = 1 and x = 2 + nodes_a, elems_a = _make_quad_grid(np.linspace(0, 1, 3), np.linspace(0, 1, 3)) + nodes_b, elems_b = _make_quad_grid(np.linspace(2, 3, 3), np.linspace(0, 1, 3)) + nodes = np.vstack([nodes_a, nodes_b]) + elements = np.vstack([elems_a, elems_b + len(nodes_a)]) + submesh = UnstructuredSubMesh(nodes, elements) + submesh.detect_box_boundaries() + + var_disc = pybamm.StateVector(slice(0, submesh.npts)) + var_disc.mesh = submesh + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts, 2)) + geometry = {"domain": {}} + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + in_second_square = pv(0.5, x=np.array([2.5]), z=np.array([0.5])) + np.testing.assert_allclose(in_second_square, 1.0, rtol=1e-10) + in_gap = pv(0.5, x=np.array([1.5]), z=np.array([0.5])) + assert np.isnan(in_gap).all() + + def test_auxiliary_domain_variable_raises(self): + geometry, submesh, _, _, _ = self._make_setup(dim=2, n=3) + # hand-build a repeated state vector, as an auxiliary-domain + # variable's discretisation would produce + repeats = 4 + var_disc = pybamm.StateVector(slice(0, submesh.npts * repeats)) + var_disc.mesh = submesh + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts * repeats, 2)) + with pytest.raises(NotImplementedError, match="auxiliary domains"): + self._make_pv(var_disc, geometry, t_sol, y_sol) + + def test_3d_dispatch_slices_and_mask(self): + geometry, submesh, _, _, var_disc = self._make_setup(dim=3, n=3) + t_sol = np.array([0.0, 1.0]) + # Spatially-constant field with linear time dependence: 7 * (1 + t) + y_sol = 7.0 * np.ones((submesh.npts, 1)) * (1 + t_sol)[np.newaxis, :] + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + assert isinstance(pv, pybamm.ProcessedVariableUnstructuredFVM) + assert pv.dimensions == 3 + + # Scalar-time spatial query inside the domain + result = pv(0.5, x=np.array([0.5]), y=np.array([0.5]), z=np.array([0.5])) + assert result.shape == (1, 1, 1) + np.testing.assert_allclose(result, 10.5, rtol=1e-10) + + # Vector time keeps time on the last axis + result_t = pv( + t_sol, x=np.array([0.3, 0.7]), y=np.array([0.5]), z=np.array([0.5]) + ) + assert result_t.shape == (2, 1, 1, 2) + np.testing.assert_allclose(result_t[..., 0], 7.0, rtol=1e-10) + np.testing.assert_allclose(result_t[..., 1], 14.0, rtol=1e-10) + + # Points outside the domain are masked (3D winding-number path) + outside = pv(0.5, x=np.array([-1.0]), y=np.array([0.5]), z=np.array([0.5])) + assert np.isnan(outside).all() + + def test_vector_field_via_solution_2d(self): + """Requesting a VectorField variable from a Solution goes through the + per-component casadi wiring and the unstructured vector-field PV.""" + geometry, submesh, disc, var, _ = self._make_setup(dim=2, n=4) + domain = "negative electrode" + + flux = pybamm.VectorField( + pybamm.PrimaryBroadcast(pybamm.Scalar(2), domain), + pybamm.PrimaryBroadcast(pybamm.Scalar(-3), domain), + ) + model = pybamm.BaseModel() + model.rhs = {var: pybamm.Scalar(0)} + model.initial_conditions = {var: pybamm.Scalar(0)} + model.variables = {"u": var, "flux": flux} + model_disc = disc.process_model(model, inplace=False) + model_disc._geometry = geometry + + centroid_x = submesh.cell_centroids[:, 0] + t_sol = np.array([0.0, 1.0]) + y_sol = centroid_x[:, np.newaxis] * (1 + t_sol)[np.newaxis, :] + solution = pybamm.Solution(t_sol, y_sol, model_disc, {}) + + scalar_pv = solution["u"] + assert isinstance(scalar_pv, pybamm.ProcessedVariableUnstructuredFVM) + + flux_pv = solution["flux"] + assert isinstance(flux_pv, pybamm.ProcessedVariableVectorFieldUnstructuredFVM) + assert flux_pv.is_vector_field + assert flux_pv.n_components == 2 + assert flux_pv.dimensions == 2 + + # entries returns one array per component + entries = flux_pv.entries + assert isinstance(entries, tuple) + np.testing.assert_allclose(entries[0], 2.0, rtol=1e-12) + np.testing.assert_allclose(entries[1], -3.0, rtol=1e-12) + + # Calling returns one array per component + comps = flux_pv(t=t_sol) + assert isinstance(comps, tuple) + assert len(comps) == 2 + np.testing.assert_allclose(comps[0], 2.0, rtol=1e-12) + np.testing.assert_allclose(comps[1], -3.0, rtol=1e-12) + + # QuickPlot has no unstructured support in this PR and must say so + with pytest.raises(NotImplementedError, match="unstructured meshes"): + pybamm.QuickPlot(solution, ["u"]) + + def test_vector_field_3d(self): + geometry, submesh, disc, _, _ = self._make_setup(dim=3, n=3) + domain = "negative electrode" + + flux = pybamm.VectorField( + pybamm.PrimaryBroadcast(pybamm.Scalar(1), domain), + pybamm.PrimaryBroadcast(pybamm.Scalar(2), domain), + pybamm.PrimaryBroadcast(pybamm.Scalar(3), domain), + ) + flux_disc = disc.process_symbol(flux) + + t_sol = np.array([0.0, 1.0]) + y_sol = np.ones((submesh.npts, 2)) + comp_casadi = [to_casadi(c, y_sol) for c in flux_disc.components] + model = pybamm.BaseModel() + model._geometry = geometry + solution = pybamm.Solution(t_sol, y_sol, model, {}) + + flux_pv = pybamm.process_variable("flux", [flux_disc], [comp_casadi], solution) + assert isinstance(flux_pv, pybamm.ProcessedVariableVectorFieldUnstructuredFVM) + assert flux_pv.dimensions == 3 + components = flux_pv( + 0.5, x=np.array([0.5]), y=np.array([0.5]), z=np.array([0.5]) + ) + assert len(components) == 3 + for value, component in zip([1.0, 2.0, 3.0], components, strict=True): + np.testing.assert_allclose(component, value, rtol=1e-10) + + def test_2d_domain_with_hole_masks_hole(self): + from pybamm.meshes.meshes import MeshGenerator + from pybamm.meshes.unstructured_submesh import ( + UnstructuredSubMesh, + _make_quad_grid, + ) + + class HoleGenerator(MeshGenerator): + """3x3 quad grid on [0,1]^2 with the centre cell removed.""" + + def __init__(self): + self.submesh_type = UnstructuredSubMesh + self.submesh_params = {} + + def __call__(self, lims, npts): + nodes, elements = _make_quad_grid( + np.linspace(0, 1, 4), np.linspace(0, 1, 4) + ) + centroids = nodes[elements].mean(axis=1) + keep = ~( + np.isclose(centroids[:, 0], 0.5) & np.isclose(centroids[:, 1], 0.5) + ) + sub = UnstructuredSubMesh(nodes, elements[keep]) + sub.detect_box_boundaries() + return sub + + domain = "negative electrode" + x = pybamm.SpatialVariable("x_n", domain=[domain], coord_sys="cartesian") + z = pybamm.SpatialVariable( + "z_2d", domain=[domain], coord_sys="cartesian", direction="tb" + ) + geometry = {domain: {x: {"min": 0.0, "max": 1.0}, z: {"min": 0.0, "max": 1.0}}} + mesh = pybamm.Mesh(geometry, {domain: HoleGenerator()}, {x: 3, z: 3}) + disc = pybamm.Discretisation(mesh, {domain: pybamm.FiniteVolumeUnstructured()}) + var = pybamm.Variable("u", domain=[domain]) + disc.set_variable_slices([var]) + var_disc = disc.process_symbol(var) + + submesh = mesh[domain] + assert submesh.npts == 8 + t_sol = np.array([0.0, 1.0]) + y_sol = 3.0 * np.ones((submesh.npts, 2)) + pv = self._make_pv(var_disc, geometry, t_sol, y_sol) + + in_hole = pv(0.5, x=np.array([0.5]), z=np.array([0.5])) + assert np.isnan(in_hole).all() + in_domain = pv(0.5, x=np.array([1 / 6]), z=np.array([1 / 6])) + np.testing.assert_allclose(in_domain, 3.0, rtol=1e-10) diff --git a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_unstructured.py b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_unstructured.py new file mode 100644 index 0000000000..be3e076f1d --- /dev/null +++ b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_unstructured.py @@ -0,0 +1,2538 @@ +""" +Unit tests for FiniteVolumeUnstructured spatial method. + +Tests cover both 2D (triangle) and 3D (tet) meshes, validating: +- TPFA Laplacian structural properties and conservation +- Green-Gauss gradient on linear fields +- Divergence (adjoint of gradient) +- Mass matrix, integrals, boundary value/flux +- Internal Neumann condition for domain coupling +""" + +import numpy as np +import pytest +from scipy.sparse import coo_matrix as sp_coo +from scipy.sparse import csr_matrix as sp_csr +from scipy.sparse.linalg import spsolve + +import pybamm +from pybamm.meshes.unstructured_submesh import ( + UnstructuredSubMesh, + _hex_grid, + _hex_to_tet, + _make_quad_grid, + _quad_to_tri, + compute_interface_data, +) +from pybamm.spatial_methods.finite_volume_unstructured import ( + FiniteVolumeUnstructured, +) + +# ====================================================================== +# Mesh helpers +# ====================================================================== + + +def _make_2d_mesh(nx=4, nz=4, x_range=(0, 1), z_range=(0, 1)): + x_edges = np.linspace(x_range[0], x_range[1], nx + 1) + z_edges = np.linspace(z_range[0], z_range[1], nz + 1) + nodes, elements = _quad_to_tri(x_edges, z_edges) + submesh = UnstructuredSubMesh(nodes, elements) + submesh.detect_box_boundaries() + return submesh + + +def _make_3d_mesh(nx=3, ny=3, nz=3, x_range=(0, 1), y_range=(0, 1), z_range=(0, 1)): + x_edges = np.linspace(x_range[0], x_range[1], nx + 1) + y_edges = np.linspace(y_range[0], y_range[1], ny + 1) + z_edges = np.linspace(z_range[0], z_range[1], nz + 1) + nodes, elements = _hex_to_tet(x_edges, y_edges, z_edges) + submesh = UnstructuredSubMesh(nodes, elements) + submesh.detect_box_boundaries() + return submesh + + +def _make_quad_mesh(nx=4, nz=4, x_range=(0, 1), z_range=(0, 1)): + """TPFA-orthogonal quadrilateral mesh (exact for linear fields).""" + x_edges = np.linspace(x_range[0], x_range[1], nx + 1) + z_edges = np.linspace(z_range[0], z_range[1], nz + 1) + nodes, elements = _make_quad_grid(x_edges, z_edges) + submesh = UnstructuredSubMesh(nodes, elements) + submesh.detect_box_boundaries() + return submesh + + +def _make_hex_mesh(nx=3, ny=3, nz=3): + """TPFA-orthogonal hexahedral mesh on the unit cube.""" + x_edges = np.linspace(0, 1, nx + 1) + y_edges = np.linspace(0, 1, ny + 1) + z_edges = np.linspace(0, 1, nz + 1) + nodes, elements = _hex_grid(x_edges, y_edges, z_edges) + submesh = UnstructuredSubMesh(nodes, elements) + submesh.detect_box_boundaries() + return submesh + + +def _make_split_2d_meshes(nx_left=3, nx_right=3, nz=3): + """Create two adjacent 2D meshes for interface testing.""" + left = _make_2d_mesh(nx_left, nz, x_range=(0, 0.5)) + right = _make_2d_mesh(nx_right, nz, x_range=(0.5, 1.0)) + compute_interface_data(left, right, left_name="left", right_name="right") + return left, right + + +def _get_internal_cells(mesh): + """Return indices of cells that do not touch any boundary face.""" + bnd_cells = set() + for indices in mesh.boundary_faces.values(): + for fi in indices: + bnd_cells.add(mesh.face_owner[fi]) + return [i for i in range(mesh.npts) if i not in bnd_cells] + + +class _MeshMap(dict): + """Minimal Mesh-like mapping that accepts PyBaMM domain lists.""" + + def __getitem__(self, key): + if isinstance(key, list): + key = tuple(key) + elif isinstance(key, str): + key = (key,) + return super().__getitem__(key) + + +def _method_with_mesh(mesh, **auxiliary_meshes): + meshes = {("test",): mesh} + meshes.update({(name,): value for name, value in auxiliary_meshes.items()}) + method = FiniteVolumeUnstructured() + method._mesh = _MeshMap(meshes) + return method + + +# ====================================================================== +# Tests: TPFA Laplacian +# ====================================================================== + + +class TestTPFALaplacian: + def test_tpfa_matrix_shape_2d(self): + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + assert L.shape == (mesh.npts, mesh.npts) + + def test_tpfa_matrix_shape_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + assert L.shape == (mesh.npts, mesh.npts) + + def test_tpfa_stiffness_symmetry_2d(self): + """The raw stiffness matrix K (before volume scaling) should be symmetric.""" + mesh = _make_2d_mesh(5, 5) + n = mesh.npts + n_int = mesh.n_internal_faces + + owner = mesh.face_owner[:n_int] + neighbor = mesh.face_neighbor[:n_int] + areas = mesh.face_areas[:n_int] + c_owner = mesh.cell_centroids[owner] + c_neighbor = mesh.cell_centroids[neighbor] + dist = np.linalg.norm(c_neighbor - c_owner, axis=1) + coeff = areas / dist + + rows = np.concatenate([owner, neighbor, owner, neighbor]) + cols = np.concatenate([neighbor, owner, owner, neighbor]) + data = np.concatenate([coeff, coeff, -coeff, -coeff]) + K = sp_csr(sp_coo((data, (rows, cols)), shape=(n, n))) + + diff = K - K.T + assert abs(diff).max() < 1e-12 + + def test_tpfa_conservation_2d(self): + """Weighted sum of L@u over all cells = 0 (internal flux conservation).""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = mesh.cell_centroids[:, 0] ** 2 + Lu = L @ u + total = np.sum(Lu * mesh.cell_volumes) + np.testing.assert_allclose(total, 0.0, atol=1e-10) + + def test_tpfa_conservation_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = mesh.cell_centroids[:, 0] ** 2 + Lu = L @ u + total = np.sum(Lu * mesh.cell_volumes) + np.testing.assert_allclose(total, 0.0, atol=1e-10) + + def test_tpfa_constant_field_2d(self): + """Laplacian of constant = 0.""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = np.ones(mesh.npts) * 7.0 + np.testing.assert_allclose(L @ u, 0.0, atol=1e-12) + + def test_tpfa_constant_field_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = np.ones(mesh.npts) * 7.0 + np.testing.assert_allclose(L @ u, 0.0, atol=1e-12) + + def test_tpfa_negative_diagonal_2d(self): + """Diagonal entries of TPFA matrix should be non-positive.""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + diag = L.diagonal() + assert np.all(diag <= 1e-15) + + +# ====================================================================== +# Tests: Neumann sign convention (PyBaMM coordinate-direction values) +# ====================================================================== + + +class TestNeumannSignConvention: + """Named sides take coordinate-direction derivatives (matching + FiniteVolume/FiniteVolume2D), so ``u = x`` needs value +1 on *both* + left and right, not the outward-normal ±1.""" + + def test_laplacian_neumann_left_right_2d(self): + mesh = _make_quad_mesh(4, 4) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + u = pybamm.Vector(mesh.cell_centroids[:, 0], domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(1), "Neumann"), + "right": (pybamm.Scalar(1), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + result = method.laplacian(variable, u, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + def test_laplacian_neumann_top_bottom_2d(self): + mesh = _make_quad_mesh(4, 4) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + u = pybamm.Vector(mesh.cell_centroids[:, 1], domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(1), "Neumann"), + "bottom": (pybamm.Scalar(1), "Neumann"), + } + } + result = method.laplacian(variable, u, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + def test_laplacian_neumann_front_back_3d(self): + mesh = _make_hex_mesh(3, 3, 3) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + u = pybamm.Vector(mesh.cell_centroids[:, 1], domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "front": (pybamm.Scalar(1), "Neumann"), + "back": (pybamm.Scalar(1), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + result = method.laplacian(variable, u, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + def test_gradient_neumann_left_2d(self): + mesh = _make_quad_mesh(4, 4) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + u = pybamm.StateVector(slice(0, mesh.npts), domains={"primary": ["test"]}) + y = mesh.cell_centroids[:, 0] + bcs = { + variable: { + "left": (pybamm.Scalar(1), "Neumann"), + "right": (pybamm.Scalar(1), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + grad = method.gradient(variable, u, bcs) + np.testing.assert_allclose(grad.components[0].evaluate(y=y), 1, atol=1e-10) + np.testing.assert_allclose(grad.components[1].evaluate(y=y), 0, atol=1e-10) + + def test_div_D_grad_neumann_left_2d(self): + mesh = _make_quad_mesh(4, 4) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + u = pybamm.Vector(mesh.cell_centroids[:, 0], domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(1), "Neumann"), + "right": (pybamm.Scalar(1), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + result = method.div_D_grad(div_symbol, variable, pybamm.Scalar(2), u, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + def test_custom_tags_use_outward_normal(self): + # Rename the axis buckets to custom tags: values are then + # outward-normal derivatives, so u = x needs -1 on the left tag. + mesh = _make_quad_mesh(4, 4) + mesh.boundary_faces = { + f"tag_{name}": faces for name, faces in mesh.boundary_faces.items() + } + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + u = pybamm.Vector(mesh.cell_centroids[:, 0], domain="test") + bcs = { + variable: { + "tag_left": (pybamm.Scalar(-1), "Neumann"), + "tag_right": (pybamm.Scalar(1), "Neumann"), + "tag_top": (pybamm.Scalar(0), "Neumann"), + "tag_bottom": (pybamm.Scalar(0), "Neumann"), + } + } + result = method.laplacian(variable, u, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + +# ====================================================================== +# Tests: operator caching +# ====================================================================== + + +class TestOperatorCaching: + def test_operators_cached_and_invalidated_on_reordering(self): + mesh = _make_2d_mesh(4, 4) + fvu = FiniteVolumeUnstructured() + assert fvu._tpfa_matrix(mesh) is fvu._tpfa_matrix(mesh) + assert fvu._green_gauss_matrices(mesh) is fvu._green_gauss_matrices(mesh) + assert fvu._divergence_matrices(mesh) is fvu._divergence_matrices(mesh) + assert fvu._div_D_grad_matrices(mesh) is fvu._div_D_grad_matrices(mesh) + + laplacian_before = fvu._tpfa_matrix(mesh) + mesh.optimize_ordering() + assert fvu._tpfa_matrix(mesh) is not laplacian_before + + def test_bc_application_does_not_mutate_cached_operators(self): + mesh = _make_quad_mesh(3, 3) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(1), "Dirichlet"), + "right": (pybamm.Scalar(0), "Dirichlet"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + laplacian_cached = method._tpfa_matrix(mesh).copy() + gauss_cached = method._green_gauss_matrices(mesh)[0].copy() + + method.laplacian(variable, values, bcs) + method.gradient(variable, values, bcs) + + assert (method._tpfa_matrix(mesh) - laplacian_cached).nnz == 0 + assert (method._green_gauss_matrices(mesh)[0] - gauss_cached).nnz == 0 + + +# ====================================================================== +# Tests: auxiliary domains (secondary/tertiary repeats) +# ====================================================================== + + +class TestAuxiliaryDomains: + def _bcs(self, variable): + return { + variable: { + "left": (pybamm.Scalar(1), "Dirichlet"), + "right": (pybamm.Scalar(2), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(3), "Dirichlet"), + } + } + + def test_laplacian_with_secondary_domain(self): + mesh = _make_quad_mesh(2, 2) + aux = _make_quad_mesh(1, 3) + method = _method_with_mesh(mesh, aux=aux) + cell_values = mesh.cell_centroids[:, 0] ** 2 + + variable = pybamm.Variable("u", domain="test") + single = method.laplacian( + variable, + pybamm.Vector(cell_values, domain="test"), + self._bcs(variable), + ) + + domains = {"primary": ["test"], "secondary": ["aux"]} + repeated_var = pybamm.Variable("u rep", domains=domains) + repeated = method.laplacian( + repeated_var, + pybamm.Vector(np.tile(cell_values, aux.npts), domains=domains), + self._bcs(repeated_var), + ) + np.testing.assert_allclose( + repeated.evaluate()[:, 0], + np.tile(single.evaluate()[:, 0], aux.npts), + atol=1e-12, + ) + + def test_gradient_with_secondary_domain(self): + mesh = _make_quad_mesh(2, 2) + aux = _make_quad_mesh(1, 3) + method = _method_with_mesh(mesh, aux=aux) + cell_values = mesh.cell_centroids[:, 0] ** 2 + + variable = pybamm.Variable("u", domain="test") + single = method.gradient( + variable, + pybamm.Vector(cell_values, domain="test"), + self._bcs(variable), + ) + + domains = {"primary": ["test"], "secondary": ["aux"]} + repeated_var = pybamm.Variable("u rep", domains=domains) + repeated = method.gradient( + repeated_var, + pybamm.Vector(np.tile(cell_values, aux.npts), domains=domains), + self._bcs(repeated_var), + ) + for single_comp, repeated_comp in zip( + single.components, repeated.components, strict=True + ): + np.testing.assert_allclose( + repeated_comp.evaluate()[:, 0], + np.tile(single_comp.evaluate()[:, 0], aux.npts), + atol=1e-12, + ) + + def test_tertiary_broadcast_size(self): + mesh = _make_quad_mesh(2, 2) + sec = _make_quad_mesh(1, 2) + ter = _make_quad_mesh(1, 5) + method = _method_with_mesh(mesh, sec=sec, ter=ter) + + child_size = mesh.npts * sec.npts + child = pybamm.Vector( + np.arange(child_size), + domains={"primary": ["test"], "secondary": ["sec"]}, + ) + domains = { + "primary": ["test"], + "secondary": ["sec"], + "tertiary": ["ter"], + } + out = method.broadcast(child, domains, "tertiary to nodes") + assert out.shape_for_testing == (child_size * ter.npts, 1) + np.testing.assert_array_equal( + out.evaluate()[:, 0], np.tile(np.arange(child_size), ter.npts) + ) + + +# ====================================================================== +# Tests: BC tag / type validation +# ====================================================================== + + +class TestBCValidation: + def _setup(self): + mesh = _make_quad_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + return mesh, method, variable, values + + def test_laplacian_unknown_side_raises(self): + _, method, variable, values = self._setup() + bcs = {variable: {"weft": (pybamm.Scalar(0), "Neumann")}} + with pytest.raises(pybamm.DiscretisationError, match="weft"): + method.laplacian(variable, values, bcs) + + def test_laplacian_unknown_bc_type_raises(self): + _, method, variable, values = self._setup() + bcs = {variable: {"left": (pybamm.Scalar(0), "Robin")}} + with pytest.raises(pybamm.DiscretisationError, match="Robin"): + method.laplacian(variable, values, bcs) + + def test_gradient_unknown_side_raises(self): + _, method, variable, values = self._setup() + bcs = {variable: {"weft": (pybamm.Scalar(0), "Neumann")}} + with pytest.raises(pybamm.DiscretisationError, match="weft"): + method.gradient(variable, values, bcs) + + def test_gradient_unknown_bc_type_raises(self): + _, method, variable, values = self._setup() + bcs = {variable: {"left": (pybamm.Scalar(0), "Dirchlet")}} + with pytest.raises(pybamm.DiscretisationError, match="Dirchlet"): + method.gradient(variable, values, bcs) + + def test_div_D_grad_unknown_side_raises(self): + _, method, variable, values = self._setup() + div_symbol = pybamm.Variable("div", domain="test") + bcs = {variable: {"weft": (pybamm.Scalar(0), "Neumann")}} + with pytest.raises(pybamm.DiscretisationError, match="weft"): + method.div_D_grad(div_symbol, variable, pybamm.Scalar(1), values, bcs) + + def test_boundary_value_unknown_side_raises(self): + _, method, variable, values = self._setup() + symbol = pybamm.BoundaryValue(variable, "missing") + with pytest.raises(pybamm.DiscretisationError, match="missing"): + method.boundary_value_or_flux(symbol, values) + + def test_boundary_integral_unknown_region_raises(self): + _, method, variable, values = self._setup() + with pytest.raises(pybamm.DiscretisationError, match="missing"): + method.boundary_integral(variable, values, "missing") + + def test_boundary_integral_entire(self): + # integral of u = 1 over the whole boundary = perimeter of unit square + mesh, method, variable, _ = self._setup() + ones = pybamm.Vector(np.ones(mesh.npts), domain="test") + result = method.boundary_integral(variable, ones, "entire") + np.testing.assert_allclose(result.evaluate().sum(), 4.0, atol=1e-12) + + def test_boundary_integral_entire_excludes_interface_faces(self): + left, _right = _make_split_2d_meshes() + # emulate interface discovery: move left mesh's right faces to an + # iface bucket + left.boundary_faces["iface_right"] = left.boundary_faces.pop("right") + method = _method_with_mesh(left) + variable = pybamm.Variable("u", domain="test") + ones = pybamm.Vector(np.ones(left.npts), domain="test") + result = method.boundary_integral(variable, ones, "entire") + # perimeter of [0, 0.5] x [0, 1] minus the shared edge (length 1) + np.testing.assert_allclose(result.evaluate().sum(), 2.0, atol=1e-12) + + def test_deleted_bucket_message_mentions_interface(self): + mesh, method, variable, values = self._setup() + mesh.boundary_faces["iface_other"] = mesh.boundary_faces.pop("right") + bcs = {variable: {"right": (pybamm.Scalar(0), "Dirichlet")}} + with pytest.raises(pybamm.DiscretisationError, match="interface"): + method.laplacian(variable, values, bcs) + + +# ====================================================================== +# Tests: Green-Gauss Gradient +# ====================================================================== + + +class TestGreenGaussGradient: + def test_gradient_constant_field_2d(self): + """Gradient of constant = 0 everywhere.""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = np.ones(mesh.npts) * 3.14 + for k in range(mesh.dimension): + np.testing.assert_allclose(G[k] @ u, 0.0, atol=1e-12) + + def test_gradient_constant_field_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = np.ones(mesh.npts) * 3.14 + for k in range(mesh.dimension): + np.testing.assert_allclose(G[k] @ u, 0.0, atol=1e-12) + + def test_gradient_linear_x_2d(self): + """Gradient of u = x should be [1, 0] on internal cells.""" + mesh = _make_2d_mesh(8, 8) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 0] + internal = _get_internal_cells(mesh) + + if internal: + np.testing.assert_allclose((G[0] @ u)[internal], 1.0, atol=1e-10) + np.testing.assert_allclose((G[1] @ u)[internal], 0.0, atol=1e-10) + + def test_gradient_linear_z_2d(self): + """Gradient of u = z should be [0, 1] on internal cells.""" + mesh = _make_2d_mesh(8, 8) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 1] + internal = _get_internal_cells(mesh) + + if internal: + np.testing.assert_allclose((G[0] @ u)[internal], 0.0, atol=1e-10) + np.testing.assert_allclose((G[1] @ u)[internal], 1.0, atol=1e-10) + + def test_gradient_linear_x_3d(self): + """Gradient of u = x on 3D tet mesh. + + On non-orthogonal tet meshes from hex splitting, the Green-Gauss + gradient with distance-weighted interpolation has O(h) error. + Boundary cells contribute a bias from zeroth-order face + extrapolation. We verify the mean is within 15% and that + internal cells are accurate. + """ + mesh = _make_3d_mesh(4, 4, 4) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 0] + grad_x = G[0] @ u + + mean_grad_x = np.sum(grad_x * mesh.cell_volumes) / mesh.cell_volumes.sum() + np.testing.assert_allclose(mean_grad_x, 1.0, atol=0.15) + + def test_gradient_linear_combo_2d(self): + """Gradient of u = 2x + 3z should be [2, 3].""" + mesh = _make_2d_mesh(8, 8) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = 2 * mesh.cell_centroids[:, 0] + 3 * mesh.cell_centroids[:, 1] + internal = _get_internal_cells(mesh) + + if internal: + np.testing.assert_allclose((G[0] @ u)[internal], 2.0, atol=1e-10) + np.testing.assert_allclose((G[1] @ u)[internal], 3.0, atol=1e-10) + + +# ====================================================================== +# Tests: Divergence +# ====================================================================== + + +class TestDivergence: + def test_divergence_matrices_shape_2d(self): + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + D = fvu._divergence_matrices(mesh) + assert len(D) == 2 + assert D[0].shape == (mesh.npts, mesh.npts) + + def test_divergence_matrices_shape_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + D = fvu._divergence_matrices(mesh) + assert len(D) == 3 + assert D[0].shape == (mesh.npts, mesh.npts) + + def test_divergence_constant_vector_field_2d(self): + """Divergence of a constant vector field = 0 on internal cells.""" + mesh = _make_2d_mesh(6, 6) + fvu = FiniteVolumeUnstructured() + D = fvu._divergence_matrices(mesh) + + Fx = np.ones(mesh.npts) * 2.0 + Fz = np.ones(mesh.npts) * 3.0 + div = D[0] @ Fx + D[1] @ Fz + + internal = _get_internal_cells(mesh) + if internal: + np.testing.assert_allclose(div[internal], 0.0, atol=1e-10) + + def test_divergence_constant_vector_field_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + D = fvu._divergence_matrices(mesh) + + Fx = np.ones(mesh.npts) * 2.0 + Fy = np.ones(mesh.npts) * 3.0 + Fz = np.ones(mesh.npts) * 4.0 + div = D[0] @ Fx + D[1] @ Fy + D[2] @ Fz + + internal = _get_internal_cells(mesh) + if internal: + np.testing.assert_allclose(div[internal], 0.0, atol=1e-10) + + +# ====================================================================== +# Tests: Mass matrix (cell volumes) +# ====================================================================== + + +class TestMassMatrix: + def test_volume_sum_2d(self): + """Sum of cell volumes = domain area.""" + mesh = _make_2d_mesh(5, 5) + np.testing.assert_allclose(mesh.cell_volumes.sum(), 1.0, atol=1e-12) + + def test_volume_sum_3d(self): + """Sum of cell volumes = domain volume.""" + mesh = _make_3d_mesh(3, 3, 3) + np.testing.assert_allclose(mesh.cell_volumes.sum(), 1.0, atol=1e-12) + + def test_volumes_positive_2d(self): + mesh = _make_2d_mesh(5, 5) + assert np.all(mesh.cell_volumes > 0) + + def test_volumes_positive_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + assert np.all(mesh.cell_volumes > 0) + + def test_volume_sum_rectangle(self): + """Non-square domain: [0,2] x [0,0.5] should have area 1.0.""" + mesh = _make_2d_mesh(6, 4, x_range=(0, 2), z_range=(0, 0.5)) + np.testing.assert_allclose(mesh.cell_volumes.sum(), 1.0, atol=1e-12) + + +# ====================================================================== +# Tests: Integral +# ====================================================================== + + +class TestIntegral: + @pytest.mark.parametrize( + ("make_mesh", "expected_linear"), + [(lambda: _make_2d_mesh(10, 10), 0.5), (lambda: _make_3d_mesh(4, 4, 4), 0.5)], + ids=["2d", "3d"], + ) + def test_definite_integral_matrix(self, make_mesh, expected_linear): + """Integral of 1 over the unit box is 1; of ``x`` is 0.5.""" + mesh = make_mesh() + method = _method_with_mesh(mesh) + child = pybamm.Variable("u", domain="test") + mat = method.definite_integral_matrix(child) + assert isinstance(mat, pybamm.Matrix) + assert mat.shape == (1, mesh.npts) + np.testing.assert_allclose(mat.entries @ np.ones(mesh.npts), 1.0, atol=1e-12) + np.testing.assert_allclose( + mat.entries @ mesh.cell_centroids[:, 0], expected_linear, atol=0.01 + ) + + def test_definite_integral_matrix_column(self): + mesh = _make_2d_mesh(3, 3) + method = _method_with_mesh(mesh) + child = pybamm.Variable("u", domain="test") + column = method.definite_integral_matrix(child, vector_type="column") + assert column.shape == (mesh.npts, 1) + np.testing.assert_array_equal(column.entries.toarray()[:, 0], mesh.cell_volumes) + with pytest.raises(pybamm.DiscretisationError, match="vector_type"): + method.definite_integral_matrix(child, vector_type="diagonal") + + def test_non_primary_integration_dimension_raises(self): + mesh = _make_2d_mesh(2, 2) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + child = pybamm.Variable( + "u", domains={"primary": ["test"], "secondary": ["aux"]} + ) + values = pybamm.Vector(np.ones(mesh.npts * aux.npts), domains=child.domains) + with pytest.raises(NotImplementedError, match="secondary"): + method.integral(child, values, "secondary") + with pytest.raises(NotImplementedError, match="secondary"): + method.definite_integral_matrix(child, integration_dimension="secondary") + + def test_definite_integral_vector_through_discretisation(self): + """``DefiniteIntegralVector`` must come back as a ``pybamm.Matrix`` so + ``process_symbol`` can shape-check it, in both orientations.""" + x = pybamm.SpatialVariable("x_n", domain=["negative electrode"]) + z = pybamm.SpatialVariable( + "z_2d", domain=["negative electrode"], coord_sys="cartesian", direction="tb" + ) + geometry = { + "negative electrode": {x: {"min": 0, "max": 1}, z: {"min": 0, "max": 1}} + } + generator = pybamm.meshes.unstructured_submesh.UnstructuredMeshGenerator() + mesh = pybamm.Mesh(geometry, {"negative electrode": generator}, {x: 3, z: 3}) + disc = pybamm.Discretisation( + mesh, {"negative electrode": FiniteVolumeUnstructured()} + ) + var = pybamm.Variable("var", domain="negative electrode") + disc.set_variable_slices([var]) + npts = mesh["negative electrode"].npts + + row = disc.process_symbol(pybamm.DefiniteIntegralVector(var)) + assert row.shape == (1, npts) + column = disc.process_symbol( + pybamm.DefiniteIntegralVector(var, vector_type="column") + ) + assert column.shape == (npts, 1) + np.testing.assert_allclose(row.evaluate() @ np.ones(npts), 1.0, atol=1e-12) + + +# ====================================================================== +# Tests: Boundary value / flux +# ====================================================================== + + +class TestBoundaryValue: + def test_boundary_faces_exist_2d(self): + mesh = _make_2d_mesh(5, 5) + assert "left" in mesh.boundary_faces + assert "right" in mesh.boundary_faces + assert "bottom" in mesh.boundary_faces + assert "top" in mesh.boundary_faces + + for tag in ["left", "right", "bottom", "top"]: + assert len(mesh.boundary_faces[tag]) > 0 + + def test_boundary_faces_exist_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + assert "left" in mesh.boundary_faces + assert "right" in mesh.boundary_faces + + def test_left_boundary_x_zero_2d(self): + """Left boundary face centroids should have x ≈ 0.""" + mesh = _make_2d_mesh(5, 5) + left_centroids = mesh.face_centroids[mesh.boundary_faces["left"]] + np.testing.assert_allclose(left_centroids[:, 0], 0.0, atol=1e-14) + + def test_right_boundary_x_one_2d(self): + """Right boundary face centroids should have x ≈ 1.""" + mesh = _make_2d_mesh(5, 5) + right_centroids = mesh.face_centroids[mesh.boundary_faces["right"]] + np.testing.assert_allclose(right_centroids[:, 0], 1.0, atol=1e-14) + + +# ====================================================================== +# Tests: Interface / internal_neumann_condition +# ====================================================================== + + +class TestInternalNeumann: + def test_interface_data_exists(self): + left, right = _make_split_2d_meshes(3, 3, 3) + assert len(left.interface_data) > 0 or len(right.interface_data) > 0 + + def test_interface_face_count(self): + """Number of interface faces should equal the number of z-boundary faces.""" + left, _right = _make_split_2d_meshes(4, 4, 4) + interface = next(iter(left.interface_data.values())) + assert len(interface["left_cells"]) > 0 + assert len(interface["right_cells"]) > 0 + assert len(interface["left_cells"]) == len(interface["right_cells"]) + + def test_interface_uniform_field(self): + """Interface gradient of uniform field = 0.""" + left, right = _make_split_2d_meshes(4, 4, 4) + interface = next(iter(left.interface_data.values())) + + left_vals = np.ones(left.npts) * 5.0 + right_vals = np.ones(right.npts) * 5.0 + + inv_dx = 1.0 / interface["cell_distances"] + grad = inv_dx * ( + right_vals[interface["right_cells"]] - left_vals[interface["left_cells"]] + ) + np.testing.assert_allclose(grad, 0.0, atol=1e-12) + + def test_interface_gradient_positive_for_increasing_x(self): + """For u = x, interface gradient should be positive.""" + left, right = _make_split_2d_meshes(4, 4, 4) + interface = next(iter(left.interface_data.values())) + + left_vals = left.cell_centroids[:, 0] + right_vals = right.cell_centroids[:, 0] + + inv_dx = 1.0 / interface["cell_distances"] + grad = inv_dx * ( + right_vals[interface["right_cells"]] - left_vals[interface["left_cells"]] + ) + assert np.all(grad > 0), "Gradient should be positive for u = x" + + def test_interface_cell_distances_positive(self): + left, _right = _make_split_2d_meshes(4, 4, 4) + interface = next(iter(left.interface_data.values())) + assert np.all(interface["cell_distances"] > 0) + + def test_interface_face_areas_positive(self): + left, _right = _make_split_2d_meshes(4, 4, 4) + interface = next(iter(left.interface_data.values())) + assert np.all(interface["face_areas"] > 0) + + +# ====================================================================== +# Tests: Conservation / divergence theorem +# ====================================================================== + + +class TestConservation: + def test_tpfa_conservation_2d(self): + """Total internal flux = 0 (conservation of Laplacian).""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = mesh.cell_centroids[:, 0] ** 2 + Lu = L @ u + total = np.sum(Lu * mesh.cell_volumes) + np.testing.assert_allclose(total, 0.0, atol=1e-10) + + def test_divergence_theorem_volume_weighted_2d(self): + """ + For F = (x, z): div(F) = 2. + Volume-weighted integral of div(F) should approach 2 * area. + The Green-Gauss divergence has boundary-cell errors, so we use + a generous tolerance. + """ + mesh = _make_2d_mesh(10, 10) + fvu = FiniteVolumeUnstructured() + D = fvu._divergence_matrices(mesh) + + Fx = mesh.cell_centroids[:, 0] + Fz = mesh.cell_centroids[:, 1] + div_F = D[0] @ Fx + D[1] @ Fz + + vol_integral = np.sum(div_F * mesh.cell_volumes) + np.testing.assert_allclose(vol_integral, 2.0, atol=0.25) + + +# ====================================================================== +# Tests: Gradient squared +# ====================================================================== + + +class TestGradientSquared: + def test_gradient_squared_linear_x_2d(self): + """|grad(x)|^2 ≈ 1 on internal cells.""" + mesh = _make_2d_mesh(8, 8) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 0] + grad_sq = sum((G[k] @ u) ** 2 for k in range(mesh.dimension)) + + internal = _get_internal_cells(mesh) + if internal: + np.testing.assert_allclose(grad_sq[internal], 1.0, atol=1e-10) + + def test_gradient_squared_constant_2d(self): + """|grad(const)|^2 = 0.""" + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = np.ones(mesh.npts) * 42.0 + grad_sq = sum((G[k] @ u) ** 2 for k in range(mesh.dimension)) + np.testing.assert_allclose(grad_sq, 0.0, atol=1e-20) + + +# ====================================================================== +# Tests: Not implemented operators +# ====================================================================== + + +class TestNotImplemented: + def test_indefinite_integral_raises(self): + fvu = FiniteVolumeUnstructured() + with pytest.raises(NotImplementedError, match="Indefinite integral"): + fvu.indefinite_integral(None, None, None) + + def test_delta_function_raises(self): + fvu = FiniteVolumeUnstructured() + with pytest.raises(NotImplementedError, match="Delta function"): + fvu.delta_function(None, None) + + +# ====================================================================== +# Tests: 3D specific +# ====================================================================== + + +class Test3D: + def test_gradient_mean_accuracy_3d(self): + """Volume-weighted mean gradient of u = x should be ~1. + + Boundary cells bias the mean via zeroth-order face extrapolation; + tolerance of 0.15 is appropriate for a 4^3 tet mesh. + """ + mesh = _make_3d_mesh(4, 4, 4) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 0] + vol = mesh.cell_volumes + total_vol = vol.sum() + + mean_gx = np.sum((G[0] @ u) * vol) / total_vol + mean_gy = np.sum((G[1] @ u) * vol) / total_vol + mean_gz = np.sum((G[2] @ u) * vol) / total_vol + + np.testing.assert_allclose(mean_gx, 1.0, atol=0.15) + np.testing.assert_allclose(mean_gy, 0.0, atol=0.15) + np.testing.assert_allclose(mean_gz, 0.0, atol=0.15) + + def test_gradient_y_mean_accuracy_3d(self): + """Volume-weighted mean gradient of u = y should be ~[0,1,0].""" + mesh = _make_3d_mesh(4, 4, 4) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 1] + vol = mesh.cell_volumes + total_vol = vol.sum() + + mean_gx = np.sum((G[0] @ u) * vol) / total_vol + mean_gy = np.sum((G[1] @ u) * vol) / total_vol + mean_gz = np.sum((G[2] @ u) * vol) / total_vol + + np.testing.assert_allclose(mean_gx, 0.0, atol=0.15) + np.testing.assert_allclose(mean_gy, 1.0, atol=0.15) + np.testing.assert_allclose(mean_gz, 0.0, atol=0.15) + + def test_gradient_z_mean_accuracy_3d(self): + """Volume-weighted mean gradient of u = z should be ~[0,0,1].""" + mesh = _make_3d_mesh(4, 4, 4) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + + u = mesh.cell_centroids[:, 2] + vol = mesh.cell_volumes + total_vol = vol.sum() + + mean_gx = np.sum((G[0] @ u) * vol) / total_vol + mean_gy = np.sum((G[1] @ u) * vol) / total_vol + mean_gz = np.sum((G[2] @ u) * vol) / total_vol + + np.testing.assert_allclose(mean_gx, 0.0, atol=0.15) + np.testing.assert_allclose(mean_gy, 0.0, atol=0.15) + np.testing.assert_allclose(mean_gz, 1.0, atol=0.15) + + def test_tpfa_constant_3d(self): + """Laplacian of constant = 0.""" + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + u = np.ones(mesh.npts) * 7.0 + np.testing.assert_allclose(L @ u, 0.0, atol=1e-12) + + def test_divergence_conservation_3d(self): + """Weighted Laplacian sum = 0 (conservation).""" + mesh = _make_3d_mesh(3, 3, 3) + fvu = FiniteVolumeUnstructured() + L = fvu._tpfa_matrix(mesh) + + u = mesh.cell_centroids[:, 0] ** 2 + Lu = L @ u + total = np.sum(Lu * mesh.cell_volumes) + np.testing.assert_allclose(total, 0.0, atol=1e-10) + + +# ====================================================================== +# Tests: Miscellaneous +# ====================================================================== + + +class TestMisc: + def test_face_count_2d(self): + """Total faces = internal + boundary.""" + mesh = _make_2d_mesh(4, 4) + n_total = len(mesh.faces) + n_bnd = sum(len(v) for v in mesh.boundary_faces.values()) + assert n_total == mesh.n_internal_faces + n_bnd + + def test_face_count_3d(self): + mesh = _make_3d_mesh(2, 2, 2) + n_total = len(mesh.faces) + n_bnd = sum(len(v) for v in mesh.boundary_faces.values()) + assert n_total == mesh.n_internal_faces + n_bnd + + def test_gradient_divergence_duality_2d(self): + """ + For the Green-Gauss method, gradient and divergence matrices are + structurally related (same interpolation weights, same normals). + Test that G_k and D_k are identical. + """ + mesh = _make_2d_mesh(5, 5) + fvu = FiniteVolumeUnstructured() + G = fvu._green_gauss_matrices(mesh) + D = fvu._divergence_matrices(mesh) + + for k in range(mesh.dimension): + diff = G[k] - D[k] + assert abs(diff).max() < 1e-14 + + def test_constructor_default_options(self): + fvu = FiniteVolumeUnstructured() + assert fvu.options is not None + assert "extrapolation" in fvu.options + + +class TestFiniteVolumeUnstructuredBehavior: + def test_build_discovers_interfaces_and_ignores_other_meshes(self): + left = _make_2d_mesh(2, 2, x_range=(0, 0.5)) + right = _make_2d_mesh(2, 2, x_range=(0.5, 1)) + structured = pybamm.SubMesh1D(np.array([0, 1]), "cartesian") + meshes = _MeshMap( + {("left",): left, ("right",): right, ("structured",): structured} + ) + + method = FiniteVolumeUnstructured() + method.build(meshes) + + assert right in [data["other_mesh"] for data in left.interface_data.values()] + assert left.npts_for_broadcast_to_nodes == left.npts + assert structured.npts_for_broadcast_to_nodes == structured.npts + + def test_build_warns_on_untagged_mesh(self, caplog): + import logging + + untagged = _make_2d_mesh(2, 2) + untagged.boundary_faces = {} + tagged = _make_2d_mesh(2, 2) + method = FiniteVolumeUnstructured() + with caplog.at_level(logging.WARNING): + method.build(_MeshMap({("untagged",): untagged, ("tagged",): tagged})) + assert "no boundary tags" in caplog.text + assert "'untagged'" in caplog.text + assert "'tagged'" not in caplog.text + + def test_interface_matching_edge_cases(self): + empty = _make_2d_mesh(1, 1) + empty.boundary_faces = {} + other = _make_2d_mesh(1, 1) + a_idx, b_idx, matched = FiniteVolumeUnstructured._interface_face_match( + empty, other + ) + assert not matched + assert a_idx.size == b_idx.size == 0 + + mesh_3d = _make_3d_mesh(1, 1, 1) + assert not FiniteVolumeUnstructured._interface_face_match(other, mesh_3d)[2] + + distant = _make_2d_mesh(1, 1, x_range=(2, 3)) + assert not FiniteVolumeUnstructured._interface_face_match(other, distant)[2] + + def test_compute_pair_interface_success_and_noops(self): + left = _make_2d_mesh(2, 2, x_range=(0, 0.5)) + right = _make_2d_mesh(2, 2, x_range=(0.5, 1)) + method = FiniteVolumeUnstructured() + + assert method._compute_pair_interface(left, right, "left", "right") + assert "iface_right" in left.boundary_faces + assert "iface_left" in right.boundary_faces + assert method._compute_pair_interface(left, right, "left", "right") is False + + far = _make_2d_mesh(1, 1, x_range=(2, 3)) + assert method._compute_pair_interface(left, far, "left", "far") is False + + shared = _make_2d_mesh(1, 1) + method._auto_compute_all_interfaces( + _MeshMap({("first",): shared, ("alias",): shared}) + ) + + def test_spatial_variable_directions_and_auxiliary_repeats(self): + mesh = _make_3d_mesh(1, 1, 1) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + domains = {"primary": ["test"], "secondary": ["aux"]} + + for name, direction, column in [ + ("x", None, 0), + ("y", None, 1), + ("z", None, 2), + ("x_n", None, 0), + ("s", "lr", 0), + ("s", "tb", 2), + ("s", "fb", 1), + ]: + symbol = pybamm.SpatialVariable(name, domains=domains, direction=direction) + actual = method.spatial_variable(symbol).evaluate().reshape(-1) + expected = np.tile(mesh.cell_centroids[:, column], aux.npts) + np.testing.assert_allclose(actual, expected) + + # ambiguous names and unknown directions raise instead of guessing x + for name, direction in [("r", None), ("zeta", None), ("s", "unknown")]: + symbol = pybamm.SpatialVariable(name, domains=domains, direction=direction) + with pytest.raises(pybamm.DomainError): + method.spatial_variable(symbol) + + def test_spatial_variable_2d_rejects_y_and_fb(self): + mesh = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh) + z = pybamm.SpatialVariable("z_2d", domain="test", direction="tb") + np.testing.assert_allclose( + method.spatial_variable(z).evaluate().reshape(-1), + mesh.cell_centroids[:, 1], + ) + # 2D meshes are x-z: y names and the front-back direction don't exist + for name, direction in [("y", None), ("s", "fb")]: + symbol = pybamm.SpatialVariable(name, domain="test", direction=direction) + with pytest.raises(pybamm.DomainError): + method.spatial_variable(symbol) + + def test_broadcast_variants(self): + mesh = _make_2d_mesh(1, 1) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + primary = {"primary": ["test"], "secondary": []} + + scalar_primary = method.broadcast(pybamm.Scalar(2), primary, "primary to nodes") + np.testing.assert_array_equal( + scalar_primary.evaluate()[:, 0], np.full(mesh.npts, 2) + ) + + vector_primary = method.broadcast( + pybamm.Vector([2, 3]), primary, "primary to nodes" + ) + np.testing.assert_array_equal( + vector_primary.evaluate()[:, 0], np.repeat([2, 3], mesh.npts) + ) + + full_domains = {"primary": ["test"], "secondary": ["aux"]} + full = method.broadcast(pybamm.Scalar(4), full_domains, "full to nodes") + np.testing.assert_array_equal( + full.evaluate()[:, 0], np.full(mesh.npts * aux.npts, 4) + ) + + secondary_child = pybamm.Vector([1, 2], domain="test") + secondary = method.broadcast(secondary_child, primary, "secondary to nodes") + np.testing.assert_array_equal(secondary.evaluate(), secondary_child.evaluate()) + assert secondary.domain == primary["primary"] + assert secondary.domains["secondary"] == primary["secondary"] + + def test_broadcast_does_not_mutate_simplified_child(self): + mesh = pybamm.SubMesh1D(np.array([0, 1]), "cartesian") + method = _method_with_mesh(mesh) + child = pybamm.StateVector(slice(0, 1)) + domains = {"primary": ["test"], "secondary": []} + + result = method.broadcast(child, domains, "full to nodes") + + assert result is not child + assert child.domain == [] + assert result.domains["primary"] == ["test"] + np.testing.assert_array_equal(result.evaluate(y=np.array([7])), [[7]]) + + def test_laplacian_and_boundary_conditions(self): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + + plain = method.laplacian(variable, values, {}) + cell_values = np.arange(mesh.npts) + expected = method._tpfa_matrix(mesh) @ cell_values + gradient_matrices, _, _ = method._least_squares_gradient(mesh, {}) + for K, G in zip( + method._cross_term_matrices(mesh), gradient_matrices, strict=True + ): + expected = expected + K @ (G @ cell_values) + np.testing.assert_allclose(plain.evaluate()[:, 0], expected) + + constant = pybamm.Vector(np.full(mesh.npts, 3), domain="test") + dirichlet_bcs = { + variable: { + side: (pybamm.Scalar(3), "Dirichlet") + for side in ["left", "right", "top", "bottom"] + } + } + np.testing.assert_allclose( + method.laplacian(variable, constant, dirichlet_bcs).evaluate(), + 0, + atol=1e-12, + ) + + neumann_bcs = { + variable: { + side: (pybamm.Scalar(0), "Neumann") + for side in ["left", "right", "top", "bottom"] + } + } + np.testing.assert_allclose( + method.laplacian(variable, constant, neumann_bcs).evaluate(), 0, atol=1e-12 + ) + + face_count = len(mesh.boundary_faces["top"]) + vector_bc = pybamm.Vector(np.arange(face_count) + 1) + _, rhs = method._apply_bcs_to_laplacian( + mesh, + method._tpfa_matrix(mesh), + pybamm.Vector(np.zeros(mesh.npts)), + {"top": (vector_bc, "Dirichlet")}, + ) + faces = mesh.boundary_faces["top"] + owners = mesh.face_owner[faces] + distance, alpha, _ = method._boundary_decomposition(mesh, faces) + coefficients = ( + mesh.face_areas[faces] * alpha / distance / mesh.cell_volumes[owners] + ) + expected_rhs = np.zeros(mesh.npts) + np.add.at(expected_rhs, owners, coefficients * (np.arange(face_count) + 1)) + np.testing.assert_allclose(rhs.evaluate()[:, 0], expected_rhs) + + def test_gradient_and_gradient_squared(self): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + constant = pybamm.Vector(np.full(mesh.npts, 3), domain="test") + dirichlet_bcs = { + variable: { + side: (pybamm.Scalar(3), "Dirichlet") + for side in ["left", "right", "top", "bottom"] + } + } + gradient = method.gradient(variable, constant, dirichlet_bcs) + for component in gradient.components: + np.testing.assert_allclose(component.evaluate(), 0, atol=1e-12) + + neumann_bcs = { + variable: { + side: (pybamm.Scalar(0), "Neumann") + for side in ["left", "right", "top", "bottom"] + } + } + for component in method.gradient(variable, constant, neumann_bcs).components: + np.testing.assert_allclose(component.evaluate(), 0, atol=1e-12) + + x_values = mesh.cell_centroids[:, 0] + values = pybamm.Vector(x_values, domain="test") + grad_squared = method.gradient_squared(variable, values, {}) + matrices, _, _ = method._least_squares_gradient(mesh, {}) + expected = sum((matrix @ x_values) ** 2 for matrix in matrices) + np.testing.assert_allclose(grad_squared.evaluate()[:, 0], expected) + + def test_divergence_input_forms_and_error(self): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + symbol = pybamm.Variable("F", domain="test") + components = [ + pybamm.Vector(np.ones(mesh.npts), domain="test"), + pybamm.Vector(np.full(mesh.npts, 2), domain="test"), + ] + + from_list = method.divergence(symbol, components, {}) + from_field = method.divergence(symbol, pybamm.VectorField(*components), {}) + np.testing.assert_allclose(from_list.evaluate(), from_field.evaluate()) + matrices = method._divergence_matrices(mesh) + expected = matrices[0] @ np.ones(mesh.npts) + expected += matrices[1] @ np.full(mesh.npts, 2) + np.testing.assert_allclose(from_list.evaluate()[:, 0], expected) + + with pytest.raises(pybamm.DiscretisationError, match="expects a VectorField"): + method.divergence(symbol, pybamm.Scalar(1), {}) + + def test_divergence_of_bc_bearing_flux_raises(self): + # div of a flux whose gradient parent has BCs is not conservative + # (the BC flux would be silently ignored), so it must raise + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + flux = -(pybamm.Scalar(2) * pybamm.grad(variable)) + components = [ + pybamm.Vector(np.ones(mesh.npts), domain="test"), + pybamm.Vector(np.ones(mesh.npts), domain="test"), + ] + vector_field = pybamm.VectorField(*components) + bcs = {variable: {"left": (pybamm.Scalar(0), "Dirichlet")}} + with pytest.raises(pybamm.DiscretisationError, match="conservative"): + method.divergence(flux, vector_field, bcs) + + # without BCs on u the same flux is fine + result = method.divergence(flux, vector_field, {}) + assert result.evaluate().shape == (mesh.npts, 1) + + def test_gradient_warns_on_bucket_without_bc(self, caplog): + import logging + + mesh = _make_quad_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(0), "Dirichlet"), + } + } + with caplog.at_level(logging.WARNING): + method.gradient(variable, values, bcs) + assert "no boundary condition" in caplog.text + assert "top" in caplog.text and "bottom" in caplog.text + + def test_div_D_grad_scalar_and_vector_coefficients(self): + mesh = _make_2d_mesh(2, 2) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + cell_values = mesh.cell_centroids[:, 0] ** 2 + values = pybamm.Vector(cell_values, domain="test") + bcs = { + variable: { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(2), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + } + } + scalar_result = method.div_D_grad( + div_symbol, variable, pybamm.Scalar(2), values, bcs + ) + + coefficient = pybamm.Vector(np.full(mesh.npts, 2), domain="test") + vector_result = method.div_D_grad( + div_symbol, variable, coefficient, values, bcs + ) + np.testing.assert_allclose( + vector_result.evaluate(), scalar_result.evaluate(), atol=1e-12 + ) + + repeated_domains = {"primary": ["test"], "secondary": ["aux"]} + repeated_div = pybamm.Variable("repeated div", domains=repeated_domains) + repeated_u = pybamm.Variable("repeated u", domains=repeated_domains) + size = mesh.npts * aux.npts + repeated_values = pybamm.Vector( + np.tile(cell_values, aux.npts), domains=repeated_domains + ) + repeated_coefficient = pybamm.Vector(np.full(size, 2), domains=repeated_domains) + repeated = method.div_D_grad( + repeated_div, + repeated_u, + repeated_coefficient, + repeated_values, + { + repeated_u: { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(2), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + } + }, + ) + np.testing.assert_allclose( + repeated.evaluate()[:, 0], + np.tile(vector_result.evaluate()[:, 0], aux.npts), + atol=1e-12, + ) + + def test_div_D_grad_per_face_bc_vectors_with_repeats(self): + # A BC value with one entry per boundary face must be shared across + # auxiliary-domain repeats, matching _bc_contribution's convention. + mesh = _make_2d_mesh(2, 2) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + cell_values = mesh.cell_centroids[:, 0] ** 2 + values = pybamm.Vector(cell_values, domain="test") + + n_left = len(mesh.boundary_faces["left"]) + n_right = len(mesh.boundary_faces["right"]) + bcs = { + variable: { + "left": (pybamm.Vector(np.linspace(1, 2, n_left)), "Dirichlet"), + "right": (pybamm.Vector(np.linspace(-1, 1, n_right)), "Neumann"), + } + } + single = method.div_D_grad(div_symbol, variable, pybamm.Scalar(2), values, bcs) + + repeated_domains = {"primary": ["test"], "secondary": ["aux"]} + repeated_div = pybamm.Variable("repeated div", domains=repeated_domains) + repeated_u = pybamm.Variable("repeated u", domains=repeated_domains) + repeated_values = pybamm.Vector( + np.tile(cell_values, aux.npts), domains=repeated_domains + ) + repeated = method.div_D_grad( + repeated_div, + repeated_u, + pybamm.Scalar(2), + repeated_values, + {repeated_u: bcs[variable]}, + ) + np.testing.assert_allclose( + repeated.evaluate()[:, 0], + np.tile(single.evaluate()[:, 0], aux.npts), + atol=1e-12, + ) + + # One entry per face per repeat passes through untiled + single_right = method.div_D_grad( + div_symbol, + variable, + pybamm.Scalar(2), + values, + {variable: {"right": bcs[variable]["right"]}}, + ) + full = method.div_D_grad( + repeated_div, + repeated_u, + pybamm.Scalar(2), + repeated_values, + { + repeated_u: { + "right": ( + pybamm.Vector(np.tile(np.linspace(-1, 1, n_right), aux.npts)), + "Neumann", + ), + } + }, + ) + np.testing.assert_allclose( + full.evaluate()[:, 0], + np.tile(single_right.evaluate()[:, 0], aux.npts), + atol=1e-12, + ) + + def test_div_D_grad_anisotropic_coefficient_raises(self): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + anisotropic = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2)) + with pytest.raises(pybamm.DiscretisationError, match="Anisotropic"): + method.div_D_grad(div_symbol, variable, anisotropic, values, {}) + + def test_integral_and_boundary_integral(self): + mesh = _make_2d_mesh(2, 2) + aux = _make_2d_mesh(1, 1) + method = _method_with_mesh(mesh, aux=aux) + domains = {"primary": ["test"], "secondary": ["aux"]} + child = pybamm.Variable("u", domains=domains) + values = pybamm.Vector(np.ones(mesh.npts * aux.npts), domains=domains) + + integral = method.integral(child, values, "primary") + np.testing.assert_allclose(integral.evaluate(), 1) + + row = method.definite_integral_matrix(child) + np.testing.assert_allclose( + row.entries.toarray()[0, : mesh.npts], mesh.cell_volumes + ) + assert row.shape == (aux.npts, mesh.npts * aux.npts) + + boundary = method.boundary_integral(child, values, "left") + np.testing.assert_allclose(boundary.evaluate(), 1) + + @pytest.mark.parametrize( + "side", + ["left", "top-right", "top-left", "bottom-right", "bottom-left"], + ) + def test_boundary_value_and_corners(self, side): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + symbol = pybamm.BoundaryValue(variable, side) + + result = method.boundary_value_or_flux(symbol, values) + assert result.domain == [] + if "-" in side: + top_bottom, left_right = side.split("-") + x = mesh.cell_centroids[:, 0] + z = mesh.cell_centroids[:, -1] + target_x = x.max() if left_right == "right" else x.min() + target_z = z.max() if top_bottom == "top" else z.min() + expected = np.argmin((x - target_x) ** 2 + (z - target_z) ** 2) + assert result.evaluate().item() == expected + else: + owners = mesh.face_owner[mesh.boundary_faces[side]] + np.testing.assert_array_equal(result.evaluate()[:, 0], owners) + + def test_boundary_gradient_raises(self): + mesh = _make_2d_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + symbol = pybamm.BoundaryGradient(variable, "left") + with pytest.raises(NotImplementedError, match="BoundaryGradient"): + method.boundary_value_or_flux(symbol, values) + + def test_corner_value_uses_boundary_cell_on_nonconvex_domain(self): + # L-shaped domain: [0,2]x[0,1] plus [0,1]x[1,2]; the top-right + # bounding-box corner (2,2) is outside the domain, and the interior + # cell nearest to it must not be picked + squares = [] + for x0 in (0, 1): + squares.append((x0, 0)) + squares.append((0, 1)) + nodes_list, elems_list = [], [] + node_ids = {} + + def nid(p): + if p not in node_ids: + node_ids[p] = len(nodes_list) + nodes_list.append(p) + return node_ids[p] + + for x0, z0 in squares: + corners = [(x0, z0), (x0 + 1, z0), (x0 + 1, z0 + 1), (x0, z0 + 1)] + elems_list.append([nid(c) for c in corners]) + mesh = UnstructuredSubMesh( + np.array(nodes_list, dtype=float), np.array(elems_list) + ) + mesh.detect_box_boundaries() + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + values = pybamm.Vector(np.arange(mesh.npts), domain="test") + symbol = pybamm.BoundaryValue(variable, "top-right") + result = method.boundary_value_or_flux(symbol, values) + chosen = int(result.evaluate().item()) + candidates = set(mesh.face_owner[mesh.boundary_faces["top"]].tolist()) | set( + mesh.face_owner[mesh.boundary_faces["right"]].tolist() + ) + assert chosen in candidates + + def test_process_binary_operators(self): + method = FiniteVolumeUnstructured() + left_components = [pybamm.StateVector(slice(0, 2)), pybamm.Vector([2, 3])] + right_components = [pybamm.Vector([4, 5]), pybamm.Vector([6, 7])] + left_field = pybamm.VectorField(*left_components) + right_field = pybamm.VectorField(*right_components) + multiplication = pybamm.Multiplication(pybamm.Scalar(1), pybamm.Scalar(2)) + + both = method.process_binary_operators( + multiplication, + None, + None, + left_field, + right_field, + ) + assert both.n_components == 2 + np.testing.assert_array_equal( + both.components[0].evaluate(y=np.array([1, 2]))[:, 0], [4, 10] + ) + np.testing.assert_array_equal(both.components[1].evaluate()[:, 0], [12, 21]) + + field_left = method.process_binary_operators( + multiplication, None, None, left_field, pybamm.Scalar(2) + ) + field_right = method.process_binary_operators( + multiplication, None, None, pybamm.Scalar(2), right_field + ) + np.testing.assert_array_equal( + field_left.components[0].evaluate(y=np.array([1, 2]))[:, 0], [2, 4] + ) + np.testing.assert_array_equal( + field_right.components[0].evaluate()[:, 0], [8, 10] + ) + + scalar = method.process_binary_operators( + multiplication, None, None, pybamm.Scalar(3), pybamm.Scalar(4) + ) + assert scalar.evaluate() == 12 + + def test_internal_neumann_unstructured_paths(self): + left, right = _make_split_2d_meshes(2, 2, 2) + method = FiniteVolumeUnstructured() + # u = x: the interface normal gradient is exactly 1 once the + # non-orthogonal cross term is included (the interface cells touch + # no external side where the fitted zero normal derivative is wrong) + left_values = pybamm.Vector(left.cell_centroids[:, 0], domain="left") + right_values = pybamm.Vector(right.cell_centroids[:, 0], domain="right") + + direct = method._internal_neumann_unstructured( + left_values, right_values, left, right, 1 + ) + np.testing.assert_allclose(direct.evaluate()[:, 0], 1.0, atol=1e-12) + + # orthogonal interface (quads): plain two-point difference + quad_left = _make_quad_mesh(2, 2, x_range=(0, 0.5)) + quad_right = _make_quad_mesh(2, 2, x_range=(0.5, 1)) + compute_interface_data(quad_left, quad_right, "left", "right") + interface = quad_left.interface_data["right"] + u_left, u_right = np.arange(quad_left.npts), np.arange(quad_right.npts) + quad_value = method._internal_neumann_unstructured( + pybamm.Vector(u_left, domain="left"), + pybamm.Vector(u_right, domain="right"), + quad_left, + quad_right, + 1, + ) + expected = ( + u_right[interface["right_cells"]] - u_left[interface["left_cells"]] + ) / interface["cell_distances"] + np.testing.assert_allclose(quad_value.evaluate()[:, 0], expected) + + left_data = left.interface_data + left.interface_data = {} + reverse = method._internal_neumann_unstructured( + left_values, right_values, left, right, 1 + ) + np.testing.assert_allclose(reverse.evaluate(), direct.evaluate()) + + # unpaired meshes raise: silently returning zeros would decouple + # the domains and solve to a wrong answer + right.interface_data = {} + with pytest.raises(pybamm.DiscretisationError, match="decoupled"): + method._internal_neumann_unstructured( + left_values, right_values, left, right, 2 + ) + left.interface_data = left_data + + def test_internal_neumann_dispatch_structured_and_mismatch(self): + method = FiniteVolumeUnstructured() + left_mesh = pybamm.SubMesh1D(np.array([0, 0.5]), "cartesian") + right_mesh = pybamm.SubMesh1D(np.array([0.5, 1]), "cartesian") + left = pybamm.Vector(np.arange(left_mesh.npts), domain="left") + right = pybamm.Vector(np.arange(right_mesh.npts), domain="right") + + structured = method.internal_neumann_condition( + left, right, left_mesh, right_mesh + ) + dx = right_mesh.nodes[0] - left_mesh.nodes[-1] + expected = (np.arange(right_mesh.npts)[0] - np.arange(left_mesh.npts)[-1]) / dx + assert structured.evaluate().item() == expected + + unstructured_left, unstructured_right = _make_split_2d_meshes(1, 1, 1) + method._mesh = _MeshMap( + { + ("aux",): _make_2d_mesh(1, 1), + ("other aux",): _make_2d_mesh(2, 1), + } + ) + left_repeated = pybamm.Vector( + np.ones(unstructured_left.npts * method.mesh["aux"].npts), + domains={"primary": ["left"], "secondary": ["aux"]}, + ) + right_repeated = pybamm.Vector( + np.ones(unstructured_right.npts * method.mesh["other aux"].npts), + domains={"primary": ["right"], "secondary": ["other aux"]}, + ) + with pytest.raises(pybamm.DomainError, match="secondary points"): + method.internal_neumann_condition( + left_repeated, + right_repeated, + unstructured_left, + unstructured_right, + ) + + def test_internal_bcs_for_concatenation(self): + left = _make_2d_mesh(1, 1, x_range=(0, 0.5)) + right = _make_2d_mesh(1, 1, x_range=(0.5, 1)) + method = FiniteVolumeUnstructured() + method._compute_pair_interface(left, right, "left", "right") + method._mesh = _MeshMap({("left",): left, ("right",): right}) + children = [ + pybamm.Variable("left temperature", domain="left"), + pybamm.Variable("right temperature", domain="right"), + ] + + class Disc: + def process_symbol(self, child): + size = method.mesh[child.domain].npts + return pybamm.Vector(np.ones(size), domains=child.domains) + + result = method.set_internal_bcs_for_concat( + Disc(), + children[0], + children, + {"left": (pybamm.Scalar(0), "Dirichlet")}, + ) + assert set(result) == set(children) + assert "iface_right" in result[children[0]] + interface_gradient, bc_type = result[children[0]]["iface_right"] + assert bc_type == "Neumann" + np.testing.assert_allclose(interface_gradient.evaluate(), 0, atol=1e-12) + + structured = pybamm.SubMesh1D(np.array([0, 1]), "cartesian") + method._mesh[("structured",)] = structured + structured_child = pybamm.Variable( + "structured temperature", domain="structured" + ) + partial = method.set_internal_bcs_for_concat( + Disc(), + children[0], + [children[0], structured_child], + {}, + ) + assert structured_child not in partial + assert partial[children[0]] == {} + + no_interface = _method_with_mesh(_make_2d_mesh(1, 1)) + assert ( + no_interface.set_internal_bcs_for_concat( + Disc(), children[0], [pybamm.Variable("u", domain="test")], {} + ) + is None + ) + + def test_concatenation_preserves_domain_order(self): + left = _make_2d_mesh(1, 1, x_range=(0, 0.5)) + right = _make_2d_mesh(1, 1, x_range=(0.5, 1)) + method = FiniteVolumeUnstructured() + method._mesh = _MeshMap({("left",): left, ("right",): right}) + left_values = pybamm.Vector([1, 2], domain="left") + right_values = pybamm.Vector([3, 4], domain="right") + + result = method.concatenation([left_values, right_values]) + + np.testing.assert_array_equal(result.evaluate()[:, 0], [1, 2, 3, 4]) + assert result.domain == ["left", "right"] + + +# ====================================================================== +# Tests: Discretisation dispatch +# ====================================================================== + + +def _get_unstructured_disc(nx=4, nz=4): + """Single-domain 2D unstructured discretisation on [0,1]^2.""" + x = pybamm.SpatialVariable( + "x_n", domain=["negative electrode"], coord_sys="cartesian" + ) + z = pybamm.SpatialVariable( + "z_2d", domain=["negative electrode"], coord_sys="cartesian", direction="tb" + ) + geometry = { + "negative electrode": {x: {"min": 0.0, "max": 1.0}, z: {"min": 0.0, "max": 1.0}} + } + mesh = pybamm.Mesh( + geometry, + { + "negative electrode": pybamm.meshes.unstructured_submesh.UnstructuredMeshGenerator() + }, + {x: nx, z: nz}, + ) + return pybamm.Discretisation( + mesh, {"negative electrode": FiniteVolumeUnstructured()} + ) + + +class TestDiscretisationDispatch: + def _disc_var_grad(self): + disc = _get_unstructured_disc() + var = pybamm.Variable("u", domain=["negative electrode"]) + disc.set_variable_slices([var]) + grad = pybamm.grad(var) + disc_grad = disc.process_symbol(grad) + u = disc.mesh["negative electrode"].cell_centroids[:, 0] + return disc, var, grad, disc_grad, u + + def test_component_of_gradient(self): + disc, _, grad, disc_grad, u = self._disc_var_grad() + comp0 = disc.process_symbol(pybamm.Component(grad, 0)) + np.testing.assert_allclose( + comp0.evaluate(y=u), + disc_grad.components[0].evaluate(y=u), + ) + comp1 = disc.process_symbol(pybamm.Component(grad, 1)) + np.testing.assert_allclose( + comp1.evaluate(y=u), + disc_grad.components[1].evaluate(y=u), + ) + + def test_component_requires_vector_field(self): + disc, var, *_ = self._disc_var_grad() + with pytest.raises( + pybamm.DiscretisationError, match="Component can only be applied" + ): + disc.process_symbol(pybamm.Component(var, 0)) + + def test_norm_of_gradient(self): + disc, _, grad, disc_grad, u = self._disc_var_grad() + norm = disc.process_symbol(pybamm.Norm(grad)) + gx = disc_grad.components[0].evaluate(y=u) + gz = disc_grad.components[1].evaluate(y=u) + np.testing.assert_allclose( + norm.evaluate(y=u), np.sqrt(gx**2 + gz**2), rtol=1e-12 + ) + + def test_norm_requires_vector_field(self): + disc, var, *_ = self._disc_var_grad() + with pytest.raises( + pybamm.DiscretisationError, match="Norm can only be applied" + ): + disc.process_symbol(pybamm.Norm(var)) + + def test_generic_unary_maps_over_components(self): + """A generic unary operator (negation) applies componentwise to a + VectorField.""" + disc, _, grad, disc_grad, u = self._disc_var_grad() + neg = disc.process_symbol(-grad) + assert isinstance(neg, pybamm.VectorField) + for k in range(2): + np.testing.assert_allclose( + neg.components[k].evaluate(y=u), + -disc_grad.components[k].evaluate(y=u), + atol=1e-12, + ) + + def test_scalar_times_gradient_lifted(self): + """Scalar * grad(u) lifts the scalar to an N-component VectorField.""" + disc, _, grad, disc_grad, u = self._disc_var_grad() + scaled = disc.process_symbol(pybamm.Scalar(2) * grad) + assert isinstance(scaled, pybamm.VectorField) + for k in range(2): + np.testing.assert_allclose( + scaled.components[k].evaluate(y=u), + 2 * disc_grad.components[k].evaluate(y=u), + atol=1e-12, + ) + + def test_gradient_times_scalar_lifted(self): + disc, _, grad, disc_grad, u = self._disc_var_grad() + scaled = disc.process_symbol(grad * pybamm.Scalar(3)) + assert isinstance(scaled, pybamm.VectorField) + for k in range(2): + np.testing.assert_allclose( + scaled.components[k].evaluate(y=u), + 3 * disc_grad.components[k].evaluate(y=u), + atol=1e-12, + ) + + def test_domainless_vector_field_binary_op(self): + """Binary ops on domainless VectorFields combine componentwise.""" + disc = _get_unstructured_disc() + vf_a = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2)) + vf_b = pybamm.VectorField(pybamm.Scalar(3), pybamm.Scalar(4)) + product = disc.process_symbol(vf_a * vf_b) + assert isinstance(product, pybamm.VectorField) + np.testing.assert_allclose(product.components[0].evaluate(), 3) + np.testing.assert_allclose(product.components[1].evaluate(), 8) + + # Scalar lifted to match the VectorField's components + scaled = disc.process_symbol(pybamm.Scalar(2) * vf_a) + assert isinstance(scaled, pybamm.VectorField) + np.testing.assert_allclose(scaled.components[0].evaluate(), 2) + np.testing.assert_allclose(scaled.components[1].evaluate(), 4) + + +class TestProcessModelConcatenation: + def test_two_domain_diffusion_steady_state(self): + """process_model on a concatenated variable dispatches internal BCs + through FiniteVolumeUnstructured; the discrete Laplacian of the exact + steady profile (linear in x) is zero.""" + x_n = pybamm.SpatialVariable( + "x_n", domain=["negative electrode"], coord_sys="cartesian" + ) + x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian") + z = pybamm.SpatialVariable( + "z_2d", + domain=["negative electrode", "separator"], + coord_sys="cartesian", + direction="tb", + ) + geometry = { + "negative electrode": { + x_n: {"min": 0.0, "max": 0.5}, + z: {"min": 0.0, "max": 1.0}, + }, + "separator": { + x_s: {"min": 0.5, "max": 1.0}, + z: {"min": 0.0, "max": 1.0}, + }, + } + gen = pybamm.meshes.unstructured_submesh.UnstructuredMeshGenerator( + element_type="quad" + ) + mesh = pybamm.Mesh( + geometry, + {"negative electrode": gen, "separator": gen}, + {x_n: 3, x_s: 3, z: 3}, + ) + disc = pybamm.Discretisation( + mesh, + { + "negative electrode": FiniteVolumeUnstructured(), + "separator": FiniteVolumeUnstructured(), + }, + ) + + var_n = pybamm.Variable("c_n", domain=["negative electrode"]) + var_s = pybamm.Variable("c_s", domain=["separator"]) + var = pybamm.concatenation(var_n, var_s) + + model = pybamm.BaseModel() + model.rhs = {var: pybamm.div(pybamm.grad(var))} + model.initial_conditions = {var: pybamm.Scalar(1)} + model.boundary_conditions = { + var: { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(1), "Dirichlet"), + } + } + model.variables = {"c": var} + model_disc = disc.process_model(model, inplace=False) + + u = np.concatenate( + [ + mesh["negative electrode"].cell_centroids[:, 0], + mesh["separator"].cell_centroids[:, 0], + ] + ) + rhs = model_disc.concatenated_rhs.evaluate(t=0, y=u).flatten() + np.testing.assert_allclose(rhs, 0.0, atol=1e-10) + + +class TestDiscretisationDispatchLifting: + def _disc_var_grad(self): + disc = _get_unstructured_disc() + var = pybamm.Variable("u", domain=["negative electrode"]) + disc.set_variable_slices([var]) + grad = pybamm.grad(var) + disc_grad = disc.process_symbol(grad) + u = disc.mesh["negative electrode"].cell_centroids[:, 0] + return disc, var, grad, disc_grad, u + + def test_gradient_minus_scalar_lifted(self): + """A right-hand Scalar is lifted to an N-component VectorField. + + A raw Subtraction node is used because operator simplification + rewrites ``x - c`` as ``-c + x``, which takes the left-Scalar path. + """ + disc, _, grad, disc_grad, u = self._disc_var_grad() + shifted = disc.process_symbol(pybamm.Subtraction(grad, pybamm.Scalar(0.5))) + assert isinstance(shifted, pybamm.VectorField) + for k in range(2): + np.testing.assert_allclose( + shifted.components[k].evaluate(y=u), + disc_grad.components[k].evaluate(y=u) - 0.5, + atol=1e-12, + ) + + def test_domainless_vector_field_minus_scalar(self): + disc = _get_unstructured_disc() + vf = pybamm.VectorField(pybamm.Scalar(3), pybamm.Scalar(4)) + shifted = disc.process_symbol(pybamm.Subtraction(vf, pybamm.Scalar(1))) + assert isinstance(shifted, pybamm.VectorField) + np.testing.assert_allclose(shifted.components[0].evaluate(), 2) + np.testing.assert_allclose(shifted.components[1].evaluate(), 3) + + def test_div_of_coefficient_times_gradient(self): + """div(D * grad(u)) is intercepted and routed to div_D_grad for both + coefficient orderings.""" + disc, var, grad, _, u = self._disc_var_grad() + base = disc.process_symbol(pybamm.div(grad)).evaluate(y=u) + scaled = disc.process_symbol(pybamm.div(pybamm.Scalar(2) * grad)).evaluate(y=u) + np.testing.assert_allclose(scaled, 2 * base, atol=1e-12) + + right_form = disc.process_symbol(pybamm.div(var * grad)).evaluate(y=u) + left_form = disc.process_symbol(pybamm.div(grad * var)).evaluate(y=u) + np.testing.assert_allclose(left_form, right_form, atol=1e-12) + + +class TestProcessModelConcatenationZStack: + def test_z_stacked_domains_use_graph_internal_bcs(self): + """Domains stacked in z: pybamm.Mesh's 1D-stack pairing fails on the + transverse mismatch, FiniteVolumeUnstructured's build() discovers the + interface by face matching, and process_model routes internal BCs + through set_internal_bcs_for_concat. The discrete Laplacian of the + exact steady profile (linear in z) is zero.""" + x_n = pybamm.SpatialVariable( + "x_n", domain=["negative electrode"], coord_sys="cartesian" + ) + z_n = pybamm.SpatialVariable( + "z_n", domain=["negative electrode"], coord_sys="cartesian" + ) + x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian") + z_s = pybamm.SpatialVariable("z_s", domain=["separator"], coord_sys="cartesian") + geometry = { + "negative electrode": { + x_n: {"min": 0.0, "max": 1.0}, + z_n: {"min": 0.0, "max": 0.5}, + }, + "separator": { + x_s: {"min": 0.0, "max": 1.0}, + z_s: {"min": 0.5, "max": 1.0}, + }, + } + gen = pybamm.meshes.unstructured_submesh.UnstructuredMeshGenerator( + element_type="quad" + ) + mesh = pybamm.Mesh( + geometry, + {"negative electrode": gen, "separator": gen}, + {x_n: 3, z_n: 3, x_s: 3, z_s: 3}, + ) + disc = pybamm.Discretisation( + mesh, + { + "negative electrode": FiniteVolumeUnstructured(), + "separator": FiniteVolumeUnstructured(), + }, + ) + # build() added graph-discovered interface buckets + assert any( + tag.startswith("iface_") + for tag in mesh["negative electrode"].boundary_faces + ) + assert any(tag.startswith("iface_") for tag in mesh["separator"].boundary_faces) + + var_n = pybamm.Variable("c_n", domain=["negative electrode"]) + var_s = pybamm.Variable("c_s", domain=["separator"]) + var = pybamm.concatenation(var_n, var_s) + + model = pybamm.BaseModel() + model.rhs = {var: pybamm.div(pybamm.grad(var))} + model.initial_conditions = {var: pybamm.Scalar(1)} + model.boundary_conditions = { + var: { + "bottom": (pybamm.Scalar(0), "Dirichlet"), + "top": (pybamm.Scalar(1), "Dirichlet"), + } + } + model.variables = {"c": var} + model_disc = disc.process_model(model, inplace=False) + + u = np.concatenate( + [ + mesh["negative electrode"].cell_centroids[:, 1], + mesh["separator"].cell_centroids[:, 1], + ] + ) + rhs = model_disc.concatenated_rhs.evaluate(t=0, y=u).flatten() + np.testing.assert_allclose(rhs, 0.0, atol=1e-10) + + +# ====================================================================== +# Tests: non-orthogonal correction +# ====================================================================== + + +def _perturb_interior_nodes(nodes, spacing, fraction=0.3, seed=0): + """Jitter interior nodes so faces are skewed as well as non-orthogonal.""" + rng = np.random.default_rng(seed) + nodes = nodes.copy() + low, high = nodes.min(axis=0), nodes.max(axis=0) + on_boundary = np.any(np.isclose(nodes, low) | np.isclose(nodes, high), axis=1) + interior = nodes[~on_boundary] + nodes[~on_boundary] = interior + rng.uniform( + -fraction * spacing, fraction * spacing, interior.shape + ) + return nodes + + +def _make_perturbed_tri_mesh(n=6): + edges = np.linspace(0, 1, n + 1) + nodes, elements = _quad_to_tri(edges, edges) + mesh = UnstructuredSubMesh(_perturb_interior_nodes(nodes, 1.0 / n), elements) + mesh.detect_box_boundaries() + return mesh + + +def _dirichlet_all_sides(mesh, u_exact): + return { + side: (pybamm.Vector(u_exact(mesh.face_centroids[faces])), "Dirichlet") + for side, faces in mesh.boundary_faces.items() + } + + +def _laplacian_system(method, mesh, bcs): + """``(L, rhs)`` with ``laplacian(u) = L @ u + rhs`` for the full operator.""" + variable = pybamm.Variable("u", domain="test") + y = pybamm.StateVector(slice(0, mesh.npts), domains={"primary": ["test"]}) + expr = method.laplacian(variable, y, {variable: bcs} if bcs else {}) + zeros = np.zeros(mesh.npts) + return sp_csr(expr.jac(y).evaluate(y=zeros)), expr.evaluate(y=zeros)[:, 0] + + +class TestNonOrthogonalCorrection: + @pytest.mark.parametrize( + "make_mesh", + [ + lambda: _make_2d_mesh(6, 6), + _make_perturbed_tri_mesh, + lambda: _make_3d_mesh(3, 3, 3), + ], + ids=["tri", "tri-perturbed", "tet"], + ) + @pytest.mark.parametrize("correction", ["over-relaxed", "minimum"]) + def test_laplacian_exact_on_linear_field(self, make_mesh, correction): + """The discrete Laplacian of a linear field vanishes on every cell; + the two-point part alone fails this on any non-orthogonal mesh.""" + mesh = make_mesh() + method = FiniteVolumeUnstructured({"non-orthogonal correction": correction}) + method._mesh = _MeshMap({("test",): mesh}) + slope = np.array([1.0, 0.7, 0.4])[: mesh.dimension] + u = mesh.cell_centroids @ slope + bcs = _dirichlet_all_sides(mesh, lambda points: points @ slope) + L, rhs = _laplacian_system(method, mesh, bcs) + np.testing.assert_allclose(L @ u + rhs, 0, atol=1e-10) + + def test_two_point_part_alone_is_not_exact(self): + mesh = _make_2d_mesh(6, 6) + L = FiniteVolumeUnstructured()._tpfa_matrix(mesh) + residual = np.abs(L @ mesh.cell_centroids[:, 0]) + assert residual[_get_internal_cells(mesh)].max() > 1 + + def test_second_order_convergence_on_triangles(self): + def u_exact(points): + return np.sin(np.pi * points[:, 0]) * np.sin(np.pi * points[:, 1]) + + errors = [] + for n in (8, 16, 32): + mesh = _make_2d_mesh(n, n) + method = _method_with_mesh(mesh) + L, rhs = _laplacian_system( + method, mesh, _dirichlet_all_sides(mesh, u_exact) + ) + source = -2 * np.pi**2 * u_exact(mesh.cell_centroids) + u = spsolve(L.tocsc(), source - rhs) + error = u - u_exact(mesh.cell_centroids) + errors.append(np.sqrt(np.sum(mesh.cell_volumes * error**2))) + rates = np.log2(np.array(errors[:-1]) / np.array(errors[1:])) + assert np.all(rates > 1.7), rates + + def test_orthogonal_mesh_has_no_cross_term(self): + mesh = _make_quad_mesh(4, 4) + method = FiniteVolumeUnstructured() + assert method._cross_term_matrices(mesh) is None + assert method._div_D_grad_matrices(mesh)[4] is None + np.testing.assert_allclose(method._decomposition(mesh)[0], 1.0) + + def test_full_operator_is_conservative(self): + mesh = _make_perturbed_tri_mesh(5) + L, _ = _laplacian_system(_method_with_mesh(mesh), mesh, None) + u = mesh.cell_centroids[:, 0] ** 2 + mesh.cell_centroids[:, 1] + np.testing.assert_allclose(np.sum((L @ u) * mesh.cell_volumes), 0, atol=1e-10) + + def test_div_D_grad_matches_laplacian_for_constant_D(self): + mesh = _make_perturbed_tri_mesh(4) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + + def u_exact(points): + return np.sin(points[:, 0]) * points[:, 1] ** 2 + + values = pybamm.Vector(u_exact(mesh.cell_centroids), domain="test") + bcs = {variable: _dirichlet_all_sides(mesh, u_exact)} + laplacian = method.laplacian(variable, values, bcs).evaluate()[:, 0] + div_grad = method.div_D_grad( + div_symbol, variable, pybamm.Scalar(2), values, bcs + ) + np.testing.assert_allclose(div_grad.evaluate()[:, 0], 2 * laplacian, atol=1e-10) + + def test_div_D_grad_exact_on_linear_field_3d(self): + mesh = _make_3d_mesh(3, 3, 3) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + slope = np.array([1.0, 0.7, 0.4]) + values = pybamm.Vector(mesh.cell_centroids @ slope, domain="test") + bcs = {variable: _dirichlet_all_sides(mesh, lambda points: points @ slope)} + result = method.div_D_grad(div_symbol, variable, pybamm.Scalar(2), values, bcs) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + def test_least_squares_gradient_exact_on_linear_field(self): + mesh = _make_perturbed_tri_mesh(5) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + slope = np.array([2.0, -3.0]) + values = pybamm.Vector(mesh.cell_centroids @ slope, domain="test") + + dirichlet = { + variable: _dirichlet_all_sides(mesh, lambda points: points @ slope) + } + components = method.gradient(variable, values, dirichlet).components + for k, component in enumerate(components): + np.testing.assert_allclose(component.evaluate()[:, 0], slope[k], atol=1e-10) + + # named sides take coordinate-direction derivatives on both ends + neumann = { + variable: { + "left": (pybamm.Scalar(2.0), "Neumann"), + "right": (pybamm.Scalar(2.0), "Neumann"), + "bottom": (pybamm.Scalar(-3.0), "Neumann"), + "top": (pybamm.Scalar(-3.0), "Neumann"), + } + } + components = method.gradient(variable, values, neumann).components + for k, component in enumerate(components): + np.testing.assert_allclose(component.evaluate()[:, 0], slope[k], atol=1e-10) + + def test_green_gauss_gradient_is_not_exact_on_skewed_mesh(self): + """Documents why the cross term cannot use the Green-Gauss gradient.""" + mesh = _make_3d_mesh(3, 3, 3) + G = FiniteVolumeUnstructured()._green_gauss_matrices(mesh) + grad_x = (G[0] @ mesh.cell_centroids[:, 0])[_get_internal_cells(mesh)] + assert np.abs(grad_x - 1).max() > 0.1 + + def test_divergence_shares_green_gauss_assembly(self): + mesh = _make_2d_mesh(3, 3) + method = FiniteVolumeUnstructured() + assert method._divergence_matrices(mesh) is method._green_gauss_matrices(mesh) + + def test_invalid_option_raises(self): + with pytest.raises(pybamm.OptionError, match="non-orthogonal correction"): + FiniteVolumeUnstructured({"non-orthogonal correction": "none"}) + + def test_option_sets_two_point_weight(self): + mesh = _make_2d_mesh(3, 3) + cos_theta = FiniteVolumeUnstructured()._face_geometry(mesh)["cos_theta"] + minimum = FiniteVolumeUnstructured({"non-orthogonal correction": "minimum"}) + over_relaxed = FiniteVolumeUnstructured() + np.testing.assert_allclose(minimum._decomposition(mesh)[0], cos_theta) + np.testing.assert_allclose(over_relaxed._decomposition(mesh)[0], 1 / cos_theta) + + def test_inverted_cell_raises(self): + mesh = _make_2d_mesh(2, 2) + mesh.face_normals[: mesh.n_internal_faces] *= -1 + with pytest.raises(pybamm.GeometryError, match="pointing away"): + FiniteVolumeUnstructured()._face_geometry(mesh) + + def test_build_warns_on_severe_non_orthogonality(self, caplog): + import logging + + # Sliver neighbour: the centroid line is ~85 degrees off the normal + nodes = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [10.0, -8.5]]) + elements = np.array([[0, 1, 2], [1, 3, 2]]) + skewed = UnstructuredSubMesh(nodes, elements) + skewed.detect_box_boundaries() + method = FiniteVolumeUnstructured() + assert method._face_geometry(skewed)["max_angle_deg"] > 70 + with caplog.at_level(logging.WARNING): + method.build(_MeshMap({("skewed",): skewed})) + assert "non-orthogonality" in caplog.text + + caplog.clear() + with caplog.at_level(logging.WARNING): + method.build(_MeshMap({("tri",): _make_2d_mesh(3, 3)})) + assert "non-orthogonality" not in caplog.text + + +# ====================================================================== +# Tests: non-orthogonal correction across domain interfaces +# ====================================================================== + + +def _two_domain_laplacian(element_type, boundary_conditions): + """Discretised ``div(grad(c))`` of a two-domain concatenation on the + unit box split at x = 0.5, returning ``(mesh, rhs_expression)``.""" + dim3 = element_type in ("tetrahedron", "hexahedron") + domains = ["negative electrode", "separator"] + x_n = pybamm.SpatialVariable("x_n", domain=[domains[0]], coord_sys="cartesian") + x_s = pybamm.SpatialVariable("x_s", domain=[domains[1]], coord_sys="cartesian") + z = pybamm.SpatialVariable( + "z_2d", domain=domains, coord_sys="cartesian", direction="tb" + ) + geometry = { + domains[0]: {x_n: {"min": 0.0, "max": 0.5}, z: {"min": 0.0, "max": 1.0}}, + domains[1]: {x_s: {"min": 0.5, "max": 1.0}, z: {"min": 0.0, "max": 1.0}}, + } + npts = {x_n: 3, x_s: 3, z: 4} + if dim3: + y = pybamm.SpatialVariable("y", domain=domains, coord_sys="cartesian") + for domain in domains: + geometry[domain][y] = {"min": 0.0, "max": 1.0} + npts[y] = 3 + generator = pybamm.meshes.unstructured_submesh.UnstructuredMeshGenerator( + element_type=element_type + ) + mesh = pybamm.Mesh(geometry, dict.fromkeys(domains, generator), npts) + disc = pybamm.Discretisation( + mesh, {domain: FiniteVolumeUnstructured() for domain in domains} + ) + var_n = pybamm.Variable("c_n", domain=[domains[0]]) + var_s = pybamm.Variable("c_s", domain=[domains[1]]) + var = pybamm.concatenation(var_n, var_s) + model = pybamm.BaseModel() + model.rhs = {var: pybamm.div(pybamm.grad(var))} + model.initial_conditions = {var: pybamm.Scalar(1)} + model.boundary_conditions = {var: boundary_conditions} + model.variables = {"c": var} + disc.process_model(model, inplace=False) + return mesh, disc.process_model(model, inplace=False).concatenated_rhs + + +class TestInterfaceCorrection: + @pytest.mark.parametrize("element_type", ["triangle", "tetrahedron"]) + def test_linear_field_exact_across_interface(self, element_type): + """u = x is the steady state of left=0, right=1 with zero flux on the + other sides; the interface flux must reproduce it on skewed pairs.""" + mesh, rhs = _two_domain_laplacian( + element_type, + { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(1), "Dirichlet"), + }, + ) + u = np.concatenate( + [ + mesh["negative electrode"].cell_centroids[:, 0], + mesh["separator"].cell_centroids[:, 0], + ] + ) + np.testing.assert_allclose(rhs.evaluate(y=u), 0, atol=1e-10) + + def test_interface_flux_is_conservative(self): + mesh, rhs = _two_domain_laplacian( + "tetrahedron", + { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + }, + ) + volumes = np.concatenate( + [mesh["negative electrode"].cell_volumes, mesh["separator"].cell_volumes] + ) + u = np.random.default_rng(1).uniform(size=len(volumes)) + np.testing.assert_allclose(volumes @ rhs.evaluate(y=u)[:, 0], 0, atol=1e-10) + + def test_interface_data_records_faces(self): + left, right = _make_split_2d_meshes(3, 3, 3) + data = left.interface_data["right"] + np.testing.assert_array_equal(data["left_faces"], left.boundary_faces["right"]) + assert set(data["right_faces"]) == set(right.boundary_faces["left"]) + np.testing.assert_array_equal( + left.face_owner[data["left_faces"]], data["left_cells"] + ) + + a = _make_2d_mesh(2, 2, x_range=(0, 0.5)) + b = _make_2d_mesh(2, 2, x_range=(0.5, 1)) + assert FiniteVolumeUnstructured()._compute_pair_interface(a, b, "a", "b") + np.testing.assert_array_equal( + a.interface_data["b"]["left_faces"], a.boundary_faces["iface_b"] + ) + np.testing.assert_array_equal( + b.interface_data["a"]["left_faces"], b.boundary_faces["iface_a"] + ) + + def test_inconsistent_face_counts_raise(self): + mesh = _make_2d_mesh(2, 2) + # give one cell an extra boundary face and another one fewer + boundary = mesh.boundary_faces["left"] + mesh.face_owner[boundary[0]] = mesh.face_owner[boundary[1]] + with pytest.raises(pybamm.DiscretisationError, match="same number of faces"): + FiniteVolumeUnstructured()._least_squares_matrices(mesh, {}) + + +class TestTimeDependentScalarBoundaryValues: + def test_neumann_value_depending_on_time(self): + """A ``pybamm.t``-dependent Neumann value evaluates for shape to + ``()``; it must still be broadcast over the side's faces.""" + mesh = _make_2d_mesh(3, 3) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + slope = 2.0 * pybamm.t + values = pybamm.Vector(mesh.cell_centroids[:, 0], domain="test") + bcs = { + variable: { + "left": (slope, "Neumann"), + "right": (slope, "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + # at t = 0.5 the prescribed slope matches u = x, so both vanish + laplacian = method.laplacian(variable, values, bcs) + np.testing.assert_allclose(laplacian.evaluate(t=0.5), 0, atol=1e-10) + div_grad = method.div_D_grad( + div_symbol, variable, pybamm.Scalar(3), values, bcs + ) + np.testing.assert_allclose(div_grad.evaluate(t=0.5), 0, atol=1e-10) + grad_x = method.gradient(variable, values, bcs).components[0] + np.testing.assert_allclose(grad_x.evaluate(t=0.5), 1, atol=1e-10) + assert np.abs(laplacian.evaluate(t=1.0)).max() > 1e-3 + + +class TestHarmonicDiffusivity: + def test_two_material_slab_is_exact(self): + """Piecewise-constant D with the jump on a face: the exact steady + profile is piecewise linear with the series-resistance flux, which the + harmonic mean reproduces and the arithmetic mean does not. Uneven + cells make the weight orientation matter.""" + x_edges = np.array([0.0, 0.3, 0.5, 0.6, 1.0]) + z_edges = np.linspace(0, 1, 3) + nodes, elements = _make_quad_grid(x_edges, z_edges) + mesh = UnstructuredSubMesh(nodes, elements) + mesh.detect_box_boundaries() + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + + x = mesh.cell_centroids[:, 0] + D1, D2 = 1.0, 25.0 + D_cells = np.where(x < 0.5, D1, D2) + flux = 1.0 / (0.5 / D1 + 0.5 / D2) + u_exact = np.where(x < 0.5, flux * x / D1, 1.0 - flux * (1.0 - x) / D2) + bcs = { + variable: { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(1), "Dirichlet"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + } + result = method.div_D_grad( + div_symbol, + variable, + pybamm.Vector(D_cells, domain="test"), + pybamm.Vector(u_exact, domain="test"), + bcs, + ) + np.testing.assert_allclose(result.evaluate(), 0, atol=1e-10) + + # the same profile is not a steady state under the arithmetic mean + _, W, S, geo, _, _ = method._div_D_grad_matrices(mesh) + G = method._div_D_grad_matrices(mesh)[0] + arithmetic = S @ ((W @ D_cells) * (G @ u_exact) * geo) + assert np.abs(arithmetic).max() > 1e-2 + + def test_face_diffusivity_is_harmonic_mean(self): + mesh = _make_quad_mesh(3, 1) + method = _method_with_mesh(mesh) + W_harmonic = method._div_D_grad_matrices(mesh)[5] + D = np.array([1.0, 4.0, 4.0]) + face_D = 1 / (W_harmonic @ (1 / D)) + # uniform cells: plain harmonic mean 2 D1 D2 / (D1 + D2) + np.testing.assert_allclose(np.sort(face_D), [1.6, 4.0]) + + def test_nonpositive_coefficient_raises(self): + mesh = _make_quad_mesh(2, 2) + method = _method_with_mesh(mesh) + variable = pybamm.Variable("u", domain="test") + div_symbol = pybamm.Variable("div", domain="test") + values = pybamm.Vector(np.ones(mesh.npts), domain="test") + with pytest.raises(pybamm.DiscretisationError, match="strictly positive"): + method.div_D_grad( + div_symbol, + variable, + pybamm.Vector(np.zeros(mesh.npts), domain="test"), + values, + {}, + ) + + +class TestOrthogonalityTolerance: + """High-aspect-ratio boxes (10 um thick, cm wide, as in a pouch cell) put + ~1e-11 of centroid rounding into the face direction; that must not switch + the wide cross-term stencil on.""" + + def test_anisotropic_boxes_stay_orthogonal(self): + quad = UnstructuredSubMesh( + *_make_quad_grid(np.linspace(0, 1e-5, 6), np.linspace(0, 0.03, 4)) + ) + hexa = UnstructuredSubMesh( + *_hex_grid( + np.linspace(0, 1e-5, 4), np.linspace(0, 0.2, 4), np.linspace(0, 0.1, 4) + ) + ) + method = FiniteVolumeUnstructured() + for mesh in (quad, hexa): + mesh.detect_box_boundaries() + _, k = method._decomposition(mesh) + assert not k.any() + assert method._cross_term_matrices(mesh) is None + assert method._div_D_grad_matrices(mesh)[4] is None + for faces in mesh.boundary_faces.values(): + assert not method._boundary_decomposition(mesh, faces)[2].any() + + def test_stencil_stays_compact_on_anisotropic_hexes(self): + mesh = UnstructuredSubMesh( + *_hex_grid( + np.linspace(0, 1e-5, 4), np.linspace(0, 0.2, 4), np.linspace(0, 0.1, 4) + ) + ) + mesh.detect_box_boundaries() + method = _method_with_mesh(mesh) + bcs = {side: (pybamm.Scalar(0), "Dirichlet") for side in mesh.boundary_faces} + L, _ = _laplacian_system(method, mesh, bcs) + # face neighbours only: self + up to 6 in 3D + assert np.diff(L.indptr).max() <= 7 + + def test_genuinely_skewed_faces_are_kept(self): + # structured split: diagonal faces are orthogonal, axis faces skewed + mesh = _make_2d_mesh(3, 3) + _, k = FiniteVolumeUnstructured()._decomposition(mesh) + skew = np.linalg.norm(k, axis=1) + assert skew.max() > 0.4 + assert (skew > 0.4).sum() >= mesh.n_internal_faces // 2