{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 5. Reproducing Previous Results\n",
    "\n",
    "This example recreates the results of\n",
    "[Bryson et al. (2015)](https://www.nature.com/articles/nature14114)\n",
    "giving the same depths of formation of two pallasite meteorites.\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": 2,
   "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",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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 are from and recreate the results\nof Bryson et al. (2015), with explanatory comments:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# These values are quoted in Bryson et al. (2015) or\n# the references therein:\n\n# material properties:\nmantle_diffusivity = 5e-7\nmantle_conductivity_value = 3.0\nmantle_density_value = 3000.0\n\nkappa_reg = 5e-8  # m^2/s\n\ncore_cp = 850.0  # J/(kg K)\ncore_density = 7800.0  # kg/m^3\n\n# geometry:\nr_planet = 200_000.0  # planetesimal radius in m\nreg_m = 8_000.0  # megaregolith thickness in m\n\n# temperatures:\ntemp_core_melting = 1200.0  # K\ntemp_init = 1600.0  # K\ntemp_surface = 250.0  # K\ncore_temp_init = 1600.0  # K\ncore_latent_heat = 270_000.0  # J/kg\n\n# discretisation:\ntimestep = 2E11  # s\ndr = 1000.0  # m\nmax_time = 400  # Myr\n\n# This value isn't explicitly listed in Bryson et al., or references\n# as Bryson et al. (2015) uses diffusivity instead\nmantle_heatcap_value = mantle_conductivity_value / (mantle_density_value * mantle_diffusivity)\n\n# Bryson et al. (2015) list regolith in km as opposed to\n# as a fraction of body radius\nreg_fraction = reg_m / r_planet  # fraction of r_planet\n\n# We don't want to incorporate the 8 km regolith when calculating core size:\n# Bryson et al. (2015) don't seem to include the regolith when calculating\n# the core size, ie the core is 50% of the non-regolith body radius.\n# We don't want to incorporate the 8 km regolith when calculating core size:\ncore_m = (r_planet - reg_m) * 0.5  # 100_000.0 # core size in m\ncore_size_factor = core_m / r_planet  # fraction of r_planet"
   ]
  },
  {
   "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": 5,
   "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 that will\n# be filled as the core freezes\nlatent = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we instantiate the core object. This will keep track of the\ntemperature of the core as the model runs, cooling as heat\nis extracted across the core-mantle boundary. This simple\neutectic core model doesn't track inner core growth, but\nthis is still a required argument to allow for future\nincorporation of more complex core objects:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "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. We want to set these values\nequal to the values used by Bryson et al. (2015):\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mantle_conductivity.setk(mantle_conductivity_value)\nmantle_heatcap.setcp(mantle_heatcap_value)\nmantle_density.setrho(mantle_density_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check that the correct values have been assigned:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0\n",
      "2000.0\n",
      "3000.0\n"
     ]
    }
   ],
   "source": [
    "print(mantle_conductivity.getk())\nprint(mantle_heatcap.getcp())\nprint(mantle_density.getrho())"
   ]
  },
  {
   "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 both the\nsurface and the core-mantle boundary.\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_mantle_bc = pytesimal.numerical_methods.surface_dirichlet_bc\nbottom_mantle_bc = pytesimal.numerical_methods.cmb_dirichlet_bc\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 for the mantle and\ncore.\n\nNow we can use the `pytesimal.analysis` module to find out more\nabout the model run. We can check when the core was freezing,\nso we can compare this to the cooling history of meteorites\nand see whether they can be expected to record magnetic remnants\nof a core dynamo:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(core_frozen,\n times_frozen,\n time_core_frozen,\n fully_frozen) = pytesimal.analysis.core_freezing(core_temperature_array,\n                                                  max_time,\n                                                  times,\n                                                  latent,\n                                                  temp_core_melting,\n                                                  timestep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can calculate arrays of cooling rates from the temperature arrays:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mantle_cooling_rates = pytesimal.analysis.cooling_rate(\n    mantle_temperature_array,\n    timestep)\ncore_cooling_rates = pytesimal.analysis.cooling_rate(core_temperature_array,\n                                                     timestep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Meteorite data (the diameter of 'cloudy-zone particles') can be used to\nestimate the rate at which the meteorites cooled through a specific\ntemperature (C. W. Yang et al., 1997). The\n`analysis.cooling_rate_cloudyzone_diameter` function calculates the cooling\nrate in K/Myr, while the `analysis.cooling_rate_to_seconds` function\nconverts this to K/s which allows comparison to our result arrays.\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d_im = 147  # cz diameter in nm\nd_esq = 158  # cz diameter in nm\n\nimilac_cooling_rate = pytesimal.analysis.cooling_rate_to_seconds(\n    pytesimal.analysis.cooling_rate_cloudyzone_diameter(d_im))\nesquel_cooling_rate = pytesimal.analysis.cooling_rate_to_seconds(\n    pytesimal.analysis.cooling_rate_cloudyzone_diameter(d_esq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this cooling rate information to find out how deep within their\nparent bodies these meteorites originally formed, and when they passed\nthrough the temperature of tetrataenite formation (when magnetism\ncan be recorded). The `analysis.meteorite_depth_and_timing()` function\nreturns the source depth of the meteorite material in the parent body\nbased on the metal cooling rates at 800 K (as a depth from surface in km and\nas a radial value from the center of the planet in m), the time that the\nmeteorite cools through the tetrataenite formation temperature in\ncomparison to the core crystallisation period, and a string defining\nthis relation between paleomagnetism recording and dynamo activity:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(im_depth,\n im_string_result,\n im_time_core_frozen,\n im_Time_of_Crossing,\n im_Critical_Radius) = pytesimal.analysis.meteorite_depth_and_timing(\n    imilac_cooling_rate,\n    mantle_temperature_array,\n    mantle_cooling_rates,\n    radii,\n    r_planet,\n    core_size_factor,\n    time_core_frozen,\n    fully_frozen,\n    dr=1000,\n    dt=timestep\n)\n\n(esq_depth,\n esq_string_result,\n esq_time_core_frozen,\n esq_Time_of_Crossing,\n esq_Critical_Radius) = pytesimal.analysis.meteorite_depth_and_timing(\n    esquel_cooling_rate,\n    mantle_temperature_array,\n    mantle_cooling_rates,\n    radii,\n    r_planet,\n    core_size_factor,\n    time_core_frozen,\n    fully_frozen,\n    dr=1000,\n    dt=timestep\n)\n\nprint(f\"Imilac depth: {im_depth}; Imilac timing: {im_string_result}\")\nprint(f\"Esquel depth: {esq_depth}; Esquel timing: {esq_string_result}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you need to save the meteorite results, they can be saved to a dictionary\nwhich can then be passed to the `load_plot_save.save_params_and_results`.\nThis allows for any number of meteorites to be analysed and only the\nrelevant data stored:\n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "meteorite_results_dict = {'Esq results':\n                              {'depth': esq_depth,\n                               'text result': esq_string_result},\n                          'Im results':\n                              {'depth': im_depth,\n                               'text result': im_string_result,\n                               'critical radius': im_Critical_Radius}}"
   ]
  },
  {
   "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,\n    timestep=timestep)"
   ]
  }
 ],
 "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
}