{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Create relaxed geodynamic 1D profile\n", "\n", "`NOTE: This notebook contains an interactive figure with sliders. it relies on python modules ipympl and mpl_interactions. `\n", "`If you want to run this notebook with static figures, restart the kernel, clear any history and set interactive = False in the first code block.`\n", "\n", "In the mantle, it is common to assume that convecting material is at chemical equilibrium; all of the reactions between phases keep pace with the changes in pressure and temperature. Because of this relaxation, physical properties such as heat capacity $C_P$, thermal expansion $\\alpha$ and compressibility $\\beta$ must be computed by numerical differentiation of the entropy $\\mathcal{S}$ and volume $\\mathcal{V}$. It is these values, rather than the unrelaxed values output as standard by BurnMan and PerpleX which should be used in geodynamic simulations.\n", "\n", "Relaxed properties can sometimes be very different from their unrelaxed counterparts. Take, for example, the univariant reaction forsterite -> Mg-wadsleyite. These transformation involves a step change in volume, and thus the relaxed compressibility at the transition is infinite. Obviously, if geodynamics software uses compressibility as an input parameter, then whichever meshing is chosen, it will completely miss the transition. There are two solutions to this problem:\n", "* Calculate the entropy and volume at the quadrature points, and calculate $\\nabla\\mathcal{S}$ and $\\nabla\\mathcal{V}$ within each cell. This method is computationally expensive and there may be convergence problems if the quadrature points are very close to the positions of near-univariant reactions.\n", "* Smooth $\\mathcal{S}(P, T)$ and $\\mathcal{V}(P, T)$ by convolution with a 2D Gaussian (in $P$ and $T$) before calculating $C_P$, $\\alpha$ and $\\beta$. A good rule of thumb is that reactions should span about 4 cells for the latent heat to be captured within a few percent.\n", "\n", "The second method is used here to create 1D material property profiles which can be directly used by $ASPECT$. The user of this notebook can vary important mineral physics parameters (rock type, potential temperature, surface gravity) and smoothing parameters (Gaussian widths)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's install some modules that we need to run interactive figures." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interactive = True\n", "if interactive:\n", " %matplotlib ipympl\n", " import ipywidgets as widgets\n", " import mpl_interactions.ipyplot as iplt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now a few imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import burnman\n", "from burnman import Layer\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import fsolve, brentq\n", "from scipy.integrate import odeint\n", "from scipy.interpolate import UnivariateSpline\n", "plt.style.use('bmh')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code block sets up some parameters for our problem,\n", "including the PerpleX model we will use to determine physical properties of the planet." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "perplex_filename = '../../burnman/data/input_perplex/in23_1.tab' # 'example23_hires.tab' # '../../burnman/data/input_perplex/in23_1.tab'\n", "potential_temperature = 1550.\n", "\n", "outer_radius = 6371.e3\n", "thickness = 550.e3\n", "n_points = 251\n", "\n", "pressure_top = 1.e5\n", "gravity_bottom = 10.\n", "\n", "depths = np.linspace(thickness, 0., n_points)\n", "\n", "\n", "rock = burnman.PerplexMaterial(perplex_filename)\n", "\n", "layer = Layer(name='Mantle', radii=outer_radius-depths)\n", "layer.set_material(rock) \n", "layer.set_temperature_mode(temperature_mode='adiabatic',\n", " temperature_top=1550.)\n", "layer.set_pressure_mode(pressure_mode='self-consistent',\n", " pressure_top=1.e5,\n", " gravity_bottom=gravity_bottom)\n", "layer.make()\n", "\n", "truncate = 4. # truncates the convolution Gaussian at 4 sigma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code block reinterpolates the entropy and volume onto a grid for smoothing purposes.\n", "As we're using a PerpleX table, we could just have read the original grid in directly, but this way, we could more easily adapt the code to use other BurnMan materials.\n", "This step takes a few seconds (10--30 s on a fast ultrabook), but we only do it once." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_gridpoints = (501, 101)\n", "\n", "min_grid_pressure = rock.bounds[0][0]\n", "max_grid_pressure = rock.bounds[0][1]\n", "min_grid_temperature = rock.bounds[1][0]\n", "max_grid_temperature = rock.bounds[1][1]\n", "\n", "grid_pressures = np.linspace(min_grid_pressure, max_grid_pressure, n_gridpoints[0])\n", "grid_temperatures = np.linspace(min_grid_temperature, max_grid_temperature, n_gridpoints[1])\n", "pp, TT = np.meshgrid(grid_pressures, grid_temperatures)\n", "mesh_shape = pp.shape\n", "pp = np.ndarray.flatten(pp)\n", "TT = np.ndarray.flatten(TT)\n", "\n", "grid_entropies = np.zeros_like(pp)\n", "grid_volumes = np.zeros_like(pp)\n", "\n", "grid_entropies, grid_volumes = layer.material.evaluate(['S', 'V'], pp, TT)\n", "\n", "grid_entropies = grid_entropies.reshape(mesh_shape)\n", "grid_volumes = grid_volumes.reshape(mesh_shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This long code block sets up three functions, which:\n", "- return the temperatures along an isentrope given a 2D S(P,T) interpolator.\n", "- return the relaxed geodynamic properties along the adiabat\n", "- save a file containing the relaxed properties" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define function to find an isentrope given a\n", "# 2D entropy interpolation function\n", "# Here we use fsolve, because we'll normally have a good starting guess\n", "# from the previous pressure\n", "def interp_isentrope(interp, pressures, entropies, T_guess):\n", " def _deltaS(args, S, P):\n", " T = args[0]\n", " return interp(P, T)[0] - S\n", " \n", " sol = [T_guess]\n", " temperatures = np.empty_like(pressures)\n", " for i in range(len(pressures)):\n", " sol = fsolve(_deltaS, sol, args=(entropies[i], pressures[i]))\n", " temperatures[i] = sol[0]\n", "\n", " return temperatures\n", "\n", "\n", "def relaxed_profile(layer, pressure_stdev, temperature_stdev,\n", " truncate):\n", "\n", " unsmoothed_T_spline = UnivariateSpline(layer.pressure[::-1], layer.temperature[::-1])\n", " unsmoothed_grid_isentrope_temperatures = unsmoothed_T_spline(grid_pressures)\n", " \n", " # Having defined the grid and calculated unsmoothed properties,\n", " # we now calculate the smoothed entropy and volume and derivatives with\n", " # respect to pressure and temperature.\n", " S_interps = burnman.tools.math.interp_smoothed_array_and_derivatives(array=grid_entropies,\n", " x_values=grid_pressures,\n", " y_values=grid_temperatures,\n", " x_stdev=pressure_stdev,\n", " y_stdev=temperature_stdev,\n", " truncate=truncate)\n", " interp_smoothed_S, interp_smoothed_dSdP, interp_smoothed_dSdT = S_interps\n", " \n", " V_interps = burnman.tools.math.interp_smoothed_array_and_derivatives(array=grid_volumes,\n", " x_values=grid_pressures,\n", " y_values=grid_temperatures,\n", " x_stdev=pressure_stdev,\n", " y_stdev=temperature_stdev,\n", " truncate=truncate)\n", " \n", " interp_smoothed_V, interp_smoothed_dVdP, interp_smoothed_dVdT = V_interps\n", " \n", " # Now we can calculate and plot the relaxed and smoothed properties along the isentrope \n", " smoothed_temperatures = interp_isentrope(interp_smoothed_S, layer.pressure[::-1], layer.S[::-1], layer.temperature[-1])[::-1]\n", " densities = layer.material.evaluate(['rho'], layer.pressure, smoothed_temperatures)[0]\n", " \n", " volumes = np.array([interp_smoothed_V(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])\n", " dSdT = np.array([interp_smoothed_dSdT(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])\n", " dVdT = np.array([interp_smoothed_dVdT(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])\n", " dVdP = np.array([interp_smoothed_dVdP(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])\n", " \n", " alphas_relaxed = dVdT / volumes\n", " compressibilities_relaxed = -dVdP / volumes\n", " specific_heats_relaxed = smoothed_temperatures * dSdT / (densities[0]*volumes[0])\n", " \n", " \n", " dT = 0.1\n", " Vpsub, Vssub = layer.material.evaluate(['p_wave_velocity', 'shear_wave_velocity'],\n", " layer.pressure, smoothed_temperatures-dT/2.)\n", " Vpadd, Vsadd = layer.material.evaluate(['p_wave_velocity', 'shear_wave_velocity'],\n", " layer.pressure, smoothed_temperatures+dT/2.)\n", "\n", " Vps = (Vpadd + Vpsub)/2.\n", " Vss = (Vsadd + Vssub)/2.\n", " dVpdT = (Vpadd - Vpsub)/dT\n", " dVsdT = (Vsadd - Vssub)/dT\n", " \n", " depths = layer.outer_radius - layer.radii\n", " return (smoothed_temperatures, layer.pressure, depths, layer.gravity, densities,\n", " alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed,\n", " Vss, Vps, dVsdT, dVpdT)\n", "\n", "flag_index = {'T': 0, 'P': 1, 'z': 2, 'g': 3, 'rho': 4,\n", " 'alpha': 5, 'beta_T': 6, 'Cp': 7,\n", " 'Vs': 8, 'Vp': 9, 'dVsdT': 10,\n", " 'dVpdT': 11}\n", "\n", "def save_relaxed_properties(layer, P_GPa_gaussian, T_K_gaussian, outfile='isentrope_properties.txt'):\n", " \"\"\"\n", " A function to output smoothed, relaxed properties for use in ASPECT\n", " depth, pressure, temperature, density, gravity, Cp (per kilo), thermal expansivity\n", " \"\"\"\n", " d = relaxed_profile(layer, P_GPa_gaussian*1.e9, T_K_gaussian, truncate)[::-1]\n", "\n", " np.savetxt(outfile, X=np.array([d[flag_index['z']],\n", " d[flag_index['P']],\n", " d[flag_index['T']],\n", " d[flag_index['rho']],\n", " d[flag_index['g']],\n", " d[flag_index['alpha']],\n", " d[flag_index['Cp']],\n", " d[flag_index['beta_T']],\n", " d[flag_index['Vs']],\n", " d[flag_index['Vp']],\n", " d[flag_index['dVsdT']],\n", " d[flag_index['dVpdT']]]).T,\n", " header=('# This ASPECT-compatible file contains material '\n", " 'properties calculated along an isentrope by the '\n", " f'BurnMan software.\\n# POINTS: {n_points}\\n'\n", " '# depth (m), pressure (Pa), temperature (K), '\n", " 'density (kg/m^3), gravity (m/s^2), '\n", " 'thermal expansivity (1/K), specific heat (J/K/kg), '\n", " 'compressibility (1/Pa), seismic Vs (m/s), '\n", " 'seismic Vp (m/s), seismic dVs/dT (m/s/K), '\n", " 'seismic dVp/dT (m/s/K)\\n'\n", " 'depth pressure '\n", " 'temperature density '\n", " 'gravity thermal_expansivity '\n", " 'specific_heat compressibility \t'\n", " 'seismic_Vs seismic_Vp '\n", " 'seismic_dVs_dT seismic_dVp_dT'),\n", " fmt='%.10e', delimiter='\\t', comments='')\n", "\n", " print('File saved to {0}'.format(outfile))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Save file" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "save_relaxed_properties(layer, P_GPa_gaussian=0.25, T_K_gaussian=0.25)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following code block attempts to make updating the multipart interactive figure as efficient as possible.\n", "It does this by calculating all of the properties in an inner function, and storing them in a global parameter (global_stored_properties). That function is then wrapped in an outer function that dictates which properties to return. If the inner function has been previously called with the same input parameters, the previously stored properties are returned, otherwise, all properties are calculated from the new input parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "global global_PT_smooth\n", "global_PT_smooth = [None, None]\n", "global_stored_properties = [None]\n", "\n", "def plot_y(flag):\n", " index = flag_index[flag]\n", " def f(x, P_GPa_gaussian, T_K_gaussian):\n", " if P_GPa_gaussian == global_PT_smooth[0] and T_K_gaussian == global_PT_smooth[1]:\n", " pass\n", " else:\n", " global_PT_smooth[0] = P_GPa_gaussian\n", " global_PT_smooth[1] = T_K_gaussian\n", " f = relaxed_profile(layer, P_GPa_gaussian*1.e9, T_K_gaussian,\n", " truncate)\n", " global_stored_properties[0] = f\n", " return global_stored_properties[0][index]\n", " return f" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = 8, 5 # inches\n", "fig = plt.figure()\n", "\n", "px, py = [2, 3]\n", "\n", "depths = layer.outer_radius - layer.radii\n", "gravity = layer.gravity\n", "x = depths/1.e3\n", "xlabel = 'Depths (km)'\n", "\n", "ax_T = fig.add_subplot(px, py, 1)\n", "ax_T.plot(x, layer.temperatures, label='unrelaxed')\n", "ax_T.set_ylabel('Temperature (K)')\n", "ax_T.set_xlabel(xlabel)\n", "\n", "ax_g = fig.add_subplot(px, py, 2)\n", "ax_g.plot(x, layer.gravity)\n", "ax_g.set_ylabel('Gravity (m/s^2)')\n", "ax_g.set_xlabel(xlabel)\n", "\n", "ax_rho = fig.add_subplot(px, py, 3)\n", "ax_rho.plot(x, layer.rho, label='$\\rho$ (kg/m$^3$)')\n", "ax_rho.plot(x, layer.v_p, label='P (km/s)')\n", "ax_rho.plot(x, layer.v_s, label='S (km/s)')\n", "ax_rho.set_ylabel('Densities/Velocities')\n", "ax_rho.set_xlabel(xlabel) \n", "\n", "ax_alpha = fig.add_subplot(px, py, 4)\n", "ax_alpha.plot(x, layer.alpha)\n", "ax_alpha.set_ylabel('alpha (/K)')\n", "ax_alpha.set_xlabel(xlabel)\n", "\n", "ax_beta = fig.add_subplot(px, py, 5)\n", "ax_beta.plot(x, layer.beta_T)\n", "ax_beta.set_ylabel('compressibilities (/Pa)')\n", "ax_beta.set_xlabel(xlabel)\n", "\n", "ax_cp = fig.add_subplot(px, py, 6)\n", "ax_cp.plot(x, layer.C_p/layer.molar_mass)\n", "ax_cp.set_ylabel('Cp (J/K/kg)')\n", "ax_cp.set_xlabel(xlabel) \n", "\n", "\n", "# Relaxed, unsmoothed properties\n", "ax_T.plot(x, plot_y('T')(x, 0., 0.), label='relaxed, unsmoothed')\n", "ax_g.plot(x, plot_y('g')(x, 0., 0.))\n", "ax_alpha.plot(x, plot_y('alpha')(x, 0., 0.))\n", "ax_beta.plot(x, plot_y('beta_T')(x, 0., 0.))\n", "ax_cp.plot(x, plot_y('Cp')(x, 0., 0.))\n", "\n", "if interactive:\n", " # Interactive smoothing\n", " P_GPa_gaussian = np.linspace(0., 3., 51)\n", " T_K_gaussian = np.linspace(0., 30, 41)\n", " controls = iplt.plot(x, plot_y('T'), P_GPa_gaussian=P_GPa_gaussian, T_K_gaussian=T_K_gaussian, ax=ax_T, label='relaxed, smoothed')\n", " _ = iplt.plot(x, plot_y('g'), controls=controls, ax=ax_g)\n", " _ = iplt.plot(x, plot_y('alpha'), controls=controls, ax=ax_alpha)\n", " _ = iplt.plot(x, plot_y('beta_T'), controls=controls, ax=ax_beta)\n", " _ = iplt.plot(x, plot_y('Cp'), controls=controls, ax=ax_cp)\n", "else:\n", " # Non-interactive smoothing\n", " P_GPa_gaussian = 0.5\n", " T_K_gaussian = 20.\n", " ax_T.plot(x, plot_y('T')(x, P_GPa_gaussian, T_K_gaussian), label='relaxed, smoothed')\n", " ax_g.plot(x, plot_y('g')(x, P_GPa_gaussian, T_K_gaussian))\n", " ax_alpha.plot(x, plot_y('alpha')(x, P_GPa_gaussian, T_K_gaussian))\n", " ax_beta.plot(x, plot_y('beta_T')(x, P_GPa_gaussian, T_K_gaussian))\n", " ax_cp.plot(x, plot_y('Cp')(x, P_GPa_gaussian, T_K_gaussian))\n", "\n", "ax_T.legend(loc='lower right',prop={'size':8})\n", "\n", "fig.set_tight_layout(True)" ] } ], "metadata": { "interpreter": { "hash": "d0e8ff0504fa29d441371c1f42f91c694e01f5e6e44698edccbd59f7213ffa15" }, "kernelspec": { "display_name": "Python 3.9.2 64-bit ('3.9.2': pyenv)", "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.9.2" }, "widgets": { "state": { "00cd58f4139c44c0afd3e5a1a8dc0ce6": { "views": [ { "cell_index": 4 } ] }, "1f5a9835b66a4a30bf5979c71220aaf5": { "views": [ { "cell_index": 5 } ] }, "293e629328764925954b515887354db7": { "views": [ { "cell_index": 3 } ] }, "2af4495b5f744137929783da0744eb4e": { "views": [ { "cell_index": 4 } ] }, "2f057a0d6f9e4f0a99d57555c2b69ad9": { "views": [ { "cell_index": 3 } ] }, "352845fdb1154917819f0fb105b827dc": { "views": [ { "cell_index": 3 } ] }, "3719d9e58af84c84964a987452c89143": { "views": [ { "cell_index": 3 } ] }, "4e7a2298c568448986e6903e7007c3ad": { "views": [ { "cell_index": 3 } ] }, "51a2351b90bc4a65b353c78dcb5b1395": { "views": [ { "cell_index": 3 } ] }, "56a7a2d4b0ce4ea7b0b69d9bf9f9b902": { "views": [ { "cell_index": 4 } ] }, "5de34cd7fe074af59b2a8e2c4990a602": { "views": [ { "cell_index": 3 } ] }, "61ca870e057b4df5b9911ac4f72c52a9": { "views": [ { "cell_index": 4 } ] }, "69653fb1aeb94ade8c7bc426d29cd1de": { "views": [ { "cell_index": 4 } ] }, "69894023bec54c3aa354c2a7e97a4b81": { "views": [ { "cell_index": 4 } ] }, "7b2df34d13664355951d3d233ac69757": { "views": [ { "cell_index": 4 } ] }, "7f10aa737d3241c9a70fc0222f75811d": { "views": [ { "cell_index": 3 } ] }, "8a0a690b97984f1da451fdf19204299a": { "views": [ { "cell_index": 4 } ] }, "8ce9093cfadb49db8d8ac1bee156c100": { "views": [ { "cell_index": 3 } ] }, "91e0694f4d164ead903a6c5a0024b973": { "views": [ { "cell_index": 3 } ] }, "943f263d581741aa92ceec77ddcb132d": { "views": [ { "cell_index": 4 } ] }, "a7e0c22d3c434e4eb5b4a9c069689a5e": { "views": [ { "cell_index": 3 } ] }, "c9e764c73b844effb3bbf512d47bf2c6": { "views": [ { "cell_index": 3 } ] }, "ccdd1a84b1f14af1832e23e9b1b1dc1f": { "views": [ { "cell_index": 6 } ] }, "d19a7f4417cb4611bff97e7d00c41802": { "views": [ { "cell_index": 4 } ] }, "df2666a8acc746229695a8d29d24df10": { "views": [ { "cell_index": 4 } ] }, "e1ed4f87995a49c7bef065bbf232e806": { "views": [ { "cell_index": 3 } ] }, "f2a1e0a22db04a3eacb620ddc63d8773": { "views": [ { "cell_index": 4 } ] } }, "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 2 }