{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 4. Planetesimal Without a Core\n",
    "\n",
    "This example shows step by step how to set up and run a model of a cooling\n",
    "planetesimal without a core. This example uses constant\n",
    "material properties in the mantle, see\n",
    "[Murphy Quinlan et\n",
    "al. (2021)](https://doi.org/10.1029/2020JE006726) and references therein.\n",
    "\n",
    "While the material properties used in this example are suitable for modelling\n",
    "the olivine mantle of a differentiated body, the model set-up may be useful\n",
    "for modelling undifferentiated meteorite parent bodies with appropriate\n",
    "thermal properties.\n",
    "\n",
    "This model set-up also allows for comparison to an analytical solution for\n",
    "conductive cooling in a sphere (see [Murphy Quinlan et\n",
    "al. (2021)](https://doi.org/10.1029/2020JE006726)).\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": [
    "As we're setting this up step-by-step instead of using the\n",
    "`pytesimal.quick_workflow` module, we need to import a\n",
    "selection of modules:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pytesimal.setup_functions\nimport pytesimal.load_plot_save\nimport pytesimal.numerical_methods\nimport pytesimal.analysis\nimport pytesimal.core_function\nimport pytesimal.mantle_properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of creating and loading a parameter file, we're going to\ndefine variables here. The values match those of the constant\nthermal properties case in Murphy Quinlan et al. (2021),\nexcept for the `core_size_factor` value:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timestep = 1E11  # s\nr_planet = 250_000.0  # m\nreg_fraction = 0.032  # fraction of r_planet\nmax_time = 400  # Myr\ntemp_core_melting = 1200.0  # K\ncore_cp = 850.0  # J/(kg K)\ncore_density = 7800.0  # kg/m^3\ntemp_init = 1600.0  # K\ntemp_surface = 250.0  # K\ncore_temp_init = 1600.0  # K\ncore_latent_heat = 270_000.0  # J/kg\nkappa_reg = 5e-8  # m^2/s\ndr = 1000.0  # m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we want to model an undifferentiated body, we don't want to include\na core in our model set up. Because we wanted it to be easy to swap\nbetween different model set-ups, the `numerical_methods.discretisation()`\nfunction still requires a core object to be instantiated. When the boundary\nconditions are set up correctly, the core object does not interact with\nthe mantle and does not influence the cooling. Setting `core_size_factor=0`\ncan lead to scalar overflow errors, which can be avoided by setting\n`core_size_factor` to a small, non-zero number that is smaller than the grid\nspacing (so the core_temperature_array has dimension 0 along radius):\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "core_size_factor = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `setup_functions.set_up()` function creates empty arrays to\nbe filled with resulting temperatures:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(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)\n\n# We define an empty list of latent heat\nlatent = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that our \"no-core\" set-up has worked:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(core_temperature_array.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we instantiate the core object. As explained previously,\nthis model set up does not include a core, but in order to make\nthe code as modular as possible, a core object still needs to be passed\nin to the main timestepping function, `numerical_methods.discretisation`.\nWe've just copied across the default arguments for convenience:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "core_values = pytesimal.core_function.IsothermalEutecticCore(\n    initial_temperature=core_temp_init,\n    melting_temperature=temp_core_melting,\n    outer_r=r_core,\n    inner_r=0,\n    rho=core_density,\n    cp=core_cp,\n    core_latent_heat=core_latent_heat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define the mantle properties. The default is to have constant\nvalues, so we don't require any arguments for this example:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(mantle_conductivity,\n mantle_heatcap,\n mantle_density) = pytesimal.mantle_properties.set_up_mantle_properties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check (or change) the value of these properties after they've been\nset up using one of the `MantleProperties` methods:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(mantle_conductivity.getk())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If temperature dependent properties are used, temperature can be passed in\nas an argument to return the value at that temperature.\n\nWe need to set up the boundary conditions for the mantle. For this example,\nwe're using fixed temperature boundary conditions at the\nsurface, and a zero-flux condition at the bottom to ensure symmetry.\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_mantle_bc = pytesimal.numerical_methods.surface_dirichlet_bc\nbottom_mantle_bc = pytesimal.numerical_methods.cmb_neumann_bc\n\n# Now we let the temperature inside the planestesimal evolve. This is the\n# slowest part of the code, because it has to iterate over all radii and\n# time.\n# This will take a minute or two!\n# Note that this function call is the exact same as in the default constant\n# and variable cases that included a core.\n\n(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": [
    "This function fills the empty arrays produced by\n`setup_functions.set_up()` with calculated temperatures.\n\nOur next step is to calculate cooling rates for the body. Usually, we\ncalculate cooling rates for the mantle and the core in the same way.\nAs the core temperature array is empty, but we still need a\n`core_cooling_rates` array to pass to our plotting function, we just set\n`core_cooling_rates = core_temperature_array` to save time computing a\nmeaningless empty array:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mantle_cooling_rates = pytesimal.analysis.cooling_rate(mantle_temperature_array, timestep)\ncore_cooling_rates = core_temperature_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an overview of the cooling history of the body, it's very useful\nto plot the temperatures and cooling rates as a heatmap through time.\nIn order to plot the results, we need to define a figure height and width,\nthen call `load_plot_save.plot_temperature_history()`,\n`load_plot_save.plot_coolingrate_history()` or `load_plot_save.two_in_one()`.\nThese functions convert the cooling rate from K/timestep to K/Myr to make\nthe results more human-readable.\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig_w = 6\nfig_h = 9\n\npytesimal.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,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results can be saved in the same way as for the constant and variable\nproperties examples.\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
}