Skip to content

Commit

Permalink
Fix upwind/downwind finite volume operators (pybamm-team#3979)
Browse files Browse the repository at this point in the history 
* pybamm-team#3526 initialise branch

* pybamm-team#3256 add UpwindDivergence and DownwindDivergence

* work in progress for fixing upwind scheme

* pybamm-team#3526 add missing argument docstring

* pybamm-team#3526 fix upwind/downwind methods

* style: pre-commit fixes

* pybamm-team#3526 get finite volume notebook from develop

* 3526 add temporary test

* ruff

* style: pre-commit fixes

* pybamm-team#3526 remove unused variable

* pybamm-team#3526 clarify wording

* pybamm-team#3526 fix tests

* pybamm-team#3526 add to CHANGELOG

* pybamm-team#3526 update coverage

* pybamm-team#3526 add integration test

* style: pre-commit fixes

* pybamm-team#3526 remove example for debugging

* Remove print for debug

Co-authored-by: Eric G. Kratz <[email protected]>

* pybamm-team#3526 add second test

* pybamm-team#3526 generalised example but did not add test as was hard to know what to measure

* style: pre-commit fixes

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric G. Kratz <[email protected]>
Co-authored-by: Valentin Sulzer <[email protected]>
  • Loading branch information
@brosaplanella @pre-commit-ci
@kratman @valentinsulzer @js1tr3
4 people authored and js1tr3 committed Aug 12, 2024
1 parent 8cb019e commit e699d39
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 37 deletions.
4 changes: 4 additions & 0 deletions pybamm/expression_tree/unary_operators.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ class UnaryOperator(pybamm.Symbol):
name of the node
child : :class:`Symbol`
child node
domains : dict
A dictionary equivalent to {'primary': domain, auxiliary_domains}.
"""

def __init__(
Expand Down Expand Up @@ -398,6 +400,8 @@ class with a :class:`Matrix`
name of the node
child : :class:`Symbol`
child node
domains : dict
A dictionary equivalent to {'primary': domain, auxiliary_domains}.
"""

def __init__(
Expand Down
54 changes: 19 additions & 35 deletions pybamm/spatial_methods/finite_volume.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -613,8 +613,13 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
n = submesh.npts
second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)

lbc_value, lbc_type = bcs["left"]
rbc_value, rbc_type = bcs["right"]
# Catch if no boundary conditions are defined
if "left" not in bcs.keys() and "right" 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))

# Add ghost node(s) to domain where necessary and count number of
# Dirichlet boundary conditions
Expand All @@ -637,7 +642,7 @@ 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 == "Neumann":
elif lbc_type in ["Neumann", None]:
lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats))
else:
raise ValueError(
Expand All @@ -656,7 +661,7 @@ 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 == "Neumann":
elif rbc_type in ["Neumann", None]:
rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats))
else:
raise ValueError(
Expand Down Expand Up @@ -1390,45 +1395,24 @@ def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction):
direction : str
Direction in which to apply the operator (upwind or downwind)
"""
submesh = self.mesh[symbol.domain]
n = submesh.npts

if symbol not in bcs:
raise pybamm.ModelError(
"Boundary conditions must be provided for " f"{direction}ing '{symbol}'"
)

if direction == "upwind":
bc, typ = bcs[symbol]["left"]
if typ != "Dirichlet":
raise pybamm.ModelError(
"Dirichlet boundary conditions must be provided for "
f"upwinding '{symbol}'"
)

concat_bc = pybamm.NumpyConcatenation(bc, discretised_symbol)

upwind_mat = vstack(
[
csr_matrix(([1], ([0], [0])), shape=(1, n + 1)),
diags([-0.5, 1.5], [0, 1], shape=(n, n + 1)),
]
)
symbol_out = pybamm.Matrix(upwind_mat) @ concat_bc
bc_side = "left"
elif direction == "downwind":
bc, typ = bcs[symbol]["right"]
if typ != "Dirichlet":
raise pybamm.ModelError(
"Dirichlet boundary conditions must be provided for "
f"downwinding '{symbol}'"
)
bc_side = "right"

concat_bc = pybamm.NumpyConcatenation(discretised_symbol, bc)
downwind_mat = vstack(
[
diags([1.5, -0.5], [0, 1], shape=(n, n + 1)),
csr_matrix(([1], ([0], [n])), shape=(1, n + 1)),
]
if bcs[symbol][bc_side][1] != "Dirichlet":
raise pybamm.ModelError(
"Dirichlet boundary conditions must be provided for "
f"{direction}ing '{symbol}'"
)
symbol_out = pybamm.Matrix(downwind_mat) @ concat_bc

# 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
49 changes: 49 additions & 0 deletions tests/integration/test_spatial_methods/test_finite_volume.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -339,6 +339,55 @@ def test_laplacian_spherical(self):
)


def solve_advection_equation(direction="upwind", source=1, bc=0):
model = pybamm.BaseModel()
x = pybamm.SpatialVariable("x", domain="domain", coord_sys="cartesian")
u = pybamm.Variable("u", domain="domain")
if direction == "upwind":
bc_side = "left"
y = x
v = pybamm.PrimaryBroadcastToEdges(1, ["domain"])
rhs = -pybamm.div(pybamm.upwind(u) * v) + source
elif direction == "downwind":
bc_side = "right"
y = 1 - x
v = pybamm.PrimaryBroadcastToEdges(-1, ["domain"])
rhs = -pybamm.div(pybamm.downwind(u) * v) + source

u_an = (bc + source * y) - (bc + source * (y - pybamm.t)) * ((y - pybamm.t) > 0)
model.boundary_conditions = {
u: {
bc_side: (pybamm.Scalar(bc), "Dirichlet"),
}
}
model.rhs = {u: rhs}
model.initial_conditions = {u: pybamm.Scalar(0)}
model.variables = {"u": u, "x": x, "analytical": u_an}
geometry = {"domain": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
submesh_types = {"domain": pybamm.Uniform1DSubMesh}
var_pts = {x: 1000}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"domain": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)
solver = pybamm.CasadiSolver()
return solver.solve(model, [0, 1])


class TestUpwindDownwind(TestCase):
def test_upwind(self):
solution = solve_advection_equation("upwind")
np.testing.assert_array_almost_equal(
solution["u"].entries, solution["analytical"].entries, decimal=2
)

def test_downwind(self):
solution = solve_advection_equation("downwind")
np.testing.assert_array_almost_equal(
solution["u"].entries, solution["analytical"].entries, decimal=2
)


if __name__ == "__main__":
print("Add -v for more debug output")
import sys
Expand Down
4 changes: 2 additions & 2 deletions tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -410,11 +410,11 @@ def test_upwind_downwind(self):
y_test = 2 * np.ones_like(nodes)
np.testing.assert_array_equal(
disc_upwind.evaluate(y=y_test),
np.concatenate([np.array([5, 0.5]), 2 * np.ones(n - 1)])[:, np.newaxis],
np.concatenate([np.array([8]), 2 * np.ones(n)])[:, np.newaxis],
)
np.testing.assert_array_equal(
disc_downwind.evaluate(y=y_test),
np.concatenate([2 * np.ones(n - 1), np.array([1.5, 3])])[:, np.newaxis],
np.concatenate([2 * np.ones(n), np.array([4])])[:, np.newaxis],
)

# Remove boundary conditions and check error is raised
Expand Down
2 changes: 2 additions & 0 deletions tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,8 @@ def test_add_ghost_nodes(self):
bcs = {"left": (pybamm.Scalar(0), "Neumann"), "right": (pybamm.Scalar(3), "x")}
with self.assertRaisesRegex(ValueError, "boundary condition must be"):
sp_meth.add_ghost_nodes(var, discretised_symbol, bcs)
with self.assertRaisesRegex(ValueError, "No boundary conditions"):
sp_meth.add_ghost_nodes(var, discretised_symbol, {})
with self.assertRaisesRegex(ValueError, "boundary condition must be"):
sp_meth.add_neumann_values(var, discretised_symbol, bcs, var.domain)

Expand Down

0 comments on commit e699d39

Please sign in to comment.