pybamm-team
Evaluate boundary values of edge-defined variables
Description
At present, if calling pybamm.BoundaryValue with a variable that is defined on cell edges, it fails with the error "Can't take the boundary value of a symbol that evaluates on edges". However, as the boundary value is by definition a cell edge, it should be straightforward to do this without requiring any extrapolation.
Motivation
Trying to solve hyperbolic PDEs, currently the finite volume method is not conservative when the flux leaving one domain is deposited in another.
Possible Implementation
The function "boundary_value" in unary_operators.py, inside the if statement "if symbol.evaluates_on_edges("primary"):"
Additional context
No response
I'm a bit unclear on what you want to do here. It sounds like you want the actual value of the symbol, not the boundary value. Say you have a 1d domain:
|-----------------|
and you have 2 variables, A, and B. A is defined on the bulk of the 1D domain, but B is only defined on the RHS boundary. Then A will have two boundary values, one on the lhs and one on the rhs. But B doesn't even have a boundary value, since it is just a 0D point. So I don't understand what boundary value you would like to get.
From the discussion I had with Will & Dharshannan, I think the main issue was that they needed to obtain the boundary value of a variable to couple with another domain, but that variable was defined on edges. My guess was that there was nothing preventing the boundary value function to work on an edge-evaluated quantity. This would just mean returning the value, no interpolation needed.
but the boundary value expression shouldn't return the value/field itself, it should return a boundary value. In the example case above there is no boundary value, so I would say an error is expected. However, if its a 2D domain:
______________________
|----------------------|
|----------------------|
|----------------------|
|----------------------|
-----------------------
If B is defined on the RHS boundary only, then it has two natural boundary values (the top right corner and the bottom right). Please excuse my awful ascii art......
Perhaps I used the wrong terminology. What we need to evaluate is not a boundary value in the sense of a boundary condition, but in the sense of an edge-defined variable evaluated on the final edge, i.e. an equivalent of PyBaMM.EvaluateAt but that simply takes the value on the final edge (which obviously coincides with the boundary) with no interpolation/extrapolation required
ok, I think it was my bad, so the edge-defined variable is defined across the whole domain (including the boundary), and yes you are right that the boundary value should just return the value at the boundary edge, no interpolation required.
Hi, any updates on this issue? Would it be possible to add an evaluate_on_edge attribute to the pybamm.BoundaryValue class, and if true, in the FiniteVolume class boundary value function just take the value of that boundary edge. (Not sure if this will work, @martinjrobins please advise if this is feasible).
I don't mind adding these changes, but would definitely like to have some supervision just to make sure nothing breaks 😅. @brosaplanella , @martinjrobins , @rtimms would anyone be free to hop on a call and discuss the changes and I can add them from my forked repo?
Hi @martinjrobins , @rtimms , @brosaplanella , would a simple addition such as the below on the boundary_value_or_flux function in finite_volume.py work (also removing the ValueError raised upstream in the unary_operators script within the boundary_value function)?:
elif isinstance(symbol, pybamm.BoundaryValue):
elif symbol.evaluates_on_edges("primary"):
if symbol.side == "left":
sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts))
elif symbol.side == "right":
sub_matrix = csr_matrix(([1], ([0], [prim_pts - 1])), shape=(1, prim_pts))
else:
raise ValueError("side must be 'left' or 'right'")
additive = pybamm.Scalar(0)
additive_multiplicative = pybamm.Scalar(1)
multiplicative = pybamm.Scalar(1)
@rtimms @martinjrobins , cool, I'll make the change on my fork and run some tests to see if 1) it works without errors and 2) conservation is achieved in the model @WillClarke25 and I are working on 🙂. Will let you know if I need any help.
@martinjrobins @rtimms I get a shape mismatch error, I think it is because when taking the edge evaluated boundary values, the matrix created follows the shape of the edges and not nodes (n+1), I am trying to manually trim the discretised_child which is used for the symbols depending on the leading side so something like:
elif symbol.evaluates_on_edges("primary"):
child_size = discretised_child.shape[0] - 1
print(discretised_child)
if symbol.side == "left":
discretised_child = discretised_child[:-1]
sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, child_size))
elif symbol.side == "right":
discretised_child = discretised_child[1:]
sub_matrix = csr_matrix(([1], ([0], [child_size - 1])), shape=(1, child_size))
else:
raise ValueError("side must be 'left' or 'right'")
additive = pybamm.Scalar(0)
additive_multiplicative = pybamm.Scalar(1)
multiplicative = pybamm.Scalar(1)
But this does not work as discretised_child is not an object that can be sliced. So how do I go about trimming it as the logic above?
@martinjrobins @rtimms I tried the below but it now gives an initial condition shape mismatch error:
elif symbol.evaluates_on_edges("primary"):
child_size = discretised_child.shape[0] # full size, no -1 here
if symbol.side == "left":
# Remove last element: keep first (N-1)
row_indices = range(child_size - 1)
col_indices = range(child_size - 1)
data = [1] * (child_size - 1)
S = csr_matrix((data, (row_indices, col_indices)), shape=(child_size - 1, child_size))
discretised_child = pybamm.Matrix(S) @ discretised_child
print(discretised_child)
# sub_matrix picks first node on sliced child (index 0)
sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, child_size - 1))
elif symbol.side == "right":
# Remove first element: keep last (N-1)
row_indices = range(child_size - 1)
col_indices = range(1, child_size)
data = [1] * (child_size - 1)
S = csr_matrix((data, (row_indices, col_indices)), shape=(child_size - 1, child_size))
discretised_child = pybamm.Matrix(S) @ discretised_child
print(discretised_child)
# sub_matrix picks last node on sliced child (index child_size - 2)
sub_matrix = csr_matrix(([1], ([0], [child_size - 2])), shape=(1, child_size - 1))
I am unsure as to how to proceed?
@rtimms @martinjrobins given that discretised_child is symbolic in nature, is there any methods defined for it to allow for slicing?