In [ ]:
%matplotlib notebook
#%matplotlib inline
import ipywidgets as widgets
from IPython.display import display
from IPython.display import clear_output
from IPython.display import FileLink, FileLinks
import burnman
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, brentq
from scipy.integrate import odeint
from scipy.interpolate import UnivariateSpline
plt.style.use('bmh')
NOTE: please wait until this page loads completely
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.
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:
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).
In [ ]:
default_file = '../../burnman/data/input_perplex/in23_1.tab' # 'example23_hires.tab' # '../../burnman/data/input_perplex/in23_1.tab'
rock_path_widget = widgets.Text(value=default_file,
placeholder='Type the *local* path here',
description='Path to PerpleX tab file:')
potlT_widget = widgets.FloatText(value=1550,
step=1, description='Potential temperature (K):')
maxP_widget = widgets.FloatText(value=25, description='Max. pressure (GPa):')
g0_widget = widgets.FloatText(value=9.81, step=0.01, description='Surface gravity (m/s$^{2}$):')
r0_widget = widgets.FloatText(value=6371., description='Outer radius (km):')
n_points_widget = widgets.BoundedIntText(value=251, min=51, max=501, description='number of profile points:')
n_P_gridpoints_widget = widgets.BoundedIntText(value=501, min=51, max=2001, description='n. P gridpoints:')
n_T_gridpoints_widget = widgets.BoundedIntText(value=101, min=11, max=1001, description='n. T gridpoints:')
max_T_gaussian_widget = widgets.FloatText(value=30., description='Max. T sigma (K):')
truncate_widget = widgets.BoundedFloatText(value=4, min=3, max=6, description='Sigma truncation')
In [ ]:
w1=widgets.VBox([rock_path_widget, potlT_widget, maxP_widget, g0_widget, r0_widget], layout=widgets.Layout(width='50%'))
w2=widgets.VBox([n_points_widget, n_P_gridpoints_widget, n_T_gridpoints_widget, max_T_gaussian_widget, truncate_widget])
w=widgets.HBox([w1, w2])
display(w)
In [ ]:
# Define fitting function to find the temperature along the isentrope
def isentrope(rock, pressures, entropy, T_guess):
def _deltaS(T, S, P, rock):
rock.set_state(P, T)
return rock.S - S
sol = T_guess
temperatures = np.empty_like(pressures)
for i, P in enumerate(pressures):
sol = brentq(_deltaS, rock.bounds[1][0], rock.bounds[1][1], args=(entropy, P, rock))
temperatures[i] = sol
return temperatures
# Define function to find an isentrope given a
# 2D entropy interpolation function
def interp_isentrope(interp, pressures, entropy, T_guess):
def _deltaS(args, S, P):
T = args[0]
return interp(P, T)[0] - S
sol = [T_guess]
temperatures = np.empty_like(pressures)
for i, P in enumerate(pressures):
sol = fsolve(_deltaS, sol, args=(entropy, P))
temperatures[i] = sol[0]
return temperatures
# Define function to self consistently calculate depth and gravity profiles
# from pressure and density profiles.
def compute_depth_gravity_profiles(pressures, densities, surface_gravity, outer_radius):
gravity = [surface_gravity] * len(pressures) # starting guess
n_gravity_iterations = 5
for i in range(n_gravity_iterations):
# Integrate the hydrostatic equation
# Make a spline fit of densities as a function of pressures
rhofunc = UnivariateSpline(pressures, densities)
# Make a spline fit of gravity as a function of depth
gfunc = UnivariateSpline(pressures, gravity)
# integrate the hydrostatic equation
depths = np.ravel(odeint((lambda p, x: 1./(gfunc(x) * rhofunc(x))), 0.0, pressures))
radii = outer_radius - depths
rhofunc = UnivariateSpline(radii[::-1], densities[::-1])
poisson = lambda p, x: 4.0 * np.pi * burnman.constants.G * rhofunc(x) * x * x
gravity = np.ravel(odeint(poisson, surface_gravity*radii[0]*radii[0], radii))
gravity = gravity / radii / radii
return depths, gravity
def relaxed_profile(rock, potential_temperature, max_pressure,
surface_gravity, outer_radius, n_points, n_gridpoints, pressure_stdev,
temperature_smoothing_factor, max_temperature_stdev, truncate,
min_grid_temperature, max_grid_temperature):
min_grid_pressure = rock.bounds[0][0]
max_grid_pressure = rock.bounds[0][1]
min_grid_temperature = rock.bounds[1][0]
max_grid_temperature = rock.bounds[1][1]
rock.set_state(1.e5, potential_temperature)
entropy = rock.S
pressures = np.linspace(1.e5, max_pressure, n_points)
temperatures = isentrope(rock, pressures, entropy, potential_temperature)
isentrope_spline = UnivariateSpline(pressures, temperatures)
grid_pressures = np.linspace(min_grid_pressure, max_grid_pressure, n_gridpoints[0])
grid_temperatures = np.linspace(min_grid_temperature, max_grid_temperature, n_gridpoints[1])
unsmoothed_grid_isentrope_temperatures = isentrope_spline(grid_pressures)
pp, TT = np.meshgrid(grid_pressures, grid_temperatures)
mesh_shape = pp.shape
pp = np.ndarray.flatten(pp)
TT = np.ndarray.flatten(TT)
grid_entropies = np.zeros_like(pp)
grid_volumes = np.zeros_like(pp)
Tdiff = np.abs(isentrope_spline(pp) - TT)
# We could compute properties over the whole grid:
# grid_entropies, grid_volumes = rock.evaluate(['S', 'V'], pp, TT)
# However, we can save some time by computing only when temperature is close enough
# to the unsmoothed isentrope to affect the smoothing.
# The maximum temperature jump for most rocks is about 50 K, so a reasonable Tmax is
# ~50 + truncate*temperature_stdev. We pad a bit more (an extra 30 K) just to be sure.
temperature_stdev = np.min([max_temperature_stdev,
temperature_smoothing_factor * pressure_stdev *
np.max(np.abs( np.gradient(unsmoothed_grid_isentrope_temperatures) )) /
(grid_pressures[1] - grid_pressures[0])])
Tdiff_max = 50 + 30 + truncate*temperature_stdev
mask = [idx for idx, Td in enumerate(Tdiff) if Td < Tdiff_max]
grid_entropies[mask], grid_volumes[mask] = rock.evaluate(['S', 'V'], pp[mask], TT[mask])
grid_entropies = grid_entropies.reshape(mesh_shape)
grid_volumes = grid_volumes.reshape(mesh_shape)
# Having defined the grid and calculated unsmoothed properties,
# we now calculate the smoothed entropy and volume and derivatives with
# respect to pressure and temperature.
S_interps = burnman.tools.interp_smoothed_array_and_derivatives(array=grid_entropies,
x_values=grid_pressures,
y_values=grid_temperatures,
x_stdev=pressure_stdev,
y_stdev=temperature_stdev,
truncate=truncate)
interp_smoothed_S, interp_smoothed_dSdP, interp_smoothed_dSdT = S_interps
V_interps = burnman.tools.interp_smoothed_array_and_derivatives(array=grid_volumes,
x_values=grid_pressures,
y_values=grid_temperatures,
x_stdev=pressure_stdev,
y_stdev=temperature_stdev,
truncate=truncate)
interp_smoothed_V, interp_smoothed_dVdP, interp_smoothed_dVdT = V_interps
# Now we can calculate and plot the relaxed and smoothed properties along the isentrope
smoothed_temperatures = interp_isentrope(interp_smoothed_S, pressures, entropy, potential_temperature)
densities = rock.evaluate(['rho'], pressures, smoothed_temperatures)[0]
depths, gravity = compute_depth_gravity_profiles(pressures, densities, surface_gravity, outer_radius)
volumes = np.array([interp_smoothed_V(p, T)[0] for (p, T) in zip(*[pressures, smoothed_temperatures])])
dSdT = np.array([interp_smoothed_dSdT(p, T)[0] for (p, T) in zip(*[pressures, smoothed_temperatures])])
dVdT = np.array([interp_smoothed_dVdT(p, T)[0] for (p, T) in zip(*[pressures, smoothed_temperatures])])
dVdP = np.array([interp_smoothed_dVdP(p, T)[0] for (p, T) in zip(*[pressures, smoothed_temperatures])])
alphas_relaxed = dVdT / volumes
compressibilities_relaxed = -dVdP / volumes
specific_heats_relaxed = smoothed_temperatures * dSdT / rock.params['molar_mass']
dT = 0.1
Vpsub, Vssub = rock.evaluate(['p_wave_velocity', 'shear_wave_velocity'],
pressures, smoothed_temperatures-dT/2.)
Vpadd, Vsadd = rock.evaluate(['p_wave_velocity', 'shear_wave_velocity'],
pressures, smoothed_temperatures+dT/2.)
Vps = (Vpadd + Vpsub)/2.
Vss = (Vsadd + Vssub)/2.
dVpdT = (Vpadd - Vpsub)/dT
dVsdT = (Vsadd - Vssub)/dT
#print('Min and max relaxed property when pressure smoothing standard deviation is {0:.2f} GPa'.format(pressure_stdev/1.e9))
#print('Specific heat: {0:.2e}, {1:.2e}'.format(np.min(specific_heats_relaxed), np.max(specific_heats_relaxed)))
#print('Thermal expansivity: {0:.2e}, {1:.2e}'.format(np.min(alphas_relaxed), np.max(alphas_relaxed)))
#print('Compressibilities: {0:.2e}, {1:.2e}\n'.format(np.min(compressibilities_relaxed), np.max(compressibilities_relaxed)))
return (smoothed_temperatures, depths, gravity, densities,
alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed,
Vss, Vps, dVsdT, dVpdT)
In [ ]:
rock = None
fig = None
potential_temperature = None
default_P_gaussian = None
default_T_smoothing = None
max_pressure = None
surface_gravity = None
outer_radius = None
n_points = None
n_gridpoints = None
max_temperature_stdev = None
truncate = None
min_grid_temperature = None
max_grid_temperature = None
T_line = None
g_line = None
alpha_line = None
beta_line = None
cp_line = None
x = None
pressures = None
def setup(button):
global rock, fig, potential_temperature, default_P_gaussian, default_T_smoothing, max_pressure, surface_gravity, outer_radius, n_points, n_gridpoints, max_temperature_stdev, truncate, min_grid_temperature, max_grid_temperature
global T_line, g_line, alpha_line, beta_line, cp_line, x, xlabel, pressures
clear_output()
rock = burnman.PerplexMaterial(rock_path_widget.value)
potential_temperature = potlT_widget.value
max_pressure = maxP_widget.value * 1.e9
surface_gravity = g0_widget.value
outer_radius = r0_widget.value * 1000.
n_points = n_points_widget.value
n_gridpoints = (n_P_gridpoints_widget.value,
n_T_gridpoints_widget.value) # number of p, T grid points for smoothing S and V
max_temperature_stdev = max_T_gaussian_widget.value # max T_stdev
truncate = truncate_widget.value # truncates the convolution Gaussian at 4 sigma
# First we calculate the isentrope at a given potential temperature
rock.set_state(1.e5, potential_temperature)
entropy = rock.S
pressures = np.linspace(1.e5, max_pressure, n_points)
temperatures = isentrope(rock, pressures, entropy, potential_temperature)
isentrope_spline = UnivariateSpline(pressures, temperatures)
# Properties can then be calculated along the isentrope
properties = rock.evaluate(['V', 'rho', 'molar_heat_capacity_p',
'thermal_expansivity', 'isothermal_compressibility',
'p_wave_velocity', 'shear_wave_velocity'],
pressures, temperatures)
volumes, densities, C_p, alphas, compressibilities, p_wave_velocities, s_wave_velocities = properties
specific_heats = C_p / rock.params['molar_mass']
depths, gravity = compute_depth_gravity_profiles(pressures, densities,
surface_gravity, outer_radius)
#x = pressures/1.e9
x = depths/1.e3
xlabel = 'Depths (km)'
plt.rcParams['figure.figsize'] = 8, 5 # inches
fig = plt.figure()
px, py = [2, 3]
ax_T = fig.add_subplot(px, py, 1)
ax_T.plot(x, temperatures, label='unrelaxed')
ax_T.set_ylabel('Temperature (K)')
ax_T.set_xlabel(xlabel)
ax_g = fig.add_subplot(px, py, 2)
ax_g.plot(x, gravity)
ax_g.set_ylabel('Gravity (m/s^2)')
ax_g.set_xlabel(xlabel)
ax_rho = fig.add_subplot(px, py, 3)
ax_rho.plot(x, densities, label='$\rho$ (kg/m$^3$)')
ax_rho.plot(x, p_wave_velocities, label='P (km/s)')
ax_rho.plot(x, s_wave_velocities, label='S (km/s)')
ax_rho.set_ylabel('Densities/Velocities')
ax_rho.set_xlabel(xlabel)
ax_alpha = fig.add_subplot(px, py, 4)
ax_alpha.plot(x, alphas)
ax_alpha.set_ylabel('alpha (/K)')
ax_alpha.set_xlabel(xlabel)
ax_beta = fig.add_subplot(px, py, 5)
ax_beta.plot(x, compressibilities)
ax_beta.set_ylabel('compressibilities (/Pa)')
ax_beta.set_xlabel(xlabel)
ax_cp = fig.add_subplot(px, py, 6)
ax_cp.plot(x, specific_heats)
ax_cp.set_ylabel('Cp (J/K/kg)')
ax_cp.set_xlabel(xlabel)
default_P_gaussian = 0.5e9
default_T_smoothing = 0.25
# Relaxed, unsmoothed properties
smoothed_temperatures, depths, gravity, densities, alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed, Vss, Vps, dVsdT, dVpdT = relaxed_profile(rock, potential_temperature, max_pressure, surface_gravity, outer_radius, n_points, n_gridpoints, 0., default_T_smoothing, max_temperature_stdev, truncate,
min_grid_temperature, max_grid_temperature)
ax_T.plot(x, smoothed_temperatures, label='relaxed')
ax_g.plot(x, gravity)
ax_alpha.plot(x, alphas_relaxed)
ax_beta.plot(x, compressibilities_relaxed)
ax_cp.plot(x, specific_heats_relaxed)
# Relaxed, smoothed properties
smoothed_temperatures, depths, gravity, densities, alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed, Vss, Vps, dVsdT, dVpdT = relaxed_profile(rock, potential_temperature, max_pressure, surface_gravity, outer_radius, n_points, n_gridpoints, default_P_gaussian, default_T_smoothing, max_temperature_stdev, truncate,
min_grid_temperature, max_grid_temperature)
T_line, = ax_T.plot(x, smoothed_temperatures, label='relaxed, smoothed')
g_line, = ax_g.plot(x, gravity)
alpha_line, = ax_alpha.plot(x, alphas_relaxed)
beta_line, = ax_beta.plot(x, compressibilities_relaxed)
cp_line, = ax_cp.plot(x, specific_heats_relaxed)
ax_T.legend(loc='lower right',prop={'size':8})
fig.tight_layout()
btn = widgets.Button(description="update")
btn.on_click(setup)
setup(None) # run setup once in the beginning
display(btn)
In [ ]:
@widgets.interact(potential_temperature=widgets.FloatText(value=potential_temperature, description='Potential temperature (K):', continuous_update=False),
P_gaussian_GPa=widgets.FloatText(value=default_P_gaussian/1.e9, min=0.0, max=5.0, description='P smoothing $\sigma$ (GPa):', continuous_update=False),
T_smoothing=widgets.FloatSlider(value = default_T_smoothing, min=0.0, max=0.5, step=0.01, description='T smoothing factor', continuous_update=False))
def update_plot(potential_temperature, P_gaussian_GPa, T_smoothing):
# Relaxed, smoothed properties
smoothed_temperatures, depths, gravity, densities, alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed, Vss, Vps, dVsdT, dVpdT = relaxed_profile(rock, potential_temperature, max_pressure, surface_gravity, outer_radius, n_points, n_gridpoints, P_gaussian_GPa*1.e9, T_smoothing, max_temperature_stdev, truncate,
min_grid_temperature, max_grid_temperature)
T_line.set_data([x, smoothed_temperatures])
g_line.set_data([x, gravity])
alpha_line.set_data([x, alphas_relaxed])
beta_line.set_data([x, compressibilities_relaxed])
cp_line.set_data([x, specific_heats_relaxed])
plt.draw()
# Convert to equal slices in depth
p_spline = UnivariateSpline(depths, pressures)
depths_eq = np.linspace(depths[0], depths[-1], n_points)
pressures_eq = p_spline(depths_eq)
smoothed_temperatures = np.interp(pressures_eq, pressures, smoothed_temperatures)
densities = np.interp(pressures_eq, pressures, densities)
gravity = np.interp(pressures_eq, pressures, gravity)
alphas_relaxed = np.interp(pressures_eq, pressures, alphas_relaxed)
specific_heats_relaxed = np.interp(pressures_eq, pressures, specific_heats_relaxed)
compressibilities_relaxed = np.interp(pressures_eq, pressures, compressibilities_relaxed)
Vss = np.interp(pressures_eq, pressures, Vss)
Vps = np.interp(pressures_eq, pressures, Vps)
dVsdT = np.interp(pressures_eq, pressures, dVsdT)
dVpdT = np.interp(pressures_eq, pressures, dVpdT)
# Finally, here's the ability to output smoothed, relaxed properties for use in ASPECT
# depth, pressure, temperature, density, gravity, Cp (per kilo), thermal expansivity
outfile = 'isentrope_properties.txt'
np.savetxt(outfile, X=np.array([depths_eq, pressures_eq, smoothed_temperatures,
densities, gravity, alphas_relaxed,
specific_heats_relaxed,
compressibilities_relaxed,
Vss, Vps, dVsdT, dVpdT]).T,
header='# This ASPECT-compatible file contains material properties calculated along an isentrope by the BurnMan software.\n# POINTS: '+str(n_points)+' \n# depth (m), pressure (Pa), temperature (K), density (kg/m^3), gravity (m/s^2), thermal expansivity (1/K), specific heat (J/K/kg), compressibility (1/Pa), seismic Vs (m/s), seismic Vp (m/s), seismic dVs/dT (m/s/K), seismic dVp/dT (m/s/K)\ndepth pressure temperature density gravity thermal_expansivity specific_heat compressibility seismic_Vs seismic_Vp seismic_dVs_dT seismic_dVp_dT',
fmt='%.10e', delimiter='\t', comments='')
print('File saved to {0}'.format(outfile))
display(FileLink(outfile))
In [ ]:
%%html
<script>
// Note:
// This html code block will
// 1. do "run-all-cells" 0.5s after the kernel is loaded
// 2. hide all code blocks (and offer a button to toggle the code)
require(
['base/js/namespace', 'jquery'],
function(jupyter, $) {
$(jupyter.events).on("kernel_ready.Kernel", function () {
// js widgets are not available immediately. Instead, trigger this a little later:
console.log("kernel_ready triggered, preparing auto run-all-cells");
setTimeout(function() {
console.log("Auto-running all cells...");
jupyter.actions.call('jupyter-notebook:run-all-cells');
}, 500);
//jupyter.actions.call('jupyter-notebook:save-notebook');
});
}
);
code_show=false;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
function init() { $('div.input').hide();}
$( document ).ready(init);
$( document ).load(init);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>
In [ ]: