Allow 3D process variables (#4804)

* #4776 initial fix

* Update src/pybamm/solvers/processed_variable_computed.py

* #4776 handle current collector case

* #4776 processed var tests

* #4776 processed var tests

* style: pre-commit fixes

* #4776 debugging processed computed var

* #4776 processed computed var shape

* #4776 add tertiary x average particle concentration

* #4776 test set initial conditions

* #4776 coverage and 3D scikit computed var

* #4776 add 3D on edge test

* #4776 add 3D on edge test processed computed

* style: pre-commit fixes

* #4776 more coverage

* #4776 add postfix

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric G. Kratz <[email protected]>

Author: 3 people

docs/source/examples/notebooks/models/coupled-degradation.ipynb

@@ -20,6 +20,9 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "\n",
+      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.2.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m25.0.1\u001b[0m\n",
+      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
       "Note: you may need to restart the kernel to use updated packages.\n"
      ]
     }
@@ -52,7 +55,7 @@
     "param = pybamm.ParameterValues(\"ORegan2022\")\n",
     "\n",
     "# Choose CasADI fast (we do a short discharge so there are no events, if events are needed choose \"fast with events\")\n",
-    "solver = pybamm.CasadiSolver(mode=\"fast\")"
+    "solver = pybamm.IDAKLUSolver()"
    ]
   },
   {
@@ -108,12 +111,12 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "7e9913de41054c7cb71be948ebf63375",
+       "model_id": "2c6b050e722f45c9921d891d4a3c1c19",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=9.999999999999998, step=0.09999999999999998)…"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=1800.0, step=18.0), Output()), _dom_classes=…"
       ]
      },
      "metadata": {},
@@ -122,7 +125,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fe0c1824eb0>"
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x17e0440d0>"
       ]
      },
      "execution_count": 4,
@@ -131,15 +134,69 @@
     }
    ],
    "source": [
+    "sim = pybamm.Simulation(\n",
+    "    model,\n",
+    "    parameter_values=param,\n",
+    "    C_rate=1,\n",
+    "    solver=solver,\n",
+    "    var_pts=var_pts,\n",
+    "    submesh_types=submesh_types,\n",
+    ")\n",
+    "sim.solve(\n",
+    "    [0, 1800]\n",
+    ")  # solving time kept short for testing purposes, feel free to extend it\n",
+    "sim.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "e86b053e",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "a79ea249a7c94d43a219067c00259084",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=3316.4079454263174, step=33.164079454263174)…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x29e78fd50>"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim = pybamm.Simulation(\n",
+    "    model,\n",
+    "    parameter_values=param,\n",
+    "    C_rate=1,\n",
+    "    solver=solver,\n",
+    "    # var_pts=var_pts,\n",
+    "    # submesh_types=submesh_types,\n",
+    ")\n",
     "sim.solve(\n",
-    "    [0, 10]\n",
+    "    [0, 3600]\n",
     ")  # solving time kept short for testing purposes, feel free to extend it\n",
     "sim.plot()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "lined-illustration",
    "metadata": {},
    "outputs": [
@@ -148,12 +205,15 @@
      "output_type": "stream",
      "text": [
       "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
-      "[2] Chang-Hui Chen, Ferran Brosa Planella, Kieran O'Regan, Dominika Gastol, W. Dhammika Widanage, and Emma Kendrick. Development of Experimental Techniques for Parameterization of Multi-scale Lithium-ion Battery Models. Journal of The Electrochemical Society, 167(8):080534, 2020. doi:10.1149/1945-7111/ab9050.\n",
-      "[3] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
-      "[4] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
-      "[5] Kieran O'Regan, Ferran Brosa Planella, W. Dhammika Widanage, and Emma Kendrick. Thermal-electrochemical parametrisation of a lithium-ion battery: mapping Li concentration and temperature dependencies. Journal of The Electrochemical Society, ():, 2021. doi:.\n",
-      "[6] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
-      "[7] Robert Timms, Scott G Marquis, Valentin Sulzer, Colin P. Please, and S Jonathan Chapman. Asymptotic Reduction of a Lithium-ion Pouch Cell Model. SIAM Journal on Applied Mathematics, 81(3):765–788, 2021. doi:10.1137/20M1336898.\n",
+      "[2] Von DAG Bruggeman. Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. i. dielektrizitätskonstanten und leitfähigkeiten der mischkörper aus isotropen substanzen. Annalen der physik, 416(7):636–664, 1935.\n",
+      "[3] Chang-Hui Chen, Ferran Brosa Planella, Kieran O'Regan, Dominika Gastol, W. Dhammika Widanage, and Emma Kendrick. Development of Experimental Techniques for Parameterization of Multi-scale Lithium-ion Battery Models. Journal of The Electrochemical Society, 167(8):080534, 2020. doi:10.1149/1945-7111/ab9050.\n",
+      "[4] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
+      "[5] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
+      "[6] Alan C. Hindmarsh. The PVODE and IDA algorithms. Technical Report, Lawrence Livermore National Lab., CA (US), 2000. doi:10.2172/802599.\n",
+      "[7] Alan C. Hindmarsh, Peter N. Brown, Keith E. Grant, Steven L. Lee, Radu Serban, Dan E. Shumaker, and Carol S. Woodward. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS), 31(3):363–396, 2005. doi:10.1145/1089014.1089020.\n",
+      "[8] Kieran O'Regan, Ferran Brosa Planella, W. Dhammika Widanage, and Emma Kendrick. Thermal-electrochemical parameters of a high energy lithium-ion cylindrical battery. Electrochimica Acta, 425:140700, 2022. doi:10.1016/j.electacta.2022.140700.\n",
+      "[9] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
+      "[10] Robert Timms, Scott G Marquis, Valentin Sulzer, Colin P. Please, and S Jonathan Chapman. Asymptotic Reduction of a Lithium-ion Pouch Cell Model. SIAM Journal on Applied Mathematics, 81(3):765–788, 2021. doi:10.1137/20M1336898.\n",
       "\n"
      ]
     }
@@ -179,7 +239,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.16"
+   "version": "3.11.5"
   },
   "toc": {
    "base_numbering": 1,

docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb

@@ -783,6 +783,13 @@ def process_symbol(self, symbol):
                 ]
             else:
                 discretised_symbol.secondary_mesh = None
+
+            # Assign tertiary mesh
+            if symbol.domains["tertiary"] != []:
+                discretised_symbol.tertiary_mesh = self.mesh[symbol.domains["tertiary"]]
+            else:
+                discretised_symbol.tertiary_mesh = None
+
             return discretised_symbol
 
     def _process_symbol(self, symbol):

src/pybamm/discretisations/discretisation.py

@@ -126,6 +126,11 @@ def generic_deserialise(cls, instance, properties):
                 else:
                     var.secondary_mesh = None
 
+                if var.domains["tertiary"] != []:
+                    var.tertiary_mesh = properties["mesh"][var.domains["tertiary"]]
+                else:
+                    var.tertiary_mesh = None
+
             if properties["geometry"]:
                 instance._geometry = pybamm.Geometry(properties["geometry"])
         else:
@@ -897,51 +902,48 @@ def set_initial_conditions_from(self, solution, inplace=True, return_type="model
         initial_conditions = {}
         if isinstance(solution, pybamm.Solution):
             solution = solution.last_state
+
+        def get_final_state_eval(final_state):
+            if isinstance(solution, pybamm.Solution):
+                final_state = final_state.data
+
+            if final_state.ndim == 0:
+                return np.array([final_state])
+            elif final_state.ndim == 1:
+                return final_state[-1:]
+            elif final_state.ndim == 2:
+                return final_state[:, -1]
+            elif final_state.ndim == 3:
+                return final_state[:, :, -1].flatten(order="F")
+            elif final_state.ndim == 4:
+                return final_state[:, :, :, -1].flatten(order="F")
+            else:
+                raise NotImplementedError("Variable must be 0D, 1D, 2D, or 3D")
+
+        def get_variable_state(var_name):
+            try:
+                return solution[var_name]
+            except KeyError as e:
+                raise pybamm.ModelError(
+                    "To update a model from a solution, each variable in "
+                    "model.initial_conditions must appear in the solution with "
+                    "the same key as the variable name. In the solution provided, "
+                    f"'{e.args[0]}' was not found."
+                ) from e
+
         for var in self.initial_conditions:
             if isinstance(var, pybamm.Variable):
-                try:
-                    final_state = solution[var.name]
-                except KeyError as e:
-                    raise pybamm.ModelError(
-                        "To update a model from a solution, each variable in "
-                        "model.initial_conditions must appear in the solution with "
-                        "the same key as the variable name. In the solution provided, "
-                        f"'{e.args[0]}' was not found."
-                    ) from e
-                if isinstance(solution, pybamm.Solution):
-                    final_state = final_state.data
-                if final_state.ndim == 0:
-                    final_state_eval = np.array([final_state])
-                elif final_state.ndim == 1:
-                    final_state_eval = final_state[-1:]
-                elif final_state.ndim == 2:
-                    final_state_eval = final_state[:, -1]
-                elif final_state.ndim == 3:
-                    final_state_eval = final_state[:, :, -1].flatten(order="F")
-                else:
-                    raise NotImplementedError("Variable must be 0D, 1D, or 2D")
+                final_state = get_variable_state(var.name)
+                final_state_eval = get_final_state_eval(final_state)
+
             elif isinstance(var, pybamm.Concatenation):
                 children = []
                 for child in var.orphans:
-                    try:
-                        final_state = solution[child.name]
-                    except KeyError as e:
-                        raise pybamm.ModelError(
-                            "To update a model from a solution, each variable in "
-                            "model.initial_conditions must appear in the solution with "
-                            "the same key as the variable name. In the solution "
-                            f"provided, {e.args[0]}"
-                        ) from e
-                    if isinstance(solution, pybamm.Solution):
-                        final_state = final_state.data
-                    if final_state.ndim == 2:
-                        final_state_eval = final_state[:, -1]
-                    else:
-                        raise NotImplementedError(
-                            "Variable in concatenation must be 1D"
-                        )
+                    final_state = get_variable_state(child.name)
+                    final_state_eval = get_final_state_eval(final_state)
                     children.append(final_state_eval)
                 final_state_eval = np.concatenate(children)
+
             else:
                 raise NotImplementedError(
                     "Variable must have type 'Variable' or 'Concatenation'"

src/pybamm/models/base_model.py

@@ -337,15 +337,7 @@ def _get_standard_concentration_distribution_variables(self, c_s):
             c_s_distribution = c_s
 
             # x-average the *tertiary* domain.
-            # NOTE: not yet implemented. Make 0.5 everywhere
-            c_s_xav_distribution = pybamm.FullBroadcast(
-                0.5,
-                [f"{domain} {phase_name}particle"],
-                {
-                    "secondary": f"{domain} {phase_name}particle size",
-                    "tertiary": "current collector",
-                },
-            )
+            c_s_xav_distribution = pybamm.x_average(c_s)
 
             # Surface concentration distribution variables
             c_s_surf_distribution = pybamm.surf(c_s)