Comparing Planet Compositions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing Planet Compositions\n",
    "\n",
    "This notebook imports `cumulative.csv`, a file produced from a grid whose sampling was *uniformly distributed* across a 1-sigma range of refractory element ratios derived from stellar abundances reported in the [Hypatia Catalog Database](https://www.hypatiacatalog.com/hypatia) described in [Hinkel et al. (2014, AJ, 148, 54)](https://doi.org/10.1088/0004-6256/148/3/54). Development of the catalog was supported by NASA's Nexus for Exoplanet System Science (NExSS) research coordination network and the Vanderbilt Initiative in Data-Intensive Astrophysics (VIDA).\n",
    "\n",
    "#### A few quick notes about cumulative.csv, and the model parameters embedded in it:\n",
    "* no Fe allowed in the mantle - only Mg, Ca, Si, Al, and O\n",
    "* stipulated 33% core by mass, core only Fe (so no light elements)\n",
    "* mass percents shown are percents of the *total* planet mass - get.adds_up() normalizes so that only the mantle is considered\n",
    "* 1700K (?) mantle potential temperature assumed - changing T changes mineralogy slightly, but not drastically\n",
    "* model details + a recent implementation: [Unterborn & Panero (2019, JGR: Planets, 124 7)](doi.org/10.1029/2018JE005844)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import external packages and internal modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import external packages\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.express as px\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "\n",
    "#Import internal packages\n",
    "import evolve\n",
    "import plot\n",
    "import fromexo\n",
    "import getall as get\n",
    "from mineralDB import minerals\n",
    "from constants import *\n",
    "\n",
    "# Convenience functions\n",
    "Pe = lambda n: format(n, '.4e')\n",
    "Pf = lambda n: format(n, '.4f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### User-input values, used across all model runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method = 'static'   # static or dynamic calculation of thermal parameters\n",
    "Pref = 5.0          # Reference pressure for Cp and alpha, GPa\n",
    "Ts = 300.0          # surface temperature, K\n",
    "Tp0 = 2000.          # starting mantle potential temperature in K        Earth = 2000.0 (initial), 1600 (present)\n",
    "Qpl = 1.0           # Relative heat production per kg mantle, vs Earth  Earth = 1.0\n",
    "tmax = 4.55         # ending time, in Ga                                Earth = 4.55\n",
    "output_file='compare_compositions_results.csv'\n",
    "output_file_2='bulk_elemental_fraction.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = fromexo.planets_from_summary()\n",
    "print(len(files.keys()), 'planets to analyze, which may contain:')\n",
    "posskeys=[]\n",
    "for f in files:\n",
    "    for i in files[f]:\n",
    "        if i not in posskeys:\n",
    "            posskeys.append(i)\n",
    "print(', '.join(sorted(posskeys)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List what aspects of a planet you'd like shown, w.r.t. thermal evolution pathways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Currently, column values in plottable array are:\n",
    "#            [t, Tp, Ra, production, loss, production / loss, \n",
    "#             Mg, Si, Ca, Al, MgSi, alpha, cp, Water]\n",
    "\n",
    "columnkeys = ['time', 'temp', 'rayleigh', 'production', 'loss', 'urey', \n",
    "             'alpha', 'cp', 'Water', 'Mg', 'Si', 'Ca', 'Al', 'MgSi']\n",
    "\n",
    "ofinterest0 = 'MgSi'\n",
    "ofinterest1 = 'Water'\n",
    "ofinterest2 = 'cp'\n",
    "\n",
    "out = open(output_file, 'w+')\n",
    "out.write('file,Mg/Si,Ca/Si,Al/Si,alpha,Cp,k,temp(K),Rayleigh,HeatLoss(W),UreyRatio,Water\\n')\n",
    "out2 = open(output_file_2, 'w+')\n",
    "out2.write('f_Mg,f_Si,f_Ca,f_Al,alpha,cp,k,tempK,Water\\n')\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate thermal evolution pathways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Hts = []  # A list of lists; column names are in columnkeys\n",
    "sep = ','  # Outputs will be comma-separated value files (*.csv)\n",
    "nplanets = 0\n",
    "print('Progress so far:')\n",
    "for file in files:\n",
    "\n",
    "    Mpl = files[file]['Mass_Me']\n",
    "    Rpl = files[file]['Radius_Re']\n",
    "    planet = {'Mpl': Mpl, 'Rpl': Rpl, 'Qpl': Qpl, 'Tp0': Tp0}\n",
    "    planet['Mp'] = files[file]['Mass_kg']\n",
    "    planet['Mc'] = planet['Mp'] * files[file]['CMF']\n",
    "    planet['Rp'] = files[file]['Radius_m']\n",
    "    planet['Rc'] = planet['Rp'] * files[file]['CRF']\n",
    "    planet['d'] = files[file]['Mantle_depth']\n",
    "    planet['Vm'] = files[file]['Mantle_vol']\n",
    "    planet['Sa'] = 4 * np.pi * planet['Rp']**2\n",
    "    planet['pm'] = files[file]['Mantle_rho']\n",
    "    planet['g'] = get.Grav * planet['Mp']/(planet['Rp']**2)\n",
    "    planet['Pcmb'] = files[file]['CMBP']\n",
    "    planet['Tcmb'] = get.CMB_T(planet['Rp'], planet['Tp0'])\n",
    "    planet['Pref'] = Pref\n",
    "    planet['Qp'] = DEFAULT['Qp']\n",
    "    planet['Ts']=300.0\n",
    "    planet['beta']=0.3\n",
    "\n",
    "\n",
    "    # Get elemental abundances for easy plotting\n",
    "    MgSi = 1 / files[file]['SiMg']\n",
    "    CaSi = files[file]['CaMg'] * MgSi\n",
    "    AlSi = files[file]['AlMg'] * MgSi\n",
    "    totmol = MgSi + CaSi + AlSi + 1\n",
    "    Mg = MgSi / totmol\n",
    "    Ca = CaSi / totmol\n",
    "    Al = AlSi / totmol\n",
    "    Si = 1 / totmol\n",
    "\n",
    "    # Given what your planet's made of, find out its water storage capacity and rheological constraints.\n",
    "    composition = files[file]['composition']\n",
    "    composition = get.adds_up(composition)\n",
    "    planet['c1'], planet['Ev'], planet['visc0'] = DEFAULT['c1'], DEFAULT['Ev'], DEFAULT['visc0']\n",
    "    thermals = get.thermals_at_P_ave(composition, planet['Pref'])\n",
    "    Water = get.average_property(composition, 'water', 0.0)\n",
    "\n",
    "    # Now start the evolutionary clock...\n",
    "    # Be sure to keep Hts, Tp, and t=0.0 here, so at the end of each run, values are reset\n",
    "    # and the clock starts again.\n",
    "    Tp = Tp0\n",
    "    dt = 0.01\n",
    "    t = 0.0\n",
    "    \n",
    "    if method == 'static':\n",
    "        #alpha, cp, k = get.Tdep_thermals(thermals, 1625)\n",
    "        planet['Cp'] = files[file]['Cp']\n",
    "        planet['alpha'] = files[file]['alpha']\n",
    "        planet['k'] = files[file]['k']\n",
    "        \n",
    "    while t <= tmax:\n",
    "\n",
    "        if method == 'dynamic': alpha, cp, k = get.Tdep_thermals(thermals, Tp)\n",
    "\n",
    "        viscT = get.viscosity(planet, Tp)\n",
    "        Ra = get.rayleigh(planet, Tp, Ts, viscT)\n",
    "\n",
    "        production = evolve.produce_heat(planet, t)\n",
    "        loss = evolve.flux_heat(planet, Tp, Ra)\n",
    "        dTp = (dt * get.seconds * (production - loss)) / (planet['Cp'] * planet['pm'] * planet['Vm'])\n",
    "\n",
    "        # This is what columnkey applies to\n",
    "        Hts.append([t, Tp, Ra, production, loss, production / loss, planet['alpha'], planet['Cp'], Water, Mg, Si, Ca, Al, MgSi])\n",
    "\n",
    "        Tp = Tp + dTp\n",
    "        t = t + dt\n",
    "\n",
    "    nplanets = nplanets + 1\n",
    "    \n",
    "    if nplanets%5==0:\n",
    "        percentcomplete=Pf(100*nplanets/len(files.keys()))[:-5]+'%...'\n",
    "        print(percentcomplete, end=\" \")\n",
    "\n",
    "    line = (file, Pf(MgSi), Pf(CaSi), Pf(AlSi), Pe(planet['alpha']), Pf(planet['Cp']), Pf(planet['k']), Pf(Tp), Pe(Ra), Pe(loss), Pf(production/loss), Pe(Water))\n",
    "    out.write(sep.join(line))\n",
    "    out.write('\\n')\n",
    "\n",
    "    line2 = (Pf(Mg), Pf(Si), Pf(Ca), Pf(Al), Pf(Tp), Pe(Water))\n",
    "    out2.write(sep.join(line2))\n",
    "    out2.write('\\n')\n",
    "\n",
    "print()\n",
    "Evolution = pd.DataFrame(np.asarray(Hts), columns=columnkeys)\n",
    "print('Done! Output file names: ' + output_file + '\\t' + output_file_2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can summarize our findings pretty easily..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes = pd.read_csv('compare_compositions_results.csv', sep=',')\n",
    "outcomes.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ...But plots can better show how chemistry seems to affect thermal pathways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "a=plot.evolution_colorcoded(Evolution, colorcolumn=ofinterest0, colortype='continuous')\n",
    "b=plot.evolution_colorcoded(Evolution, colorcolumn=ofinterest1, colortype='continuous')\n",
    "c=plot.evolution_colorcoded(Evolution, colorcolumn=ofinterest2, colortype='continuous')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots can also show how factors correlate over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cortime = px.scatter_matrix(Evolution, \n",
    "    dimensions=['time', 'temp', 'rayleigh', 'loss'],\n",
    "    color='MgSi',\n",
    "    color_discrete_sequence=px.colors.qualitative.Vivid,\n",
    "    title=\"Scatter matrix of planet data set\",\n",
    "    labels=columnkeys)\n",
    "cortime.update_layout(width=1000, height=1000)\n",
    "cortime.update_traces(diagonal_visible=False, showupperhalf=False, marker=dict(size=2))\n",
    "cortime.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ...Or across parameters of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes = pd.read_csv('compare_compositions_results.csv', sep=',')\n",
    "corparams = px.scatter_matrix(outcomes, \n",
    "    dimensions=['Mg/Si', 'temp(K)', 'alpha', 'Cp', 'k', 'Water'],\n",
    "    color=\"temp(K)\",\n",
    "    color_discrete_sequence=px.colors.qualitative.Vivid,\n",
    "    title=\"Scatter matrix of planet data set\",\n",
    "    labels=columnkeys)\n",
    "corparams.update_layout(width=1000, height=1000)\n",
    "corparams.update_traces(diagonal_visible=False, showupperhalf=False, marker=dict(size=8))\n",
    "corparams.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corparallel = px.parallel_coordinates(outcomes, color='Mg/Si', color_continuous_scale='Rainbow', color_continuous_midpoint=np.mean(outcomes['Mg/Si']),\n",
    "                              dimensions=['alpha', 'Cp', 'k', 'Water', 'temp(K)']).update_layout(template='ggplot2+presentation')\n",
    "corparallel\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Want to save one of these interactive figures?\n",
    "Copy, paste, and customize this in a \"code\" cell below (i.e. with figurename changed to indicate which you want, e.g. cortime, corparams, corparallel, etc.):\n",
    "\n",
    "\n",
    "``import plotly.io as pio\n",
    "pio.write_html(figurename, file='myfigure.html', auto_open=False)``"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}