From 5a9456e0112602949dd324d5644b900f9e703f5c Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Fri, 31 Jul 2026 14:39:04 -0700 Subject: [PATCH 1/2] feat: add unstructured 2D/3D DFN battery models Add BasicDFN2DUnstructured and BasicDFN3DUnstructured on top of the meshing, spatial-method, and plotting stack, plus Max/Min reduction and domain-size helpers needed by those models. Co-authored-by: Cursor --- CHANGELOG.md | 1 + docs/source/api/models/lithium_ion/dfn.rst | 6 + .../src/pybamm/expression_tree/functions.py | 32 +- .../src/pybamm/expression_tree/symbol.py | 8 + .../lithium_ion/__init__.py | 2 + .../lithium_ion/basic_dfn_2d_unstructured.py | 419 ++++++++++++++++ .../lithium_ion/basic_dfn_3d_unstructured.py | 450 ++++++++++++++++++ .../test_lithium_ion/test_basic_models.py | 47 ++ packages/pybamm/tests/strategies/symbols.py | 1 + .../unit/test_expression_tree/test_symbol.py | 8 + .../test_lithium_ion/test_basic_models.py | 8 + 11 files changed, 978 insertions(+), 4 deletions(-) create mode 100644 packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py create mode 100644 packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_3d_unstructured.py diff --git a/CHANGELOG.md b/CHANGELOG.md index fbf31ef93b..73af719825 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ ## Features +- Added basic unstructured 2D/3D DFN models (`BasicDFN2DUnstructured`, `BasicDFN3DUnstructured`). ([#5690](https://github.com/pybamm-team/PyBaMM/pull/5690)) - Added VTK-based plotting (`VTKQuickPlot`) for unstructured mesh solutions, including headless CI OpenGL setup. ([#5689](https://github.com/pybamm-team/PyBaMM/pull/5689)) - Added `FiniteVolumeUnstructured` spatial method and unstructured processed-variable support for cell-centered data on arbitrary meshes. The TPFA Laplacian carries an implicit non-orthogonal correction (`"non-orthogonal correction"` option: `"over-relaxed"` or `"minimum"`) and gradients use a least-squares reconstruction, so both are exact on linear fields and second-order on skewed triangle and tetrahedral meshes. Diffusion coefficients reach faces through the distance-weighted harmonic mean, as in `FiniteVolume`, so material interfaces carry the exact series flux. ([#5688](https://github.com/pybamm-team/PyBaMM/pull/5688)) - Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. Hexahedra must have planar faces (warped hexes raise a `GeometryError`), and `UserSuppliedUnstructuredMesh` accepts tetrahedral, triangular, and quadrilateral cells only. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687)) diff --git a/docs/source/api/models/lithium_ion/dfn.rst b/docs/source/api/models/lithium_ion/dfn.rst index 4213f4ccf9..6231e88c55 100644 --- a/docs/source/api/models/lithium_ion/dfn.rst +++ b/docs/source/api/models/lithium_ion/dfn.rst @@ -13,4 +13,10 @@ Doyle-Fuller-Newman (DFN) .. autoclass:: pybamm.lithium_ion.BasicDFNHalfCell :members: +.. autoclass:: pybamm.lithium_ion.BasicDFN2DUnstructured + :members: + +.. autoclass:: pybamm.lithium_ion.BasicDFN3DUnstructured + :members: + .. footbibliography:: diff --git a/packages/pybamm/src/pybamm/expression_tree/functions.py b/packages/pybamm/src/pybamm/expression_tree/functions.py index dca18046be..794ca7af38 100644 --- a/packages/pybamm/src/pybamm/expression_tree/functions.py +++ b/packages/pybamm/src/pybamm/expression_tree/functions.py @@ -719,8 +719,32 @@ def log10(child: pybamm.Symbol): return log(child, base=10) -class Max(SpecificFunction): - """Max function.""" +class Reduction(SpecificFunction): + """Base class for reduction operations that collapse a spatial + field to a scalar (e.g. max, min). Automatically clears domains + and returns scalar shape.""" + + def __init__(self, function: Callable, child: pybamm.Symbol): + super().__init__(function, child) + self.clear_domains() + + @classmethod + def _from_json(cls, snippet: dict): + """See :meth:`pybamm.SpecificFunction._from_json()`. + + ``SpecificFunction._from_json`` bypasses ``__init__``, so the domains + inherited from the child have to be cleared again here. + """ + instance = super()._from_json(snippet) + instance.clear_domains() + return instance + + def _evaluate_for_shape(self): + return np.nan * np.ones((1, 1)) + + +class Max(Reduction): + """Max function (reduction to scalar).""" def __init__(self, child): super().__init__(np.max, child) @@ -750,8 +774,8 @@ def max(child: pybamm.Symbol): return pybamm.simplify_if_constant(Max(child)) -class Min(SpecificFunction): - """Min function.""" +class Min(Reduction): + """Min function (reduction to scalar).""" def __init__(self, child): super().__init__(np.min, child) diff --git a/packages/pybamm/src/pybamm/expression_tree/symbol.py b/packages/pybamm/src/pybamm/expression_tree/symbol.py index f2cff8a1c4..733be05796 100644 --- a/packages/pybamm/src/pybamm/expression_tree/symbol.py +++ b/packages/pybamm/src/pybamm/expression_tree/symbol.py @@ -45,12 +45,20 @@ def domain_size(domain: list[str] | str): fixed_domain_sizes = { "current collector": 3, "negative particle": 5, + "negative primary particle": 5, + "negative secondary particle": 5, "positive particle": 7, + "positive primary particle": 7, + "positive secondary particle": 7, "negative electrode": 11, "separator": 13, "positive electrode": 17, "negative particle size": 19, + "negative primary particle size": 19, + "negative secondary particle size": 19, "positive particle size": 23, + "positive primary particle size": 23, + "positive secondary particle size": 23, } if domain in [[], None]: size = 1 diff --git a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/__init__.py b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/__init__.py index 79ce665a55..2c31fcb2e6 100644 --- a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/__init__.py +++ b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/__init__.py @@ -24,6 +24,8 @@ from .newman_tobias import NewmanTobias from .basic_dfn import BasicDFN from .basic_dfn_2d import BasicDFN2D +from .basic_dfn_2d_unstructured import BasicDFN2DUnstructured +from .basic_dfn_3d_unstructured import BasicDFN3DUnstructured from .basic_spm import BasicSPM from .basic_spm_with_3d_thermal import Basic3DThermalSPM from .basic_dfn_half_cell import BasicDFNHalfCell diff --git a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py new file mode 100644 index 0000000000..7d8af2378f --- /dev/null +++ b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py @@ -0,0 +1,419 @@ +# +# Basic Doyle-Fuller-Newman (DFN) Model — 2D Unstructured FVM +# +import pybamm +from pybamm.models.full_battery_models.lithium_ion.base_lithium_ion_model import ( + BaseModel, +) + + +class BasicDFN2DUnstructured(BaseModel): + """Doyle-Fuller-Newman (DFN) model on a 2D unstructured mesh. + + Identical physics to :class:`BasicDFN2D` but uses + :class:`~pybamm.FiniteVolumeUnstructured` on triangle or quad elements + instead of the structured tensor-product grid. + + Parameters + ---------- + name : str, optional + The name of the model. + element_type : str, optional + Element type for the built-in mesh generator: ``"quad"`` (default, + TPFA-orthogonal) or ``"triangle"``. + """ + + def __init__( + self, + name="Doyle-Fuller-Newman model (2D unstructured)", + element_type="quad", + ): + super().__init__(name=name) + self._element_type = element_type + pybamm.citations.register("Marquis2019") + + ###################### + # Variables + ###################### + Q = pybamm.Variable("Discharge capacity [A.h]") + + x = pybamm.SpatialVariable( + "x", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="lr", + ) + x_n = pybamm.SpatialVariable( + "x_n", domain="negative electrode", coord_sys="cartesian", direction="lr" + ) + x_s = pybamm.SpatialVariable( + "x_s", domain="separator", coord_sys="cartesian", direction="lr" + ) + x_p = pybamm.SpatialVariable( + "x_p", domain="positive electrode", coord_sys="cartesian", direction="lr" + ) + z_n = pybamm.SpatialVariable( + "z_n", domain="negative electrode", coord_sys="cartesian", direction="tb" + ) + z_s = pybamm.SpatialVariable( + "z_s", domain="separator", coord_sys="cartesian", direction="tb" + ) + z_p = pybamm.SpatialVariable( + "z_p", domain="positive electrode", coord_sys="cartesian", direction="tb" + ) + z = pybamm.SpatialVariable( + "z", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + + c_e_n = pybamm.Variable( + "Negative electrolyte concentration [mol.m-3]", + domain="negative electrode", + ) + c_e_s = pybamm.Variable( + "Separator electrolyte concentration [mol.m-3]", + domain="separator", + ) + c_e_p = pybamm.Variable( + "Positive electrolyte concentration [mol.m-3]", + domain="positive electrode", + ) + c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p) + + phi_e_n = pybamm.Variable( + "Negative electrolyte potential [V]", + domain="negative electrode", + ) + phi_e_s = pybamm.Variable( + "Separator electrolyte potential [V]", + domain="separator", + ) + phi_e_p = pybamm.Variable( + "Positive electrolyte potential [V]", + domain="positive electrode", + ) + phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p) + + phi_s_n = pybamm.Variable( + "Negative electrode potential [V]", domain="negative electrode" + ) + phi_s_p = pybamm.Variable( + "Positive electrode potential [V]", + domain="positive electrode", + ) + c_s_n = pybamm.Variable( + "Negative particle concentration [mol.m-3]", + domain="negative particle", + auxiliary_domains={"secondary": "negative electrode"}, + ) + c_s_p = pybamm.Variable( + "Positive particle concentration [mol.m-3]", + domain="positive particle", + auxiliary_domains={"secondary": "positive electrode"}, + ) + + T = self.param.T_init + + ###################### + # Other set-up + ###################### + i_cell = self.param.current_density_with_time + + eps_n = pybamm.FunctionParameter( + "Negative electrode porosity", + {"Through-cell distance (x) [m]": x_n, "Vertical distance (z) [m]": z_n}, + ) + eps_s = pybamm.FunctionParameter( + "Separator porosity", + {"Through-cell distance (x) [m]": x_s, "Vertical distance (z) [m]": z_s}, + ) + eps_p = pybamm.FunctionParameter( + "Positive electrode porosity", + {"Through-cell distance (x) [m]": x_p, "Vertical distance (z) [m]": z_p}, + ) + eps = pybamm.concatenation(eps_n, eps_s, eps_p) + + eps_s_n = pybamm.FunctionParameter( + "Negative electrode active material volume fraction", + {"Through-cell distance (x) [m]": x_n, "Vertical distance (z) [m]": z_n}, + ) + eps_s_p = pybamm.FunctionParameter( + "Positive electrode active material volume fraction", + {"Through-cell distance (x) [m]": x_p, "Vertical distance (z) [m]": z_p}, + ) + + tor = pybamm.concatenation( + eps_n**self.param.n.b_e, eps_s**self.param.s.b_e, eps_p**self.param.p.b_e + ) + a_n = 3 * self.param.n.prim.epsilon_s_av / self.param.n.prim.R_typ + a_p = 3 * self.param.p.prim.epsilon_s_av / self.param.p.prim.R_typ + + # Interfacial reactions + c_s_surf_n = pybamm.surf(c_s_n) + sto_surf_n = c_s_surf_n / self.param.n.prim.c_max + j0_n = self.param.n.prim.j0(c_e_n, c_s_surf_n, T) + delta_phi_n = phi_s_n - phi_e_n + eta_n = delta_phi_n - self.param.n.prim.U(sto_surf_n, T) + Feta_RT_n = self.param.F * eta_n / (self.param.R * T) + j_n = 2 * j0_n * pybamm.sinh(self.param.n.prim.ne / 2 * Feta_RT_n) + + c_s_surf_p = pybamm.surf(c_s_p) + sto_surf_p = c_s_surf_p / self.param.p.prim.c_max + j0_p = self.param.p.prim.j0(c_e_p, c_s_surf_p, T) + delta_phi_p = phi_s_p - phi_e_p + eta_p = delta_phi_p - self.param.p.prim.U(sto_surf_p, T) + Feta_RT_p = self.param.F * eta_p / (self.param.R * T) + j_s = pybamm.PrimaryBroadcast(0, "separator") + j_p = 2 * j0_p * pybamm.sinh(self.param.p.prim.ne / 2 * Feta_RT_p) + + a_j_n = a_n * j_n + a_j_p = a_p * j_p + a_j = pybamm.concatenation(a_j_n, j_s, a_j_p) + + ###################### + # State of Charge + ###################### + current = self.param.current_with_time + self.rhs[Q] = current / 3600 + self.initial_conditions[Q] = pybamm.Scalar(0) + + ###################### + # Particles + ###################### + N_s_n = -self.param.n.prim.D(c_s_n, T) * pybamm.grad(c_s_n) + N_s_p = -self.param.p.prim.D(c_s_p, T) * pybamm.grad(c_s_p) + self.rhs[c_s_n] = -pybamm.div(N_s_n) + self.rhs[c_s_p] = -pybamm.div(N_s_p) + self.boundary_conditions[c_s_n] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": ( + -j_n / (self.param.F * pybamm.surf(self.param.n.prim.D(c_s_n, T))), + "Neumann", + ), + } + self.boundary_conditions[c_s_p] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": ( + -j_p / (self.param.F * pybamm.surf(self.param.p.prim.D(c_s_p, T))), + "Neumann", + ), + } + self.initial_conditions[c_s_n] = self.param.n.prim.c_init + self.initial_conditions[c_s_p] = self.param.p.prim.c_init + + c_s_n_av = pybamm.RAverage(c_s_n) + c_s_p_av = pybamm.RAverage(c_s_p) + solid_lithium_negative = pybamm.Integral(c_s_n_av * eps_s_n, [x_n, z_n]) + solid_lithium_positive = pybamm.Integral(c_s_p_av * eps_s_p, [x_p, z_p]) + total_solid_lithium = solid_lithium_negative + solid_lithium_positive + + ###################### + # Current in the solid + ###################### + sigma_eff_n = self.param.n.sigma(None, T) * eps_s_n**self.param.n.b_s + sigma_eff_p = self.param.p.sigma(None, T) * eps_s_p**self.param.p.b_s + self.algebraic[phi_s_n] = ( + self.param.L_x**2 + * self.param.L_z**2 + * (pybamm.div(-sigma_eff_n * pybamm.grad(phi_s_n)) + a_j_n) + ) + self.algebraic[phi_s_p] = ( + self.param.L_x**2 + * self.param.L_z**2 + * (pybamm.div(-sigma_eff_p * pybamm.grad(phi_s_p)) + a_j_p) + ) + self.boundary_conditions[phi_s_n] = { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + self.boundary_conditions[phi_s_p] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[phi_s_n] = pybamm.Scalar(0) + self.initial_conditions[phi_s_p] = self.param.ocv_init + ###################### + # Current in the electrolyte + ###################### + kappa_eff = self.param.kappa_e(c_e, T) * tor + kappa_D_eff = kappa_eff * self.param.chiRT_over_Fc(c_e, T) + self.algebraic[phi_e] = ( + self.param.L_x**2 + * self.param.L_z**2 + * ( + pybamm.div(kappa_D_eff * pybamm.grad(c_e)) + - pybamm.div(kappa_eff * pybamm.grad(phi_e)) + - a_j + ) + ) + self.boundary_conditions[phi_e] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[phi_e] = -self.param.n.prim.U_init + + ###################### + # Electrolyte concentration + ###################### + D_e_eff = tor * self.param.D_e(c_e, T) + self.rhs[c_e] = (1 / eps) * ( + pybamm.div(D_e_eff * pybamm.grad(c_e)) + + (1 - self.param.t_plus(c_e, T)) * a_j / self.param.F + ) + self.boundary_conditions[c_e] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[c_e] = self.param.c_e_init + + ###################### + # (Some) variables + ###################### + voltage = pybamm.boundary_value(phi_s_p, "top-right") + num_cells = pybamm.Parameter( + "Number of cells connected in series to make a battery" + ) + total_lithium = pybamm.Integral(c_e * eps, [x, z]) + self.variables = { + "Negative particle concentration [mol.m-3]": c_s_n, + "Total lithium [mol]": total_lithium, + "Negative particle surface concentration [mol.m-3]": c_s_surf_n, + "Electrolyte concentration [mol.m-3]": c_e, + "Negative electrolyte concentration [mol.m-3]": c_e_n, + "Separator electrolyte concentration [mol.m-3]": c_e_s, + "Positive electrolyte concentration [mol.m-3]": c_e_p, + "Positive particle concentration [mol.m-3]": c_s_p, + "Positive particle surface concentration [mol.m-3]": c_s_surf_p, + "Current [A]": current, + "Current variable [A]": current, + "Negative electrode potential [V]": phi_s_n, + "Electrolyte potential [V]": phi_e, + "Negative electrolyte potential [V]": phi_e_n, + "Separator electrolyte potential [V]": phi_e_s, + "Positive electrolyte potential [V]": phi_e_p, + "Positive electrode potential [V]": phi_s_p, + "Voltage [V]": voltage, + "Battery voltage [V]": voltage * num_cells, + "Time [s]": pybamm.t, + "Discharge capacity [A.h]": Q, + "x": x, + "z": z, + "Current density [A.m-2]": a_j, + "Electrolyte current density [A.m-2]": a_j, + "x_n": x_n, + "x_s": x_s, + "x_p": x_p, + "z_n": z_n, + "z_s": z_s, + "z_p": z_p, + "Negative electrode surface concentration [mol.m-3]": c_s_surf_n, + "Negative electrode surface stoichiometry": sto_surf_n, + "Positive electrode surface concentration [mol.m-3]": c_s_surf_p, + "Positive electrode surface stoichiometry": sto_surf_p, + "Positive electrode surface potential difference [V]": delta_phi_p, + "Negative electrode surface potential difference [V]": delta_phi_n, + "Positive electrode overpotential [V]": eta_p, + "Negative electrode overpotential [V]": eta_n, + "Positive electrode ocp [V]": self.param.p.prim.U(sto_surf_p, T), + "Negative electrode ocp [V]": self.param.n.prim.U(sto_surf_n, T), + "Positive electrode current density [A.m-2]": j_p, + "Negative electrode current density [A.m-2]": j_n, + "Electrolyte flux [mol.m-2.s-1]": D_e_eff, + "Positive solid lithium [mol]": solid_lithium_positive, + "Negative solid lithium [mol]": solid_lithium_negative, + "Total solid lithium [mol]": total_solid_lithium, + } + self.events += [ + pybamm.Event("Minimum voltage [V]", voltage - self.param.voltage_low_cut), + pybamm.Event("Maximum voltage [V]", self.param.voltage_high_cut - voltage), + ] + + @property + def default_geometry(self): + z_2d = pybamm.SpatialVariable( + "z_2d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + return { + "negative electrode": { + "x_n": {"min": 0, "max": self.param.n.L}, + z_2d: {"min": 0, "max": self.param.L_z}, + }, + "separator": { + "x_s": {"min": self.param.n.L, "max": self.param.n.L + self.param.s.L}, + z_2d: {"min": 0, "max": self.param.L_z}, + }, + "positive electrode": { + "x_p": { + "min": self.param.n.L + self.param.s.L, + "max": self.param.n.L + self.param.s.L + self.param.p.L, + }, + z_2d: {"min": 0, "max": self.param.L_z}, + }, + "positive particle": { + "r_p": {"min": 0, "max": self.param.p.prim.R_typ}, + }, + "negative particle": { + "r_n": {"min": 0, "max": self.param.n.prim.R_typ}, + }, + "current collector": { + "z": {"position": 0}, + }, + } + + @property + def default_spatial_methods(self): + return { + "negative electrode": pybamm.FiniteVolumeUnstructured(), + "separator": pybamm.FiniteVolumeUnstructured(), + "positive electrode": pybamm.FiniteVolumeUnstructured(), + "positive particle": pybamm.FiniteVolume(), + "negative particle": pybamm.FiniteVolume(), + "current collector": pybamm.ZeroDimensionalSpatialMethod(), + } + + @property + def default_submesh_types(self): + return { + "negative electrode": pybamm.UnstructuredMeshGenerator( + element_type=self._element_type + ), + "separator": pybamm.UnstructuredMeshGenerator( + element_type=self._element_type + ), + "positive electrode": pybamm.UnstructuredMeshGenerator( + element_type=self._element_type + ), + "positive particle": pybamm.Uniform1DSubMesh, + "negative particle": pybamm.Uniform1DSubMesh, + "current collector": pybamm.SubMesh0D, + } + + @property + def default_var_pts(self): + z_2d = pybamm.SpatialVariable( + "z_2d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + return { + "x_n": 20, + "x_s": 30, + "x_p": 20, + "r_p": 20, + "r_n": 20, + z_2d: 10, + } diff --git a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_3d_unstructured.py b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_3d_unstructured.py new file mode 100644 index 0000000000..18f04a733e --- /dev/null +++ b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_3d_unstructured.py @@ -0,0 +1,450 @@ +# +# Basic Doyle-Fuller-Newman (DFN) Model — 3D Unstructured FVM +# +import pybamm +from pybamm.models.full_battery_models.lithium_ion.base_lithium_ion_model import ( + BaseModel, +) + + +class BasicDFN3DUnstructured(BaseModel): + """Doyle-Fuller-Newman (DFN) model on a 3D unstructured mesh. + + Extends :class:`BasicDFN2DUnstructured` to three spatial dimensions + (x, y, z) using tetrahedral elements. The through-cell direction is + *x*, the width direction is *y*, and the height direction is *z*. + + Parameters + ---------- + name : str, optional + The name of the model. + """ + + def __init__( + self, + name="Doyle-Fuller-Newman model (3D unstructured)", + ): + super().__init__(name=name) + pybamm.citations.register("Marquis2019") + + ###################### + # Variables + ###################### + Q = pybamm.Variable("Discharge capacity [A.h]") + + x = pybamm.SpatialVariable( + "x", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="lr", + ) + x_n = pybamm.SpatialVariable( + "x_n", domain="negative electrode", coord_sys="cartesian", direction="lr" + ) + x_s = pybamm.SpatialVariable( + "x_s", domain="separator", coord_sys="cartesian", direction="lr" + ) + x_p = pybamm.SpatialVariable( + "x_p", domain="positive electrode", coord_sys="cartesian", direction="lr" + ) + y_n = pybamm.SpatialVariable( + "y_n", domain="negative electrode", coord_sys="cartesian", direction="fb" + ) + y_s = pybamm.SpatialVariable( + "y_s", domain="separator", coord_sys="cartesian", direction="fb" + ) + y_p = pybamm.SpatialVariable( + "y_p", domain="positive electrode", coord_sys="cartesian", direction="fb" + ) + y = pybamm.SpatialVariable( + "y", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="fb", + ) + z_n = pybamm.SpatialVariable( + "z_n", domain="negative electrode", coord_sys="cartesian", direction="tb" + ) + z_s = pybamm.SpatialVariable( + "z_s", domain="separator", coord_sys="cartesian", direction="tb" + ) + z_p = pybamm.SpatialVariable( + "z_p", domain="positive electrode", coord_sys="cartesian", direction="tb" + ) + z = pybamm.SpatialVariable( + "z", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + + c_e_n = pybamm.Variable( + "Negative electrolyte concentration [mol.m-3]", + domain="negative electrode", + ) + c_e_s = pybamm.Variable( + "Separator electrolyte concentration [mol.m-3]", + domain="separator", + ) + c_e_p = pybamm.Variable( + "Positive electrolyte concentration [mol.m-3]", + domain="positive electrode", + ) + c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p) + + phi_e_n = pybamm.Variable( + "Negative electrolyte potential [V]", + domain="negative electrode", + ) + phi_e_s = pybamm.Variable( + "Separator electrolyte potential [V]", + domain="separator", + ) + phi_e_p = pybamm.Variable( + "Positive electrolyte potential [V]", + domain="positive electrode", + ) + phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p) + + phi_s_n = pybamm.Variable( + "Negative electrode potential [V]", domain="negative electrode" + ) + phi_s_p = pybamm.Variable( + "Positive electrode potential [V]", + domain="positive electrode", + ) + c_s_n = pybamm.Variable( + "Negative particle concentration [mol.m-3]", + domain="negative particle", + auxiliary_domains={"secondary": "negative electrode"}, + ) + c_s_p = pybamm.Variable( + "Positive particle concentration [mol.m-3]", + domain="positive particle", + auxiliary_domains={"secondary": "positive electrode"}, + ) + + T = self.param.T_init + + ###################### + # Other set-up + ###################### + i_cell = self.param.current_density_with_time + + eps_n = pybamm.FunctionParameter( + "Negative electrode porosity", + {"Through-cell distance (x) [m]": x_n, "Vertical distance (z) [m]": z_n}, + ) + eps_s = pybamm.FunctionParameter( + "Separator porosity", + {"Through-cell distance (x) [m]": x_s, "Vertical distance (z) [m]": z_s}, + ) + eps_p = pybamm.FunctionParameter( + "Positive electrode porosity", + {"Through-cell distance (x) [m]": x_p, "Vertical distance (z) [m]": z_p}, + ) + eps = pybamm.concatenation(eps_n, eps_s, eps_p) + + eps_s_n = pybamm.FunctionParameter( + "Negative electrode active material volume fraction", + {"Through-cell distance (x) [m]": x_n, "Vertical distance (z) [m]": z_n}, + ) + eps_s_p = pybamm.FunctionParameter( + "Positive electrode active material volume fraction", + {"Through-cell distance (x) [m]": x_p, "Vertical distance (z) [m]": z_p}, + ) + + tor = pybamm.concatenation( + eps_n**self.param.n.b_e, eps_s**self.param.s.b_e, eps_p**self.param.p.b_e + ) + a_n = 3 * self.param.n.prim.epsilon_s_av / self.param.n.prim.R_typ + a_p = 3 * self.param.p.prim.epsilon_s_av / self.param.p.prim.R_typ + + # Interfacial reactions + c_s_surf_n = pybamm.surf(c_s_n) + sto_surf_n = c_s_surf_n / self.param.n.prim.c_max + j0_n = self.param.n.prim.j0(c_e_n, c_s_surf_n, T) + delta_phi_n = phi_s_n - phi_e_n + eta_n = delta_phi_n - self.param.n.prim.U(sto_surf_n, T) + Feta_RT_n = self.param.F * eta_n / (self.param.R * T) + j_n = 2 * j0_n * pybamm.sinh(self.param.n.prim.ne / 2 * Feta_RT_n) + + c_s_surf_p = pybamm.surf(c_s_p) + sto_surf_p = c_s_surf_p / self.param.p.prim.c_max + j0_p = self.param.p.prim.j0(c_e_p, c_s_surf_p, T) + delta_phi_p = phi_s_p - phi_e_p + eta_p = delta_phi_p - self.param.p.prim.U(sto_surf_p, T) + Feta_RT_p = self.param.F * eta_p / (self.param.R * T) + j_s = pybamm.PrimaryBroadcast(0, "separator") + j_p = 2 * j0_p * pybamm.sinh(self.param.p.prim.ne / 2 * Feta_RT_p) + + a_j_n = a_n * j_n + a_j_p = a_p * j_p + a_j = pybamm.concatenation(a_j_n, j_s, a_j_p) + + ###################### + # State of Charge + ###################### + current = self.param.current_with_time + self.rhs[Q] = current / 3600 + self.initial_conditions[Q] = pybamm.Scalar(0) + + ###################### + # Particles + ###################### + N_s_n = -self.param.n.prim.D(c_s_n, T) * pybamm.grad(c_s_n) + N_s_p = -self.param.p.prim.D(c_s_p, T) * pybamm.grad(c_s_p) + self.rhs[c_s_n] = -pybamm.div(N_s_n) + self.rhs[c_s_p] = -pybamm.div(N_s_p) + self.boundary_conditions[c_s_n] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": ( + -j_n / (self.param.F * pybamm.surf(self.param.n.prim.D(c_s_n, T))), + "Neumann", + ), + } + self.boundary_conditions[c_s_p] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": ( + -j_p / (self.param.F * pybamm.surf(self.param.p.prim.D(c_s_p, T))), + "Neumann", + ), + } + self.initial_conditions[c_s_n] = self.param.n.prim.c_init + self.initial_conditions[c_s_p] = self.param.p.prim.c_init + + c_s_n_av = pybamm.RAverage(c_s_n) + c_s_p_av = pybamm.RAverage(c_s_p) + solid_lithium_negative = pybamm.Integral(c_s_n_av * eps_s_n, [x_n, y_n, z_n]) + solid_lithium_positive = pybamm.Integral(c_s_p_av * eps_s_p, [x_p, y_p, z_p]) + total_solid_lithium = solid_lithium_negative + solid_lithium_positive + + ###################### + # Current in the solid + ###################### + sigma_eff_n = self.param.n.sigma(None, T) * eps_s_n**self.param.n.b_s + sigma_eff_p = self.param.p.sigma(None, T) * eps_s_p**self.param.p.b_s + L_scale = self.param.L_x**2 * self.param.L_z**2 + self.algebraic[phi_s_n] = L_scale * ( + pybamm.div(-sigma_eff_n * pybamm.grad(phi_s_n)) + a_j_n + ) + self.algebraic[phi_s_p] = L_scale * ( + pybamm.div(-sigma_eff_p * pybamm.grad(phi_s_p)) + a_j_p + ) + self.boundary_conditions[phi_s_n] = { + "left": (pybamm.Scalar(0), "Dirichlet"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + "front": (pybamm.Scalar(0), "Neumann"), + "back": (pybamm.Scalar(0), "Neumann"), + } + self.boundary_conditions[phi_s_p] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": ( + i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), + "Neumann", + ), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + "front": (pybamm.Scalar(0), "Neumann"), + "back": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[phi_s_n] = pybamm.Scalar(0) + self.initial_conditions[phi_s_p] = self.param.ocv_init + ###################### + # Current in the electrolyte + ###################### + kappa_eff = self.param.kappa_e(c_e, T) * tor + kappa_D_eff = kappa_eff * self.param.chiRT_over_Fc(c_e, T) + self.algebraic[phi_e] = L_scale * ( + pybamm.div(kappa_D_eff * pybamm.grad(c_e)) + - pybamm.div(kappa_eff * pybamm.grad(phi_e)) + - a_j + ) + self.boundary_conditions[phi_e] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + "front": (pybamm.Scalar(0), "Neumann"), + "back": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[phi_e] = -self.param.n.prim.U_init + + ###################### + # Electrolyte concentration + ###################### + D_e_eff = tor * self.param.D_e(c_e, T) + self.rhs[c_e] = (1 / eps) * ( + pybamm.div(D_e_eff * pybamm.grad(c_e)) + + (1 - self.param.t_plus(c_e, T)) * a_j / self.param.F + ) + self.boundary_conditions[c_e] = { + "left": (pybamm.Scalar(0), "Neumann"), + "right": (pybamm.Scalar(0), "Neumann"), + "top": (pybamm.Scalar(0), "Neumann"), + "bottom": (pybamm.Scalar(0), "Neumann"), + "front": (pybamm.Scalar(0), "Neumann"), + "back": (pybamm.Scalar(0), "Neumann"), + } + self.initial_conditions[c_e] = self.param.c_e_init + + ###################### + # (Some) variables + ###################### + voltage = pybamm.boundary_value(phi_s_p, "top-right") + num_cells = pybamm.Parameter( + "Number of cells connected in series to make a battery" + ) + total_lithium = pybamm.Integral(c_e * eps, [x, y, z]) + self.variables = { + "Negative particle concentration [mol.m-3]": c_s_n, + "Total lithium [mol]": total_lithium, + "Negative particle surface concentration [mol.m-3]": c_s_surf_n, + "Electrolyte concentration [mol.m-3]": c_e, + "Negative electrolyte concentration [mol.m-3]": c_e_n, + "Separator electrolyte concentration [mol.m-3]": c_e_s, + "Positive electrolyte concentration [mol.m-3]": c_e_p, + "Positive particle concentration [mol.m-3]": c_s_p, + "Positive particle surface concentration [mol.m-3]": c_s_surf_p, + "Current [A]": current, + "Current variable [A]": current, + "Negative electrode potential [V]": phi_s_n, + "Electrolyte potential [V]": phi_e, + "Negative electrolyte potential [V]": phi_e_n, + "Separator electrolyte potential [V]": phi_e_s, + "Positive electrolyte potential [V]": phi_e_p, + "Positive electrode potential [V]": phi_s_p, + "Voltage [V]": voltage, + "Battery voltage [V]": voltage * num_cells, + "Time [s]": pybamm.t, + "Discharge capacity [A.h]": Q, + "x": x, + "y": y, + "z": z, + "Current density [A.m-2]": a_j, + "Electrolyte current density [A.m-2]": a_j, + "x_n": x_n, + "x_s": x_s, + "x_p": x_p, + "y_n": y_n, + "y_s": y_s, + "y_p": y_p, + "z_n": z_n, + "z_s": z_s, + "z_p": z_p, + "Negative electrode surface concentration [mol.m-3]": c_s_surf_n, + "Negative electrode surface stoichiometry": sto_surf_n, + "Positive electrode surface concentration [mol.m-3]": c_s_surf_p, + "Positive electrode surface stoichiometry": sto_surf_p, + "Positive electrode surface potential difference [V]": delta_phi_p, + "Negative electrode surface potential difference [V]": delta_phi_n, + "Positive electrode overpotential [V]": eta_p, + "Negative electrode overpotential [V]": eta_n, + "Positive electrode ocp [V]": self.param.p.prim.U(sto_surf_p, T), + "Negative electrode ocp [V]": self.param.n.prim.U(sto_surf_n, T), + "Positive electrode current density [A.m-2]": j_p, + "Negative electrode current density [A.m-2]": j_n, + "Electrolyte flux [mol.m-2.s-1]": D_e_eff, + "Positive solid lithium [mol]": solid_lithium_positive, + "Negative solid lithium [mol]": solid_lithium_negative, + "Total solid lithium [mol]": total_solid_lithium, + } + self.events += [ + pybamm.Event("Minimum voltage [V]", voltage - self.param.voltage_low_cut), + pybamm.Event("Maximum voltage [V]", self.param.voltage_high_cut - voltage), + ] + + @property + def default_geometry(self): + y_3d = pybamm.SpatialVariable( + "y_3d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="fb", + ) + z_3d = pybamm.SpatialVariable( + "z_3d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + return { + "negative electrode": { + "x_n": {"min": 0, "max": self.param.n.L}, + y_3d: {"min": 0, "max": self.param.L_y}, + z_3d: {"min": 0, "max": self.param.L_z}, + }, + "separator": { + "x_s": { + "min": self.param.n.L, + "max": self.param.n.L + self.param.s.L, + }, + y_3d: {"min": 0, "max": self.param.L_y}, + z_3d: {"min": 0, "max": self.param.L_z}, + }, + "positive electrode": { + "x_p": { + "min": self.param.n.L + self.param.s.L, + "max": self.param.n.L + self.param.s.L + self.param.p.L, + }, + y_3d: {"min": 0, "max": self.param.L_y}, + z_3d: {"min": 0, "max": self.param.L_z}, + }, + "positive particle": { + "r_p": {"min": 0, "max": self.param.p.prim.R_typ}, + }, + "negative particle": { + "r_n": {"min": 0, "max": self.param.n.prim.R_typ}, + }, + "current collector": { + "z": {"position": 0}, + }, + } + + @property + def default_spatial_methods(self): + return { + "negative electrode": pybamm.FiniteVolumeUnstructured(), + "separator": pybamm.FiniteVolumeUnstructured(), + "positive electrode": pybamm.FiniteVolumeUnstructured(), + "positive particle": pybamm.FiniteVolume(), + "negative particle": pybamm.FiniteVolume(), + "current collector": pybamm.ZeroDimensionalSpatialMethod(), + } + + @property + def default_submesh_types(self): + return { + "negative electrode": pybamm.UnstructuredMeshGenerator(), + "separator": pybamm.UnstructuredMeshGenerator(), + "positive electrode": pybamm.UnstructuredMeshGenerator(), + "positive particle": pybamm.Uniform1DSubMesh, + "negative particle": pybamm.Uniform1DSubMesh, + "current collector": pybamm.SubMesh0D, + } + + @property + def default_var_pts(self): + y_3d = pybamm.SpatialVariable( + "y_3d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="fb", + ) + z_3d = pybamm.SpatialVariable( + "z_3d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + return { + "x_n": 5, + "x_s": 5, + "x_p": 5, + "r_p": 10, + "r_n": 10, + y_3d: 3, + z_3d: 3, + } diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py index 837682b9d6..3a73e7a0cb 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py @@ -52,3 +52,50 @@ class TestBasicDFNHalfCell(BaseBasicModelTest): def setup(self): options = {"working electrode": "positive"} self.model = pybamm.lithium_ion.BasicDFNHalfCell(options) + + +class TestBasicDFN2DUnstructured: + def test_solves_and_matches_structured(self): + import numpy as np + + z_2d = pybamm.SpatialVariable( + "z_2d", + domain=["negative electrode", "separator", "positive electrode"], + coord_sys="cartesian", + direction="tb", + ) + var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_p": 10, "r_n": 10, z_2d: 3} + t_eval = np.linspace(0, 3600, 20) + + model_s = pybamm.lithium_ion.BasicDFN2D() + sim_s = pybamm.Simulation(model_s, var_pts=var_pts) + sol_s = sim_s.solve(t_eval) + + model_u = pybamm.lithium_ion.BasicDFN2DUnstructured(element_type="quad") + sim_u = pybamm.Simulation(model_u, var_pts=var_pts) + sol_u = sim_u.solve(t_eval) + + V_s = sol_s["Voltage [V]"](t=t_eval) + V_u = sol_u["Voltage [V]"](t=t_eval) + + np.testing.assert_allclose(V_u, V_s, atol=5e-3) + + +class TestBasicDFN3DUnstructured: + def test_solves_and_matches_2d(self): + import numpy as np + + t_eval = np.linspace(0, 3600, 20) + + model_2d = pybamm.lithium_ion.BasicDFN2DUnstructured(element_type="quad") + sim_2d = pybamm.Simulation(model_2d) + sol_2d = sim_2d.solve(t_eval) + + model_3d = pybamm.lithium_ion.BasicDFN3DUnstructured() + sim_3d = pybamm.Simulation(model_3d) + sol_3d = sim_3d.solve(t_eval) + + V_2d = sol_2d["Voltage [V]"](t=t_eval) + V_3d = sol_3d["Voltage [V]"](t=t_eval) + + np.testing.assert_allclose(V_3d, V_2d, atol=5e-3) diff --git a/packages/pybamm/tests/strategies/symbols.py b/packages/pybamm/tests/strategies/symbols.py index 4ad7a89627..6a4afa08cb 100644 --- a/packages/pybamm/tests/strategies/symbols.py +++ b/packages/pybamm/tests/strategies/symbols.py @@ -1037,6 +1037,7 @@ def _vector_branch( pybamm.StateVectorBase, # abstract base; StateVector + StateVectorDot cover it pybamm.Function, # to_json() raises NotImplementedError — only SpecificFunction subclasses round-trip pybamm.SpecificFunction, # base for named funcs; direct instantiation not useful + pybamm.Reduction, # abstract base for scalar reductions; Max and Min cover it pybamm.Broadcast, # abstract base; PrimaryBroadcast/Secondary/Full covered pybamm.Integral, # base for domain-constrained integrals; heavyweight constructor pybamm.IndependentVariable, # abstract base; Time + SpatialVariable covered diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_symbol.py b/packages/pybamm/tests/unit/test_expression_tree/test_symbol.py index 9e4255dff3..a506db6a6d 100644 --- a/packages/pybamm/tests/unit/test_expression_tree/test_symbol.py +++ b/packages/pybamm/tests/unit/test_expression_tree/test_symbol.py @@ -24,6 +24,14 @@ def test_fixed_domains(self): assert domain_size(["negative electrode"]) == 11 assert domain_size(["separator"]) == 13 assert domain_size(["positive electrode"]) == 17 + assert domain_size(["negative primary particle"]) == 5 + assert domain_size(["negative secondary particle"]) == 5 + assert domain_size(["positive primary particle"]) == 7 + assert domain_size(["positive secondary particle"]) == 7 + assert domain_size(["negative primary particle size"]) == 19 + assert domain_size(["negative secondary particle size"]) == 19 + assert domain_size(["positive primary particle size"]) == 23 + assert domain_size(["positive secondary particle size"]) == 23 def test_fixed_domains_are_additive(self): assert domain_size(["negative electrode", "separator"]) == 11 + 13 diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py index a8b2e7395e..c9f5fd69a4 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py @@ -47,3 +47,11 @@ def test_dfn_2d_vector_field_variable(self): component = solution["Electrolyte current density x [A.m-2]"] assert np.all(np.isfinite(component(t=5))) + + def test_dfn_2d_unstructured(self): + model = pybamm.lithium_ion.BasicDFN2DUnstructured(element_type="quad") + model.check_well_posedness() + + def test_dfn_3d_unstructured(self): + model = pybamm.lithium_ion.BasicDFN3DUnstructured() + model.check_well_posedness() From b6dcc1a22d5f190f69fd32a8706bfe1d1ff7e987 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Fri, 31 Jul 2026 14:39:25 -0700 Subject: [PATCH 2/2] style: match export/strategy whitespace from full branch Co-authored-by: Cursor --- packages/pybamm/src/pybamm/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/packages/pybamm/src/pybamm/__init__.py b/packages/pybamm/src/pybamm/__init__.py index f06bea5525..cd5152282c 100644 --- a/packages/pybamm/src/pybamm/__init__.py +++ b/packages/pybamm/src/pybamm/__init__.py @@ -168,7 +168,6 @@ UserSuppliedSubmesh3D, ) - from .meshes.unstructured_submesh import ( UnstructuredSubMesh, UnstructuredMeshGenerator,