diff --git a/.lycheeignore b/.lycheeignore index 5b74eceeae..8fa7101cd8 100644 --- a/.lycheeignore +++ b/.lycheeignore @@ -16,3 +16,7 @@ file:///home/runner/work/PyBaMM/PyBaMM/docs/source/user_guide/fundamentals/pybam # Telemetry https://us.i.posthog.com + +# Live site behind a LiteSpeed anti-bot WAF that returns 415 to CI/datacenter +# IPs (works fine from browsers and locally) — false positive, not a dead link +https://bpxstandard.com/ diff --git a/CHANGELOG.md b/CHANGELOG.md index a080f18316..e477d3c6b7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +# [v26.6.1.2](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.2) - 2026-06-15 + +## Bug fixes + +- `RegulariseSqrtAndPower` no longer regularises state-independent bases, fixing corrupted small rate constants in exchange-current density functions. ([#5600](https://github.com/pybamm-team/PyBaMM/pull/5600)) +- Fixed `Serialise.load_custom_model` leaving a loaded lithium-ion model's cached `param` built from default options instead of the restored ones, so a model loaded with non-default options (e.g. composite `"particle phases"`) had parameters inconsistent with its options. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599)) +- Fixed `BatteryModelOptions` letting two invalid MSMR option sets pass validation and then crash during parameter construction with `invalid literal for int() with base 10: 'none'`; both now raise a clear `OptionError`. The all-or-nothing MSMR check now detects MSMR inside a per-electrode tuple (e.g. `"open-circuit potential": ("MSMR", "single")`), and an MSMR model must set `"number of MSMR reactions"` to a positive integer rather than the default `"none"`. ([#5599](https://github.com/pybamm-team/PyBaMM/pull/5599)) +- Fixed `Serialise.serialise_custom_model` not recording a custom model's discretisation recipe (`geometry`, `var_pts`, `submesh_types`, `spatial_methods`), so a custom model reloaded via the `BaseModel` fallback (when its defining package is unavailable) was no longer discretisable. ([#5609](https://github.com/pybamm-team/PyBaMM/pull/5609)) + # [v26.6.1.1](https://github.com/pybamm-team/PyBaMM/tree/v26.6.1.1) - 2026-06-11 ## Bug fixes diff --git a/CITATION.cff b/CITATION.cff index cacb1087a9..f5c599bafc 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -24,6 +24,6 @@ keywords: - "expression tree" - "python" - "symbolic differentiation" -version: "26.6.1.1" +version: "26.6.1.2" repository-code: "https://github.com/pybamm-team/PyBaMM" title: "Python Battery Mathematical Modelling (PyBaMM)" diff --git a/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb b/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb index ca561750fb..33bca4759c 100644 --- a/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb +++ b/docs/source/examples/notebooks/performance/01-simulation-pipeline.ipynb @@ -48,7 +48,7 @@ "4. **Solver solve** - The next step is to solve the equations. This is done by the numerical solver, which in most cases is one of the solvers provided by the [Sundials](https://sundials.readthedocs.io/en/latest/#) suite of solvers. The solver takes the CasADi functions generated in the previous step and uses them to solve the equations numerically.\n", "5. **Post-processing** - The output of the previous step is the solution to the equations at a set of time points. The final step is to post-process this solution to get the variables of interest. This is done by the [`pybamm.Solution`](https://docs.pybamm.org/en/stable/source/api/solvers/solution.html) and [`pybamm.ParameterValues`](https://docs.pybamm.org/en/stable/source/api/parameters/parameter_values.html#pybamm.ParameterValues) classes, which take the solution and the parameter values and return the variables of interest interpolated at the desired time points.\n", "\n", - "For every simulation you run with PyBaMM, these steps are always performed in order to generate each solution, even if they are hidden from you by various convenience classes like [`pybamm.Simulation`](https://docs.pybamm.org/en/stable/source/api/simulation.html) . Understanding the pipeline is important because it can help you identify where the performance bottlenecks are in your simulations, and how you can optimise them to get the best performance.\n", + "For every simulation you run with PyBaMM, these steps are always performed in order to generate each solution, even if they are hidden from you by various convenience classes like [`pybamm.Simulation`](https://docs.pybamm.org/en/stable/source/api/simulation/simulation.html) . Understanding the pipeline is important because it can help you identify where the performance bottlenecks are in your simulations, and how you can optimise them to get the best performance.\n", "\n", "Below is a simple example that demonstrates the different parts of the pipeline, using the DFN model. While you probably have seen something similar to this script before, the one difference is that we have split the solver setup (step 3) and solve (step 4) parts of the pipeline into two separate steps. Often these are combined into the first call to the solver, but we have split them here to make it clear that they are separate steps." ] diff --git a/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb b/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb index c77b1596a2..101081f749 100644 --- a/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb +++ b/docs/source/examples/notebooks/solvers/idaklu-jax-interface.ipynb @@ -8,7 +8,7 @@ "source": [ "# IDAKLU-JAX interface\n", "\n", - "The IDAKLU-JAX interface requires that PyBaMM is installed with the [optional JAX solver enabled](https://docs.pybamm.org/en/stable/source/user_guide/installation/gnu-linux-mac.html#optional-jaxsolver) (`pip install pybamm[jax]`) and requires at least Python 3.10.\n", + "The IDAKLU-JAX interface requires that PyBaMM is installed with the [optional JAX solver enabled](https://docs.pybamm.org/en/stable/source/user_guide/installation/index.html#optional-jaxsolver) (`pip install pybamm[jax]`) and requires at least Python 3.10.\n", "\n", "PyBaMM provides two mechanisms to interface battery models with JAX. The first (JaxSolver) implements PyBaMM models directly in native JAX, and as such provides the greatest flexibility. However, these models can be very slow to compile, especially during their initial run, and can require large amounts of memory.\n", "\n", diff --git a/src/pybamm/expression_tree/operations/regularise.py b/src/pybamm/expression_tree/operations/regularise.py index 98bc7b239c..aeac2621df 100644 --- a/src/pybamm/expression_tree/operations/regularise.py +++ b/src/pybamm/expression_tree/operations/regularise.py @@ -11,9 +11,15 @@ class RegulariseSqrtAndPower: """ Callable that replaces Sqrt and Power nodes with RegPower nodes. - All Sqrt and Power nodes are replaced with RegPower. If the base of the - operation matches a symbol in the scales map, that scale is used; - otherwise the default scale (None) is used. + A Sqrt or Power node is replaced with RegPower when the base of the + operation matches a symbol in the scales map (using that scale) or when + the base is state-dependent, i.e. contains a variable, state vector or + time (using the default scale of 1). Bases that are independent of the + state (e.g. ``pybamm.Parameter`` or ``pybamm.InputParameter`` rate + constants) are left unchanged: they have no derivative to regularise, and + RegPower with the default scale of 1 corrupts values smaller than the + regularisation tolerance delta (e.g. a rate constant of ~1e-9 raised to + the power 0.5 evaluates ~470x too small). Parameters ---------- @@ -100,20 +106,25 @@ def _process(self, sym, resolved_scales): new_children = [self._process(child, resolved_scales) for child in sym.children] - if isinstance(sym, pybamm.Sqrt): - child = new_children[0] - scale = self._get_scale(child, resolved_scales) - return pybamm.RegPower(child, 0.5, scale=scale) - - if isinstance(sym, pybamm.Power): - base, exponent = new_children + if isinstance(sym, pybamm.Sqrt | pybamm.Power): + base = new_children[0] + exponent = 0.5 if isinstance(sym, pybamm.Sqrt) else new_children[1] scale = self._get_scale(base, resolved_scales) - return pybamm.RegPower(base, exponent, scale=scale) + if scale is not None or self._depends_on_state(base): + return pybamm.RegPower(base, exponent, scale=scale) if any(n is not o for n, o in zip(new_children, sym.children, strict=True)): return sym.create_copy(new_children=new_children) return sym + @staticmethod + def _depends_on_state(expr): + """Whether the expression depends on the state (or time).""" + return any( + isinstance(node, pybamm.VariableBase | pybamm.StateVectorBase | pybamm.Time) + for node in expr.pre_order() + ) + def _get_scale(self, expr, resolved_scales): """Get scale for an expression, defaulting to None. diff --git a/src/pybamm/expression_tree/operations/serialise.py b/src/pybamm/expression_tree/operations/serialise.py index b11dac109c..3c83580d3a 100644 --- a/src/pybamm/expression_tree/operations/serialise.py +++ b/src/pybamm/expression_tree/operations/serialise.py @@ -534,6 +534,15 @@ def serialise_custom_model(model: pybamm.BaseModel, compress: bool = False) -> d }, } + # Capture the discretisation recipe (geometry / var_pts / submesh / + # spatial methods). These live on subclass ``default_*`` properties and + # are otherwise lost when the model is reloaded in an environment where + # the defining package can't be imported (the base-class fallback yields + # empty defaults), leaving the model undiscretisable. + discretisation = Serialise._serialise_discretisation(model) + if discretisation: + model_content["discretisation"] = discretisation + SCHEMA_VERSION = "1.1" model_json = { "schema_version": SCHEMA_VERSION, @@ -628,6 +637,35 @@ def save_custom_model( except Exception as e: raise ValueError(f"Failed to save custom model: {e}") from e + @staticmethod + def _encode_geometry_value(value): + """Recursively convert any :class:`pybamm.Symbol` in a geometry value to + its JSON form, descending through nested dicts/lists. + + Geometry limits can nest symbols arbitrarily deep (e.g. current-collector + ``tabs -> negative -> z_centre`` is a :class:`pybamm.Parameter`), which a + single-level pass leaves as raw, non-JSON-serialisable objects. + """ + if isinstance(value, pybamm.Symbol): + return convert_symbol_to_json(value) + if isinstance(value, dict): + return {k: Serialise._encode_geometry_value(v) for k, v in value.items()} + if isinstance(value, (list, tuple)): + return [Serialise._encode_geometry_value(v) for v in value] + return value + + @staticmethod + def _decode_geometry_value(value): + """Inverse of :meth:`_encode_geometry_value`: rebuild :class:`pybamm.Symbol` + objects from their JSON form at any nesting depth.""" + if isinstance(value, dict): + if "$type" in value or "type" in value: + return convert_symbol_from_json(value) + return {k: Serialise._decode_geometry_value(v) for k, v in value.items()} + if isinstance(value, list): + return [Serialise._decode_geometry_value(v) for v in value] + return value + @staticmethod def serialise_custom_geometry(geometry: pybamm.Geometry) -> dict: """ @@ -654,26 +692,15 @@ def serialise_custom_geometry(geometry: pybamm.Geometry) -> dict: geometry_dict_serialized[domain]["symbol_" + key_str] = ( convert_symbol_to_json(key) ) - # Serialize the value dict - serialized_value = {} - for k, v in value.items(): - if isinstance(v, pybamm.Symbol): - serialized_value[k] = convert_symbol_to_json(v) - else: - serialized_value[k] = v - geometry_dict_serialized[domain][key_str] = serialized_value + geometry_dict_serialized[domain][key_str] = ( + Serialise._encode_geometry_value(value) + ) elif isinstance(key, str): - # String keys (like 'tabs') - keep as is - if isinstance(value, dict): - serialized_value = {} - for k, v in value.items(): - if isinstance(v, pybamm.Symbol): - serialized_value[k] = convert_symbol_to_json(v) - else: - serialized_value[k] = v - geometry_dict_serialized[domain][key] = serialized_value - else: - geometry_dict_serialized[domain][key] = value + # String keys (like 'tabs'), whose values may nest Symbols + # arbitrarily deep (e.g. tabs -> negative -> z_centre). + geometry_dict_serialized[domain][key] = ( + Serialise._encode_geometry_value(value) + ) SCHEMA_VERSION = "1.1" geometry_json = { @@ -794,29 +821,71 @@ def load_custom_geometry(filename: str | dict) -> pybamm.Geometry: if key in symbol_keys: # Use the reconstructed SpatialVariable as key spatial_var = symbol_keys[key] - reconstructed_value = {} - for k, v in value.items(): - if isinstance(v, dict) and ("$type" in v or "type" in v): - # Reconstruct PyBaMM Symbol using convert_symbol_from_json - reconstructed_value[k] = convert_symbol_from_json(v) - else: - reconstructed_value[k] = v - reconstructed_geometry[domain][spatial_var] = reconstructed_value + reconstructed_geometry[domain][spatial_var] = ( + Serialise._decode_geometry_value(value) + ) else: - # String key (like 'tabs') - if isinstance(value, dict): - reconstructed_value = {} - for k, v in value.items(): - if isinstance(v, dict) and ("$type" in v or "type" in v): - reconstructed_value[k] = convert_symbol_from_json(v) - else: - reconstructed_value[k] = v - reconstructed_geometry[domain][key] = reconstructed_value - else: - reconstructed_geometry[domain][key] = value + # String key (like 'tabs'), values may nest Symbols deep + reconstructed_geometry[domain][key] = ( + Serialise._decode_geometry_value(value) + ) return pybamm.Geometry(reconstructed_geometry) + @staticmethod + def _serialise_discretisation(model) -> dict: + """Serialise a model's discretisation recipe (geometry, var_pts, submesh + types, spatial methods) into a JSON-serialisable dict. + + Returns an empty dict when the model exposes no discretisation defaults + (e.g. a plain ``pybamm.BaseModel``), so models without custom defaults + round-trip exactly as before. + """ + out: dict = {} + geometry = getattr(model, "default_geometry", None) or {} + if geometry: + out["geometry"] = Serialise.serialise_custom_geometry( + pybamm.Geometry(geometry) + ) + var_pts = getattr(model, "default_var_pts", None) or {} + if var_pts: + out["var_pts"] = Serialise.serialise_var_pts(var_pts) + submesh_types = getattr(model, "default_submesh_types", None) or {} + if submesh_types: + out["submesh_types"] = Serialise.serialise_submesh_types(submesh_types) + spatial_methods = getattr(model, "default_spatial_methods", None) or {} + if spatial_methods: + out["spatial_methods"] = Serialise.serialise_spatial_methods( + spatial_methods + ) + return out + + @staticmethod + def _restore_discretisation(model, discretisation: dict | None) -> None: + """Restore a discretisation recipe (from :meth:`_serialise_discretisation`) + onto ``model``, so its ``default_*`` properties return the saved values + even when the defining subclass could not be imported on load. + + No-op when ``discretisation`` is falsy (older payloads or models without + custom discretisation defaults). + """ + if not discretisation: + return + if "geometry" in discretisation: + model._default_geometry = dict( + Serialise.load_custom_geometry(discretisation["geometry"]) + ) + if "var_pts" in discretisation: + model._default_var_pts = Serialise.load_var_pts(discretisation["var_pts"]) + if "submesh_types" in discretisation: + model._default_submesh_types = Serialise.load_submesh_types( + discretisation["submesh_types"] + ) + if "spatial_methods" in discretisation: + model._default_spatial_methods = Serialise.load_spatial_methods( + discretisation["spatial_methods"] + ) + @staticmethod def serialise_spatial_method_item(method) -> dict: """Serialise a spatial method. The class is encoded via the kernel's @@ -1386,10 +1455,9 @@ def load_custom_model(filename: str | dict) -> pybamm.BaseModel: model = base_cls() model.name = model_data["name"] model.schema_version = schema_version - # Restore options so round-trip serialisation produces an equivalent - # model. A JSON round-trip turns tuple-valued options (e.g. "particle - # phases": ("2", "1")) into lists, which pybamm's options validation - # rejects, so convert lists back to tuples before assigning. + # JSON turns tuple-valued options into lists (rejected by validation); + # convert back. Assigning options rebuilds the options-derived param via + # the setter, so it matches the restored options. opts = model_data.get("options", {}) if opts is not None: model.options = Serialise._convert_options(opts) @@ -1495,6 +1563,10 @@ def load_custom_model(filename: str | dict) -> pybamm.BaseModel: # Restore observable state model._solution_observable = False + # Restore the discretisation recipe onto the (possibly fallback) model so + # it stays discretisable even when the defining subclass wasn't importable. + Serialise._restore_discretisation(model, model_data.get("discretisation")) + return model @staticmethod diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py index 333ec48e4f..6432cc39bd 100644 --- a/src/pybamm/models/base_model.py +++ b/src/pybamm/models/base_model.py @@ -641,6 +641,12 @@ def param(self): def param(self, values): self._param = values + def _rebuild_param(self): + """Rebuild ``param`` from ``self.options``; no-op unless ``param`` is + options-derived. Called by the ``options`` setter whenever options are + (re)assigned, so subclasses keep ``param`` in sync automatically. + """ + @property def options(self): """Returns the model options dictionary that allows customization of the model's behavior.""" @@ -682,22 +688,22 @@ def geometry(self): @property def default_var_pts(self): """Returns a dictionary of the default variable points for the model, which is empty by default.""" - return {} + return getattr(self, "_default_var_pts", None) or {} @property def default_geometry(self): """Returns a dictionary of the default geometry for the model, which is empty by default.""" - return {} + return getattr(self, "_default_geometry", None) or {} @property def default_submesh_types(self): """Returns a dictionary of the default submesh types for the model, which is empty by default.""" - return {} + return getattr(self, "_default_submesh_types", None) or {} @property def default_spatial_methods(self): """Returns a dictionary of the default spatial methods for the model, which is empty by default.""" - return {} + return getattr(self, "_default_spatial_methods", None) or {} @property def default_solver(self): diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py index 8bf74ba715..162a7f6beb 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -18,6 +18,24 @@ def represents_positive_integer(s): return val > 0 +def _is_msmr(value): + """True if an option requests MSMR, including inside a per-electrode tuple.""" + if isinstance(value, (tuple, list)): + return "MSMR" in value + return value == "MSMR" + + +def _per_electrode(value, index): + """Resolve a possibly per-electrode option to one electrode. + + Per-electrode tuples are ``(negative, positive)``; a scalar applies to + both. ``index`` is 0 for the negative electrode, 1 for the positive. + """ + if isinstance(value, (tuple, list)): + return value[index] + return value + + def _rename_option(options_dict, option_name, old_name, new_name): if option_name not in options_dict: return @@ -621,13 +639,11 @@ def __init__(self, extra_options): f"Option '{name}' not recognised. Best matches are {options.get_best_matches(name)}" ) - # If any of "open-circuit potential", "particle" or "intercalation kinetics" is - # "MSMR" then all of them must be "MSMR". - # Note: this check is currently performed on full cells, but is loosened for - # half-cells where you must pass a tuple of options to only set MSMR models in - # the working electrode + # All-or-nothing on full cells: if any of OCP/particle/intercalation + # kinetics requests MSMR (incl. inside a per-electrode tuple), all must. + # Half-cells are loosened -- a tuple sets MSMR in the working electrode. msmr_check_list = [ - options[opt] == "MSMR" + _is_msmr(options[opt]) for opt in ["open-circuit potential", "particle", "intercalation kinetics"] ] if ( @@ -640,6 +656,29 @@ def __init__(self, extra_options): "'intercalation kinetics' is 'MSMR' then all of them must be 'MSMR'" ) + # Validate per electrode so a mixed full cell (MSMR in one electrode, + # conventional in the other) is accepted. + if options["working electrode"] == "both": + electrode_domains = [("negative", 0), ("positive", 1)] + else: + electrode_domains = [("positive", 1)] + for domain, index in electrode_domains: + domain_uses_msmr = any( + _per_electrode(options[opt], index) == "MSMR" + for opt in [ + "open-circuit potential", + "particle", + "intercalation kinetics", + ] + ) + domain_count = _per_electrode(options["number of MSMR reactions"], index) + if domain_uses_msmr and not represents_positive_integer(domain_count): + raise pybamm.OptionError( + "'number of MSMR reactions' must be a positive integer for " + f"the {domain} electrode when it uses 'MSMR' " + f"(got {domain_count!r})" + ) + # If "SEI film resistance" is "distributed" then "total interfacial current # density as a state" must be "true" if options["SEI film resistance"] == "distributed": @@ -1209,6 +1248,9 @@ def options(self, extra_options): f"must use surface formulation to solve {self!s} with hydrolysis" ) self._options = options + # rebuild whenever options are (re)assigned. + # No-op unless the subclass overrides ``_rebuild_param``. + self._rebuild_param() def set_standard_output_variables(self): # Time diff --git a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py index 3f48532fda..207b49115f 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py @@ -29,13 +29,17 @@ class BaseModel(pybamm.BaseBatteryModel): def __init__(self, options=None, name="Unnamed lithium-ion model", build=False): super().__init__(options, name) - self.param = pybamm.LithiumIonParameters(self.options) self.set_standard_output_variables() # Li models should calculate esoh by default self._calc_esoh = True + def _rebuild_param(self): + # param is options-derived, so rebuild it whenever options are + # (re)assigned: at construction and on deserialisation. + self.param = pybamm.LithiumIonParameters(self.options) + def set_submodels(self, build): self.set_external_circuit_submodel() self.set_porosity_submodel() diff --git a/tests/unit/test_expression_tree/test_operations/test_regularise.py b/tests/unit/test_expression_tree/test_operations/test_regularise.py index c17ea9555e..12b20ae568 100644 --- a/tests/unit/test_expression_tree/test_operations/test_regularise.py +++ b/tests/unit/test_expression_tree/test_operations/test_regularise.py @@ -2,6 +2,8 @@ # Tests for the RegulariseSqrtAndPower class # +import pytest + import pybamm @@ -40,16 +42,16 @@ def test_basic_power_replacement(self): expr = c_e**0.5 result = regulariser(expr, inputs=inputs) - # The original c_e**0.5 Power node should be replaced - # (result will have different Power nodes from reg_power formula) - assert result != expr + # Matched base: RegPower with the registered scale + assert isinstance(result, pybamm.RegPower) + assert result.children[2] == pybamm.Scalar(1000.0) - def test_scale_default_to_one(self): - """Test that unrecognized expressions get scale=None.""" + def test_state_dependent_base_default_scale(self): + """State-dependent base with no registered scale: regularised with + the default scale of 1.""" c_e = pybamm.Variable("c_e") c_s = pybamm.Variable("c_s") - # Only c_e has a scale, c_s should get scale=None inputs = {"c_e": c_e, "c_s": c_s} regulariser = pybamm.RegulariseSqrtAndPower( {c_e: pybamm.Scalar(1000.0)}, @@ -59,14 +61,162 @@ def test_scale_default_to_one(self): expr = pybamm.sqrt(c_s) result = regulariser(expr, inputs=inputs) - # Should be replaced with RegPower (no Sqrt) - has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result.pre_order()) - assert not has_sqrt - # Check it's a RegPower with scale=1 (default) assert isinstance(result, pybamm.RegPower) - # Scale is the third child assert result.children[2] == pybamm.Scalar(1) + def test_normalised_state_dependent_base_still_regularised(self): + """Normalised base (ORegan2022-style (c_e / c_e_ref) ** (1 - alpha)) + matches no scale but is still regularised with the default scale.""" + c_e = pybamm.Variable("c_e") + inputs = {"c_e": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + + c_e_ref = pybamm.Parameter("Reference concentration") + expr = (c_e / c_e_ref) ** 0.208 + result = regulariser(expr, inputs=inputs) + + assert isinstance(result, pybamm.RegPower) + assert result.children[2] == pybamm.Scalar(1) + + def test_parameter_power_not_corrupted(self): + """State-independent Parameter base is not regularised: RegPower with + the default scale would corrupt a small rate constant.""" + c_e = pybamm.Variable("c_e") + inputs = {"c_e": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + + rate_constant = pybamm.Parameter("Rate constant") + expr = rate_constant**0.5 + result = regulariser(expr, inputs=inputs) + + assert isinstance(result, pybamm.Power) + assert not any(isinstance(n, pybamm.RegPower) for n in result.pre_order()) + + parameter_values = pybamm.ParameterValues({"Rate constant": 5e-9}) + value = parameter_values.process_symbol(result).evaluate() + assert value == pytest.approx(5e-9**0.5, rel=1e-12) + + def test_input_parameter_power_not_corrupted(self): + """State-independent InputParameter base is not regularised.""" + c_e = pybamm.Variable("c_e") + inputs = {"c_e": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + + rate_constant = pybamm.InputParameter("Rate constant") + expr = rate_constant**0.5 + result = regulariser(expr, inputs=inputs) + + assert isinstance(result, pybamm.Power) + value = result.evaluate(inputs={"Rate constant": 5e-9}) + assert value == pytest.approx(5e-9**0.5, rel=1e-12) + + def test_constant_base_state_dependent_exponent(self): + """Constant base, state-dependent exponent: outer Power kept, sqrt in + the exponent regularised.""" + c_e = pybamm.Variable("c_e") + inputs = {"c_e": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + + base = pybamm.Parameter("Rate constant") + expr = base ** pybamm.sqrt(c_e) + result = regulariser(expr, inputs=inputs) + + assert isinstance(result, pybamm.Power) + assert isinstance(result.children[0], pybamm.Parameter) + assert any(isinstance(n, pybamm.RegPower) for n in result.pre_order()) + assert not any(isinstance(n, pybamm.Sqrt) for n in result.pre_order()) + + def test_existing_reg_power_preserved(self): + """A manual RegPower node (e.g. MSMR j0_j) is a Function, not a + Sqrt/Power, so it is left untouched, even when nested.""" + c_e = pybamm.Variable("c_e") + c_e_ref = pybamm.Parameter("c_e_ref") + aj = pybamm.Parameter("aj") + inputs = {"c_e": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + + expr = pybamm.reg_power(c_e / c_e_ref, 1 - aj) + result = regulariser(expr, inputs=inputs) + assert result is expr + + # Nested: reg_power preserved, sibling Power regularised + nested = pybamm.reg_power(c_e / c_e_ref, 1 - aj) * c_e**0.5 + result = regulariser(nested, inputs=inputs) + n_reg_power = sum( + 1 for n in result.pre_order() if isinstance(n, pybamm.RegPower) + ) + assert n_reg_power == 2 + assert not any(isinstance(n, pybamm.Power) for n in result.pre_order()) + assert not any(isinstance(n, pybamm.Sqrt) for n in result.pre_order()) + + @pytest.mark.parametrize("set_name", ["Chen2020", "ORegan2022", "Ecker2015"]) + def test_j0_pipeline_produces_regpower(self, set_name): + """Built-in j0 pipeline regularises all concentration powers, + including ORegan2022's unmatched normalised bases.""" + param = pybamm.LithiumIonParameters() + c_e = pybamm.Variable("c_e") + c_s_surf = pybamm.Variable("c_s_surf") + T = pybamm.Variable("T") + j0 = param.n.prim.j0(c_e, c_s_surf, T) + + parameter_values = pybamm.ParameterValues(set_name) + processed = parameter_values.process_symbol(j0) + + n_reg_power = sum( + 1 for n in processed.pre_order() if isinstance(n, pybamm.RegPower) + ) + n_plain_power = sum( + 1 for n in processed.pre_order() if isinstance(n, pybamm.Power) + ) + assert n_reg_power == 3, f"{set_name}: expected 3 RegPower, got {n_reg_power}" + assert n_plain_power == 0, ( + f"{set_name}: expected no plain Power nodes, got {n_plain_power}" + ) + + def test_function_parameter_post_processor_small_parameter(self): + """Full j0-style path: a FunctionParameter raising a Parameter to a + fractional power must evaluate exactly.""" + + def exchange_current(c_e): + rate_constant = pybamm.Parameter("Rate constant") + return rate_constant**0.5 * c_e**0.5 + + c_e = pybamm.Variable("c_e") + inputs = {"Electrolyte concentration [mol.m-3]": c_e} + regulariser = pybamm.RegulariseSqrtAndPower( + {c_e: pybamm.Scalar(1000.0)}, + inputs=inputs, + ) + function_parameter = pybamm.FunctionParameter( + "Exchange-current density [A.m-2]", + {"Electrolyte concentration [mol.m-3]": pybamm.Scalar(1000.0)}, + post_processor=regulariser, + ) + parameter_values = pybamm.ParameterValues( + { + "Exchange-current density [A.m-2]": exchange_current, + "Rate constant": 5e-9, + } + ) + value = parameter_values.process_symbol(function_parameter).evaluate() + expected = 5e-9**0.5 * 1000.0**0.5 + assert value == pytest.approx(expected, rel=1e-6) + def test_exact_match_only(self): """Test that only exact matches are used for scales.""" c_s = pybamm.Variable("c_s") @@ -85,7 +235,7 @@ def test_exact_match_only(self): has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result1.pre_order()) assert not has_sqrt - # sqrt(c_s / c_s_max) should also be replaced + # sqrt(c_s / c_s_max) should also be replaced (state-dependent base) expr2 = pybamm.sqrt(c_s / c_s_max) result2 = regulariser(expr2, inputs=inputs) has_sqrt = any(isinstance(n, pybamm.Sqrt) for n in result2.pre_order()) diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py index cf25799524..cd14f8b638 100644 --- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py +++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py @@ -444,6 +444,60 @@ def test_options(self): "number of MSMR reactions": "1.5", } ) + # MSMR inside a per-electrode tuple must be rejected cleanly, not slip + # through to crash later as int('none') in parameter construction. + with pytest.raises(pybamm.OptionError, match=r"MSMR"): + pybamm.BaseBatteryModel({"open-circuit potential": ("MSMR", "single")}) + with pytest.raises(pybamm.OptionError, match=r"MSMR"): + pybamm.BaseBatteryModel({"particle": ("MSMR", "Fickian diffusion")}) + # MSMR with the default "none" reaction count must also be rejected cleanly. + with pytest.raises(pybamm.OptionError, match=r"MSMR"): + pybamm.BaseBatteryModel( + { + "open-circuit potential": "MSMR", + "particle": "MSMR", + "intercalation kinetics": "MSMR", + } + ) + + def test_msmr_mixed_electrode_accepted(self): + # Verify non-MSMR electrode's "none" reaction count is legitimate + options = { + "open-circuit potential": ("MSMR", "single"), + "particle": ("MSMR", "Fickian diffusion"), + "intercalation kinetics": ("MSMR", "symmetric Butler-Volmer"), + "number of MSMR reactions": ("3", "none"), + } + model = pybamm.BaseBatteryModel(options) + assert model.options["number of MSMR reactions"] == ("3", "none") + + def test_msmr_half_cell_does_not_validate_counter_electrode(self): + # On a half cell the negative electrode is lithium metal, not a porous + # MSMR electrode, so a scalar "MSMR" option must not force a reaction + # count on it. + options = { + "working electrode": "positive", + "open-circuit potential": "MSMR", + "particle": "MSMR", + "intercalation kinetics": "MSMR", + "number of MSMR reactions": ("none", "6"), + } + model = pybamm.BaseBatteryModel(options) + assert model.options["number of MSMR reactions"] == ("none", "6") + + def test_msmr_half_cell_still_validates_working_electrode(self): + # The working electrode is still validated: a "none" count + # there must be rejected cleanly rather than slipping through. + with pytest.raises(pybamm.OptionError, match=r"positive electrode"): + pybamm.BaseBatteryModel( + { + "working electrode": "positive", + "open-circuit potential": "MSMR", + "particle": "MSMR", + "intercalation kinetics": "MSMR", + "number of MSMR reactions": ("6", "none"), + } + ) def test_build_twice(self): model = pybamm.lithium_ion.SPM() # need to pick a model to set vars and build diff --git a/tests/unit/test_serialisation/test_round_trip.py b/tests/unit/test_serialisation/test_round_trip.py index 99dc8ffbcc..d8f1d58fca 100644 --- a/tests/unit/test_serialisation/test_round_trip.py +++ b/tests/unit/test_serialisation/test_round_trip.py @@ -11,12 +11,17 @@ from __future__ import annotations +import ast +import inspect import json +import textwrap import warnings from datetime import datetime import numpy as np import pytest +from hypothesis import given, reject, settings +from hypothesis import strategies as st import pybamm from pybamm.expression_tree.operations.serialise import ( @@ -24,6 +29,7 @@ convert_symbol_from_json, convert_symbol_to_json, ) +from pybamm.models.full_battery_models.base_battery_model import BatteryModelOptions from tests.unit.test_serialisation._helpers import ( _experiments_equal, _solver_init_args_equal, @@ -500,3 +506,342 @@ def test_preserves_structure(self, symbol): f"{type(symbol).__name__} structural id changed across round-trip: " f"{symbol!r} -> {restored!r}" ) + + +# Generative option round-trip fuzzing. JSON has no tuple type, so a tuple option +# (e.g. "particle phases": ("2", "1")) deserialises as a list, which validation +# rejects. Fuzzes the build-free to_json/from_json path; the deterministic +# fixtures below cover to_config/from_config (which builds the model). + +# Option registry, sourced from pybamm so the strategy never hardcodes values. +_POSSIBLE_OPTIONS = BatteryModelOptions({}).possible_options + +_STRING_VALUES = { + key: [v for v in values if isinstance(v, str)] + for key, values in _POSSIBLE_OPTIONS.items() +} + +# Drawable keys. "dimensionality" (the only int option) is excluded: a non-zero +# value needs companion geometry options the random strategy can't assemble; its +# int round-trip is covered by test_int_option_survives_json_round_trip_as_int. +_STRING_OPTION_KEYS = sorted( + key for key, values in _STRING_VALUES.items() if key != "dimensionality" and values +) + +# Options accepting a 2-tuple of strings; mirrors the list in BatteryModelOptions +# (test_tuple_capable_keys_match_validation_source guards drift). +_TUPLE_CAPABLE_KEYS = frozenset( + { + "diffusivity", + "exchange-current density", + "intercalation kinetics", + "interface utilisation", + "lithium plating", + "loss of active material", + "number of MSMR reactions", + "open-circuit potential", + "particle", + "particle mechanics", + "particle phases", + "particle size", + "SEI", + "SEI on cracks", + "stress-induced diffusion", + } +) + +# SPM matches the original regression; DFN has a broader option surface. +_OPTION_MODEL_CLASSES = [pybamm.lithium_ion.SPM, pybamm.lithium_ion.DFN] + + +def _assert_param_matches_options(model): + """Assert a loaded lithium-ion model's ``param`` matches its restored options. + + Guards #5598; lead-acid ``param`` is options-independent, so skip it. + """ + if isinstance(model, pybamm.lithium_ion.BaseModel): + assert dict(model.param.options) == dict(model.options), ( + f"param built from {dict(model.param.options)!r} != restored options " + f"{dict(model.options)!r}; load_custom_model must rebuild param." + ) + # Symptom-level guard: a two-phase negative electrode names its primary + # active-material fraction "Primary: ..."; a param left built from + # single-phase options would not. Only meaningful on a full cell -- a + # half cell's negative electrode is lithium metal with no porous + # particle, so ``param.n.prim`` has no ``epsilon_s``. + if ( + model.options["working electrode"] == "both" + and model.options.negative["particle phases"] == "2" + ): + epsilon_s_name = model.param.n.prim.epsilon_s.name + assert epsilon_s_name.startswith("Primary:"), ( + f"negative primary active material fraction is {epsilon_s_name!r}; " + "expected a 'Primary:' prefix for a two-phase electrode -- param " + "was not rebuilt from the restored composite options." + ) + + +def test_options_setter_rebuilds_param(): + """Assigning ``options`` post-construction rebuilds the options-derived + ``param`` via the setter, keeping it in sync without an explicit hook call. + """ + # A complete, valid two-phase options set (as a load path would restore). + full_opts = dict(pybamm.lithium_ion.SPM({"particle phases": ("2", "1")}).options) + + model = pybamm.lithium_ion.SPM(build=False) # default: single particle phase + model.options = full_opts + # The setter rebuilds param, so it tracks the reassigned options. + assert dict(model.param.options) == dict(model.options) + + +@st.composite +def valid_option_dicts(draw): + """Draw a small registry-valid option dict, biased to exercise tuples.""" + keys = draw( + st.lists( + st.sampled_from(_STRING_OPTION_KEYS), + min_size=1, + max_size=3, + unique=True, + ) + ) + options = {} + for key in keys: + choices = _STRING_VALUES[key] + if key in _TUPLE_CAPABLE_KEYS and draw(st.booleans()): + options[key] = ( + draw(st.sampled_from(choices)), + draw(st.sampled_from(choices)), + ) + else: + options[key] = draw(st.sampled_from(choices)) + return options + + +def test_tuple_capable_keys_match_validation_source(): + """``_TUPLE_CAPABLE_KEYS`` must match the tuple-capable list in + BatteryModelOptions -- AST-extracted so it can't silently drift. + """ + tree = ast.parse(textwrap.dedent(inspect.getsource(BatteryModelOptions))) + candidates = [ + { + elt.value + for elt in node.elts + if isinstance(elt, ast.Constant) and isinstance(elt.value, str) + } + for node in ast.walk(tree) + if isinstance(node, ast.List) + ] + # The tuple-capable list is the one literal containing these marker keys. + matches = [c for c in candidates if {"particle phases", "SEI on cracks"} <= c] + assert len(matches) == 1, ( + "Could not uniquely locate the tuple-capable option list in " + "BatteryModelOptions source; update the extraction in this test." + ) + assert matches[0] == set(_TUPLE_CAPABLE_KEYS), ( + "Tuple-capable options in pybamm have changed. Update " + "_TUPLE_CAPABLE_KEYS so option fuzzing keeps generating tuple values " + f"for every such key. Source: {sorted(matches[0])}, " + f"test constant: {sorted(_TUPLE_CAPABLE_KEYS)}." + ) + assert _TUPLE_CAPABLE_KEYS <= set(_POSSIBLE_OPTIONS), ( + "_TUPLE_CAPABLE_KEYS contains keys absent from the option registry." + ) + + +class TestModelOptionsRoundTrip: + """Any valid option dict must survive a JSON round-trip with exact types.""" + + @pytest.mark.parametrize( + "model_cls", _OPTION_MODEL_CLASSES, ids=lambda c: c.__name__ + ) + @settings(max_examples=100, deadline=None) + @given(options=valid_option_dicts()) + def test_options_survive_json_round_trip(self, model_cls, options): + try: + # Construction is the precondition, not the SUT; pass a copy + # (pybamm mutates the dict). + model = model_cls(dict(options), build=False) + except (pybamm.OptionError, NotImplementedError): + # Invalid/unsupported option set -> discard. Any other exception is a + # real bug that must surface (a broad except here once hid an MSMR + # int('none') crash). + reject() + + expected = dict(model.options) + + via_json = pybamm.BaseModel.from_json( + json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder)) + ) + assert _values_equal(expected, dict(via_json.options)), ( + f"{model_cls.__name__} options lost or mutated via " + f"to_json/from_json for input {options!r}: " + f"{expected!r} -> {dict(via_json.options)!r}" + ) + # Options surviving isn't enough -- the derived param must match too. + _assert_param_matches_options(via_json) + + +def test_int_option_survives_json_round_trip_as_int(): + """``dimensionality`` (the only int option) must round-trip as an int, not be + coerced to/from bool (the #5495 shape, ``True == 1``). Bool-like options are + strings "true"/"false" (pybamm rejects Python bools), already fuzzed; a + non-zero value needs companion geometry options to construct. + """ + options = { + "dimensionality": 1, + "current collector": "potential pair", + "cell geometry": "pouch", + } + model = pybamm.lithium_ion.DFN(dict(options), build=False) + + via_json = pybamm.BaseModel.from_json( + json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder)) + ) + restored = via_json.options["dimensionality"] + assert restored == 1 + assert isinstance(restored, int) and not isinstance(restored, bool) + + +# Buildable tuple options that must round-trip through both load paths -- +# deterministic, always-run coverage of the bug shape. +_TUPLE_OPTION_FIXTURES = [ + (pybamm.lithium_ion.SPM, "particle phases", ("2", "1")), + (pybamm.lithium_ion.SPM, "SEI", ("none", "constant")), + (pybamm.lithium_ion.DFN, "particle size", ("single", "distribution")), + (pybamm.lithium_ion.DFN, "particle", ("Fickian diffusion", "uniform profile")), +] + + +class TestTupleOptionRoundTripRegression: + @pytest.mark.parametrize( + "model_cls,key,value", + _TUPLE_OPTION_FIXTURES, + ids=["SPM-particle phases", "SPM-SEI", "DFN-particle size", "DFN-particle"], + ) + def test_tuple_option_survives_both_paths(self, model_cls, key, value): + model = model_cls({key: value}, build=False) + + via_json = pybamm.BaseModel.from_json( + json.loads(json.dumps(model.to_json(), default=Serialise._json_encoder)) + ) + assert via_json.options[key] == value + assert isinstance(via_json.options[key], tuple) + # Composite "particle phases" is the case (param stuck single-phase). + _assert_param_matches_options(via_json) + + via_config = pybamm.BaseModel.from_config( + json.loads(json.dumps(model.to_config())) + ) + assert via_config.options[key] == value + assert isinstance(via_config.options[key], tuple) + + +# Discretised-model round-trip (save_model/load_model): previously only checked +# "load and solve". These verify structural identity and solution equivalence +# (mesh transitively: same equations + solution => same mesh). + +_DISCRETISED_MODEL_CLASSES = [pybamm.lithium_ion.SPM, pybamm.lithium_ion.DFN] + + +def _discretised_model(model_cls, options=None): + """Build, parameterise and discretise a model; return it and its mesh + (the mesh is needed by ``serialise_model`` for named variables). + """ + model = model_cls(options) + geometry = model.default_geometry + param = model.default_parameter_values + param.process_model(model) + param.process_geometry(geometry) + mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) + disc = pybamm.Discretisation(mesh, model.default_spatial_methods) + disc.process_model(model) + return model, mesh + + +def _round_trip_discretised(model, mesh): + """Serialise a discretised model and load it back, via dict (no temp file).""" + serialised = Serialise().serialise_model(model, mesh) + return Serialise().load_model(json.loads(json.dumps(serialised))) + + +@pytest.mark.parametrize( + "model_cls", _DISCRETISED_MODEL_CLASSES, ids=lambda c: c.__name__ +) +class TestDiscretisedModelRoundTrip: + def test_derived_state_matches_original(self, model_cls): + """Concatenated equations, mass matrix, bounds and options must match.""" + model, mesh = _discretised_model(model_cls) + loaded = _round_trip_discretised(model, mesh) + + for attr in ( + "concatenated_rhs", + "concatenated_algebraic", + "concatenated_initial_conditions", + "mass_matrix", + ): + original = getattr(model, attr) + restored = getattr(loaded, attr) + assert original.id == restored.id, ( + f"{model_cls.__name__}.{attr} changed across the discretised " + f"round-trip: {original.id} != {restored.id}" + ) + + assert all( + np.array_equal(a, b) + for a, b in zip(model.bounds, loaded.bounds, strict=True) + ), f"{model_cls.__name__} bounds changed across the discretised round-trip" + assert dict(model.options) == dict(loaded.options) + + # Events are serialised too -- verify them structurally rather than + # relying on the solution check. + assert len(model.events) == len(loaded.events) + for original_event, restored_event in zip( + model.events, loaded.events, strict=True + ): + assert original_event.name == restored_event.name + assert original_event.event_type == restored_event.event_type + assert original_event.expression.id == restored_event.expression.id, ( + f"{model_cls.__name__} event {original_event.name!r} changed " + "across the discretised round-trip" + ) + + def test_solution_matches_original(self, model_cls): + """The reloaded model must solve to the same answer as the original.""" + model, mesh = _discretised_model(model_cls) + # Serialise before solving so any solve-time caching cannot leak in. + loaded = _round_trip_discretised(model, mesh) + + t_eval = [0, 3600] + original_solution = model.default_solver.solve(model, t_eval) + loaded_solution = loaded.default_solver.solve(loaded, t_eval) + + assert np.allclose( + original_solution.y, loaded_solution.y, rtol=1e-6, atol=1e-8 + ), ( + f"{model_cls.__name__} state trajectory diverged after a " + "discretised round-trip" + ) + assert np.allclose( + original_solution["Voltage [V]"].entries, + loaded_solution["Voltage [V]"].entries, + rtol=1e-6, + atol=1e-8, + ), f"{model_cls.__name__} Voltage [V] diverged after a discretised round-trip" + + +def test_nondefault_tuple_option_survives_discretised_round_trip(): + """A tuple-valued option must survive the discretised save/load path too. + + The parametrised tests above only use default options; this guards + ``load_model``'s list->tuple conversion. ``particle`` is used because it + discretises with the default parameter set (composite ``particle phases`` + would need a phase-prefixed one). + """ + options = {"particle": ("Fickian diffusion", "uniform profile")} + model, mesh = _discretised_model(pybamm.lithium_ion.DFN, options) + loaded = _round_trip_discretised(model, mesh) + + assert loaded.options["particle"] == ("Fickian diffusion", "uniform profile") + assert isinstance(loaded.options["particle"], tuple) + _assert_param_matches_options(loaded) diff --git a/tests/unit/test_serialisation/test_serialisation.py b/tests/unit/test_serialisation/test_serialisation.py index c1f7ba4455..f7a966ab21 100644 --- a/tests/unit/test_serialisation/test_serialisation.py +++ b/tests/unit/test_serialisation/test_serialisation.py @@ -844,6 +844,112 @@ def test_load_legacy_payload_without_base_class_mro(self, tmp_path): assert type(loaded) is pybamm.BaseModel + def test_load_restores_discretisation_after_base_model_fallback(self, tmp_path): + """A custom model's discretisation recipe (geometry / var_pts / submesh + types / spatial methods) survives a load that falls back to + ``pybamm.BaseModel`` because the defining subclass can't be imported.""" + + class _DiscModel(pybamm.BaseModel): + @property + def default_geometry(self): + return { + "positive particle": { + "r_p": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1e-5)} + }, + "current collector": {"z": {"position": 1}}, + } + + @property + def default_var_pts(self): + return {"r_p": 20} + + @property + def default_submesh_types(self): + return { + "positive particle": pybamm.Uniform1DSubMesh, + "current collector": pybamm.SubMesh0D, + } + + @property + def default_spatial_methods(self): + return { + "positive particle": pybamm.FiniteVolume(), + "current collector": pybamm.ZeroDimensionalSpatialMethod(), + } + + model = _DiscModel(name="DiscModel") + a = pybamm.Variable("a") + model.rhs = {a: a} + model.initial_conditions = {a: pybamm.Scalar(1)} + model.variables = {"a": a} + + file_path = tmp_path / "model.json" + Serialise.save_custom_model(model, filename=str(file_path)) + + # The discretisation recipe is captured at serialise time... + with open(file_path) as f: + data = json.load(f) + assert "discretisation" in data["model"] + # ...and force the import-failed fallback to a bare pybamm.BaseModel. + data["model"]["base_class"] = "nonexistent_module.DiscModel" + data["model"]["base_class_mro"] = ["nonexistent_module.DiscModel"] + with open(file_path, "w") as f: + json.dump(data, f) + + with pytest.warns(UserWarning, match=r"Falling back to pybamm\.BaseModel"): + loaded = Serialise.load_custom_model(str(file_path)) + + assert type(loaded) is pybamm.BaseModel + # The fallback model is still discretisable: defaults are restored. + assert list(loaded.default_geometry.keys()) == [ + "positive particle", + "current collector", + ] + assert loaded.default_var_pts == {"r_p": 20} + assert set(loaded.default_submesh_types) == { + "positive particle", + "current collector", + } + assert set(loaded.default_spatial_methods) == { + "positive particle", + "current collector", + } + + def test_plain_custom_model_has_no_discretisation_section(self, tmp_path): + """A model with no custom discretisation defaults emits no + ``discretisation`` section, so existing payloads are unaffected.""" + model = pybamm.BaseModel(name="Plain") + a = pybamm.Variable("a") + model.rhs = {a: -a} + model.initial_conditions = {a: pybamm.Scalar(1)} + model.variables = {"a": a} + + payload = Serialise.serialise_custom_model(model) + assert "discretisation" not in payload["model"] + # Round-trips without error. + Serialise.load_custom_model(payload) + + def test_custom_geometry_round_trip_with_nested_tab_symbols(self): + """Current-collector geometry nests symbols arbitrarily deep + (``tabs -> negative -> z_centre`` is a Parameter). The geometry must + JSON-round-trip without leaking raw, non-serialisable Symbol objects.""" + options = { + "dimensionality": 1, + "current collector": "potential pair", + "cell geometry": "pouch", + } + model = pybamm.lithium_ion.DFN(dict(options), build=False) + geometry = pybamm.Geometry(model.default_geometry) + + serialised = Serialise.serialise_custom_geometry(geometry) + # The whole payload must be JSON-serialisable (no raw Symbols left). + serialised = json.loads(json.dumps(serialised, default=Serialise._json_encoder)) + loaded = Serialise.load_custom_geometry(serialised) + + z_centre = loaded["current collector"]["tabs"]["negative"]["z_centre"] + assert isinstance(z_centre, pybamm.Parameter) + assert z_centre.name == "Negative tab centre z-coordinate [m]" + def test_function_parameter_with_diff_variable_serialisation(self): x = pybamm.Variable("x") diff_var = pybamm.Variable("r")