diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py index c26dfed0fa..7c7421c845 100644 --- a/src/pybamm/__init__.py +++ b/src/pybamm/__init__.py @@ -43,6 +43,7 @@ from .expression_tree.independent_variable import * from .expression_tree.independent_variable import t from .expression_tree.vector import Vector +from .expression_tree.vector_field import VectorField from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot from .expression_tree.exceptions import * @@ -139,6 +140,10 @@ SpectralVolume1DSubMesh, SymbolicUniform1DSubMesh, ) +from .meshes.two_dimensional_submeshes import ( + SubMesh2D, + Uniform2DSubMesh, +) from .meshes.scikit_fem_submeshes import ( ScikitSubMesh2D, ScikitUniform2DSubMesh, @@ -154,6 +159,7 @@ from .spatial_methods.spatial_method import SpatialMethod from .spatial_methods.zero_dimensional_method import ZeroDimensionalSpatialMethod from .spatial_methods.finite_volume import FiniteVolume +from .spatial_methods.finite_volume_2d import FiniteVolume2D from .spatial_methods.spectral_volume import SpectralVolume from .spatial_methods.scikit_finite_element import ScikitFiniteElement diff --git a/src/pybamm/discretisations/discretisation.py b/src/pybamm/discretisations/discretisation.py index 5efb67921f..5e5f2c1d2b 100644 --- a/src/pybamm/discretisations/discretisation.py +++ b/src/pybamm/discretisations/discretisation.py @@ -270,27 +270,65 @@ def set_variable_slices(self, variables): end = 0 lower_bounds = [] upper_bounds = [] + + # Make sure we hit concatenation variables first. + # concat_variables = [var for var in variables if isinstance(var, pybamm.ConcatenationVariable)] + # not_concat_variables = [var for var in variables if not isinstance(var, pybamm.ConcatenationVariable)] + # variables = concat_variables + not_concat_variables + # Iterate through unpacked variables, adding appropriate slices to y_slices for variable in variables: + if variable in y_slices: + continue # Add up the size of all the domains in variable.domain if isinstance(variable, pybamm.ConcatenationVariable): - start_ = start spatial_method = self.spatial_methods[variable.domain[0]] + dimension = spatial_method.mesh[variable.domain[0]].dimension + start_ = start children = variable.children meshes = OrderedDict() + lr_points = OrderedDict() + tb_points = OrderedDict() for child in children: meshes[child] = [spatial_method.mesh[dom] for dom in child.domain] + if dimension == 2: + lr_points[child] = sum( + spatial_method.mesh[dom].npts_lr for dom in child.domain + ) + tb_points[child] = sum( + spatial_method.mesh[dom].npts_tb for dom in child.domain + ) sec_points = spatial_method._get_auxiliary_domain_repeats( variable.domains ) for _ in range(sec_points): + start_this_child = start_ for child, mesh in meshes.items(): for domain_mesh in mesh: end += domain_mesh.npts_for_broadcast_to_nodes # Add to slices - y_slices[child].append(slice(start_, end)) - y_slices_explicit[child].append(slice(start_, end)) - # Increment start_ + if dimension == 2: + other_children = set(meshes.keys()) - {child} + num_pts_to_skip = sum( + lr_points[other_child] for other_child in other_children + ) + for row in range(tb_points[child]): + start_this_row = ( + start_this_child + + (lr_points[child] + num_pts_to_skip) * row + ) + end_this_row = start_this_row + lr_points[child] + y_slices[child].append( + slice(start_this_row, end_this_row) + ) + y_slices_explicit[child].append( + slice(start_this_row, end_this_row) + ) + start_this_child += lr_points[child] + else: + y_slices[child].append(slice(start_, end)) + y_slices_explicit[child].append(slice(start_, end)) + # Increment start_ start_ = end else: end += self._get_variable_size(variable) @@ -808,6 +846,22 @@ def _process_symbol(self, symbol): if isinstance(symbol, pybamm.BinaryOperator): # Pre-process children left, right = symbol.children + # Catch case where diffusion is a scalar and turn it into an identity matrix vector field. + if len(symbol.domain) != 0: + spatial_method = self.spatial_methods[symbol.domain[0]] + else: + spatial_method = None + if isinstance(spatial_method, pybamm.FiniteVolume2D): + if isinstance(left, pybamm.Scalar) and ( + isinstance(right, pybamm.VectorField) + or isinstance(right, pybamm.Gradient) + ): + left = pybamm.VectorField(left, left) + elif isinstance(right, pybamm.Scalar) and ( + isinstance(left, pybamm.VectorField) + or isinstance(left, pybamm.Gradient) + ): + right = pybamm.VectorField(right, right) disc_left = self.process_symbol(left) disc_right = self.process_symbol(right) if symbol.domain == []: @@ -875,7 +929,10 @@ def _process_symbol(self, symbol): symbol.integration_variable[0].domain[0] ] out = integral_spatial_method.integral( - child, disc_child, symbol._integration_dimension + child, + disc_child, + symbol._integration_dimension, + symbol.integration_variable, ) out.copy_domains(symbol) return out @@ -915,6 +972,14 @@ def _process_symbol(self, symbol): return child_spatial_method.evaluate_at( symbol, disc_child, symbol.position ) + elif isinstance(symbol, pybamm.UpwindDownwind2D): + return spatial_method.upwind_or_downwind( + child, + disc_child, + self.bcs, + symbol.lr_direction, + symbol.tb_direction, + ) elif isinstance(symbol, pybamm.UpwindDownwind): direction = symbol.name # upwind or downwind return spatial_method.upwind_or_downwind( @@ -923,6 +988,16 @@ def _process_symbol(self, symbol): elif isinstance(symbol, pybamm.NotConstant): # After discretisation, we can make the symbol constant return disc_child + elif isinstance(symbol, pybamm.Magnitude): + if not isinstance(disc_child, pybamm.VectorField): + raise ValueError("Magnitude can only be applied to a vector field") + direction = symbol.direction + if direction == "lr": + return disc_child.lr_field + elif direction == "tb": + return disc_child.tb_field + else: + raise ValueError("Invalid direction") else: return symbol.create_copy(new_children=[disc_child]) @@ -994,6 +1069,11 @@ def _process_symbol(self, symbol): new_symbol = self.process_symbol(symbol.children[0]) return new_symbol + elif isinstance(symbol, pybamm.VectorField): + left_symbol = self.process_symbol(symbol.lr_field) + right_symbol = self.process_symbol(symbol.tb_field) + return symbol.create_copy(new_children=[left_symbol, right_symbol]) + else: # Backup option: return the object return symbol diff --git a/src/pybamm/expression_tree/__init__.py b/src/pybamm/expression_tree/__init__.py index 7ac80e5353..20a68e53a4 100644 --- a/src/pybamm/expression_tree/__init__.py +++ b/src/pybamm/expression_tree/__init__.py @@ -2,4 +2,4 @@ 'concatenations', 'exceptions', 'functions', 'independent_variable', 'input_parameter', 'interpolant', 'matrix', 'operations', 'parameter', 'printing', 'scalar', 'state_vector', 'symbol', - 'unary_operators', 'variable', 'vector', 'discrete_time_sum' ] + 'unary_operators', 'variable', 'vector', 'discrete_time_sum', 'vector_field', 'magnitude' ] diff --git a/src/pybamm/expression_tree/concatenations.py b/src/pybamm/expression_tree/concatenations.py index 5ad32493ca..77556247c0 100644 --- a/src/pybamm/expression_tree/concatenations.py +++ b/src/pybamm/expression_tree/concatenations.py @@ -451,7 +451,7 @@ class SparseStack(Concatenation): The equations to concatenate """ - def __init__(self, *children): + def __init__(self, *children, name="sparse_stack"): children = list(children) if not any(issparse(child.evaluate_for_shape()) for child in children): concatenation_function = np.vstack @@ -459,7 +459,7 @@ def __init__(self, *children): concatenation_function = vstack super().__init__( *children, - name="sparse_stack", + name=name, check_domain=False, concat_fun=concatenation_function, ) diff --git a/src/pybamm/expression_tree/independent_variable.py b/src/pybamm/expression_tree/independent_variable.py index 672b89d322..ce2e0fc673 100644 --- a/src/pybamm/expression_tree/independent_variable.py +++ b/src/pybamm/expression_tree/independent_variable.py @@ -142,12 +142,14 @@ def __init__( auxiliary_domains: AuxiliaryDomainType = None, domains: DomainsType = None, coord_sys=None, + direction=None, ) -> None: self.coord_sys = coord_sys super().__init__( name, domain=domain, auxiliary_domains=auxiliary_domains, domains=domains ) domain = self.domain + self.direction = direction if domain == []: raise ValueError("domain must be provided") diff --git a/src/pybamm/expression_tree/unary_operators.py b/src/pybamm/expression_tree/unary_operators.py index bc7a9168de..f5ac22cdae 100644 --- a/src/pybamm/expression_tree/unary_operators.py +++ b/src/pybamm/expression_tree/unary_operators.py @@ -904,8 +904,17 @@ def __init__(self, child, region="entire"): name += "negative tab" elif region == "positive tab": name += "positive tab" + elif region == "top": + name += "top" + elif region == "bottom": + name += "bottom" + elif region == "left": + name += "left" + elif region == "right": + name += "right" self.region = region super().__init__(name, child, domains) + self.domains = {} def set_id(self): """See :meth:`pybamm.Symbol.set_id()`""" @@ -1187,6 +1196,10 @@ class UpwindDownwind(SpatialOperator): """ def __init__(self, name, child): + self._perform_checks(child) + super().__init__(name, child) + + def _perform_checks(self, child): if child.domain == []: raise pybamm.DomainError( f"Cannot upwind '{child}' since its domain is empty. " @@ -1197,13 +1210,42 @@ def __init__(self, name, child): raise TypeError( f"Cannot upwind '{child}' since it does not " + "evaluate on nodes." ) - super().__init__(name, child) def _evaluates_on_edges(self, dimension: str) -> bool: """See :meth:`pybamm.Symbol._evaluates_on_edges()`.""" return True +class UpwindDownwind2D(UpwindDownwind): + """ + A node in the expression tree representing an upwinding or downwinding operator. + Usually to be used for better stability in convection-dominated equations. + """ + + def __init__(self, child, lr_direction, tb_direction): + super().__init__("upwind_downwind_2d", child) + self.lr_direction = lr_direction + self.tb_direction = tb_direction + + def _unary_new_copy(self, child, perform_simplifications=True): + """See :meth:`UnaryOperator._unary_new_copy()`.""" + return self.__class__(child, self.lr_direction, self.tb_direction) + + +class Magnitude(UnaryOperator): + """ + A node in the expression tree representing the magnitude of a vector field. + """ + + def __init__(self, child, direction): + super().__init__("magnitude" + f"({direction})", child) + self.direction = direction + + def _unary_new_copy(self, child, perform_simplifications=True): + """See :meth:`UnaryOperator._unary_new_copy()`.""" + return self.__class__(child, self.direction) + + class Upwind(UpwindDownwind): """ Upwinding operator. To be used if flow velocity is positive (left to right). diff --git a/src/pybamm/expression_tree/vector_field.py b/src/pybamm/expression_tree/vector_field.py new file mode 100644 index 0000000000..f5b040aabb --- /dev/null +++ b/src/pybamm/expression_tree/vector_field.py @@ -0,0 +1,32 @@ +import pybamm +from typing import Optional + + +class VectorField(pybamm.Symbol): + """ + A node in the expression tree representing a vector field. + """ + + def __init__(self, lr_field, tb_field): + children = [lr_field, tb_field] + if lr_field.domain != tb_field.domain: + raise ValueError("lr_field and tb_field must have the same domain") + super().__init__(name="vector_field", children=children, domain=lr_field.domain) + self.lr_field = lr_field + self.tb_field = tb_field + + def create_copy(self, new_children: Optional[list[pybamm.Symbol]] = None): + if new_children is None: + new_children = [self.lr_field, self.tb_field] + return VectorField(*new_children) + + def _evaluate_for_shape(self): + return self.lr_field.evaluate_for_shape() + + def evaluates_on_edges(self, dimension: str) -> bool: + left_evaluates_on_edges = self.lr_field.evaluates_on_edges(dimension) + right_evaluates_on_edges = self.tb_field.evaluates_on_edges(dimension) + if left_evaluates_on_edges == right_evaluates_on_edges: + return left_evaluates_on_edges + else: + raise ValueError("lr_field and tb_field must have the same domain") diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py index 2f0977d9a6..7212ab6649 100644 --- a/src/pybamm/meshes/meshes.py +++ b/src/pybamm/meshes/meshes.py @@ -4,6 +4,7 @@ import numbers import numpy as np import pybamm +import warnings class Mesh(dict): @@ -74,7 +75,7 @@ def __init__(self, geometry, submesh_types, var_pts): and var.domain[0] in geometry.keys() ): raise KeyError( - f"Points not given for a variable in domain '{domain}'" + f"Points not given for variable '{var.name}' in domain '{domain}'" ) # Otherwise add to the dictionary of submesh points submesh_pts[domain][var.name] = var_name_pts[var.name] @@ -170,8 +171,44 @@ def combine_submeshes(self, *submeshnames): raise ValueError("Submesh domains being combined cannot be empty") # Check that the final edge of each submesh is the same as the first edge of the # next submesh + # TODO: We need a more robust way to check whether the submeshes are being combined lr or tb for i in range(len(submeshnames) - 1): - if self[submeshnames[i]].edges[-1] == self[submeshnames[i + 1]].edges[0]: + if self[submeshnames[i]].dimension != self[submeshnames[i + 1]].dimension: + raise pybamm.DomainError( + "Cannot combine submeshes of different dimensions" + ) + elif self[submeshnames[i]].dimension == 2: + if "left" in submeshnames[i] or "right" in submeshnames[i + 1]: + # Make sure that the lr edges are aligned + if ( + self[submeshnames[i]].edges_lr[-1] + != self[submeshnames[i + 1]].edges_lr[0] + ): + raise pybamm.DomainError("lr edges are not aligned") + elif ( + self[submeshnames[i]].edges_tb + != self[submeshnames[i + 1]].edges_tb + ).any(): + raise pybamm.DomainError("tb edges are not aligned") + else: + pass + + elif "bottom" in submeshnames[i] or "top" in submeshnames[i + 1]: + # Make sure that the tb edges are aligned + if ( + self[submeshnames[i]].edges_tb[-1] + != self[submeshnames[i + 1]].edges_tb[0] + ): + raise pybamm.DomainError("tb edges are not aligned") + elif ( + self[submeshnames[i]].edges_lr + != self[submeshnames[i + 1]].edges_lr + ).any(): + raise pybamm.DomainError("lr edges are not aligned") + else: + pass + pass + elif self[submeshnames[i]].edges[-1] == self[submeshnames[i + 1]].edges[0]: # submeshes are aligned, all good pass elif hasattr(self[submeshnames[i]], "min") or hasattr( @@ -189,12 +226,51 @@ def combine_submeshes(self, *submeshnames): raise pybamm.DomainError( "trying to combine two meshes in different coordinate systems" ) - combined_submesh_edges = np.concatenate( - [self[submeshnames[0]].edges] - + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]] - ) + coord_sys = self[submeshnames[0]].coord_sys - submesh = pybamm.SubMesh1D(combined_submesh_edges, coord_sys) + if self[submeshnames[0]].dimension == 1: + combined_submesh_edges = np.concatenate( + [self[submeshnames[0]].edges] + + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]] + ) + submesh = pybamm.SubMesh1D(combined_submesh_edges, coord_sys) + elif self[submeshnames[0]].dimension == 2: + # If it's an lr concatenation, then we only need to concatenate the edges_lr + if "left" in submeshnames[0] or "right" in submeshnames[-1]: + combined_submesh_edges_lr = np.concatenate( + [self[submeshnames[0]].edges_lr] + + [ + self[submeshname].edges_lr[1:] + for submeshname in submeshnames[1:] + ] + ) + combined_submesh_edges_tb = self[submeshnames[0]].edges_tb + elif "top" in submeshnames[0] or "bottom" in submeshnames[-1]: + combined_submesh_edges_tb = np.concatenate( + [self[submeshnames[0]].edges_tb] + + [ + self[submeshname].edges_tb[1:] + for submeshname in submeshnames[1:] + ] + ) + combined_submesh_edges_lr = self[submeshnames[0]].edges_lr + else: + warnings.warn( + "Could not determine how to combine submeshes. Assuming lr concatenation.", + stacklevel=2, + ) + combined_submesh_edges_lr = np.concatenate( + [self[submeshnames[0]].edges_lr] + + [ + self[submeshname].edges_lr[1:] + for submeshname in submeshnames[1:] + ] + ) + combined_submesh_edges_tb = self[submeshnames[0]].edges_tb + submesh = pybamm.SubMesh2D( + combined_submesh_edges_lr, combined_submesh_edges_tb, coord_sys + ) + if getattr(self[submeshnames[0]], "length", None) is not None: # Assume that the ghost cells have the same length as the first submesh if any("ghost" in submeshname for submeshname in submeshnames): @@ -209,12 +285,32 @@ def combine_submeshes(self, *submeshnames): submesh.length = submesh_length submesh.min = submesh_min # add in internal boundaries - for submeshname in submeshnames[1:]: + for i, submeshname in enumerate(submeshnames[1:]): + i = i + 1 if getattr(self[submeshname], "length", None) is not None: min = self[submeshname].min else: min = 0 - submesh.internal_boundaries.append(self[submeshname].edges[0] + min) + if submesh.dimension == 1: + submesh.internal_boundaries.append(self[submeshname].edges[0] + min) + elif ( + "left" in submeshname + or "left" in submeshnames[i - 1] + or "right" in submeshname + ): + submesh.internal_boundaries.append(self[submeshname].edges_lr[0] + min) + elif ( + "top" in submeshname + or "top" in submeshnames[i - 1] + or "bottom" in submeshname + ): + submesh.internal_boundaries.append(self[submeshname].edges_tb[0] + min) + else: + warnings.warn( + "Could not determine how to combine submeshes. Assuming lr concatenation.", + stacklevel=2, + ) + submesh.internal_boundaries.append(self[submeshname].edges_lr[0] + min) return submesh def add_ghost_meshes(self): @@ -224,6 +320,7 @@ def add_ghost_meshes(self): This will be useful for calculating the gradient with Dirichlet BCs. """ # Get all submeshes relating to space (i.e. exclude time) + # Ignore 2D submeshes for now... submeshes = [ (domain, submesh) for domain, submesh in self.items() @@ -233,24 +330,25 @@ def add_ghost_meshes(self): ) ] for domain, submesh in submeshes: - edges = submesh.edges - - # left ghost cell: two edges, one node, to the left of existing submesh - lgs_edges = np.array([2 * edges[0] - edges[1], edges[0]]) - lgs_submesh = pybamm.SubMesh1D(lgs_edges, submesh.coord_sys) - if getattr(submesh, "length", None) is not None: - lgs_submesh.length = submesh.length - lgs_submesh.min = submesh.min - self[domain[0] + "_left ghost cell"] = lgs_submesh - - # right ghost cell: two edges, one node, to the right of - # existing submesh - rgs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]]) - rgs_submesh = pybamm.SubMesh1D(rgs_edges, submesh.coord_sys) - if getattr(submesh, "length", None) is not None: - rgs_submesh.length = submesh.length - rgs_submesh.min = submesh.min - self[domain[0] + "_right ghost cell"] = rgs_submesh + if submesh.dimension == 1: + self[domain[0] + "_left ghost cell"] = submesh.create_ghost_cell("left") + self[domain[0] + "_right ghost cell"] = submesh.create_ghost_cell( + "right" + ) + elif submesh.dimension == 2: + self[domain[0] + "_left ghost cell"] = submesh.create_ghost_cell("left") + self[domain[0] + "_right ghost cell"] = submesh.create_ghost_cell( + "right" + ) + self[domain[0] + "_bottom ghost cell"] = submesh.create_ghost_cell( + "bottom" + ) + self[domain[0] + "_top ghost cell"] = submesh.create_ghost_cell("top") + else: + raise NotImplementedError( + "ghost cells not implemented for submeshes of dimension " + + str(submesh.dimension) + ) @property def geometry(self): diff --git a/src/pybamm/meshes/one_dimensional_submeshes.py b/src/pybamm/meshes/one_dimensional_submeshes.py index 26d8f1293e..2cf1ddfb6f 100644 --- a/src/pybamm/meshes/one_dimensional_submeshes.py +++ b/src/pybamm/meshes/one_dimensional_submeshes.py @@ -32,6 +32,7 @@ def __init__(self, edges, coord_sys, tabs=None): self.npts = self.nodes.size self.coord_sys = coord_sys self.internal_boundaries = [] + self.dimension = 1 # Add tab locations in terms of "left" and "right" if tabs and "negative tab" not in tabs.keys(): @@ -84,6 +85,22 @@ def to_json(self): return json_dict + def create_ghost_cell(self, side): + edges = self.edges + if side == "left": + gs_edges = np.array([2 * edges[0] - edges[1], edges[0]]) + elif side == "right": + gs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]]) + else: + raise NotImplementedError( + "Only left and right ghost cells are implemented for 1D submeshes" + ) + gs_submesh = pybamm.SubMesh1D(gs_edges, self.coord_sys) + if hasattr(self, "length") and getattr(self, "length", None) is not None: + gs_submesh.length = self.length + gs_submesh.min = self.min + return gs_submesh + class SymbolicUniform1DSubMesh(SubMesh1D): def __init__(self, lims, npts, tabs=None): @@ -107,6 +124,7 @@ def __init__(self, lims, npts, tabs=None): self.npts = self.nodes.size self.coord_sys = coord_sys self.internal_boundaries = [] + self.dimension = 1 class Uniform1DSubMesh(SubMesh1D): diff --git a/src/pybamm/meshes/two_dimensional_submeshes.py b/src/pybamm/meshes/two_dimensional_submeshes.py new file mode 100644 index 0000000000..61d694bd46 --- /dev/null +++ b/src/pybamm/meshes/two_dimensional_submeshes.py @@ -0,0 +1,156 @@ +# +# Two-dimensional submeshes +# +import pybamm +from .meshes import SubMesh + +import numpy as np + + +class SubMesh2D(SubMesh): + """ + 2D submesh class. + Contains the position of the nodes, the number of mesh points, and + (optionally) information about the tab locations. + + Parameters + ---------- + edges_lr : array_like + An array containing the points corresponding to the edges of the submesh + in the left-right direction + edges_tb : array_like + An array containing the points corresponding to the edges of the submesh + in the top-bottom direction + coord_sys_lr : string + The coordinate system of the submesh in the left-right direction + coord_sys_tb : string + The coordinate system of the submesh in the top-bottom direction + tabs : dict, optional + A dictionary that contains information about the size and location of + the tabs + """ + + def __init__(self, edges_lr, edges_tb, coord_sys, tabs=None): + self.edges_lr = edges_lr + self.edges_tb = edges_tb + self.nodes_lr = (self.edges_lr[1:] + self.edges_lr[:-1]) / 2 + self.nodes_tb = (self.edges_tb[1:] + self.edges_tb[:-1]) / 2 + self.d_edges_lr = np.diff(self.edges_lr) + self.d_edges_tb = np.diff(self.edges_tb) + self.d_nodes_lr = np.diff(self.nodes_lr) + self.d_nodes_tb = np.diff(self.nodes_tb) + self.npts_lr = self.nodes_lr.size + self.npts_tb = self.nodes_tb.size + self.npts = self.npts_lr * self.npts_tb + self.coord_sys = coord_sys + self.internal_boundaries = [] + self.dimension = 2 + + # Add tab locations in terms of "left" and "right" + self.tabs = tabs + + def read_lims(self, lims): + # Separate limits and tabs + # Read and remove tabs. If "tabs" is not a key in "lims", then tabs is set to + # "None" and nothing is removed from lims + tabs = lims.pop("tabs", None) + + # check that only one variable passed in + if len(lims) != 2: + raise pybamm.GeometryError("lims should only contain two variables") + + ((spatial_var_lr, spatial_lims_lr), (spatial_var_tb, spatial_lims_tb)) = ( + lims.items() + ) + + if isinstance(spatial_var_lr, str): + spatial_var_lr = getattr(pybamm.standard_spatial_vars, spatial_var_lr) + if isinstance(spatial_var_tb, str): + spatial_var_tb = getattr(pybamm.standard_spatial_vars, spatial_var_tb) + + return spatial_var_lr, spatial_lims_lr, spatial_var_tb, spatial_lims_tb, tabs + + def to_json(self): + json_dict = { + "edges_lr": self.edges_lr.tolist(), + "edges_tb": self.edges_tb.tolist(), + "coord_sys_lr": self.coord_sys_lr, + "coord_sys_tb": self.coord_sys_tb, + } + + if hasattr(self, "tabs"): + json_dict["tabs"] = self.tabs + + return json_dict + + def create_ghost_cell(self, side): + if side == "left": + gs_edges_lr = np.array( + [2 * self.edges_lr[0] - self.edges_lr[1], self.edges_lr[0]] + ) + gs_edges_tb = self.edges_tb + elif side == "right": + gs_edges_lr = np.array( + [self.edges_lr[-1], 2 * self.edges_lr[-1] - self.edges_lr[-2]] + ) + gs_edges_tb = self.edges_tb + elif side == "top": + gs_edges_lr = self.edges_lr + gs_edges_tb = np.array( + [2 * self.edges_tb[0] - self.edges_tb[1], self.edges_tb[0]] + ) + elif side == "bottom": + gs_edges_lr = self.edges_lr + gs_edges_tb = np.array( + [self.edges_tb[-1], 2 * self.edges_tb[-1] - self.edges_tb[-2]] + ) + else: + raise ValueError(f"Invalid side: {side}") + gs_submesh = pybamm.SubMesh2D(gs_edges_lr, gs_edges_tb, self.coord_sys) + return gs_submesh + + +class Uniform2DSubMesh(SubMesh2D): + """ + A 2D submesh with uniform spacing in both dimensions + """ + + def __init__(self, lims, npts): + spatial_var_lr, spatial_lims_lr, spatial_var_tb, spatial_lims_tb, tabs = ( + self.read_lims(lims) + ) + npts_lr = npts[spatial_var_lr.name] + npts_tb = npts[spatial_var_tb.name] + + edges_lr = np.linspace( + spatial_lims_lr["min"], spatial_lims_lr["max"], npts_lr + 1 + ) + edges_tb = np.linspace( + spatial_lims_tb["min"], spatial_lims_tb["max"], npts_tb + 1 + ) + if spatial_var_lr.coord_sys != spatial_var_tb.coord_sys: + raise pybamm.GeometryError( + "Coordinate systems must be the same for 2D submeshes" + ) + coord_sys = spatial_var_lr.coord_sys + + super().__init__(edges_lr, edges_tb, coord_sys, tabs=tabs) + + def read_lims(self, lims): + # Separate limits and tabs + # Read and remove tabs. If "tabs" is not a key in "lims", then tabs is set to + # "None" and nothing is removed from lims + tabs = lims.pop("tabs", None) + + # check that only one variable passed in + if len(lims) != 2: + raise pybamm.GeometryError("lims should only contain two variables") + + ((spatial_var1, spatial_lims1), (spatial_var2, spatial_lims2)) = lims.items() + + if isinstance(spatial_var1, str): + spatial_var1 = getattr(pybamm.standard_spatial_vars, spatial_var1) + if isinstance(spatial_var2, str): + spatial_var2 = getattr(pybamm.standard_spatial_vars, spatial_var2) + + return spatial_var1, spatial_lims1, spatial_var2, spatial_lims2, tabs diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index ba8c6c324e..f86ca14d97 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/src/pybamm/parameters/parameter_values.py @@ -569,7 +569,15 @@ def process_boundary_conditions(self, model): "positive tab": pos. tab bc "no tab": no tab bc}. """ new_boundary_conditions = {} - sides = ["left", "right", "negative tab", "positive tab", "no tab"] + sides = [ + "left", + "right", + "negative tab", + "positive tab", + "no tab", + "top", + "bottom", + ] for variable, bcs in model.boundary_conditions.items(): processed_variable = self.process_symbol(variable) new_boundary_conditions[processed_variable] = {} @@ -787,6 +795,11 @@ def _process_symbol(self, symbol): new_children = [self.process_symbol(child) for child in symbol.children] return symbol.create_copy(new_children) + elif isinstance(symbol, pybamm.VectorField): + left_symbol = self.process_symbol(symbol.lr_field) + right_symbol = self.process_symbol(symbol.tb_field) + return symbol.create_copy(new_children=[left_symbol, right_symbol]) + # Variables: update scale elif isinstance(symbol, pybamm.Variable): new_symbol = symbol.create_copy() diff --git a/src/pybamm/solvers/processed_variable.py b/src/pybamm/solvers/processed_variable.py index b49a5fec64..d9c55a1eca 100644 --- a/src/pybamm/solvers/processed_variable.py +++ b/src/pybamm/solvers/processed_variable.py @@ -807,6 +807,23 @@ def _observe_postfix(self, entries, t): return entries +class ProcessedVariable2DReal(ProcessedVariable): + def _shape(self, t): + return [self.mesh.npts, len(t)] + + def initialise(self): + if self.entries_raw_initialized: + return + + entries = self.observe_raw() + + # entries_for_interp, coords = self._interp_setup(entries, t) + + self._entries_raw = entries + # self._entries_for_interp_raw = entries_for_interp + # self._coords_raw = coords + + class ProcessedVariable3D(ProcessedVariable): """ An object that can be evaluated at arbitrary (scalars or vectors) t and x, and @@ -1100,10 +1117,15 @@ def process_variable(base_variables, *args, **kwargs): ): return ProcessedVariable3DSciKitFEM(base_variables, *args, **kwargs) + if mesh and hasattr(mesh, "edges_lr") and hasattr(mesh, "edges_tb"): + return ProcessedVariable2DReal(base_variables, *args, **kwargs) + # check variable shape if len(base_eval_shape) == 0 or base_eval_shape[0] == 1: return ProcessedVariable0D(base_variables, *args, **kwargs) + if mesh is None: + return ProcessedVariable2DReal(base_variables, *args, **kwargs) n = mesh.npts base_shape = base_eval_shape[0] # Try some shapes that could make the variable a 1D variable diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py index 98ec1b3467..066b2a7b19 100644 --- a/src/pybamm/solvers/solution.py +++ b/src/pybamm/solvers/solution.py @@ -459,6 +459,21 @@ def _update_variable(self, variable): vars_pybamm[i] = var_pybamm elif variable in model._variables_casadi: var_casadi = model._variables_casadi[variable] + elif isinstance(var_pybamm, pybamm.VectorField): + var_casadi_lr = self.process_casadi_var( + var_pybamm.lr_field, + inputs, + ys.shape, + ) + var_casadi_tb = self.process_casadi_var( + var_pybamm.tb_field, + inputs, + ys.shape, + ) + var_casadi = { + "lr": var_casadi_lr, + "tb": var_casadi_tb, + } else: var_casadi = self.process_casadi_var( var_pybamm, diff --git a/src/pybamm/spatial_methods/__init__.py b/src/pybamm/spatial_methods/__init__.py index 56b9988536..fc070c315c 100644 --- a/src/pybamm/spatial_methods/__init__.py +++ b/src/pybamm/spatial_methods/__init__.py @@ -1,2 +1,2 @@ __all__ = ['finite_volume', 'scikit_finite_element', 'spatial_method', - 'spectral_volume', 'zero_dimensional_method'] + 'spectral_volume', 'zero_dimensional_method', 'finite_volume_2d'] diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index dff5728f25..f54094dd8d 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -226,7 +226,9 @@ def laplacian(self, symbol, discretised_symbol, boundary_conditions): grad = self.gradient(symbol, discretised_symbol, boundary_conditions) return self.divergence(grad, grad, boundary_conditions) - def integral(self, child, discretised_child, integration_dimension): + def integral( + self, child, discretised_child, integration_dimension, integration_variable + ): """Vector-vector dot product to implement the integral operator.""" integration_vector = self.definite_integral_matrix( child, integration_dimension=integration_dimension diff --git a/src/pybamm/spatial_methods/finite_volume_2d.py b/src/pybamm/spatial_methods/finite_volume_2d.py new file mode 100644 index 0000000000..4ae540ec76 --- /dev/null +++ b/src/pybamm/spatial_methods/finite_volume_2d.py @@ -0,0 +1,2066 @@ +# +# Finite Volume discretisation class +# +import pybamm + +from scipy.sparse import ( + diags, + spdiags, + eye, + kron, + csr_matrix, + vstack, + hstack, + coo_matrix, + block_diag, +) +import numpy as np + + +class FiniteVolume2D(pybamm.SpatialMethod): + """ + A class which implements the steps specific to the finite volume method during + discretisation. + + For broadcast and mass_matrix, we follow the default behaviour from SpatialMethod. + + Parameters + ---------- + options : dict-like, optional + A dictionary of options to be passed to the spatial method. The only option + currently available is "extrapolation", which has options for "order" and "use_bcs". + It sets the order separately for `pybamm.BoundaryValue` and `pybamm.BoundaryGradient`. + Default is "linear" for the value and quadratic for the gradient. + """ + + def __init__(self, options=None): + super().__init__(options) + + def build(self, mesh): + super().build(mesh) + + # add npts_for_broadcast to mesh domains for this particular discretisation + for dom in mesh.keys(): + mesh[dom].npts_for_broadcast_to_nodes = mesh[dom].npts + + def spatial_variable(self, symbol): + """ + Creates a discretised spatial variable compatible with + the FiniteVolume method. + + Parameters + ---------- + symbol : :class:`pybamm.SpatialVariable` + The spatial variable to be discretised. + + Returns + ------- + :class:`pybamm.Vector` + Contains the discretised spatial variable + """ + symbol_mesh = self.mesh[symbol.domain] + symbol_direction = symbol.direction + if symbol_mesh.dimension != 2: + raise ValueError(f"Spatial variable {symbol} is not in 2D") + repeats = self._get_auxiliary_domain_repeats(symbol.domains) + # Vector should be size npts_lr x npts_tb + # Do LR first, then TB + if symbol.evaluates_on_edges("primary"): + LR, TB = np.meshgrid(symbol_mesh.edges_lr, symbol_mesh.edges_tb) + lr = LR.flatten() + tb = TB.flatten() + else: + LR, TB = np.meshgrid(symbol_mesh.nodes_lr, symbol_mesh.nodes_tb) + lr = LR.flatten() + tb = TB.flatten() + if symbol_direction == "lr": + entries = np.tile(lr, repeats) + elif symbol_direction == "tb": + entries = np.tile(tb, repeats) + return pybamm.Vector(entries, domains=symbol.domains) + + def _gradient(self, symbol, discretised_symbol, boundary_conditions, direction): + """ + Gradient with a specific direction (lr or tb) + """ + domain = symbol.domain + + # Add Dirichlet boundary conditions, if defined + if direction == "lr": + relevant_bcs = ["left", "right"] + elif direction == "tb": + relevant_bcs = ["top", "bottom"] + else: + raise ValueError(f"Direction {direction} not supported") + + if symbol in boundary_conditions: + bcs = { + key: boundary_conditions[symbol][key] + for key in relevant_bcs + if key in boundary_conditions[symbol] + } + if any(bc[1] == "Dirichlet" for bc in bcs.values()): + # add ghost nodes and update domain + discretised_symbol, domain = self.add_ghost_nodes( + symbol, discretised_symbol, bcs + ) + gradient_matrix = self.gradient_matrix(domain, symbol.domains, direction) + + grad = gradient_matrix @ discretised_symbol + grad.copy_domains(symbol) + + # Add Neumann boundary conditions, if defined + if symbol in boundary_conditions: + bcs = { + key: boundary_conditions[symbol][key] + for key in relevant_bcs + if key in boundary_conditions[symbol] + } + if any(bc[1] == "Neumann" for bc in bcs.values()): + grad = self.add_neumann_values(symbol, grad, bcs, domain) + + return grad + + def gradient(self, symbol, discretised_symbol, boundary_conditions): + """Matrix-vector multiplication to implement the gradient operator. + See :meth:`pybamm.SpatialMethod.gradient` + """ + # Multiply by gradient matrix + grad_lr = self._gradient(symbol, discretised_symbol, boundary_conditions, "lr") + grad_tb = self._gradient(symbol, discretised_symbol, boundary_conditions, "tb") + grad = pybamm.VectorField(grad_lr, grad_tb) + return grad + + def gradient_matrix(self, domain, domains, direction): + """ + Gradient matrix for finite volumes in the appropriate domain. + Equivalent to grad(y) = (y[1:] - y[:-1])/dx + + Parameters + ---------- + domains : list + The domain in which to compute the gradient matrix, including ghost nodes + + Returns + ------- + :class:`pybamm.Matrix` + The (sparse) finite volume gradient matrix for the domain + """ + # Create appropriate submesh by combining submeshes in primary domain + submesh = self.mesh[domain] + + # Create matrix using submesh + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + e_lr = 1 / submesh.d_nodes_lr + e_tb = 1 / submesh.d_nodes_tb + if direction == "lr": + sub_matrix = diags([-e_lr, e_lr], [0, 1], shape=((n_lr - 1), n_lr)) + sub_matrix = block_diag((sub_matrix,) * n_tb) + elif direction == "tb": + e_tb = np.repeat(e_tb, n_lr) + sub_matrix = diags( + [-e_tb, e_tb], [0, n_lr], shape=(n_lr * (n_tb - 1), n_lr * n_tb) + ) + + # number of repeats + second_dim_repeats = self._get_auxiliary_domain_repeats(domains) + + # generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), which is + # not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should not be an + # issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + return pybamm.Matrix(matrix) + + def divergence(self, symbol, discretised_symbol, boundary_conditions): + """Matrix-vector multiplication to implement the divergence operator. + See :meth:`pybamm.SpatialMethod.divergence` + """ + + divergence_matrix_lr = self.divergence_matrix(symbol.domains, "lr") + divergence_matrix_tb = self.divergence_matrix(symbol.domains, "tb") + + grad_lr = discretised_symbol.lr_field + grad_tb = discretised_symbol.tb_field + + div_lr = divergence_matrix_lr @ grad_lr + div_tb = divergence_matrix_tb @ grad_tb + + out = div_lr + div_tb + + return out + + def divergence_matrix(self, domains, direction): + """ + Divergence with a specific direction (lr or tb) + """ + submesh = self.mesh[domains["primary"]] + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + e_lr = 1 / submesh.d_edges_lr + e_tb = 1 / submesh.d_edges_tb + if direction == "lr": + sub_matrix = diags([-e_lr, e_lr], [0, 1], shape=(n_lr, n_lr + 1)) + sub_matrix = block_diag((sub_matrix,) * n_tb) + elif direction == "tb": + e_tb = np.repeat(e_tb, n_lr + 1) + sub_matrix = diags( + [-e_tb, e_tb], [0, n_lr], shape=(n_lr * n_tb, n_lr * (n_tb + 1)) + ) + return pybamm.Matrix(sub_matrix) + + def integral( + self, child, discretised_child, integration_dimension, integration_variable + ): + """matrix-vector product to implement the integral operator.""" + integration_matrix = self.definite_integral_matrix( + child, + integration_dimension=integration_dimension, + integration_variable=integration_variable, + ) + if len(integration_variable) > 1: + dir_1 = integration_variable[0].direction + dir_2 = integration_variable[1].direction + if dir_1 == dir_2: + raise ValueError( + "Integration variables must be in different directions" + ) + else: + one_dimensional_matrix = self.one_dimensional_integral_matrix( + child, dir_2 + ) + integration_matrix = one_dimensional_matrix @ integration_matrix + else: + pass + domains = child.domains + second_dim_repeats = self._get_auxiliary_domain_repeats(domains) + integration_matrix = kron(eye(second_dim_repeats), integration_matrix) + return pybamm.Matrix(integration_matrix) @ discretised_child + + def laplacian(self, symbol, discretised_symbol, boundary_conditions): + """ + Laplacian operator, implemented as div(grad(.)) + See :meth:`pybamm.SpatialMethod.laplacian` + """ + grad = self.gradient(symbol, discretised_symbol, boundary_conditions) + return self.divergence(grad, grad, boundary_conditions) + + def definite_integral_matrix( + self, + child, + vector_type="row", + integration_dimension="primary", + integration_variable=None, + ): + if integration_variable is None: + raise ValueError("Integration variable must be provided for 2D integration") + else: + integration_direction = integration_variable[0].direction + + domains = child.domains + domain = child.domains[integration_dimension] + submesh = self.mesh[domain] + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + + if integration_dimension == "primary": + # Create appropriate submesh by combining submeshes in domain + submesh = self.mesh[domains["primary"]] + + # Create vector of ones for primary domain submesh + if integration_direction == "lr": + d_edges = submesh.d_edges_lr + cols_list = [] + rows_list = [] + for n in range(n_tb): + cols = np.arange(n * n_lr, (n + 1) * n_lr) + rows = np.ones(n_lr) * n + cols_list.append(cols) + rows_list.append(rows) + cols = np.concatenate(cols_list) + rows = np.concatenate(rows_list) + sub_matrix = csr_matrix( + (np.tile(d_edges, n_tb), (rows, cols)), shape=(n_tb, n_lr * n_tb) + ) + elif integration_direction == "tb": + d_edges = submesh.d_edges_tb + cols_list = [] + rows_list = [] + for n in range(n_tb): + rows = np.arange(0, n_lr) + cols = np.arange(n * n_lr, (n + 1) * n_lr) + cols_list.append(cols) + rows_list.append(rows) + cols = np.concatenate(cols_list) + rows = np.concatenate(rows_list) + sub_matrix = csr_matrix( + (np.tile(d_edges, n_lr), (rows, cols)), shape=(n_lr, n_lr * n_tb) + ) + + # repeat matrix for each node in secondary dimensions + else: + raise NotImplementedError( + "Only primary integration dimension is implemented for 2D integration" + ) + + return sub_matrix + + def one_dimensional_integral_matrix(self, child, direction): + """ + One-dimensional integral matrix for finite volumes in the appropriate domain. + Equivalent to int(y) = sum(y[i] * dx[i]) + """ + submesh = self.mesh[child.domain] + domains = child.domains + + # Create vector of ones for primary domain submesh + if direction == "lr": + d_edges = submesh.d_edges_lr + elif direction == "tb": + d_edges = submesh.d_edges_tb + + # repeat matrix for each node in secondary dimensions + second_dim_repeats = self._get_auxiliary_domain_repeats(domains) + # generate full matrix from the submatrix + matrix = kron(eye(second_dim_repeats), d_edges) + + return matrix + + def boundary_integral(self, child, discretised_child, region): + """ + Boundary integral operator, implemented as int(grad(.)) + See :meth:`pybamm.SpatialMethod.boundary_integral` + """ + symbol = pybamm.BoundaryValue(child, region) + boundary_value = self.boundary_value_or_flux(symbol, discretised_child) + if region == "left" or region == "right": + direction = "lr" + elif region == "top" or region == "bottom": + direction = "tb" + else: + raise ValueError(f"Region {region} not supported") + + integral_matrix = self.one_dimensional_integral_matrix(child, direction) + domains = child.domains + second_dim_repeats = self._get_auxiliary_domain_repeats(domains) + integral_matrix = kron(eye(second_dim_repeats), integral_matrix) + return pybamm.Matrix(integral_matrix) @ boundary_value + + def indefinite_integral(self, child, discretised_child, direction): + raise NotImplementedError + + def indefinite_integral_matrix_edges(self, domains, direction): + raise NotImplementedError + + def indefinite_integral_matrix_nodes(self, domains, direction): + raise NotImplementedError + + def delta_function(self, symbol, discretised_symbol): + raise NotImplementedError + + def internal_neumann_condition( + self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh + ): + """ + A method to find the internal Neumann conditions between two symbols + on adjacent subdomains. + + Parameters + ---------- + left_symbol_disc : :class:`pybamm.Symbol` + The discretised symbol on the left subdomain + right_symbol_disc : :class:`pybamm.Symbol` + The discretised symbol on the right subdomain + left_mesh : list + The mesh on the left subdomain + right_mesh : list + The mesh on the right subdomain + """ + + left_npts = left_mesh.npts + left_npts_lr = left_mesh.npts_lr + right_npts = right_mesh.npts + right_npts_lr = right_mesh.npts_lr + + second_dim_repeats = self._get_auxiliary_domain_repeats( + left_symbol_disc.domains + ) + + if second_dim_repeats != self._get_auxiliary_domain_repeats( + right_symbol_disc.domains + ): + raise pybamm.DomainError( + """Number of secondary points in subdomains do not match""" + ) + + left_sub_matrix = np.zeros((1, left_npts)) + left_sub_matrix[0][left_npts_lr - 1 : left_npts_lr : left_npts] = 1 + + left_matrix = pybamm.Matrix( + csr_matrix(kron(eye(second_dim_repeats), left_sub_matrix)) + ) + + right_sub_matrix = np.zeros((1, right_npts)) + right_sub_matrix[0][0:right_npts_lr:right_npts] = 1 + right_matrix = pybamm.Matrix( + csr_matrix(kron(eye(second_dim_repeats), right_sub_matrix)) + ) + + # Finite volume derivative + # Remove domains to avoid clash + right_mesh_x = right_mesh.nodes_lr[0] + left_mesh_x = left_mesh.nodes_lr[-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 + + def add_ghost_nodes(self, symbol, discretised_symbol, bcs): + """ + Add ghost nodes to a symbol. + + For Dirichlet bcs, for a boundary condition "y = a at the left-hand boundary", + we concatenate a ghost node to the start of the vector y with value "2*a - y1" + where y1 is the value of the first node. + Similarly for the right-hand boundary condition and top and bottom boundaries. + + For Neumann bcs no ghost nodes are added. Instead, the exact value provided + by the boundary condition is used at the cell edge when calculating the + gradient (see :meth:`pybamm.FiniteVolume2D.add_neumann_values`). + + Parameters + ---------- + symbol : :class:`pybamm.SpatialVariable` + The variable to be discretised + discretised_symbol : :class:`pybamm.Vector` + Contains the discretised variable + bcs : dict of tuples (:class:`pybamm.Scalar`, str) + Dictionary (with keys "left", "right", "top", "bottom") of boundary conditions. Each + boundary condition consists of a value and a flag indicating its type + (e.g. "Dirichlet") + + Returns + ------- + :class:`pybamm.Symbol` + `Matrix @ discretised_symbol + bcs_vector`. When evaluated, this gives the + discretised_symbol, with appropriate ghost nodes concatenated at each end. + + """ + domain = symbol.domain + submesh = self.mesh[domain] + + # Prepare sizes and empty bcs_vector + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + n = submesh.npts + second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + # Catch if no boundary conditions are defined + if ( + "left" not in bcs.keys() + and "right" not in bcs.keys() + and "top" not in bcs.keys() + and "bottom" not in bcs.keys() + ): + raise ValueError(f"No boundary conditions have been provided for {symbol}") + + # Allow to only pass one boundary condition (for upwind/downwind) + lbc_value, lbc_type = bcs.get("left", (None, None)) + rbc_value, rbc_type = bcs.get("right", (None, None)) + tbc_value, tbc_type = bcs.get("top", (None, None)) + bbc_value, bbc_type = bcs.get("bottom", (None, None)) + + # Add ghost node(s) to domain where necessary and count number of + # Dirichlet boundary conditions + # [left, top, n, bottom, right] + n_bcs = 0 + + if tbc_type == "Dirichlet" and bbc_type != "Dirichlet": + if isinstance(domain, list) or isinstance(domain, tuple): + domain = [(d + "_top ghost cell", d) for d in domain] + else: + domain = [domain + "_top ghost cell", domain] + n_bcs += 1 + elif tbc_type != "Dirichlet" and bbc_type == "Dirichlet": + if isinstance(domain, list) or isinstance(domain, tuple): + domain = [(d, d + "_bottom ghost cell") for d in domain] + else: + domain = [domain, domain + "_bottom ghost cell"] + n_bcs += 1 + elif tbc_type == "Dirichlet" and bbc_type == "Dirichlet": + if isinstance(domain, list) or isinstance(domain, tuple): + domain = [ + (d + "_top ghost cell", d, d + "_bottom ghost cell") for d in domain + ] + else: + domain = [ + domain + "_top ghost cell", + domain, + domain + "_bottom ghost cell", + ] + n_bcs += 2 + + if lbc_type == "Dirichlet": + domain = [domain[0] + "_left ghost cell", *domain] + n_bcs += 1 + if rbc_type == "Dirichlet": + domain = [*domain, domain[-1] + "_right ghost cell"] + n_bcs += 1 + + # Calculate values for ghost nodes for any Dirichlet boundary conditions + if lbc_type == "Dirichlet": + # Create matrix to extract the leftmost column of values + lbc_sub_matrix = coo_matrix(([1], ([0], [0])), shape=(n_lr + n_bcs, 1)) + lbc_matrix = csr_matrix(kron(eye(second_dim_repeats), lbc_sub_matrix)) + lbc_matrix = vstack( + [ + lbc_matrix, + ] + * n_tb + ) + if lbc_value.evaluates_to_number(): + left_ghost_constant = ( + 2 * lbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + ) + else: + left_ghost_constant = 2 * lbc_value + + lbc_vector = pybamm.Matrix(lbc_matrix) @ left_ghost_constant + elif lbc_type == "Neumann": + lbc_vector = pybamm.Vector( + np.zeros((n_lr + n_bcs) * second_dim_repeats * n_tb) + ) + else: + lbc_vector = pybamm.Vector( + np.zeros((n_tb + n_bcs) * second_dim_repeats * n_lr) + ) + + # Calculate values for ghost nodes for any Dirichlet boundary conditions + if rbc_type == "Dirichlet": + # Create matrix to extract the leftmost column of values + rbc_sub_matrix = coo_matrix( + ([1], ([n_lr + n_bcs - 1], [0])), shape=(n_lr + n_bcs, 1) + ) + rbc_matrix = csr_matrix(kron(eye(second_dim_repeats), rbc_sub_matrix)) + rbc_matrix = vstack( + [ + rbc_matrix, + ] + * n_tb + ) + if rbc_value.evaluates_to_number(): + right_ghost_constant = ( + 2 * rbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + ) + else: + right_ghost_constant = 2 * rbc_value + rbc_vector = pybamm.Matrix(rbc_matrix) @ right_ghost_constant + elif rbc_type == "Neumann": + rbc_vector = pybamm.Vector( + np.zeros((n_lr + n_bcs) * second_dim_repeats * n_tb) + ) + else: + rbc_vector = pybamm.Vector( + np.zeros((n_tb + n_bcs) * second_dim_repeats * n_lr) + ) + + # Calculate values for ghost nodes for any Dirichlet boundary conditions + if tbc_type == "Dirichlet": + # Create matrix to extract the leftmost column of values + row_indices = np.arange(0, n_lr) + col_indices = np.zeros(len(row_indices)) + vals = np.ones(len(row_indices)) + tbc_sub_matrix = coo_matrix( + (vals, (row_indices, col_indices)), shape=((n_tb + n_bcs) * n_lr, 1) + ) + tbc_matrix = csr_matrix(kron(eye(second_dim_repeats), tbc_sub_matrix)) + + if tbc_value.evaluates_to_number(): + top_ghost_constant = ( + 2 * tbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + ) + else: + top_ghost_constant = 2 * tbc_value + tbc_vector = pybamm.Matrix(tbc_matrix) @ top_ghost_constant + elif tbc_type == "Neumann": + tbc_vector = pybamm.Vector( + np.zeros((n_tb + n_bcs) * second_dim_repeats * n_lr) + ) + else: + tbc_vector = pybamm.Vector( + np.zeros((n_lr + n_bcs) * second_dim_repeats * n_tb) + ) + + # Calculate values for ghost nodes for any Dirichlet boundary conditions + if bbc_type == "Dirichlet": + # Create matrix to extract the leftmost column of values + row_indices = np.arange( + (n_lr * (n_tb + n_bcs)) - n_lr, n_lr * (n_tb + n_bcs) + ) + col_indices = np.zeros(len(row_indices)) + vals = np.ones(len(row_indices)) + bbc_sub_matrix = coo_matrix( + (vals, (row_indices, col_indices)), shape=((n_tb + n_bcs) * n_lr, 1) + ) + bbc_matrix = csr_matrix(kron(eye(second_dim_repeats), bbc_sub_matrix)) + + if bbc_value.evaluates_to_number(): + bottom_ghost_constant = ( + 2 * bbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + ) + else: + bottom_ghost_constant = 2 * bbc_value + bbc_vector = pybamm.Matrix(bbc_matrix) @ bottom_ghost_constant + elif bbc_type == "Neumann": + bbc_vector = pybamm.Vector( + np.zeros((n_tb + n_bcs) * second_dim_repeats * n_lr) + ) + else: + bbc_vector = pybamm.Vector( + np.zeros((n_lr + n_bcs) * second_dim_repeats * n_tb) + ) + + bcs_vector = lbc_vector + rbc_vector + tbc_vector + bbc_vector + # Need to match the domain. E.g. in the case of the boundary condition + # on the particle, the gradient has domain particle but the bcs_vector + # has domain electrode, since it is a function of the macroscopic variables + bcs_vector.copy_domains(discretised_symbol) + + # Make matrix to calculate ghost nodes + # coo_matrix takes inputs (data, (row, col)) and puts data[i] at the point + # (row[i], col[i]) for each index of data. + if lbc_type == "Dirichlet": + left_ghost_vector = coo_matrix(([-1], ([0], [0])), shape=(1, n_lr)) + else: + left_ghost_vector = None + if rbc_type == "Dirichlet": + right_ghost_vector = coo_matrix(([-1], ([0], [n_lr - 1])), shape=(1, n_lr)) + else: + right_ghost_vector = None + + if tbc_type == "Dirichlet": + row_indices = np.arange(0, n_lr) + col_indices = np.arange(0, n_lr) + top_ghost_vector = coo_matrix( + (-np.ones(n_lr), (row_indices, col_indices)), shape=(n_lr, n) + ) + else: + top_ghost_vector = None + if bbc_type == "Dirichlet": + row_indices = np.arange(0, n_lr) + col_indices = np.arange(n - n_lr, n) + bottom_ghost_vector = coo_matrix( + (-np.ones(n_lr), (row_indices, col_indices)), shape=(n_lr, n) + ) + else: + bottom_ghost_vector = None + + if lbc_type == "Dirichlet" or rbc_type == "Dirichlet": + sub_matrix = vstack([left_ghost_vector, eye(n_lr), right_ghost_vector]) + sub_matrix = block_diag((sub_matrix,) * n_tb) + else: + sub_matrix = vstack( + [ + left_ghost_vector, + top_ghost_vector, + eye(n), + bottom_ghost_vector, + right_ghost_vector, + ] + ) + + # repeat matrix for secondary dimensions + # Convert to csr_matrix so that we can take the index (row-slicing), which is + # not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should not be an + # issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + + new_symbol = pybamm.Matrix(matrix) @ discretised_symbol + bcs_vector + + return new_symbol, domain + + def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): + """ + Add the known values of the gradient from Neumann boundary conditions to + the discretised gradient. + + Dirichlet bcs are implemented using ghost nodes, see + :meth:`pybamm.FiniteVolume.add_ghost_nodes`. + + Parameters + ---------- + symbol : :class:`pybamm.SpatialVariable` + The variable to be discretised + discretised_gradient : :class:`pybamm.Vector` + Contains the discretised gradient of symbol + bcs : dict of tuples (:class:`pybamm.Scalar`, str) + Dictionary (with keys "left" and "right") of boundary conditions. Each + boundary condition consists of a value and a flag indicating its type + (e.g. "Dirichlet") + domain : list of strings + The domain of the gradient of the symbol (may include ghost nodes) + + Returns + ------- + :class:`pybamm.Symbol` + `Matrix @ discretised_gradient + bcs_vector`. When evaluated, this gives the + discretised_gradient, with the values of the Neumann boundary conditions + concatenated at each end (if given). + + """ + # get relevant grid points + submesh = self.mesh[domain] + + # Prepare sizes and empty bcs_vector + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + lbc_value, lbc_type = bcs.get("left", (None, None)) + rbc_value, rbc_type = bcs.get("right", (None, None)) + tbc_value, tbc_type = bcs.get("top", (None, None)) + bbc_value, bbc_type = bcs.get("bottom", (None, None)) + + # Count number of Neumann boundary conditions + n_bcs = 0 + if lbc_type == "Neumann": + n_bcs += 1 + if rbc_type == "Neumann": + n_bcs += 1 + if tbc_type == "Neumann": + n_bcs += 1 + if bbc_type == "Neumann": + n_bcs += 1 + + # Add any values from Neumann boundary conditions to the bcs vector + if lbc_type == "Neumann" and lbc_value != 0: + lbc_sub_matrix = coo_matrix(([1], ([0], [0])), shape=(n_lr - 1 + n_bcs, 1)) + lbc_matrix = csr_matrix(kron(eye(second_dim_repeats), lbc_sub_matrix)) + lbc_matrix = vstack([lbc_matrix] * n_tb) + if lbc_value.evaluates_to_number(): + left_bc = lbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + else: + left_bc = lbc_value + lbc_vector = pybamm.Matrix(lbc_matrix) @ left_bc + elif lbc_type == "Dirichlet" or (lbc_type == "Neumann" and lbc_value == 0): + lbc_vector = pybamm.Vector( + np.zeros((n_lr - 1 + n_bcs) * second_dim_repeats * n_tb) + ) + elif lbc_type is None and rbc_type is None: + lbc_vector = pybamm.Vector( + np.zeros((n_tb - 1 + n_bcs) * second_dim_repeats * n_lr) + ) + else: + raise ValueError( + f"boundary condition must be Dirichlet or Neumann, not '{rbc_type}'" + ) + + if rbc_type == "Neumann" and rbc_value != 0: + rbc_sub_matrix = coo_matrix( + ([1], ([n_lr + n_bcs - 2], [0])), shape=(n_lr - 1 + n_bcs, 1) + ) + rbc_matrix = csr_matrix(kron(eye(second_dim_repeats), rbc_sub_matrix)) + rbc_matrix = vstack([rbc_matrix] * n_tb) + if rbc_value.evaluates_to_number(): + right_bc = rbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + else: + right_bc = rbc_value + rbc_vector = pybamm.Matrix(rbc_matrix) @ right_bc + elif rbc_type == "Dirichlet" or (rbc_type == "Neumann" and rbc_value == 0): + rbc_vector = pybamm.Vector( + np.zeros((n_lr - 1 + n_bcs) * second_dim_repeats * n_tb) + ) + elif rbc_type is None and lbc_type is None: + rbc_vector = pybamm.Vector( + np.zeros((n_tb - 1 + n_bcs) * second_dim_repeats * n_lr) + ) + else: + raise ValueError( + f"boundary condition must be Dirichlet or Neumann, not '{rbc_type}'" + ) + + # Add any values from Neumann boundary conditions to the bcs vector + if tbc_type == "Neumann" and tbc_value != 0: + row_indices = np.arange(0, n_lr) + col_indices = np.zeros(len(row_indices)) + vals = np.ones(len(row_indices)) + tbc_sub_matrix = coo_matrix( + (vals, (row_indices, col_indices)), shape=((n_tb - 1 + n_bcs) * n_lr, 1) + ) + tbc_matrix = csr_matrix(kron(eye(second_dim_repeats), tbc_sub_matrix)) + if tbc_value.evaluates_to_number(): + top_bc = tbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + else: + top_bc = tbc_value + tbc_vector = pybamm.Matrix(tbc_matrix) @ top_bc + elif tbc_type == "Dirichlet" or (tbc_type == "Neumann" and tbc_value == 0): + tbc_vector = pybamm.Vector( + np.zeros((n_tb - 1 + n_bcs) * second_dim_repeats * n_lr) + ) + elif tbc_type is None and bbc_type is None: + tbc_vector = pybamm.Vector( + np.zeros((n_lr - 1 + n_bcs) * second_dim_repeats * n_tb) + ) + else: + raise ValueError( + f"boundary condition must be Dirichlet or Neumann, not '{tbc_type}'" + ) + + # Add any values from Neumann boundary conditions to the bcs vector + if bbc_type == "Neumann" and bbc_value != 0: + row_indices = np.arange( + (n_lr * (n_tb - 1 + n_bcs)) - n_lr, n_lr * (n_tb - 1 + n_bcs) + ) + col_indices = np.zeros(len(row_indices)) + vals = np.ones(len(row_indices)) + bbc_sub_matrix = coo_matrix( + (vals, (row_indices, col_indices)), shape=((n_tb - 1 + n_bcs) * n_lr, 1) + ) + bbc_matrix = csr_matrix(kron(eye(second_dim_repeats), bbc_sub_matrix)) + if bbc_value.evaluates_to_number(): + bottom_bc = bbc_value * pybamm.Vector(np.ones(second_dim_repeats)) + else: + bottom_bc = bbc_value + bbc_vector = pybamm.Matrix(bbc_matrix) @ bottom_bc + elif bbc_type == "Dirichlet" or (bbc_type == "Neumann" and bbc_value == 0): + bbc_vector = pybamm.Vector( + np.zeros((n_tb - 1 + n_bcs) * second_dim_repeats * n_lr) + ) + elif bbc_type is None and tbc_type is None: + bbc_vector = pybamm.Vector( + np.zeros((n_lr - 1 + n_bcs) * second_dim_repeats * n_tb) + ) + else: + raise ValueError( + f"boundary condition must be Dirichlet or Neumann, not '{bbc_type}'" + ) + + bcs_vector = lbc_vector + rbc_vector + tbc_vector + bbc_vector + # Need to match the domain. E.g. in the case of the boundary condition + # on the particle, the gradient has domain particle but the bcs_vector + # has domain electrode, since it is a function of the macroscopic variables + bcs_vector.copy_domains(discretised_gradient) + + # Make matrix which makes "gaps" in the the discretised gradient into + # which the known Neumann values will be added. E.g. in 1D if the left + # boundary condition is Dirichlet and the right Neumann, this matrix will + # act to append a zero to the end of the discretised gradient + if lbc_type == "Neumann": + left_vector = csr_matrix((1, n_lr - 1)) + else: + left_vector = None + if rbc_type == "Neumann": + right_vector = csr_matrix((1, n_lr - 1)) + else: + right_vector = None + if tbc_type == "Neumann": + top_vector = csr_matrix((n_lr, (n_tb - 1) * n_lr)) + else: + top_vector = None + if bbc_type == "Neumann": + bottom_vector = csr_matrix((n_lr, (n_tb - 1) * n_lr)) + else: + bottom_vector = None + + if lbc_type == "Neumann" or rbc_type == "Neumann": + sub_matrix = vstack([left_vector, eye(n_lr - 1), right_vector]) + sub_matrix = block_diag((sub_matrix,) * n_tb) + elif tbc_type == "Neumann" or bbc_type == "Neumann": + sub_matrix = vstack([top_vector, eye((n_tb - 1) * n_lr), bottom_vector]) + + # repeat matrix for secondary dimensions + # Convert to csr_matrix so that we can take the index (row-slicing), which is + # not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should not be an + # issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + + new_gradient = pybamm.Matrix(matrix) @ discretised_gradient + bcs_vector + + return new_gradient + + def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): + """ + Uses extrapolation to get the boundary value or flux of a variable in the + Finite Volume Method. + + See :meth:`pybamm.SpatialMethod.boundary_value` + """ + # Find the number of submeshes + submesh = self.mesh[discretised_child.domain] + + if "-" in symbol.side: + side_first, side_second = symbol.side.split("-") + if ("top" in side_first or "bottom" in side_first) and ( + "right" in side_second or "left" in side_second + ): + side_first, side_second = side_second, side_first + else: + side_first = symbol.side + side_second = None + + repeats = self._get_auxiliary_domain_repeats(discretised_child.domains) + + if bcs is None: + bcs = {} + + extrap_order_gradient = self.options["extrapolation"]["order"]["gradient"] + extrap_order_value = self.options["extrapolation"]["order"]["value"] + use_bcs = self.options["extrapolation"]["use bcs"] + + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + + nodes_lr = submesh.nodes_lr + edges_lr = submesh.edges_lr + nodes_tb = submesh.nodes_tb + edges_tb = submesh.edges_tb + + dx0_lr = nodes_lr[0] - edges_lr[0] + dx1_lr = submesh.d_nodes_lr[0] + # dx2_lr = submesh.d_nodes_lr[1] + + dxN_lr = edges_lr[-1] - nodes_lr[-1] + dxNm1_lr = submesh.d_nodes_lr[-1] + # dxNm2_lr = submesh.d_nodes_lr[-2] + + dx0_tb = nodes_tb[0] - edges_tb[0] + dx1_tb = submesh.d_nodes_tb[0] + # dx2_tb = submesh.d_nodes_tb[1] + + dxN_tb = edges_tb[-1] - nodes_tb[-1] + dxNm1_tb = submesh.d_nodes_tb[-1] + # dxNm2_tb = submesh.d_nodes_tb[-2] + + child = symbol.child + + # Create submatrix to compute boundary values or fluxes + # Derivation of extrapolation formula can be found at: + # https://github.com/Scottmar93/extrapolation-coefficents/tree/master + if isinstance(symbol, pybamm.BoundaryValue): + if use_bcs and pybamm.has_bc_of_form(child, side_first, bcs, "Dirichlet"): + if side_first == "left" or side_first == "right": + # just use the value from the bc: f(x*) + sub_matrix = csr_matrix((n_tb, n_tb * n_lr)) + additive = bcs[child][side_first][0] + elif side_first == "top" or side_first == "bottom": + sub_matrix = csr_matrix((n_lr, n_lr * n_tb)) + additive = bcs[child][side_first][0] + else: + raise ValueError( + "side_first must be 'left', 'right', 'top', or 'bottom'" + ) + + elif side_first == "left": + if extrap_order_value == "linear": + # to find value at x* use formula: + # f(x*) = f_1 - (dx0 / dx1) (f_2 - f_1) + + if use_bcs and pybamm.has_bc_of_form( + child, side_first, bcs, "Neumann" + ): + dx0 = dx0_lr + row_indices = np.arange(0, n_tb) + col_indices_0 = np.arange(0, n_tb * n_lr, n_lr) + vals_0 = np.ones(n_tb) + sub_matrix = csr_matrix( + ( + vals_0, + ( + row_indices, + col_indices_0, + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = -dx0 * bcs[child][side_first][0] + + else: + dx0 = dx0_lr + dx1 = dx1_lr + row_indices = np.arange(0, n_tb) + col_indices_0 = np.arange(0, n_tb * n_lr, n_lr) + col_indices_1 = col_indices_0 + 1 + vals_0 = np.ones(n_tb) * (1 + (dx0 / dx1)) + vals_1 = np.ones(n_tb) * (-(dx0 / dx1)) + sub_matrix = csr_matrix( + ( + np.hstack([vals_0, vals_1]), + ( + np.hstack([row_indices, row_indices]), + np.hstack([col_indices_0, col_indices_1]), + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = pybamm.Scalar(0) + + elif extrap_order_value == "quadratic": + if use_bcs and pybamm.has_bc_of_form( + child, symbol.side, bcs, "Neumann" + ): + raise NotImplementedError + + else: + raise NotImplementedError + + else: + raise NotImplementedError + + elif side_first == "right": + if extrap_order_value == "linear": + if use_bcs and pybamm.has_bc_of_form( + child, side_first, bcs, "Neumann" + ): + dxN = dxN_lr + row_indices = np.arange(0, n_tb) + col_indices_N = np.arange(n_lr - 1, n_lr * n_tb, n_lr) + vals_N = np.ones(n_tb) + sub_matrix = csr_matrix( + ( + vals_0, + ( + row_indices, + col_indices_N, + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = dxN * bcs[child][side_first][0] + + else: + # to find value at x* use formula: + # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1) + dxN = dxN_lr + dxNm1 = dxNm1_lr + row_indices = np.arange(0, n_tb) + col_indices_Nm1 = np.arange(n_lr - 2, n_lr * n_tb, n_lr) + col_indices_N = col_indices_Nm1 + 1 + vals_Nm1 = np.ones(n_tb) * (-(dxN / dxNm1)) + vals_N = np.ones(n_tb) * (1 + (dxN / dxNm1)) + sub_matrix = csr_matrix( + ( + np.hstack([vals_Nm1, vals_N]), + ( + np.hstack([row_indices, row_indices]), + np.hstack([col_indices_Nm1, col_indices_N]), + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = pybamm.Scalar(0) + elif extrap_order_value == "quadratic": + if use_bcs and pybamm.has_bc_of_form( + child, symbol.side, bcs, "Neumann" + ): + raise NotImplementedError + else: + raise NotImplementedError + else: + raise NotImplementedError + + elif side_first == "top": + if extrap_order_value == "linear": + if use_bcs and pybamm.has_bc_of_form( + child, side_first, bcs, "Neumann" + ): + dx0 = dx0_tb + first_val = np.ones(n_lr) + sub_matrix = spdiags( + first_val, + [ + 0, + ], + n_lr, + n_lr * n_tb, + ) + additive = -dx0 * bcs[child][side_first][0] + else: + dx0 = dx0_tb + dx1 = dx1_tb + first_val = (1 + (dx0 / dx1)) * np.ones(n_lr) + second_val = -(dx0 / dx1) * np.ones(n_lr) + rows_first = np.arange(0, n_lr) + rows_second = rows_first + cols_first = np.arange(0, n_lr) + cols_second = np.arange(n_lr, 2 * n_lr) + rows = np.concatenate([rows_first, rows_second]) + cols = np.concatenate([cols_first, cols_second]) + vals = np.concatenate([first_val, second_val]) + sub_matrix = csr_matrix( + (vals, (rows, cols)), + shape=(n_lr, n_lr * n_tb), + ) + additive = pybamm.Scalar(0) + + elif extrap_order_value == "quadratic": + if use_bcs and pybamm.has_bc_of_form( + child, symbol.side, bcs, "Neumann" + ): + raise NotImplementedError + else: + raise NotImplementedError + else: + raise NotImplementedError + + elif side_first == "bottom": + if extrap_order_value == "linear": + if use_bcs and pybamm.has_bc_of_form( + child, side_first, bcs, "Neumann" + ): + dxNm1 = dxNm1_tb + dxN = dxN_tb + val_N = np.ones(n_lr) + sub_matrix = spdiags( + val_N, + [ + (n_tb - 1) * n_lr, + ], + n_lr, + n_lr * n_tb, + ) + additive = dxN * bcs[child][side_first][0] + else: + dx0 = dxNm1_tb + dx1 = dxN_tb + first_val = -(dx0 / dx1) * np.ones(n_lr) + second_val = (1 + (dx0 / dx1)) * np.ones(n_lr) + rows_first = np.arange(0, n_lr) + rows_second = np.arange(0, n_lr) + cols_first = np.arange((n_tb - 2) * n_lr, (n_tb - 1) * n_lr) + cols_second = np.arange((n_tb - 1) * n_lr, n_tb * n_lr) + rows = np.concatenate([rows_first, rows_second]) + cols = np.concatenate([cols_first, cols_second]) + vals = np.concatenate([first_val, second_val]) + sub_matrix = csr_matrix( + (vals, (rows, cols)), + shape=(n_lr, n_lr * n_tb), + ) + additive = pybamm.Scalar(0) + + elif extrap_order_value == "quadratic": + if use_bcs and pybamm.has_bc_of_form( + child, symbol.side, bcs, "Neumann" + ): + raise NotImplementedError + else: + raise NotImplementedError + else: + raise NotImplementedError + + if side_second == "top": + if use_bcs and pybamm.has_bc_of_form( + child, side_second, bcs, "Neumann" + ): + dx0 = dx0_tb + additive = -dx0 * bcs[child][side_second][0] + sub_matrix_second = csr_matrix( + ( + [ + 1, + ], + ( + [ + 0, + ], + [ + 0, + ], + ), + ), + shape=(1, n_tb), + ) + sub_matrix = sub_matrix_second @ sub_matrix + else: + dx0 = dx0_tb + dx1 = dx1_tb + row_indices = [0, 0] + col_indices = [0, 1] + vals = [1 + (dx0 / dx1), -(dx0 / dx1)] + sub_matrix_second = csr_matrix( + (vals, (row_indices, col_indices)), shape=(1, n_tb) + ) + sub_matrix = sub_matrix_second @ sub_matrix + elif side_second == "bottom": + if use_bcs and pybamm.has_bc_of_form( + child, side_second, bcs, "Neumann" + ): + dxN = dxN_tb + additive = dxN * bcs[child][side_second][0] + sub_matrix_second = csr_matrix( + ( + [ + 1, + ], + ( + [ + 0, + ], + [ + n_tb - 1, + ], + ), + ), + shape=(1, n_tb), + ) + sub_matrix = sub_matrix_second @ sub_matrix + else: + dxN = dxN_tb + dxNm1 = dxNm1_tb + row_indices = [0, 0] + col_indices = [n_tb - 2, n_tb - 1] + vals = [-(dxN / dxNm1), 1 + (dxN / dxNm1)] + sub_matrix_second = csr_matrix( + (vals, (row_indices, col_indices)), shape=(1, n_tb) + ) + sub_matrix = sub_matrix_second @ sub_matrix + elif side_second is None: + pass + else: + raise ValueError("side_second must be 'top' or 'bottom'") + + elif isinstance(symbol, pybamm.BoundaryGradient): + if use_bcs and pybamm.has_bc_of_form(child, side_first, bcs, "Neumann"): + if side_first == "left" or side_first == "right": + # just use the value from the bc: f(x*) + sub_matrix = csr_matrix((n_tb, n_tb * n_lr)) + additive = bcs[child][side_first][0] + elif side_first == "top" or side_first == "bottom": + sub_matrix = csr_matrix((n_lr, n_lr * n_tb)) + additive = bcs[child][side_first][0] + else: + raise ValueError( + "side_first must be 'left', 'right', 'top', or 'bottom'" + ) + + elif side_first == "left": + if extrap_order_gradient == "linear": + # f'(x*) = (f_2 - f_1) / dx1 + dx1 = dx1_lr + row_indices = np.arange(0, n_tb) + col_indices_0 = np.arange(0, n_tb * n_lr, n_lr) + col_indices_1 = col_indices_0 + 1 + vals_0 = np.ones(n_tb) * (-1 / dx1) + vals_1 = np.ones(n_tb) * (1 / dx1) + sub_matrix = csr_matrix( + ( + np.hstack([vals_0, vals_1]), + ( + np.hstack([row_indices, row_indices]), + np.hstack([col_indices_0, col_indices_1]), + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = pybamm.Scalar(0) + + elif extrap_order_gradient == "quadratic": + raise NotImplementedError + + else: + raise NotImplementedError + + elif side_first == "right": + if extrap_order_gradient == "linear": + # use formula: + # f'(x*) = (f_N - f_Nm1) / dxNm1 + dxN = dxN_lr + dxNm1 = dxNm1_lr + row_indices = np.arange(0, n_tb) + col_indices_Nm1 = np.arange(n_lr - 2, n_lr * n_tb, n_lr) + col_indices_N = col_indices_Nm1 + 1 + vals_Nm1 = np.ones(n_tb) * (-1 / dxNm1) + vals_N = np.ones(n_tb) * (1 / dxNm1) + sub_matrix = csr_matrix( + ( + np.hstack([vals_Nm1, vals_N]), + ( + np.hstack([row_indices, row_indices]), + np.hstack([col_indices_Nm1, col_indices_N]), + ), + ), + shape=(n_tb, n_tb * n_lr), + ) + additive = pybamm.Scalar(0) + + elif extrap_order_gradient == "quadratic": + raise NotImplementedError + else: + raise NotImplementedError + + elif side_first == "top": + dx1 = dx1_tb + first_val = (-1 / dx1) * np.ones(n_lr) + second_val = (1 / dx1) * np.ones(n_lr) + sub_matrix = spdiags( + [first_val, second_val], [0, n_lr], n_lr, n_lr * n_tb + ) + additive = pybamm.Scalar(0) + + elif side_first == "bottom": + dx0 = dxNm1_tb + dx1 = dxN_tb + first_val = -(dx0 / dx1) * np.ones(n_lr) + second_val = (1 + (dx0 / dx1)) * np.ones(n_lr) + sub_matrix = spdiags( + [first_val, second_val], + [(n_tb - 2) * n_lr, (n_tb - 1) * n_lr], + n_lr, + n_lr * n_tb, + ) + additive = pybamm.Scalar(0) + + if side_second == "top": + dx0 = dx0_tb + dx1 = dx1_tb + row_indices = [0, 0] + col_indices = [0, 1] + vals = [-1 / dx1, 1 / dx1] + sub_matrix_second = csr_matrix( + (vals, (row_indices, col_indices)), shape=(1, n_tb) + ) + sub_matrix = sub_matrix_second @ sub_matrix + elif side_second == "bottom": + dxN = dxN_tb + dxNm1 = dxNm1_tb + row_indices = [0, 0] + col_indices = [n_tb - 2, n_tb - 1] + vals = [-1 / dxNm1, 1 / dxNm1] + sub_matrix_second = csr_matrix( + (vals, (row_indices, col_indices)), shape=(1, n_tb) + ) + sub_matrix = sub_matrix_second @ sub_matrix + elif side_second is None: + pass + else: + raise ValueError("side_second must be 'top' or 'bottom'") + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), which is + # not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should not be an + # issue + matrix = csr_matrix(kron(eye(repeats), sub_matrix)) + + # Return boundary value with domain given by symbol + matrix = pybamm.Matrix(matrix) + boundary_value = matrix @ discretised_child + boundary_value.copy_domains(symbol) + + additive.copy_domains(symbol) + boundary_value += additive + + return boundary_value + + def evaluate_at(self, symbol, discretised_child, position): + """ + Returns the symbol evaluated at a given position in space. + + Parameters + ---------- + symbol: :class:`pybamm.Symbol` + The boundary value or flux symbol + discretised_child : :class:`pybamm.StateVector` + The discretised variable from which to calculate the boundary value + position : :class:`pybamm.Scalar` + The point in one-dimensional space at which to evaluate the symbol. + + Returns + ------- + :class:`pybamm.MatrixMultiplication` + The variable representing the value at the given point. + """ + raise NotImplementedError + + def process_binary_operators(self, bin_op, left, right, disc_left, disc_right): + """Discretise binary operators in model equations. Performs appropriate + averaging of diffusivities if one of the children is a gradient operator, so + that discretised sizes match up. For this averaging we use the harmonic + mean [1]. + + [1] Recktenwald, Gerald. "The control-volume finite-difference approximation to + the diffusion equation." (2012). + + Parameters + ---------- + bin_op : :class:`pybamm.BinaryOperator` + Binary operator to discretise + left : :class:`pybamm.Symbol` + The left child of `bin_op` + right : :class:`pybamm.Symbol` + The right child of `bin_op` + disc_left : :class:`pybamm.Symbol` + The discretised left child of `bin_op` + disc_right : :class:`pybamm.Symbol` + The discretised right child of `bin_op` + Returns + ------- + :class:`pybamm.BinaryOperator` + Discretised binary operator + + """ + """Discretise binary operators in model equations. Performs appropriate + averaging of diffusivities if one of the children is a gradient operator, so + that discretised sizes match up. For this averaging we use the harmonic + mean [1]. + + [1] Recktenwald, Gerald. "The control-volume finite-difference approximation to + the diffusion equation." (2012). + + Parameters + ---------- + bin_op : :class:`pybamm.BinaryOperator` + Binary operator to discretise + left : :class:`pybamm.Symbol` + The left child of `bin_op` + right : :class:`pybamm.Symbol` + The right child of `bin_op` + disc_left : :class:`pybamm.Symbol` + The discretised left child of `bin_op` + disc_right : :class:`pybamm.Symbol` + The discretised right child of `bin_op` + Returns + ------- + :class:`pybamm.BinaryOperator` + Discretised binary operator + + """ + # Post-processing to make sure discretised dimensions match + left_evaluates_on_edges = left.evaluates_on_edges("primary") + right_evaluates_on_edges = right.evaluates_on_edges("primary") + + # This could be cleaned up a bit, but it works for now. + if hasattr(disc_left, "lr_field") and hasattr(disc_right, "lr_field"): + if right_evaluates_on_edges and not left_evaluates_on_edges: + if isinstance(right, pybamm.Gradient): + method = "harmonic" + disc_left_lr = self.node_to_edge( + disc_left.lr_field, method=method, direction="lr" + ) + disc_left_tb = self.node_to_edge( + disc_left.tb_field, method=method, direction="tb" + ) + disc_left = pybamm.VectorField(disc_left_lr, disc_left_tb) + else: + method = "arithmetic" + disc_left_lr = self.node_to_edge( + disc_left.lr_field, method=method, direction="lr" + ) + disc_left_tb = self.node_to_edge( + disc_left.tb_field, method=method, direction="tb" + ) + disc_left = pybamm.VectorField(disc_left_lr, disc_left_tb) + elif left_evaluates_on_edges and not right_evaluates_on_edges: + if isinstance(left, pybamm.Gradient): + method = "harmonic" + disc_right_lr = self.node_to_edge( + disc_right.lr_field, method=method, direction="lr" + ) + disc_right_tb = self.node_to_edge( + disc_right.tb_field, method=method, direction="tb" + ) + disc_right = pybamm.VectorField(disc_right_lr, disc_right_tb) + else: + method = "arithmetic" + disc_right_lr = self.node_to_edge( + disc_right.lr_field, method=method, direction="lr" + ) + disc_right_tb = self.node_to_edge( + disc_right.tb_field, method=method, direction="tb" + ) + disc_right = pybamm.VectorField(disc_right_lr, disc_right_tb) + # both are vector fields, so we need make a new vector field. + lr_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left.lr_field, disc_right.lr_field]) + ) + tb_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left.tb_field, disc_right.tb_field]) + ) + return pybamm.VectorField(lr_field, tb_field) + elif hasattr(disc_left, "lr_field") and not hasattr(disc_right, "lr_field"): + # one is a vector field, so we need to make a new vector field. + if left_evaluates_on_edges and not right_evaluates_on_edges: + if isinstance(left, pybamm.Gradient): + method = "harmonic" + disc_right_lr = self.node_to_edge( + disc_right, method=method, direction="lr" + ) + disc_right_tb = self.node_to_edge( + disc_right, method=method, direction="tb" + ) + disc_right = pybamm.VectorField(disc_right_lr, disc_right_tb) + else: + method = "arithmetic" + lr_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left.lr_field, disc_right]) + ) + tb_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left.tb_field, disc_right]) + ) + return pybamm.VectorField(lr_field, tb_field) + elif not hasattr(disc_left, "lr_field") and hasattr(disc_right, "lr_field"): + # one is a vector field, so we need to make a new vector field. + if right_evaluates_on_edges and not left_evaluates_on_edges: + if isinstance(right, pybamm.Gradient): + method = "harmonic" + disc_left_lr = self.node_to_edge( + disc_left, method=method, direction="lr" + ) + disc_left_tb = self.node_to_edge( + disc_left, method=method, direction="tb" + ) + disc_left = pybamm.VectorField(disc_left_lr, disc_left_tb) + else: + method = "arithmetic" + lr_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left, disc_right.lr_field]) + ) + tb_field = pybamm.simplify_if_constant( + bin_op.create_copy([disc_left, disc_right.tb_field]) + ) + return pybamm.VectorField(lr_field, tb_field) + else: + pass + # inner product takes fluxes from edges to nodes + if isinstance(bin_op, pybamm.Inner): + if left_evaluates_on_edges: + disc_left = self.edge_to_node(disc_left) + if right_evaluates_on_edges: + disc_right = self.edge_to_node(disc_right) + + # If neither child evaluates on edges, or both children have gradients, + # no need to do any averaging + elif left_evaluates_on_edges == right_evaluates_on_edges: + pass + # If only left child evaluates on edges, map right child onto edges + # using the harmonic mean if the left child is a gradient (i.e. this + # binary operator represents a flux) + elif left_evaluates_on_edges and not right_evaluates_on_edges: + if isinstance(left, pybamm.Gradient): + method = "harmonic" + disc_right_lr = self.node_to_edge( + disc_right, method=method, direction="lr" + ) + disc_right_tb = self.node_to_edge( + disc_right, method=method, direction="tb" + ) + disc_right = pybamm.VectorField(disc_right_lr, disc_right_tb) + else: + method = "arithmetic" + + # If only right child evaluates on edges, map left child onto edges + # using the harmonic mean if the right child is a gradient (i.e. this + # binary operator represents a flux) + elif right_evaluates_on_edges and not left_evaluates_on_edges: + if isinstance(right, pybamm.Gradient): + method = "harmonic" + disc_left_lr = self.node_to_edge( + disc_left, method=method, direction="lr" + ) + disc_left_tb = self.node_to_edge( + disc_left, method=method, direction="tb" + ) + disc_left = pybamm.VectorField(disc_left_lr, disc_left_tb) + else: + method = "arithmetic" + + # Return new binary operator with appropriate class + out = pybamm.simplify_if_constant(bin_op.create_copy([disc_left, disc_right])) + + return out + + def concatenation(self, disc_children): + """Discrete concatenation, taking `edge_to_node` for children that evaluate on + edges. + See :meth:`pybamm.SpatialMethod.concatenation` + """ + for idx, child in enumerate(disc_children): + submesh = self.mesh[child.domain] + repeats = self._get_auxiliary_domain_repeats(child.domains) + n_nodes = len(submesh.nodes_lr) * len(submesh.nodes_tb) * repeats + n_edges = len(submesh.edges_lr) * len(submesh.edges_tb) * repeats + child_size = child.size + if child_size != n_nodes: + # Average any children that evaluate on the edges (size n_edges) to + # evaluate on nodes instead, so that concatenation works properly + if child_size == n_edges: + disc_children[idx] = self.edge_to_node(child) + else: + raise pybamm.ShapeError( + "child must have size n_nodes (number of nodes in the mesh) " + "or n_edges (number of edges in the mesh)" + ) + # EXPERIMENTAL: Need to reorder things for 2D + if not all(isinstance(child, pybamm.StateVector) for child in disc_children): + # All will have the same number of points in the tb direction, so we just need to get the lr points + lr_mesh_points = [ + self.mesh[child.domain[0]].npts_lr for child in disc_children + ] + tb_mesh_points = self.mesh[disc_children[0].domain[0]].npts_tb + num_children = len(disc_children) + eyes = [np.eye(lr_mesh_points[i]) for i in range(num_children)] + num_blocks = num_children * tb_mesh_points + block_mat = [ + [ + np.zeros((lr_mesh_points_this, lr_mesh_points_this)) + for lr_mesh_points_this in lr_mesh_points + for _ in range(tb_mesh_points) + ] + for _ in range(num_blocks) + ] + i = 0 + for tb_idx in range(tb_mesh_points): + for child_idx in range(num_children): + if isinstance(disc_children[child_idx], pybamm.StateVector): + block_mat[i][tb_idx + child_idx * tb_mesh_points] = eyes[ + child_idx + ] + i += 1 + else: + block_mat[i][tb_idx + child_idx * tb_mesh_points] = eyes[ + child_idx + ] + i += 1 + block_mat = csr_matrix(np.block(block_mat)) + block_mat = kron(eye(repeats), block_mat) + matrix = pybamm.Matrix(block_mat) + repeats = self._get_auxiliary_domain_repeats(disc_children[0].domains) + return matrix @ pybamm.domain_concatenation(disc_children, self.mesh) + else: + return pybamm.domain_concatenation(disc_children, self.mesh) + + def edge_to_node(self, discretised_symbol, method="arithmetic"): + """ + Convert a discretised symbol evaluated on the cell edges to a discretised symbol + evaluated on the cell nodes. + See :meth:`pybamm.FiniteVolume.shift` + """ + return self.shift(discretised_symbol, "edge to node", method) + + def node_to_edge(self, discretised_symbol, method="arithmetic", direction="lr"): + """ + Convert a discretised symbol evaluated on the cell nodes to a discretised symbol + evaluated on the cell edges. + See :meth:`pybamm.FiniteVolume.shift` + """ + return self.shift(discretised_symbol, "node to edge", method, direction) + + def shift(self, discretised_symbol, shift_key, method, direction="lr"): + """ + Convert a discretised symbol evaluated at edges/nodes, to a discretised symbol + evaluated at nodes/edges. Can be the arithmetic mean or the harmonic mean. + + Note: when computing fluxes at cell edges it is better to take the + harmonic mean based on [1]. + + [1] Recktenwald, Gerald. "The control-volume finite-difference approximation to + the diffusion equation." (2012). + + Parameters + ---------- + discretised_symbol : :class:`pybamm.Symbol` + Symbol to be averaged. When evaluated, this symbol returns either a scalar + or an array of shape (n,) or (n+1,), where n is the number of points in the + mesh for the symbol's domain (n = self.mesh[symbol.domain].npts) + shift_key : str + Whether to shift from nodes to edges ("node to edge"), or from edges to + nodes ("edge to node") + method : str + Whether to use the "arithmetic" or "harmonic" mean + + Returns + ------- + :class:`pybamm.Symbol` + Averaged symbol. When evaluated, this returns either a scalar or an array of + shape (n+1,) (if `shift_key = "node to edge"`) or (n,) (if + `shift_key = "edge to node"`) + """ + + def arithmetic_mean(array, direction): + """Calculate the arithmetic mean of an array using matrix multiplication""" + # Create appropriate submesh by combining submeshes in domain + submesh = self.mesh[array.domain] + + # Create 1D matrix using submesh + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + + if shift_key == "node to edge": + if direction == "lr": + sub_matrix_left = csr_matrix( + ([1.5, -0.5], ([0, 0], [0, 1])), shape=(1, n_lr) + ) + sub_matrix_center = diags( + [0.5, 0.5], [0, 1], shape=(n_lr - 1, n_lr) + ) + sub_matrix_right = csr_matrix( + ([-0.5, 1.5], ([0, 0], [n_lr - 2, n_lr - 1])), shape=(1, n_lr) + ) + sub_matrix = vstack( + [sub_matrix_left, sub_matrix_center, sub_matrix_right] + ) + sub_matrix = block_diag((sub_matrix,) * n_tb) + elif direction == "tb": + # Matrix to compute values at the exterior edges + one_fives_top = np.ones(n_lr) * 1.5 + neg_zero_fives_top = np.ones(n_lr) * -0.5 + rows = np.arange(0, n_lr) + cols_first = np.arange(0, n_lr) + cols_second = np.arange(n_lr, 2 * n_lr) + data = np.hstack([one_fives_top, neg_zero_fives_top]) + cols = np.hstack([cols_first, cols_second]) + rows = np.hstack([rows, rows]) + sub_matrix_top = csr_matrix( + (data, (rows, cols)), shape=(n_lr, n_lr * n_tb) + ) + cols_first = np.arange((n_tb - 2) * (n_lr), (n_tb - 1) * (n_lr)) + cols_second = np.arange((n_tb - 1) * (n_lr), (n_tb) * (n_lr)) + data = np.hstack([neg_zero_fives_top, one_fives_top]) + cols = np.hstack([cols_first, cols_second]) + sub_matrix_bottom = csr_matrix( + (data, (rows, cols)), shape=(n_lr, n_lr * n_tb) + ) + data = np.ones((n_tb - 1) * n_lr) * 0.5 + data = np.hstack([data, data]) + rows = np.arange(0, (n_tb - 1) * n_lr) + rows = np.hstack([rows, rows]) + cols_first = np.arange(0, (n_tb - 1) * n_lr) + cols_second = np.arange(n_lr, n_lr * n_tb) + cols = np.hstack([cols_first, cols_second]) + sub_matrix_center = csr_matrix( + (data, (rows, cols)), shape=(n_lr * (n_tb - 1), n_lr * n_tb) + ) + sub_matrix = vstack( + [ + sub_matrix_top, + sub_matrix_center, + sub_matrix_bottom, + ] + ) + + elif shift_key == "edge to node": + raise NotImplementedError + else: + raise ValueError(f"shift key '{shift_key}' not recognised") + # Second dimension length + second_dim_repeats = self._get_auxiliary_domain_repeats( + discretised_symbol.domains + ) + + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), which + # is not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should not be an + # issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + + return pybamm.Matrix(matrix) @ array + + def harmonic_mean(array, direction): + """ + Calculate the harmonic mean of an array using matrix multiplication. + The harmonic mean is computed as + + .. math:: + D_{eff} = \\frac{D_1 D_2}{\\beta D_2 + (1 - \\beta) D_1}, + + where + + .. math:: + \\beta = \\frac{\\Delta x_1}{\\Delta x_2 + \\Delta x_1} + + accounts for the difference in the control volume widths. This is the + definiton from [1], which is the same as that in [2] but with slightly + different notation. + + [1] Torchio, M et al. "LIONSIMBA: A Matlab Framework Based on a Finite + Volume Model Suitable for Li-Ion Battery Design, Simulation, and Control." + (2016). + [2] Recktenwald, Gerald. "The control-volume finite-difference + approximation to the diffusion equation." (2012). + """ + # Create appropriate submesh by combining submeshes in domain + submesh = self.mesh[array.domain] + + # Get second dimension length for use later + second_dim_repeats = self._get_auxiliary_domain_repeats( + discretised_symbol.domains + ) + + # Create 1D matrix using submesh + n = submesh.npts + n_lr = submesh.npts_lr + n_tb = submesh.npts_tb + + if shift_key == "node to edge": + if direction == "lr": + # Matrix to compute values at the exterior edges + edges_sub_matrix_left = csr_matrix( + ([1.5, -0.5], ([0, 0], [0, 1])), shape=(1, n_lr) + ) + edges_sub_matrix_center = csr_matrix((n_lr - 1, n_lr)) + edges_sub_matrix_right = csr_matrix( + ([-0.5, 1.5], ([0, 0], [n_lr - 2, n_lr - 1])), shape=(1, n_lr) + ) + edges_sub_matrix = vstack( + [ + edges_sub_matrix_left, + edges_sub_matrix_center, + edges_sub_matrix_right, + ] + ) + edges_sub_matrix = block_diag((edges_sub_matrix,) * n_tb) + + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), + # which is not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should + # not be an issue + edges_matrix = csr_matrix( + kron(eye(second_dim_repeats), edges_sub_matrix) + ) + + # Matrix to extract the node values running from the first node + # to the penultimate node in the primary dimension (D_1 in the + # definiton of the harmonic mean) + sub_matrix_D1 = hstack([eye(n_lr - 1), csr_matrix((n_lr - 1, 1))]) + sub_matrix_D1 = block_diag((sub_matrix_D1,) * n_tb) + matrix_D1 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D1)) + D1 = pybamm.Matrix(matrix_D1) @ array + + # Matrix to extract the node values running from the second node + # to the final node in the primary dimension (D_2 in the + # definiton of the harmonic mean) + sub_matrix_D2 = hstack([csr_matrix((n_lr - 1, 1)), eye(n_lr - 1)]) + sub_matrix_D2 = block_diag((sub_matrix_D2,) * n_tb) + matrix_D2 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D2)) + D2 = pybamm.Matrix(matrix_D2) @ array + + # Compute weight beta + dx = submesh.d_edges_lr + sub_beta = (dx[:-1] / (dx[1:] + dx[:-1]))[:, np.newaxis] + sub_beta = np.repeat(sub_beta, n_tb, axis=0) + beta = pybamm.Array( + np.kron(np.ones((second_dim_repeats, 1)), sub_beta) + ) + + # dx_real = dx * length, therefore, beta is unchanged + # Compute harmonic mean on internal edges + D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta)) + + # Matrix to pad zeros at the beginning and end of the array where + # the exterior edge values will be added + sub_matrix = vstack( + [ + csr_matrix((1, n_lr - 1)), + eye(n_lr - 1), + csr_matrix((1, n_lr - 1)), + ] + ) + sub_matrix = block_diag((sub_matrix,) * n_tb) + + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), + # which is not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should + # not be an issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + + return ( + pybamm.Matrix(edges_matrix) @ array + + pybamm.Matrix(matrix) @ D_eff + ) + elif direction == "tb": + # Matrix to compute values at the exterior edges + # Matrix to compute values at the exterior edges + one_fives_top = np.ones(n_lr) * 1.5 + neg_zero_fives_top = np.ones(n_lr) * -0.5 + rows = np.arange(0, n_lr) + cols_first = np.arange(0, n_lr) + cols_second = np.arange(n_lr, 2 * n_lr) + data = np.hstack([one_fives_top, neg_zero_fives_top]) + cols = np.hstack([cols_first, cols_second]) + rows = np.hstack([rows, rows]) + edges_sub_matrix_top = csr_matrix( + (data, (rows, cols)), shape=(n_lr, n_lr * n_tb) + ) + cols_first = np.arange((n_tb - 2) * (n_lr), (n_tb - 1) * (n_lr)) + cols_second = np.arange((n_tb - 1) * (n_lr), (n_tb) * (n_lr)) + data = np.hstack([neg_zero_fives_top, one_fives_top]) + cols = np.hstack([cols_first, cols_second]) + edges_sub_matrix_bottom = csr_matrix( + (data, (rows, cols)), shape=(n_lr, n_lr * n_tb) + ) + edges_sub_matrix_center = csr_matrix( + ((n_tb - 1) * n_lr, n_tb * n_lr) + ) + edges_sub_matrix = vstack( + [ + edges_sub_matrix_top, + edges_sub_matrix_center, + edges_sub_matrix_bottom, + ] + ) + + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), + # which is not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should + # not be an issue + edges_matrix = csr_matrix( + kron(eye(second_dim_repeats), edges_sub_matrix) + ) + + # Matrix to extract the node values running from the first node + # to the penultimate node in the primary dimension (D_1 in the + # definiton of the harmonic mean) + sub_matrix_D1 = hstack( + [eye(n_lr * (n_tb - 1)), csr_matrix((n_lr * (n_tb - 1), n_lr))] + ) + matrix_D1 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D1)) + D1 = pybamm.Matrix(matrix_D1) @ array + + # Matrix to extract the node values running from the second node + # to the final node in the primary dimension (D_2 in the + # definiton of the harmonic mean) + sub_matrix_D2 = hstack( + [csr_matrix((n_lr * (n_tb - 1), n_lr)), eye(n_lr * (n_tb - 1))] + ) + matrix_D2 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D2)) + D2 = pybamm.Matrix(matrix_D2) @ array + + # Compute weight beta + dx = submesh.d_edges_tb + sub_beta = (dx[:-1] / (dx[1:] + dx[:-1]))[:, np.newaxis] + sub_beta = np.repeat(sub_beta, n_lr, axis=0) + beta = pybamm.Array( + np.kron(np.ones((second_dim_repeats, 1)), sub_beta) + ) + + # dx_real = dx * length, therefore, beta is unchanged + # Compute harmonic mean on internal edges + # Note: add small number to denominator to regularise D_eff + D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta)) + + # Matrix to pad zeros at the beginning and end of the array where + # the exterior edge values will be added + sub_matrix = vstack( + [ + csr_matrix((n_lr, n_lr * (n_tb - 1))), + eye(n_lr * (n_tb - 1)), + csr_matrix((n_lr, n_lr * (n_tb - 1))), + ] + ) + + # Generate full matrix from the submatrix + # Convert to csr_matrix so that we can take the index (row-slicing), + # which is not supported by the default kron format + # Note that this makes column-slicing inefficient, but this should + # not be an issue + matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix)) + + return ( + pybamm.Matrix(edges_matrix) @ array + + pybamm.Matrix(matrix) @ D_eff + ) + + elif shift_key == "edge to node": + # Matrix to extract the edge values running from the first edge + # to the penultimate edge in the primary dimension (D_1 in the + # definiton of the harmonic mean) + raise NotImplementedError + sub_matrix_D1 = hstack([eye(n), csr_matrix((n, 1))]) + matrix_D1 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D1)) + D1 = pybamm.Matrix(matrix_D1) @ array + + # Matrix to extract the edge values running from the second edge + # to the final edge in the primary dimension (D_2 in the + # definiton of the harmonic mean) + sub_matrix_D2 = hstack([csr_matrix((n, 1)), eye(n)]) + matrix_D2 = csr_matrix(kron(eye(second_dim_repeats), sub_matrix_D2)) + D2 = pybamm.Matrix(matrix_D2) @ array + + # Compute weight beta + dx0 = submesh.nodes[0] - submesh.edges[0] # first edge to node + dxN = submesh.edges[-1] - submesh.nodes[-1] # last node to edge + dx = np.concatenate(([dx0], submesh.d_nodes, [dxN])) + sub_beta = (dx[:-1] / (dx[1:] + dx[:-1]))[:, np.newaxis] + beta = pybamm.Array(np.kron(np.ones((second_dim_repeats, 1)), sub_beta)) + + # Compute harmonic mean on nodes + # Note: add small number to denominator to regularise D_eff + D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta) + 1e-16) + + return D_eff + + else: + raise ValueError(f"shift key '{shift_key}' not recognised") + + # If discretised_symbol evaluates to number there is no need to average + if discretised_symbol.size == 1: + out = discretised_symbol + elif method == "arithmetic": + out = arithmetic_mean(discretised_symbol, direction) + elif method == "harmonic": + out = harmonic_mean(discretised_symbol, direction) + else: + raise ValueError(f"method '{method}' not recognised") + return out + + def upwind_or_downwind( + self, symbol, discretised_symbol, bcs, lr_direction, tb_direction + ): + """ + Implement an upwinding operator. Currently, this requires the symbol to have + a Dirichlet boundary condition on the left side or top side (for upwinding) or right side + or bottom side (for downwinding). + + Parameters + ---------- + symbol : :class:`pybamm.SpatialVariable` + The variable to be discretised + discretised_gradient : :class:`pybamm.Vector` + Contains the discretised gradient of symbol + bcs : dict of tuples (:class:`pybamm.Scalar`, str) + Dictionary (with keys "left" and "right") of boundary conditions. Each + boundary condition consists of a value and a flag indicating its type + (e.g. "Dirichlet") + lr_direction : str + Direction in which to apply the operator (upwind or downwind) in the lr direction + tb_direction : str + Direction in which to apply the operator (upwind or downwind) in the tb direction + """ + if symbol not in bcs: + raise pybamm.ModelError( + f"Boundary conditions must be provided for {lr_direction}ing '{symbol}' and {tb_direction}ing '{symbol}'" + ) + + if lr_direction == "upwind": + bc_side = "left" + elif lr_direction == "downwind": + bc_side = "right" + else: + raise ValueError(f"direction '{lr_direction}' not recognised") + + if tb_direction == "upwind": + bc_side = "top" + elif tb_direction == "downwind": + bc_side = "bottom" + else: + raise ValueError(f"direction '{tb_direction}' not recognised") + + if bcs[symbol][bc_side][1] != "Dirichlet": + raise pybamm.ModelError( + "Dirichlet boundary conditions must be provided for " + f"{lr_direction}ing '{symbol}' and {tb_direction}ing '{symbol}'" + ) + + # Extract only the relevant boundary condition as the model might have both + bc_subset = {bc_side: bcs[symbol][bc_side]} + symbol_out, _ = self.add_ghost_nodes(symbol, discretised_symbol, bc_subset) + return symbol_out diff --git a/src/pybamm/spatial_methods/scikit_finite_element.py b/src/pybamm/spatial_methods/scikit_finite_element.py index 23d30e93eb..d5a26fe7ed 100644 --- a/src/pybamm/spatial_methods/scikit_finite_element.py +++ b/src/pybamm/spatial_methods/scikit_finite_element.py @@ -283,7 +283,9 @@ def stiffness_form(u, v, w): return pybamm.Matrix(stiffness) - def integral(self, child, discretised_child, integration_dimension): + def integral( + self, child, discretised_child, integration_dimension, integration_variable + ): """Vector-vector dot product to implement the integral operator. See :meth:`pybamm.SpatialMethod.integral` """ diff --git a/tests/__init__.py b/tests/__init__.py index 8227cd0650..445f60a360 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -24,6 +24,7 @@ from .shared import ( get_mesh_for_testing, + get_mesh_for_testing_2d, get_mesh_for_testing_symbolic, get_p2d_mesh_for_testing, get_size_distribution_mesh_for_testing, diff --git a/tests/shared.py b/tests/shared.py index 1af9152cd3..77a1223c1d 100644 --- a/tests/shared.py +++ b/tests/shared.py @@ -119,6 +119,86 @@ def get_mesh_for_testing( return pybamm.Mesh(geometry, submesh_types, var_pts) +def get_mesh_for_testing_2d( + xpts=None, + rpts=10, + Rpts=10, + ypts=15, + zpts=15, +): + param = pybamm.ParameterValues( + values={ + "Electrode height [m]": 0.5, + "Negative electrode thickness [m]": 1 / 3, + "Separator thickness [m]": 1 / 3, + "Positive electrode thickness [m]": 1 / 3, + "Negative particle radius [m]": 0.5, + "Positive particle radius [m]": 0.5, + } + ) + + x = pybamm.SpatialVariable( + "x", ["negative electrode", "separator", "positive electrode"], direction="lr" + ) + z = pybamm.SpatialVariable( + "z", ["negative electrode", "separator", "positive electrode"], direction="tb" + ) + + geometry = { + "negative electrode": { + x: { + "min": pybamm.Scalar(0), + "max": pybamm.Parameter("Negative electrode thickness [m]"), + }, + z: { + "min": pybamm.Scalar(0), + "max": pybamm.Parameter("Electrode height [m]"), + }, + }, + "separator": { + x: { + "min": pybamm.Parameter("Negative electrode thickness [m]"), + "max": pybamm.Parameter("Separator thickness [m]") + + pybamm.Parameter("Negative electrode thickness [m]"), + }, + z: { + "min": pybamm.Scalar(0), + "max": pybamm.Parameter("Electrode height [m]"), + }, + }, + "positive electrode": { + x: { + "min": pybamm.Parameter("Separator thickness [m]") + + pybamm.Parameter("Negative electrode thickness [m]"), + "max": pybamm.Parameter("Positive electrode thickness [m]") + + pybamm.Parameter("Separator thickness [m]") + + pybamm.Parameter("Negative electrode thickness [m]"), + }, + z: { + "min": pybamm.Scalar(0), + "max": pybamm.Parameter("Electrode height [m]"), + }, + }, + } + param.process_geometry(geometry) + + submesh_types = { + "negative electrode": pybamm.Uniform2DSubMesh, + "separator": pybamm.Uniform2DSubMesh, + "positive electrode": pybamm.Uniform2DSubMesh, + } + + if xpts is None: + xn_pts = 40 + else: + xn_pts = xpts + var_pts = { + x: xn_pts, + z: zpts, + } + return pybamm.Mesh(geometry, submesh_types, var_pts) + + def get_p2d_mesh_for_testing(xpts=None, rpts=10): geometry = pybamm.battery_geometry() return get_mesh_for_testing(xpts=xpts, rpts=rpts, geometry=geometry) diff --git a/tests/unit/test_meshes/test_two_dimensional_submesh.py b/tests/unit/test_meshes/test_two_dimensional_submesh.py new file mode 100644 index 0000000000..2cf011065d --- /dev/null +++ b/tests/unit/test_meshes/test_two_dimensional_submesh.py @@ -0,0 +1,54 @@ +import pytest +import pybamm + + +@pytest.fixture() +def x(): + x = pybamm.SpatialVariable("x", domain=["my 2d domain"], coord_sys="cartesian") + return x + + +@pytest.fixture() +def y(): + return pybamm.SpatialVariable("y", domain=["my 2d domain"], coord_sys="cartesian") + + +@pytest.fixture() +def geometry(x, y): + geometry = { + "my 2d domain": { + x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}, + y: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}, + } + } + return geometry + + +class TestUniform2DSubMesh: + def test_exceptions(self): + lims = {"a": 1} + with pytest.raises(pybamm.GeometryError): + pybamm.Uniform2DSubMesh(lims, None) + + def test_symmetric_mesh_creation_no_parameters(self, x, y, geometry): + submesh_types = {"my 2d domain": pybamm.Uniform2DSubMesh} + var_pts = {x: 20, y: 20} + + # create mesh + mesh = pybamm.Mesh(geometry, submesh_types, var_pts) + + # check boundary locations + assert mesh["my 2d domain"].edges_lr[0] == 0 + assert mesh["my 2d domain"].edges_lr[-1] == 1 + assert mesh["my 2d domain"].edges_tb[0] == 0 + assert mesh["my 2d domain"].edges_tb[-1] == 1 + + # check number of edges and nodes + assert len(mesh["my 2d domain"].nodes_lr) == var_pts[x] + assert len(mesh["my 2d domain"].nodes_tb) == var_pts[y] + assert ( + len(mesh["my 2d domain"].edges_lr) == len(mesh["my 2d domain"].nodes_lr) + 1 + ) + assert ( + len(mesh["my 2d domain"].edges_tb) == len(mesh["my 2d domain"].nodes_tb) + 1 + ) diff --git a/tests/unit/test_spatial_methods/test_finite_volume_2d/test_2d_finite_volume_integration.py b/tests/unit/test_spatial_methods/test_finite_volume_2d/test_2d_finite_volume_integration.py new file mode 100644 index 0000000000..ae632f7d98 --- /dev/null +++ b/tests/unit/test_spatial_methods/test_finite_volume_2d/test_2d_finite_volume_integration.py @@ -0,0 +1,64 @@ +# +# Tests for integration using Finite Volume method +# + +import pybamm +import numpy as np + +from tests import get_mesh_for_testing_2d + + +class TestFiniteVolumeIntegration: + def test_definite_integral(self): + # create discretisation + mesh = get_mesh_for_testing_2d() + spatial_methods = { + "negative electrode": pybamm.FiniteVolume2D(), + "separator": pybamm.FiniteVolume2D(), + "positive electrode": pybamm.FiniteVolume2D(), + } + disc = pybamm.Discretisation(mesh, spatial_methods) + # lengths + ln = mesh["negative electrode"].edges_lr[-1] + ls = mesh["separator"].edges_lr[-1] - ln + l_tb = mesh["negative electrode"].edges_tb[-1] + + # macroscale variable (lr integration) + var = pybamm.Variable("var", domain=["negative electrode", "separator"]) + x = pybamm.SpatialVariable( + "x", ["negative electrode", "separator"], direction="lr" + ) + integral_eqn = pybamm.Integral(var, x) + disc.set_variable_slices([var]) + integral_eqn_disc = disc.process_symbol(integral_eqn) + submesh = mesh[("negative electrode", "separator")] + constant_y = np.ones(submesh.npts_lr * submesh.npts_tb) + assert (integral_eqn_disc.evaluate(None, constant_y) == ln + ls).all() + + var = pybamm.Variable("var", domain=["negative electrode", "separator"]) + z = pybamm.SpatialVariable( + "z", ["negative electrode", "separator"], direction="tb" + ) + integral_eqn_tb = pybamm.Integral(var, z) + disc.set_variable_slices([var]) + integral_eqn_disc_tb = disc.process_symbol(integral_eqn_tb) + submesh = mesh[("negative electrode", "separator")] + constant_y = np.ones(submesh.npts_lr * submesh.npts_tb) + np.testing.assert_allclose( + integral_eqn_disc_tb.evaluate(None, constant_y), l_tb + ) + + var = pybamm.Variable("var", domain=["negative electrode"]) + x = pybamm.SpatialVariable("x", ["negative electrode"], direction="lr") + integral_eqn = pybamm.Integral(var, x) + disc.set_variable_slices([var]) + integral_eqn_disc = disc.process_symbol(integral_eqn) + submesh = mesh["negative electrode"] + LR, TB = np.meshgrid(submesh.edges_lr, submesh.edges_tb) + lr = LR.flatten() + np.testing.assert_allclose( + integral_eqn_disc.evaluate(None, lr), + (ln + ls) ** 2 / 2, + rtol=1e-7, + atol=1e-6, + ) diff --git a/tests/unit/test_spatial_methods/test_finite_volume_2d/test_finite_volume_2d.py b/tests/unit/test_spatial_methods/test_finite_volume_2d/test_finite_volume_2d.py new file mode 100644 index 0000000000..2948c74702 --- /dev/null +++ b/tests/unit/test_spatial_methods/test_finite_volume_2d/test_finite_volume_2d.py @@ -0,0 +1,58 @@ +import pybamm +import numpy as np +import pytest + +from tests import get_mesh_for_testing_2d + + +class TestFiniteVolume2D: + def test_node_to_edge_to_node(self): + # Create discretisation + mesh = get_mesh_for_testing_2d() + fin_vol = pybamm.FiniteVolume2D() + fin_vol.build(mesh) + n_lr = mesh["negative electrode"].npts_lr + n_tb = mesh["negative electrode"].npts_tb + + # node to edge + c = pybamm.StateVector(slice(0, n_lr * n_tb), domain=["negative electrode"]) + y_test = np.ones(n_lr * n_tb) * 2 + diffusivity_c_ari_lr = fin_vol.node_to_edge( + c, method="arithmetic", direction="lr" + ) + np.testing.assert_array_equal( + diffusivity_c_ari_lr.evaluate(None, y_test), + np.ones(((n_lr + 1) * (n_tb), 1)) * 2, + ) + diffusivity_c_har_lr = fin_vol.node_to_edge( + c, method="harmonic", direction="lr" + ) + np.testing.assert_array_equal( + diffusivity_c_har_lr.evaluate(None, y_test), + np.ones(((n_lr + 1) * (n_tb), 1)) * 2, + ) + diffusivity_c_ari_tb = fin_vol.node_to_edge( + c, method="arithmetic", direction="tb" + ) + np.testing.assert_array_equal( + diffusivity_c_ari_tb.evaluate(None, y_test), + np.ones(((n_lr) * (n_tb + 1), 1)) * 2, + ) + diffusivity_c_har_tb = fin_vol.node_to_edge( + c, method="harmonic", direction="tb" + ) + np.testing.assert_array_equal( + diffusivity_c_har_tb.evaluate(None, y_test), + np.ones(((n_lr) * (n_tb + 1), 1)) * 2, + ) + + # bad shift key + with pytest.raises(ValueError, match="shift key"): + fin_vol.shift(c, "bad shift key", "arithmetic") + + with pytest.raises(ValueError, match="shift key"): + fin_vol.shift(c, "bad shift key", "harmonic") + + # bad method + with pytest.raises(ValueError, match="method"): + fin_vol.shift(c, "shift key", "bad method")