{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 3. Variable Properties\n", "\n", "This example shows step by step how to set up and run a model of a cooling\n", "planetesimal without using a parameters file. This example uses temperature\n", "dependent material properties in the mantle and reproduces the results of\n", "[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", "\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):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timestep = 1E11 # s\nr_planet = 250_000.0 # m\ncore_size_factor = 0.5 # fraction of r_planet\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": [ "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 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": 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. We want to use temperature-dependent\nproperties, so we need to specify this:\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(\n cond_constant='n',\n density_constant='n',\n heat_cap_constant='n'\n)" ] }, { "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(f' Mantle conductivity {mantle_conductivity.getk()} W/(m K)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When temperature dependent properties are used, temperature can be passed in\nas an argument to return the value at that temperature. The default\ntemperature is 295 K, so if no temperature argument is passed, the\nfunction is evaluated at `T=295`.\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": null, "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# 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# The mantle property objects are passed in in the same way as in the\n# example with constant thermal properties.\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(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)\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)\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, )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a few formats or ways to save the results. The temperature and\ncooling rate arrays can be saved as compressed `.npz` arrays, to be loaded\nat a later time and replotted/new meteorite formation depths calculated etc.\nThe input parameters can be saved as a `.json` file, which allows the run to\nbe documented and provides all the metadata needed to reproduce the results.\nFor either of these, a results folder and results filename is needed. The\nfolder can be defined relative to the working directory, or with an absolute\npath. An absolute path usually results in less confusion!\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# define folder and check it exists:\nfolder = 'workflow'\npytesimal.load_plot_save.check_folder_exists(folder)\n# define a results filename prefix:\nresult_filename = 'variable_workflow_results'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result arrays can now be saved:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pytesimal.load_plot_save.save_result_arrays(result_filename,\n folder,\n mantle_temperature_array,\n core_temperature_array,\n mantle_cooling_rates,\n core_cooling_rates)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to save the result parameter file, we also need to define a\n`run_ID`, a descriptive string to identify the model run, and clarify\nwhether we used constant or variable thermal properties:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "run_ID = 'Example run with default variable properties'\ncond_constant = 'n'\ndensity_constant = 'n'\nheat_cap_constant = 'n'\n\npytesimal.load_plot_save.save_params_and_results(\n result_filename, run_ID, folder, timestep, r_planet, core_size_factor,\n reg_fraction, max_time, temp_core_melting, mantle_heatcap.getcp(),\n mantle_density.getrho(), mantle_conductivity.getk(), core_cp, 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, time_core_frozen, fully_frozen,\n meteorite_results=meteorite_results_dict,\n latent_list_len=len(latent))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This results file can then be loaded as a parameter file if you want to\nrepeat the same set up.\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 }