feat(zoo): add the multilayer pouch stack with a 3D temperature field per zone

[BradyPlanden]

Description

This brings the multilayer pouch model from #5489 (by mleot) into the model zoo, as the zoo's first contributed model, and resolves #5505. #5489's seven commits are rebased onto main with their authorship, and one commit on top moves the models out of pybamm.lithium_ion into packages/pybamm-model-zoo/src/pybamm_model_zoo/multilayer_3d_thermal/. The net diff touches nothing under packages/pybamm/.

The stack's num_physical_layers unit cells are lumped into num_subdivisions zones through its thickness. Each zone runs an SPM, SPMe, or DFN for its unit cells in parallel and carries its own 3D temperature field on a scikit-fem mesh, whose volume average its kinetics and transport see. Adjacent zones exchange heat through a contact resistance, and the zones connect in parallel, with their current fractions solved for, or in series. A cooling plate on one face or a hot layer can then move current between layers, which Basic3DThermalSPM's single temperature and electrochemistry cannot.

Checking each variant against PyBaMM's own model, held isothermal, turned up four bugs in #5489:

	#5489	This PR
SPMe electrolyte	A net source term: c_e rose from 1000 to 4160 mol.m-3 in 30 minutes, 41 mV from SPMe	Resolves c_e(x) with the composite SPMe's concentration overpotential and electrolyte and solid ohmic drops: 0.34 mV from SPMe, lithium conserved
DFN	Does not build on main, where electrode conductivity takes the stoichiometry	Builds, 0.15 mV from DFN, migration kept inside the electrolyte flux as in BasicDFN
apply_stack_scaling in series	Scaled the capacity by every unit cell, so 1C drove each cell N times too hard	Scales by the unit cells that share the current
Default parameters	No face heat transfer coefficients, so the model could not build from its own defaults	Supplied where a parameter set lacks them

The SPM matched SPM to 0.03 mV before and after.

• pybamm.source refuses every domain but "cell" and "current collector". Add Model for Multilayer Stacked Pouch Cells with thru-stack heat generation & diffusion. #5489 widened that check in core; here _compat.source assembles the same mass-matrix source on the per-zone domains instead, and names the upstream change that would retire it.
• MultiLayer3DThermalSPM is the registered model. MultiLayer3DThermalSPMe and MultiLayer3DThermalDFN share its folder, stacking, thermal coupling and geometry, and override only each zone's electrochemistry.
• The DFN computes its electrolyte ohmic heat from its own fields, where Add Model for Multilayer Stacked Pouch Cells with thru-stack heat generation & diffusion. #5489 used the initial concentration.
• Tests pin the isothermal reduction to all three core models, the adiabatic energy balance (heat stored equals heat generated, which pins the interface coupling), heat equal to current times lost voltage, electrolyte lithium conservation, even current sharing in a symmetric stack, what a C-rate means in parallel and in series, and how one-sided cooling moves current to the warm side and back. They run in about 8 s.
• Two print-only examples replace Add Model for Multilayer Stacked Pouch Cells with thru-stack heat generation & diffusion. #5489's three plotting scripts, since the zoo runs examples with warnings as errors: one-sided cooling of a 24-cell stack, and SPM against SPMe against DFN in the same stack.
• The unreleased num_layers/layers_per_zone keywords are dropped, invalid arguments raise pybamm.OptionError, and "Layer i current [A]" is now the zone's current, with "Layer i per-unit-cell current [A]" beside it.
• mleot is the maintainer in model.toml, in the community tier, so zoo.info() and the docs table send users to them for support. Code owners need write access, which mleot does not have, so @pybamm-team/maintainers own the folder in CODEOWNERS and approve its pull requests as needed.

The README lists what is not validated. The largest gap is the thermal conductivity: it is PyBaMM's lambda_eff, a thickness-weighted in-plane average that includes the current collectors, applied in every direction, so conduction through the stack is overstated.

Also found, for a separate core PR: pybamm.x_average(eps * c_e) returns x_average(eps) * x_average(c_e) when eps is a concatenation of per-region broadcasts, because XAverage.symbol_is_constant treats it as uniform. This entry sums its electrolyte lithium region by region to avoid it.

Type of change

A model zoo entry, so its bullet is in packages/pybamm-model-zoo/CHANGELOG.md and PyBaMM's changelog is untouched.

Important checks:

Please confirm the following before marking the PR as ready for review:

• No style issues: nox -s pre-commit
• All tests pass: nox -s tests
• The documentation builds: nox -s doctests
• Code is commented for hard-to-understand areas
• Tests added that prove fix is effective or that feature works

[mleot]

Hi @BradyPlanden thanks for putting this together. I did some more thinking and I believe that there may be a cleaner implementation which is more DRY in terms of using the existing pybamm models and not defining from scratch. I am working on some changes I will commit them soon.

[mleot]

@BradyPlanden I opened #5840 against this branch, as I couldn't commit to this PR directly. The main difference is that each zone uses PyBaMM's own SPM, SPMe or DFN under the stack's options, so options like particle phases, hysteresis and intercalation kinetics now work in every zone. That replaces the hand-written zone equations, including your SPMe and DFN fixes, which the new zones no longer need. We achieve this with pybamm.replace.

It also fixes four issues I found checking against PyBaMM's lumped model, all present since #5489:

• zones spanned L_x which doesn't include foils, but carried the current collectors' heat capacity, so a stack heated 20% too fast

• through-stack conduction used the in-plane lambda_eff; now they use the proper formualtion:
  $\lambda_\perp = L_{\text{cell}} \big/ \sum_k \left( L_k / \lambda_k \right)$

• the parallel current split was singular at rest; zone currents are now solved directly, so layers can equilibrate at rest while their sum still equals the applied current.

• each unit cell counted whole foils, twice over for double-sided electrodes. Double-sided is now the default.

Layer-by-layer tab resistance is ready too, but I've kept it for a separate PR once it has more evidence behind it. Happy to adjust anything to fit the zoo.