Merge pull request #1278 from pybamm-team/issue-1194-initial-conditions

Issue 1194 initial conditions

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added functionality to initialize a model using the solution from another model ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278))
 -   Added submodels for active material ([#1262](https://github.com/pybamm-team/PyBaMM/pull/1262))
 -   Added composite surface form electrolyte models: `CompositeDifferential` and `CompositeAlgebraic` ([#1207](https://github.com/pybamm-team/PyBaMM/issues/1207))
 
@@ -11,6 +12,7 @@
 
 ## Bug fixes
 
+-   Fixed `Simulation` and `model.new_copy` to fix a bug where changes to the model were overwritten ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278))
 ## Breaking changes
 
 -   Operations such as `1*x` and `0+x` now directly return `x`. This can be bypassed by explicitly creating the binary operators, e.g. `pybamm.Multiplication(1, x)` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252))

docs/source/models/submodels/external_circuit/index.rst

@@ -13,4 +13,3 @@ variable to be constant.
 
   current_control_external_circuit
   function_control_external_circuit
-  experiment_events

examples/notebooks/initialize-model-with-solution.ipynb

@@ -0,0 +1,316 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Initializing a model"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Example showing how to initialize a model with another model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[33mWARNING: You are using pip version 20.2.4; however, version 20.3 is available.\n",
+      "You should consider upgrading via the '/Users/vsulzer/Documents/Energy_storage/PyBaMM/.tox/dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n",
+      "Note: you may need to restart the kernel to use updated packages.\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pip install pybamm -q\n",
+    "\n",
+    "import pybamm\n",
+    "import pandas as pd\n",
+    "import os\n",
+    "\n",
+    "os.chdir(pybamm.__path__[0]+'/..')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Solve a model with a drive cycle"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = pybamm.lithium_ion.DFN()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Set up drive cycle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import drive cycle from file\n",
+    "drive_cycle = pd.read_csv(\n",
+    "    \"pybamm/input/drive_cycles/US06.csv\", comment=\"#\", header=None\n",
+    ").to_numpy()\n",
+    "# create interpolant\n",
+    "param = model.default_parameter_values\n",
+    "timescale = param.evaluate(model.timescale)\n",
+    "current_interpolant = pybamm.Interpolant(drive_cycle, timescale * pybamm.t)\n",
+    "# set drive cycle\n",
+    "param[\"Current function [A]\"] = current_interpolant"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Create and run simulation using the CasadiSolver in \"fast\" mode, remembering to pass in the updated parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim_US06_1 = pybamm.Simulation(\n",
+    "    model, parameter_values=param, solver=pybamm.CasadiSolver(mode=\"fast\")\n",
+    ")\n",
+    "sol_US06_1 = sim_US06_1.solve()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Update initial conditions based on a solution and solve again"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now pre-charge with CCCV, update the initial conditions, and solve again with the US06 drive cycle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "experiment = pybamm.Experiment(\n",
+    "    [\"Charge at 1 A until 4.1 V\", \"Hold at 4.1 V until 50 mA\"]\n",
+    ")\n",
+    "sim_cccv = pybamm.Simulation(model, experiment=experiment)\n",
+    "sol_cccv = sim_cccv.solve()\n",
+    "\n",
+    "# MODEL RE-INITIALIZATION: #############################################################\n",
+    "# Now initialize the model with the solution of the charge, and then discharge with\n",
+    "# the US06 drive cycle\n",
+    "# We could also do this inplace by setting inplace to True, which modifies the original\n",
+    "# model in place\n",
+    "new_model = model.set_initial_conditions_from(sol_cccv, inplace=False)\n",
+    "########################################################################################\n",
+    "\n",
+    "sim_US06_2 = pybamm.Simulation(\n",
+    "    new_model, parameter_values=param, solver=pybamm.CasadiSolver(mode=\"fast\")\n",
+    ")\n",
+    "sol_US06_2 = sim_US06_2.solve()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot both solutions, we can clearly see the difference now that initial conditions have been updated"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "1e429e6661e24d43bb04d9e4e43f7732",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=600.0, step=6.0), Output()), _dom_classes=('…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x13eb1c0d0>"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pybamm.dynamic_plot(\n",
+    "    [sol_US06_1, sol_US06_2], labels=[\"Default initial conditions\", \"Fully charged\"]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Initialize using a different model"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also initialize the model using the solution of a different model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spm = pybamm.lithium_ion.SPM()\n",
+    "sim_spm_cccv = pybamm.Simulation(spm, experiment=experiment)\n",
+    "sol_spm_cccv = sim_spm_cccv.solve()\n",
+    "\n",
+    "# MODEL RE-INITIALIZATION: #############################################################\n",
+    "# Now initialize the model with the solution of the charge, and then discharge with\n",
+    "# the US06 drive cycle\n",
+    "# We could also do this inplace by setting inplace to True, which modifies the original\n",
+    "# model in place\n",
+    "new_dfn = model.set_initial_conditions_from(sol_spm_cccv, inplace=False)\n",
+    "########################################################################################\n",
+    "\n",
+    "sim_US06_3 = pybamm.Simulation(\n",
+    "    new_dfn, parameter_values=param, solver=pybamm.CasadiSolver(mode=\"fast\")\n",
+    ")\n",
+    "sol_US06_3 = sim_US06_3.solve()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now the model initialized by the DFN and the model initialized by the SPM give the same solution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d33aa8dd8ac54446968d608c46eb24b9",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=600.0, step=6.0), Output()), _dom_classes=('…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x1407ae940>"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "pybamm.dynamic_plot(\n",
+    "    [sol_US06_1, sol_US06_2, sol_US06_3], \n",
+    "    labels=[\"Default initial conditions\", \"Fully charged (from DFN)\", \"Fully charged (from SPM)\"]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.6"
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": true
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

pybamm/__init__.py

@@ -116,6 +116,7 @@ def version(formatted=False):
 from .expression_tree.operations.jacobian import Jacobian
 from .expression_tree.operations.convert_to_casadi import CasadiConverter
 from .expression_tree.operations.unpack_symbols import SymbolUnpacker
+from .expression_tree.operations.replace_symbols import SymbolReplacer
 
 #
 # Model classes

pybamm/discretisations/discretisation.py

@@ -179,7 +179,7 @@ def process_model(self, model, inplace=True, check_model=True):
             model_disc = model
         else:
             # create an empty copy of the original model
-            model_disc = model.new_copy()
+            model_disc = model.new_empty_copy()
 
         # Keep a record of y_slices in the model
         model_disc.y_slices = self.y_slices_explicit
@@ -216,6 +216,11 @@ def process_model(self, model, inplace=True, check_model=True):
             processed_events.append(processed_event)
         model_disc.events = processed_events
 
+        # Set external variables
+        model_disc.external_variables = [
+            self.process_symbol(var) for var in model.external_variables
+        ]
+
         # Create mass matrix
         pybamm.logger.info("Create mass matrix for {}".format(model.name))
         model_disc.mass_matrix, model_disc.mass_matrix_inv = self.create_mass_matrix(
@@ -1166,7 +1171,7 @@ def check_initial_conditions_rhs(self, model):
             == y0.shape[0]
         ), pybamm.ModelError(
             """
-            Concatenation of (rhs, algebraic) and initial_conditions must have the
+            'Concatenation of (rhs, algebraic)' and 'initial_conditions' must have the
             same shape after discretisation but rhs.shape = {}, algebraic.shape = {},
             and initial_conditions.shape = {}.
             """.format(