diff --git a/examples/scripts/flux_bc.py b/examples/scripts/flux_bc.py new file mode 100644 index 0000000000..f98becbb1b --- /dev/null +++ b/examples/scripts/flux_bc.py @@ -0,0 +1,124 @@ +from copy import deepcopy + +import matplotlib.pyplot as plt +import numpy as np + +import pybamm + +pybamm.set_logging_level("DEBUG") + +model = pybamm.BaseModel() + +R = pybamm.Parameter("Particle radius") +D = pybamm.Parameter("Diffusion coefficient") +c0 = pybamm.Parameter("Initial concentration") + +c = pybamm.Variable("Concentration", domain="negative particle") + +# governing equations +N = -D * pybamm.grad(c) # flux +N.set_do_not_simplify() + +dcdt = -pybamm.div(N) +print("automatic sign handling", dcdt) +model.rhs = {c: dcdt} + + +# initial conditions +model.initial_conditions = {c: c0} + +model.variables = { + "Concentration": c, + "Flux": N, + "Surface concentration": pybamm.surf(c), +} + +model2 = deepcopy(model) + +# boundary conditions +lbc = pybamm.Scalar(0) +rbc = pybamm.Scalar(1) +model.boundary_conditions = { + c: { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(-2), ("Flux", N)), + } +} + +param = pybamm.ParameterValues( + { + "Particle radius": 10e-6, + "Diffusion coefficient": 3.9e-14, + "Initial concentration": 2.5e4, + } +) + +r = pybamm.SpatialVariable( + "r", domain=["negative particle"], coord_sys="spherical polar" +) +geometry = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": R}}} + + +param.process_model(model) +param.process_geometry(geometry) + +submesh_types = {"negative particle": pybamm.Uniform1DSubMesh} +var_pts = {r: 20} +mesh = pybamm.Mesh(geometry, submesh_types, var_pts) + +spatial_methods = {"negative particle": pybamm.FiniteVolume()} +disc = pybamm.Discretisation(mesh, spatial_methods) +disc.process_model(model) + +# solve +solver = pybamm.ScipySolver() +t = np.linspace(0, 3600, 600) +solution = solver.solve(model, t) + +# post-process, so that the solution can be called at any time t or spaceå r +# (using interpolation) +c_sol = solution["Concentration"] +c_surf = solution["Surface concentration"] + +# plot +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4)) + +ax1.plot(solution.t, c_surf(solution.t)) +ax1.set_xlabel("Time [s]") +ax1.set_ylabel("Surface concentration") + +r = mesh["negative particle"].nodes # radial position +time = 1000 # time in seconds +ax2.plot(r * 1e6, c_sol(t=time, r=r), label=f"t={time}[s]") +ax2.set_xlabel("Particle radius") +ax2.set_ylabel("Concentration") + + +# Comparison with Neumann +model2.boundary_conditions = { + c: { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(2) / D, "Neumann"), + } +} + +param.process_model(model2) +disc.process_model(model2) + +# solve +solver = pybamm.ScipySolver() +t = np.linspace(0, 3600, 600) +solution = solver.solve(model2, t) + +# post-process, so that the solution can be called at any time t or spaceå r +# (using interpolation) +c = solution["Concentration"] +c_surf = solution["Surface concentration"] + +# Plot +ax1.plot(solution.t, c_surf(solution.t), "--") +ax2.plot(r * 1e6, c(t=time, r=r), "--", label=f"t={time}[s]") +ax2.legend() + +plt.tight_layout() +plt.show() diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py index e7fb211457..45c7593137 100644 --- a/src/pybamm/__init__.py +++ b/src/pybamm/__init__.py @@ -7,6 +7,7 @@ root_dir, load, is_constant_and_can_evaluate, + is_flux_boundary_condition, ) from .util import ( get_parameters_filepath, diff --git a/src/pybamm/discretisations/discretisation.py b/src/pybamm/discretisations/discretisation.py index 4f5562225b..7ed5e7746b 100644 --- a/src/pybamm/discretisations/discretisation.py +++ b/src/pybamm/discretisations/discretisation.py @@ -951,13 +951,13 @@ def _process_symbol(self, symbol): new_children=[disc_left.tb_field, disc_right.tb_field] ) ) - return pybamm.VectorField(disc_lr, disc_tb) - - return pybamm.simplify_if_constant( - symbol.create_copy(new_children=[disc_left, disc_right]) - ) + processed_symbol = pybamm.VectorField(disc_lr, disc_tb) + else: + processed_symbol = pybamm.simplify_if_constant( + symbol.create_copy(new_children=[disc_left, disc_right]) + ) else: - return spatial_method.process_binary_operators( + processed_symbol = spatial_method.process_binary_operators( symbol, left, right, disc_left, disc_right ) elif isinstance(symbol, pybamm._BaseAverage): @@ -974,7 +974,7 @@ def _process_symbol(self, symbol): x = symbol.integration_variable v = pybamm.ones_like(child) average = pybamm.Integral(child, x) / pybamm.Integral(v, x) - return self.process_symbol(average) + processed_symbol = self.process_symbol(average) elif isinstance(symbol, pybamm.UnaryOperator): child = symbol.child @@ -984,31 +984,39 @@ def _process_symbol(self, symbol): child_spatial_method = self.spatial_methods[child.domain[0]] if isinstance(symbol, pybamm.Gradient): - return child_spatial_method.gradient(child, disc_child, self.bcs) + processed_symbol = child_spatial_method.gradient( + child, disc_child, self.bcs + ) elif isinstance(symbol, pybamm.Divergence): - return child_spatial_method.divergence(child, disc_child, self.bcs) + processed_symbol = child_spatial_method.divergence( + child, disc_child, self.bcs + ) elif isinstance(symbol, pybamm.Laplacian): - return child_spatial_method.laplacian(child, disc_child, self.bcs) + processed_symbol = child_spatial_method.laplacian( + child, disc_child, self.bcs + ) elif isinstance(symbol, pybamm.GradientSquared): - return child_spatial_method.gradient_squared( + processed_symbol = child_spatial_method.gradient_squared( child, disc_child, self.bcs ) elif isinstance(symbol, pybamm.Mass): - return child_spatial_method.mass_matrix(child, self.bcs) + processed_symbol = child_spatial_method.mass_matrix(child, self.bcs) elif isinstance(symbol, pybamm.BoundaryMass): - return child_spatial_method.boundary_mass_matrix(child, self.bcs) + processed_symbol = child_spatial_method.boundary_mass_matrix( + child, self.bcs + ) elif isinstance(symbol, pybamm.IndefiniteIntegral): - return child_spatial_method.indefinite_integral( + processed_symbol = child_spatial_method.indefinite_integral( child, disc_child, "forward" ) elif isinstance(symbol, pybamm.BackwardIndefiniteIntegral): - return child_spatial_method.indefinite_integral( + processed_symbol = child_spatial_method.indefinite_integral( child, disc_child, "backward" ) @@ -1023,10 +1031,10 @@ def _process_symbol(self, symbol): symbol.integration_variable, ) out.copy_domains(symbol) - return out + processed_symbol = out elif isinstance(symbol, pybamm.DefiniteIntegralVector): - return child_spatial_method.definite_integral_matrix( + processed_symbol = child_spatial_method.definite_integral_matrix( child, vector_type=symbol.vector_type ) @@ -1034,7 +1042,7 @@ def _process_symbol(self, symbol): child_spatial_method = self.spatial_methods[ symbol.integration_domain[0] ] - return child_spatial_method.one_dimensional_integral( + processed_symbol = child_spatial_method.one_dimensional_integral( symbol, child, disc_child, @@ -1042,19 +1050,19 @@ def _process_symbol(self, symbol): symbol.direction, ) elif isinstance(symbol, pybamm.BoundaryIntegral): - return child_spatial_method.boundary_integral( + processed_symbol = child_spatial_method.boundary_integral( child, disc_child, symbol.region ) elif isinstance(symbol, pybamm.Broadcast): # Broadcast new_child to the domain specified by symbol.domain # Different discretisations may broadcast differently - return spatial_method.broadcast( + processed_symbol = spatial_method.broadcast( disc_child, symbol.domains, symbol.broadcast_type ) elif isinstance(symbol, pybamm.DeltaFunction): - return spatial_method.delta_function(symbol, disc_child) + processed_symbol = spatial_method.delta_function(symbol, disc_child) elif isinstance(symbol, pybamm.BoundaryOperator): # if boundary operator applied on "negative tab" or @@ -1064,15 +1072,15 @@ def _process_symbol(self, symbol): mesh = self.mesh[symbol.children[0].domain[0]] if isinstance(mesh, pybamm.SubMesh1D): symbol.side = mesh.tabs[symbol.side] - return child_spatial_method.boundary_value_or_flux( + processed_symbol = child_spatial_method.boundary_value_or_flux( symbol, disc_child, self.bcs ) elif isinstance(symbol, pybamm.EvaluateAt): - return child_spatial_method.evaluate_at( + processed_symbol = child_spatial_method.evaluate_at( symbol, disc_child, symbol.position ) elif isinstance(symbol, pybamm.UpwindDownwind2D): - return spatial_method.upwind_or_downwind( + processed_symbol = spatial_method.upwind_or_downwind( child, disc_child, self.bcs, @@ -1080,46 +1088,46 @@ def _process_symbol(self, symbol): symbol.tb_direction, ) elif isinstance(symbol, pybamm.NodeToEdge2D): - return spatial_method.node_to_edge( + processed_symbol = spatial_method.node_to_edge( disc_child, method="arithmetic", direction=symbol.direction, ) elif isinstance(symbol, pybamm.UpwindDownwind): direction = symbol.name # upwind or downwind - return spatial_method.upwind_or_downwind( + processed_symbol = spatial_method.upwind_or_downwind( child, disc_child, self.bcs, direction ) elif isinstance(symbol, pybamm.NotConstant): # After discretisation, we can make the symbol constant - return disc_child + processed_symbol = 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 + processed_symbol = disc_child.lr_field elif direction == "tb": - return disc_child.tb_field + processed_symbol = disc_child.tb_field else: raise ValueError("Invalid direction") else: if isinstance(disc_child, pybamm.VectorField): - return pybamm.VectorField( + processed_symbol = pybamm.VectorField( symbol.create_copy(new_children=[disc_child.lr_field]), symbol.create_copy(new_children=[disc_child.tb_field]), ) else: - return symbol.create_copy(new_children=[disc_child]) + processed_symbol = symbol.create_copy(new_children=[disc_child]) elif isinstance(symbol, (pybamm.Function, pybamm.Conditional)): disc_children = [self.process_symbol(child) for child in symbol.children] - return symbol.create_copy(disc_children) + processed_symbol = symbol.create_copy(disc_children) elif isinstance(symbol, pybamm.VariableDot): # Add symbol's reference and multiply by the symbol's scale # so that the state vector is of order 1 - return symbol.reference + symbol.scale * pybamm.StateVectorDot( + processed_symbol = symbol.reference + symbol.scale * pybamm.StateVectorDot( *self.y_slices[symbol.get_variable()], domains=symbol.domains, ) @@ -1138,12 +1146,12 @@ def _process_symbol(self, symbol): ) from error # Add symbol's reference and multiply by the symbol's scale # so that the state vector is of order 1 - return symbol.reference + symbol.scale * pybamm.StateVector( + processed_symbol = symbol.reference + symbol.scale * pybamm.StateVector( *y_slices, domains=symbol.domains ) elif isinstance(symbol, pybamm.SpatialVariable): - return spatial_method.spatial_variable(symbol) + processed_symbol = spatial_method.spatial_variable(symbol) elif isinstance(symbol, pybamm.ConcatenationVariable): # create new children without scale and reference @@ -1159,12 +1167,12 @@ def _process_symbol(self, symbol): self.y_slices = old_y_slices new_symbol = spatial_method.concatenation(new_children) # apply scale to the whole concatenation - return symbol.reference + symbol.scale * new_symbol + processed_symbol = symbol.reference + symbol.scale * new_symbol elif isinstance(symbol, pybamm.Concatenation): new_children = [self.process_symbol(child) for child in symbol.children] new_symbol = spatial_method.concatenation(new_children) - return new_symbol + processed_symbol = new_symbol elif isinstance(symbol, pybamm.InputParameter): if symbol.domain != []: @@ -1173,7 +1181,7 @@ def _process_symbol(self, symbol): expected_size = None if symbol._expected_size is None: symbol._expected_size = expected_size - return symbol.create_copy() + processed_symbol = symbol.create_copy() elif isinstance(symbol, pybamm.CoupledVariable): raise pybamm.DiscretisationError( @@ -1185,19 +1193,37 @@ def _process_symbol(self, symbol): # VectorField is a subclass of TensorField, handle it first for specificity left_symbol = self.process_symbol(symbol.lr_field) right_symbol = self.process_symbol(symbol.tb_field) - return symbol.create_copy(new_children=[left_symbol, right_symbol]) + processed_symbol = symbol.create_copy( + new_children=[left_symbol, right_symbol] + ) elif isinstance(symbol, pybamm.TensorField): # General TensorField handling (rank-2 tensors) processed_children = [self.process_symbol(c) for c in symbol.children] - return symbol.create_copy(new_children=processed_children) + processed_symbol = symbol.create_copy(new_children=processed_children) elif isinstance(symbol, pybamm.Constant): # after discretisation we just care about the value, not the name - return self.process_symbol(pybamm.Scalar(symbol.value)) + processed_symbol = self.process_symbol(pybamm.Scalar(symbol.value)) else: # Backup option: return the object - return symbol + processed_symbol = symbol + + if self.bcs: + key_id = next(iter(self.bcs.keys())) + for bc in self.bcs[key_id].values(): + if pybamm.is_flux_boundary_condition(bc[1]) and bc[1][1] == symbol: + if not isinstance(spatial_method, pybamm.FiniteVolume): + raise ValueError( + "Flux boundary conditions are only implemented for 1D finite volumes." + ) + else: + processed_symbol = spatial_method.add_flux_values( + symbol, processed_symbol, self.bcs[key_id] + ) + return processed_symbol + + return processed_symbol def concatenate(self, *symbols, sparse=False): if sparse: diff --git a/src/pybamm/expression_tree/symbol.py b/src/pybamm/expression_tree/symbol.py index e3c39db6b9..d2b1d122a0 100644 --- a/src/pybamm/expression_tree/symbol.py +++ b/src/pybamm/expression_tree/symbol.py @@ -167,6 +167,9 @@ def simplify_if_constant(symbol: pybamm.Symbol): Utility function to simplify an expression tree if it evaluates to a constant scalar, vector or matrix. Also handles TensorField by simplifying each component. """ + if symbol.do_not_simplify: + return symbol + # Handle TensorField by simplifying each component if isinstance(symbol, pybamm.TensorField): simplified_children = [simplify_if_constant(child) for child in symbol.children] @@ -264,6 +267,8 @@ def __init__( ): self.test_shape() + self._do_not_simplify = False + @classmethod def _from_json(cls, snippet: dict): """ @@ -445,6 +450,13 @@ def read_domain_or_domains( ) return domains + @property + def do_not_simplify(self): + return self._do_not_simplify + + def set_do_not_simplify(self): + self._do_not_simplify = True + @property def id(self): return self._id diff --git a/src/pybamm/expression_tree/unary_operators.py b/src/pybamm/expression_tree/unary_operators.py index d11bca032d..3a7ab75198 100644 --- a/src/pybamm/expression_tree/unary_operators.py +++ b/src/pybamm/expression_tree/unary_operators.py @@ -1498,12 +1498,13 @@ def div(symbol): new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain) return pybamm.PrimaryBroadcast(new_child, symbol.domain) # Divergence commutes with Negate operator - if isinstance(symbol, pybamm.Negate): - return -div(symbol.orphans[0]) - elif isinstance(symbol, pybamm.Multiplication | pybamm.Division): - left, right = symbol.orphans - if isinstance(left, pybamm.Negate): - return -div(symbol._binary_new_copy(left.orphans[0], right)) + if not symbol.do_not_simplify: + if isinstance(symbol, pybamm.Negate): + return -div(symbol.orphans[0]) + elif isinstance(symbol, pybamm.Multiplication | pybamm.Division): + left, right = symbol.orphans + if isinstance(left, pybamm.Negate): + return -div(symbol._binary_new_copy(left.orphans[0], right)) # Last resort return Divergence(symbol) diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py index ba710211e1..06ace7f971 100644 --- a/src/pybamm/models/base_model.py +++ b/src/pybamm/models/base_model.py @@ -2462,9 +2462,12 @@ def check_and_convert_bcs(self, boundary_conditions): eqn, typ = boundary_conditions[var][side] boundary_conditions[var][side] = (pybamm.Scalar(eqn), typ) # Check types - if bc[1] not in ["Dirichlet", "Neumann"]: + if bc[1] not in [ + "Dirichlet", + "Neumann", + ] and not pybamm.is_flux_boundary_condition(bc[1]): raise pybamm.ModelError( - f"boundary condition types must be Dirichlet or Neumann, not '{bc[1]}'" + f"boundary condition types must be Dirichlet, Flux or Neumann, not '{bc[1]}'" ) return boundary_conditions diff --git a/src/pybamm/parameters/parameter_substitutor.py b/src/pybamm/parameters/parameter_substitutor.py index 796309c6df..0e75d7fa24 100644 --- a/src/pybamm/parameters/parameter_substitutor.py +++ b/src/pybamm/parameters/parameter_substitutor.py @@ -571,6 +571,8 @@ def process_boundary_conditions( pybamm.logger.verbose( f"Processing parameters for {variable!r} ({side} bc)" ) + if pybamm.is_flux_boundary_condition(typ): + typ = (typ[0], self.process_symbol(typ[1])) processed_bc = (self.process_symbol(bc), typ) new_boundary_conditions[processed_variable][side] = processed_bc except KeyError as err: diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index e9053d5f15..cfa5ffba4d 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -110,7 +110,10 @@ def gradient(self, symbol, discretised_symbol, boundary_conditions): # Add Neumann boundary conditions, if defined if symbol in boundary_conditions: bcs = boundary_conditions[symbol] - if any(bc[1] == "Neumann" for bc in bcs.values()): + if any( + bc[1] == "Neumann" or pybamm.is_flux_boundary_condition(bc[1]) + for bc in bcs.values() + ): out = self.add_neumann_values(symbol, out, bcs, domain) return out @@ -854,7 +857,9 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs): else: left_ghost_constant = 2 * lbc_value lbc_vector = pybamm.Matrix(lbc_matrix) @ left_ghost_constant - elif lbc_type in ["Neumann", None]: + elif lbc_type in ["Neumann", None] or pybamm.is_flux_boundary_condition( + lbc_type + ): lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -875,7 +880,9 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs): else: right_ghost_constant = 2 * rbc_value rbc_vector = pybamm.Matrix(rbc_matrix) @ right_ghost_constant - elif rbc_type in ["Neumann", None]: + elif rbc_type in ["Neumann", None] or pybamm.is_flux_boundary_condition( + rbc_type + ): rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -955,9 +962,9 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): # Count number of Neumann boundary conditions n_bcs = 0 - if lbc_type == "Neumann": + if lbc_type == "Neumann" or pybamm.is_flux_boundary_condition(lbc_type): n_bcs += 1 - if rbc_type == "Neumann": + if rbc_type == "Neumann" or pybamm.is_flux_boundary_condition(rbc_type): n_bcs += 1 # Add any values from Neumann boundary conditions to the bcs vector @@ -971,7 +978,11 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): 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): + elif ( + lbc_type == "Dirichlet" + or (lbc_type == "Neumann" and lbc_value == 0) + or pybamm.is_flux_boundary_condition(lbc_type) + ): lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -989,7 +1000,11 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): 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): + elif ( + rbc_type == "Dirichlet" + or (rbc_type == "Neumann" and rbc_value == 0) + or pybamm.is_flux_boundary_condition(rbc_type) + ): rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -1006,11 +1021,11 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): # 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": + if lbc_type == "Neumann" or pybamm.is_flux_boundary_condition(lbc_type): left_vector = csr_matrix((1, n)) else: left_vector = None - if rbc_type == "Neumann": + if rbc_type == "Neumann" or pybamm.is_flux_boundary_condition(rbc_type): right_vector = csr_matrix((1, n)) else: right_vector = None @@ -1027,6 +1042,99 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain): return new_gradient + def add_flux_values(self, symbol, processed_symbol, bcs): + """ + Add the known values of the flux boundary conditions to + the discretised symbol. + + Flux boundary are a variation of Neumann boundary conditions that + instead of being directly implemented for the gradient, + are given terms of an expression(gradient). Add the gradient level, + homogeneous Neumann boundary conditions are added to ensure correct + dimensions. + + Dirichlet bcs are implemented using ghost nodes, see + :meth:`pybamm.FiniteVolume.add_ghost_nodes`. + + Parameters + ---------- + symbol : :class:`pybamm.Symbol` + processed_symbol : :class:`pybamm.Symbol` + Contains the discretised symbol without boundary conditions + 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") + + Returns + ------- + :class:`pybamm.Symbol` + `processed_symbol + bcs_vector`. + + """ + domain = symbol.domain + n_bcs = 0 + lbc_value, lbc_type = bcs["left"] + rbc_value, rbc_type = bcs["right"] + + # Count number of Neumann and flux boundary conditions and adjust domain for dirichlet conditions + if lbc_type == "Neumann" or pybamm.is_flux_boundary_condition(lbc_type): + n_bcs += 1 + elif lbc_type == "Dirichlet": + domain = [domain[0] + "_left ghost cell", *domain] + + if rbc_type == "Neumann" or pybamm.is_flux_boundary_condition(rbc_type): + n_bcs += 1 + elif rbc_type == "Dirichlet": + domain = [*domain, domain[-1] + "_right ghost cell"] + + # get relevant grid points + submesh = self.mesh[domain] + + # Prepare sizes and empty bcs_vector + n = submesh.npts - 1 + second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) + + # Add any values from flux boundary conditions to the bcs vector + if pybamm.is_flux_boundary_condition(lbc_type) and lbc_value != 0: + lbc_sub_matrix = coo_matrix(([1.0], ([0], [0])), shape=(n + n_bcs, 1)) + lbc_matrix = csr_matrix( + kron(eye(second_dim_repeats, dtype=np.float64), lbc_sub_matrix) + ) + 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 + else: + lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) + + if pybamm.is_flux_boundary_condition(rbc_type) and rbc_value != 0: + rbc_sub_matrix = coo_matrix( + ([1.0], ([n + n_bcs - 1], [0])), shape=(n + n_bcs, 1) + ) + rbc_matrix = csr_matrix( + kron(eye(second_dim_repeats, dtype=np.float64), rbc_sub_matrix) + ) + 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 + else: + rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) + + bcs_vector = lbc_vector + rbc_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(processed_symbol) + + new_processed_symbol = processed_symbol + bcs_vector + + return new_processed_symbol + def _get_boundary_submesh_length(self, side: str, domains: list[str]): if side == "left": return self.mesh[domains[0]].length @@ -1057,7 +1165,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): See :meth:`pybamm.SpatialMethod.boundary_value` """ - + pybamm.logger.debug("TODO") # Find the number of submeshes submesh = self.mesh[discretised_child.domain] diff --git a/src/pybamm/util.py b/src/pybamm/util.py index 70d8849729..f7dca78b4c 100644 --- a/src/pybamm/util.py +++ b/src/pybamm/util.py @@ -418,3 +418,11 @@ def import_optional_dependency(module_name, attribute=None): except ModuleNotFoundError as error: # Raise an ModuleNotFoundError if the module or attribute is not available raise ModuleNotFoundError(err_msg) from error + + +def is_flux_boundary_condition(bc_type): + return ( + isinstance(bc_type, tuple) + and bc_type[0] == "Flux" + and isinstance(bc_type[1], pybamm.Symbol) + )