{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exchange-current density from resistance\n",
    "\n",
    "In this example, we show how to fit the exchange-current density to resistance data calculated from synthetic pulse data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pybamm\n",
    "import ionworkspipeline as iwp\n",
    "import pandas as pd"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by generating synthetic pulse data for a half-cell. We use the standard symmetric exchange-current density function and will fit the reference value. We set up the model and parameters,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.lithium_ion.SPMe(\n",
    "    {\"working electrode\": \"positive\", \"particle\": \"uniform profile\"}\n",
    ")\n",
    "model.events = []\n",
    "parameter_values = iwp.ParameterValues(\"Xu2019\")\n",
    "\n",
    "\n",
    "def j0(c_e, c_s_surf, c_s_max, T):\n",
    "    j0_ref = pybamm.Parameter(\n",
    "        \"Positive electrode reference exchange-current density [A.m-2]\"\n",
    "    )\n",
    "    c_e_init = pybamm.Parameter(\"Initial concentration in electrolyte [mol.m-3]\")\n",
    "\n",
    "    return (\n",
    "        j0_ref\n",
    "        * (c_e / c_e_init) ** 0.5\n",
    "        * (c_s_surf / c_s_max) ** 0.5\n",
    "        * (1 - c_s_surf / c_s_max) ** 0.5\n",
    "    )\n",
    "\n",
    "\n",
    "parameter_values.update(\n",
    "    {\n",
    "        \"Positive electrode reference exchange-current density [A.m-2]\": 5,\n",
    "        \"Positive electrode exchange-current density [A.m-2]\": j0,\n",
    "    },\n",
    "    check_already_exists=False,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then simulate a GITT experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_rate = 1 / 15\n",
    "pause_s = 10  # 10s\n",
    "step_s = 5 * 60  # 5m\n",
    "step_h = step_s / 3600\n",
    "rest_s = 1 * 3600  # 1h\n",
    "max_cycles = int(1 / C_rate / step_h * 2)\n",
    "V_min = parameter_values[\"Lower voltage cut-off [V]\"]\n",
    "\n",
    "gitt_experiment = pybamm.Experiment(\n",
    "    [\n",
    "        (\n",
    "            pybamm.step.rest(pause_s, period=1),\n",
    "            pybamm.step.c_rate(\n",
    "                C_rate,\n",
    "                duration=step_s,\n",
    "                period=step_s / 100,\n",
    "                termination=f\"{V_min} V\",\n",
    "            ),\n",
    "            pybamm.step.rest(rest_s, period=rest_s / 100),\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "sim = iwp.Simulation(\n",
    "    model, parameter_values=parameter_values, experiment=gitt_experiment\n",
    ")\n",
    "sol = None\n",
    "i = 0\n",
    "\n",
    "# Run until the voltage drops below the lower cut-off\n",
    "# during step 1 (the discharge step)\n",
    "while sol is None or sol.cycles[-1].steps[1].termination == \"final time\":\n",
    "    sol = sim.solve(starting_solution=sol)\n",
    "    i += 1\n",
    "    if i > max_cycles:\n",
    "        raise ValueError(\"Reached maximum number of cycles\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We export the solution data to a dataframe. \n",
    "\n",
    "> Note: for synthetic data we can directly store the stoichiometry, whereas in practice we would need to compute this from the half-cell capacity and the charge passed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = [\n",
    "    \"Time [s]\",\n",
    "    \"Current [A]\",\n",
    "    \"Voltage [V]\",\n",
    "    \"Positive electrode stoichiometry\",\n",
    "]\n",
    "df = pd.DataFrame(sim.solution.get_data_dict(variables))\n",
    "df = df.rename(columns={\"Cycle\": \"Cycle number\", \"Step\": \"Step number\"})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we calculate the 1s resistance from the solution data using the function `iwp.objectives.pulse.calculate_pulse_resistance`. We pass in the step number \"1\" which is the \"on\" step of our pulse data (steps 0 and 2 are the \"off\" steps). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = [1]\n",
    "dts = [1]\n",
    "data = iwp.objectives.pulse.calculate_pulse_resistance(\"positive\", df, steps, dts)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then plot the calculated resistance as a function of stoichiometry. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = data.plot(x=\"SOC\", y=\"Resistance [Ohm]\", style=\"o-\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we set up the objective that will be used to perform the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = iwp.objectives.HalfCellResistanceModel(\"positive\")\n",
    "objective = iwp.objectives.Resistance(data, {\"model\": model})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We fit the exchange-current density and an extra scalar \"contact resistance\" parameter that accounts for any other resistances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = {\n",
    "    \"Positive electrode reference exchange-current density [A.m-2]\": iwp.Parameter(\n",
    "        \"j0_ref\", initial_value=1, bounds=(0.1, 10)\n",
    "    ),\n",
    "    \"Contact resistance [Ohm]\": iwp.Parameter(\"R_contact\"),\n",
    "}\n",
    "data_fit = iwp.DataFit(objective, parameters=parameters)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We pass in the known parameter values and run the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure we're not accidentally initializing with the correct values by passing\n",
    "# them in\n",
    "params_for_pipeline = {k: v for k, v in parameter_values.items() if k not in parameters}\n",
    "data_fit.run(params_for_pipeline)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = data_fit.plot_fit_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "6b12f1f625627d6c4a17ef696ddbbbd9bd4b12881121180de40e09e7956aa05c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}