{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 0. Introduction\n",
    "\n",
    "Pytesimal is a finite difference code to perform numerical models of a\n",
    "conductively cooling planetesimal, both with constant and temperature-dependent\n",
    "properties.\n",
    "\n",
    "In this example, we walk through the theoretical background of the model and\n",
    "explain step-by-step how the code works.\n",
    "\n",
    "In order for this to run without installing the package, we put\n",
    "the directory on the PYTHONPATH. This should be removed if you\n",
    "have installed the package and are running this from your\n",
    "machine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Setup PYTHONPATH to allow pytesimal\n",
    "# import without install of package\n",
    "import setup_path"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model set-up\n",
    "\n",
    "Pytesimal allows the modelling of a conductively cooling body, with the\n",
    "following different regions:\n",
    "\n",
    "* An isothermal convecting core: this can be replaced with a more complex core\n",
    "  model or can be switched off to make a core-less body.\n",
    "* A discretised, conductive mantle: this is the region of focus in Pytesimal.\n",
    "  The material properties for this region can be temperature-dependent or\n",
    "  constant.\n",
    "* A discretised megaregolith: this region is also conductively cooling;\n",
    "  material properties can only be constant in this region (constant\n",
    "  diffusivity). This region can be switched off.\n",
    "\n",
    "These different configurations, along with different material properties,\n",
    "can be set up to replicate a wide variety of different planetesimal\n",
    "geometries.\n",
    "\n",
    "In order to set up our model, we first import the required packages:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pytesimal.setup_functions\n",
    "import pytesimal.load_plot_save\n",
    "import pytesimal.numerical_methods\n",
    "import pytesimal.analysis\n",
    "import pytesimal.core_function\n",
    "import pytesimal.mantle_properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of setting up a model run is to use a parameter file.\n",
    "The parameter file is essentially a dictionary holding values\n",
    "for different variables, including the planetesimal radius,\n",
    "core size and regolith thickness, material properties for\n",
    "the body, and values to define the numerical discretisation.\n",
    "\n",
    "We can generate a\n",
    "default parameter file by calling the `make_default_param_file` function.\n",
    "This function can be called with a `filepath` argument, specifying where\n",
    "the file should be saved and what it should be called:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filepath = 'parameters.txt'\n",
    "pytesimal.load_plot_save.make_default_param_file(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This provides a template file for you to edit to suit your specific\n",
    "model set up (see documentation on\n",
    "[parameter files](https://pytesimal.readthedocs.io/en/latest/README.html#parameter-files)\n",
    "for more information on the content and layout of the parameter\n",
    "file). This file can be edited and then loaded in - in practise,\n",
    "you wouldn't do this all in one script like we have here - you\n",
    "would create and edit a parameter file (or copy and edit a\n",
    "pre-existing one), and then in a separate script, would\n",
    "load the parameter file and run the model.\n",
    "\n",
    "As we're just going to use the default values from the\n",
    "parameter file, we'll just load it straight back in without\n",
    "editing it:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(run_ID, folder, timestep, r_planet, core_size_factor,\n",
    " reg_fraction, max_time, temp_core_melting, mantle_heatcap_value,\n",
    " mantle_density_value, mantle_cond_value, core_heatcap, core_density,\n",
    " temp_init, temp_surface, core_temp_init, core_latent_heat,\n",
    " kappa_reg, dr, cond_constant, density_constant,\n",
    " heat_cap_constant) = pytesimal.load_plot_save.load_params_from_file(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This big collection of parameters will be fed in to our model!\n",
    "\n",
    "To make this example run a bit faster, we're going to change the\n",
    "timestep from $1 \\times 10^{11}$ to $2 \\times 10^{11}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timestep = 2e11\n",
    "\n",
    "#\n",
    "# In order to discretise the different regions and record temperatures\n",
    "# at radii points at each timestep, we need to set up some arrays\n",
    "# that match the geometry we've passed in using the parameter file.\n",
    "# These arrays will be placeholders until the numerical method fills\n",
    "# them with values:\n",
    "\n",
    "(r_core,\n",
    " radii,\n",
    " core_radii,\n",
    " reg_thickness,\n",
    " where_regolith,\n",
    " times,\n",
    " mantle_temperature_array,\n",
    " core_temperature_array) = pytesimal.setup_functions.set_up(timestep,\n",
    "                                                            r_planet,\n",
    "                                                            core_size_factor,\n",
    "                                                            reg_fraction,\n",
    "                                                            max_time,\n",
    "                                                            dr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The planetesimal core\n",
    "\n",
    "Before we set up the numerical scheme for the discretised regions of\n",
    "the planetesimal, we need to instantiate the core object.\n",
    "This will keep track of the\n",
    "temperature of the core as the model runs, cooling as heat\n",
    "is extracted across the core-mantle boundary. The heat extracted in one timestep\n",
    "($P_{\\mathrm{CMB}}$) is:\n",
    "\n",
    "$$ P_{\\mathrm{CMB}} = - {A}_{\\mathrm{c}} k_{\\mathrm{m}} \\frac{\\partial T}{\\partial r}\\bigg\\vert _{r = r_\\mathrm{c}} $$\n",
    "\n",
    "where $A_\\mathrm{c}$ is the core surface area, $r_\\mathrm{c}$ is the\n",
    "core radius, and $k_\\mathrm{m}$ is the thermal conductivity at the base of\n",
    "the mantle or discretised region. The corresponding change in core\n",
    "boundary temperature $\\Delta T$ is:\n",
    "\n",
    "$$\\Delta T = - \\frac{P_{\\mathrm{CMB}}}{\\rho_{\\mathrm{c}} C_{\\mathrm{c}} V_{\\mathrm{c}}} \\delta t$$\n",
    "\n",
    "where $\\rho_{\\mathrm{c}}$ and $C_{\\mathrm{c}}$ are the density\n",
    "and heat capacity of the core, and $V_{\\mathrm{c}}$ is the volume of the core.\n",
    "\n",
    "Once the core reaches its freezing temperature, the temperature is pinned.\n",
    "Latent heat is extracted until the total latent heat associated with core\n",
    "crystallisation has been removed. We need to set up an empty list\n",
    "to keep track of this latent heat:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "latent = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simple\n",
    "eutectic core model doesn't track inner core growth, but\n",
    "this is still a required argument to allow for future\n",
    "incorporation of more complex core objects:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "core_values = pytesimal.core_function.IsothermalEutecticCore(initial_temperature=core_temp_init,\n",
    "                                                             melting_temperature=temp_core_melting,\n",
    "                                                             outer_r=r_core, inner_r=0, rho=core_density,\n",
    "                                                             cp=core_heatcap,\n",
    "                                                             core_latent_heat=core_latent_heat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conductive cooling\n",
    "\n",
    "The conductively cooling regions in the planetesimal can be described in\n",
    "1D by the heat equation in spherical geometry:\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t} \\rho C=\n",
    "    \\frac{1}{r^2} \\frac{\\partial}{\\partial r}\\left(k r^2 \\frac{\\partial T}{\\partial r} \\right),$$\n",
    "\n",
    "where $T$ is temperature, $t$ is time, $\\rho$ is density,\n",
    "$C$ is heat capacity, $k$ is thermal conductivity, and\n",
    "$r$ is radius.\n",
    "\n",
    "We need to define the thermal properties for this region\n",
    "($\\rho$, $C$, and $k$). In our parameter\n",
    "file, we've already defined these as constant in temperature\n",
    "and have listed values. We just need to pass those arguments\n",
    "to the `set_up_mantle_properties` function:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(mantle_conductivity, mantle_heatcap, mantle_density) = pytesimal.mantle_properties.set_up_mantle_properties(\n",
    "    cond_constant,\n",
    "    density_constant,\n",
    "    heat_cap_constant,\n",
    "    mantle_density_value,\n",
    "    mantle_heatcap_value,\n",
    "    mantle_cond_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conduction equation is solved numerically using an explicit finite\n",
    "difference scheme, Forward-Time Central-Space (FTCS). FTCS gives first-order\n",
    "convergence in time and second-order in space, and is conditionally stable when\n",
    "applied to the heat equation. We can quickly calculate the diffusivity in the\n",
    "mantle and then check we meet Von Neumann stability criteria for the\n",
    "mantle and the megaregolith:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mantle_diffusivity = pytesimal.numerical_methods.calculate_diffusivity(mantle_cond_value, mantle_heatcap_value,\n",
    "                                                                       mantle_density_value)\n",
    "\n",
    "mantle_stability = pytesimal.numerical_methods.check_stability(mantle_diffusivity, timestep, dr)\n",
    "\n",
    "reg_stability = pytesimal.numerical_methods.check_stability(kappa_reg, timestep, dr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the boundary conditions for the mantle. For this example,\n",
    "we're using fixed temperature boundary conditions at both the\n",
    "surface and the core-mantle boundary: at the planetesimal's surface,\n",
    "the temperature is held at a fixed temperature specified in the\n",
    "parameter file (250 K), while at the core-mantle boundary, the\n",
    "temperature is updated by the core.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_mantle_bc = pytesimal.numerical_methods.surface_dirichlet_bc\n",
    "bottom_mantle_bc = pytesimal.numerical_methods.cmb_dirichlet_bc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we pass our boundary conditions, core object, initial\n",
    "temperature, material properties, geometry, and arrays to\n",
    "the `discretisation` function, which returns arrays of\n",
    "temperatures in the mantle and the core, and a list of latent\n",
    "heat values during core crystallisation:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(mantle_temperature_array,\n",
    " core_temperature_array,\n",
    " latent,\n",
    " ) = pytesimal.numerical_methods.discretisation(\n",
    "    core_values,\n",
    "    latent,\n",
    "    temp_init,\n",
    "    core_temp_init,\n",
    "    top_mantle_bc,\n",
    "    bottom_mantle_bc,\n",
    "    temp_surface,\n",
    "    mantle_temperature_array,\n",
    "    dr,\n",
    "    core_temperature_array,\n",
    "    timestep,\n",
    "    r_core,\n",
    "    radii,\n",
    "    times,\n",
    "    where_regolith,\n",
    "    kappa_reg,\n",
    "    mantle_conductivity,\n",
    "    mantle_heatcap,\n",
    "    mantle_density)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing results\n",
    "\n",
    "We can calculate the cooling rates in the planetesimal\n",
    "with the `analysis` module:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mantle_cooling_rates = pytesimal.analysis.cooling_rate(mantle_temperature_array, timestep)\n",
    "core_cooling_rates = pytesimal.analysis.cooling_rate(core_temperature_array, timestep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then plot the results using the `load_plot_save`\n",
    "module. The `two_in_one` function can be called to quickly\n",
    "plot both the temperature and cooling rates:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig_w = 6\n",
    "fig_h = 9\n",
    "\n",
    "pytesimal.load_plot_save.two_in_one(\n",
    "    fig_w,\n",
    "    fig_h,\n",
    "    mantle_temperature_array,\n",
    "    core_temperature_array,\n",
    "    mantle_cooling_rates,\n",
    "    core_cooling_rates,\n",
    "    timestep=timestep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further information\n",
    "\n",
    "Other examples are available in the\n",
    "[gallery](https://pytesimal.readthedocs.io/en/latest/examples/index.html),\n",
    "and further information on the theoretical background can be found in\n",
    "[Murphy Quinlan et\n",
    "al. (2021)](https://doi.org/10.1029/2020JE006726).\n",
    "\n"
   ]
  }
 ],
 "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}