Source code for burnman.classes.planet

from __future__ import print_function
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU
# GPL v2 or later.

import numpy as np
import warnings
from .material import material_property


[docs]class Planet(object): """ A class to build (self-consistent) Planets made out of Layers (``burnman.Layer``). By default the planet is set to be self-consistent (with zero pressure at the surface and zero gravity at the center), but this can be overwritte using the set_pressure_mode(). Pressure_modes defined in the individual layers will be ignored. If temperature modes are already set for each of the layers, when the planet is initialized, the planet will be built immediately. """ def __init__(self, name, layers, n_max_iterations=50, max_delta=1.e-5, verbose=False): """ Parameters ---------- name : string Name of planet layers : list of ``burnman.Layer`` Layers to build the planet out of (layers are sorted within the planet) n_max_iterations : int Maximum number of iterations to reach self-consistent planet (default = 50) max_delta : float Relative update to the center pressure of the planet between iterations to stop iterations (default = 1.e-5) """ # sort layers self.layers = sorted(layers, key=lambda x: x.inner_radius) # assert layers attach to one another if len(self.layers) > 1: for i in range(1, len(self.layers)): assert(self.layers[i].inner_radius == self.layers[i - 1].outer_radius) self.name = name self.radii = self.evaluate(['radii']) self.n_slices = len(self.radii) self.radius_planet = max(self.radii) self.volume = 4./3.*np.pi*np.power(self.radius_planet, 3.) for layer in self.layers: layer.n_start = np.where(self.radii == layer.inner_radius)[0][-1] layer.n_end = np.where(self.radii == layer.outer_radius)[0][0] + 1 self._cached = {} self.verbose = verbose self.set_pressure_mode( n_max_iterations=n_max_iterations, max_delta=max_delta) def __iter__(self): """ Planet object will iterate over Layers. """ return list(self.layers).__iter__() def __str__(self): """ Prints details of the planet """ writing = '{0} consists of {1} layers:\n'.format( self.name, len(self.layers)) for layer in self: writing = writing + layer.__str__() return writing
[docs] def reset(self): """ Resets all cached material properties. It is typically not required for the user to call this function. """ self._cached = {}
[docs] def get_layer(self, name): """ Returns a layer with a given name Parameters ---------- name : string Given name of a layer """ for layer in self.layers: if layer.name == name: return layer raise LookupError()
[docs] def get_layer_by_radius(self, radius): """ Returns a layer in which this radius lies Parameters ---------- radius : float radius at which to evaluate the layer """ for layer in self.layers: if layer.outer_radius >= radius: return layer raise LookupError()
[docs] def evaluate(self, properties, radlist=None): """ Function that is generally used to evaluate properties of the different layers and stitch them together. If asking for different radii than the internal radlist, pressure and temperature values are interpolated and the layer material evaluated at those pressures and temperatures. Parameters ---------- properties : list of strings List of properties to evaluate radlist : array of floats Radii to evaluate properties at. If left empty, internal radius lists are used. Returns ------- properties_array : numpy array 1D or 2D array of requested properties (1D if only one property was requested) """ if radlist is None: values = np.empty([len(properties), np.sum([len(layer.radii) for layer in self.layers])]) for i, prop in enumerate(properties): if prop == 'depth': values[i] = np.array( [self.radius_planet - r for layer in self.layers for r in layer.radii]) else: j = 0 for layer in self.layers: vals = getattr(layer, prop) values[i, j:j + len(vals)] = vals j += len(vals) else: values = np.empty([len(properties), len(radlist)]) l_idx = [i for i, layer in enumerate(self.layers) for r in radlist if r >= layer.inner_radius and r <= layer.outer_radius] for j, r in enumerate(radlist): values[:, j] = self.layers[l_idx[j]].evaluate( properties, [r], self.radius_planet).T[0] if values.shape[0] == 1: values = values[0] return values
[docs] def set_pressure_mode(self, pressure_mode='self-consistent', pressures=None, pressure_top=0., gravity_bottom=0., n_max_iterations=50, max_delta=1.e-5): """ Sets the pressure mode of the planet by user-defined values are in a self-consistent fashion. pressure_mode is 'user-defined' or 'self-consistent'. The default for the planet is self-consistent, with zero pressure at the surface and zero pressure at the center. Parameters ---------- pressure_mode : string This can be set to 'user-defined' or 'self-consistent' pressures : array of floats Pressures (Pa) to set layer to ('user-defined'). This should be the same length as defined radius array for the layer pressure_top : float Pressure (Pa) at the top of the layer. gravity_bottom : float gravity (m/s^2) at the bottom the layer n_max_iterations : int Maximum number of iterations to reach self-consistent pressures (default = 50) """ self.reset() assert(pressure_mode == 'user-defined' or pressure_mode == 'self-consistent') self.pressure_mode = pressure_mode self.gravity_bottom = gravity_bottom if pressure_mode == 'user-defined': assert(len(pressures) == len(self.radii)) self._pressures = pressures warnings.warn('User-defined pressures mean that the planet is ' 'unlikely to be self-consistent') if pressure_mode == 'self-consistent': self.pressure_top = pressure_top self.n_max_iterations = n_max_iterations self.max_delta = max_delta
[docs] def make(self): """ This routine needs to be called before evaluating any properties. If pressures and temperatures are self-consistent, they are computed across the planet here. Also initializes an array of materials in each Layer to compute properties from. """ self.reset() for layer in self.layers: assert(layer.temperature_mode is not None) if self.pressure_mode == 'user-defined': self._temperatures = self._evaluate_temperature(self._pressures) if self.pressure_mode == 'self-consistent': new_press = self.pressure_top + \ (-self.radii + max(self.radii)) * \ 1.e3 # initial pressure curve guess temperatures = self._evaluate_temperature(new_press) # Make it self-consistent!!! i = 0 while i < self.n_max_iterations: i += 1 ref_press = new_press new_grav, new_press = self._evaluate_eos(new_press, temperatures, self.gravity_bottom, self.pressure_top) temperatures = self._evaluate_temperature(new_press) rel_err = abs((max(ref_press) - max(new_press)) / max(new_press)) if self.verbose: print(f'Iteration {i:0d} maximum relative pressure error: ' f'{rel_err:.1e}') if rel_err < self.max_delta: break self.pressures = new_press self.temperatures = temperatures self._gravity = new_grav for layer in self.layers: layer.sublayers = [] layer.pressures = self.pressures[layer.n_start: layer.n_end] layer.temperatures = self.temperatures[layer.n_start: layer.n_end] layer.gravity_bottom = self._gravity[layer.n_start - 1] layer.pressure_mode = 'set-in-planet' for i in range(len(layer.radii)): layer.sublayers.append(layer.material.copy()) layer.sublayers[i].set_state( layer.pressures[i], layer.temperatures[i])
def _evaluate_eos(self, pressures, temperatures, gravity_bottom, pressure_top): """ Used to update the pressure profile in set_state() """ density = self._evaluate_density(pressures, temperatures) grav = self._compute_gravity(density, gravity_bottom) press = self._compute_pressure(density, grav, pressure_top) return grav, press def _evaluate_density(self, pressures, temperatures): """ Used to update the density profile in _evaluate_eos() """ density = [] for layer in self.layers: density.append(layer.material.evaluate(['density'], pressures[layer.n_start:layer.n_end], temperatures[layer.n_start:layer.n_end])) return np.squeeze(np.hstack(density)) def _evaluate_temperature(self, pressures): """ Returns the temperatures of different layers for given pressures. Used by set_state() """ temps = [] temperature_top = None for layer in self.layers[::-1]: if temperature_top is None or layer.temperature_top is not None: temperature_top = layer.temperature_top temps.extend(layer._evaluate_temperature( (pressures[layer.n_start:layer.n_end]), temperature_top)[::-1]) temperature_top = temps[-1] return np.hstack(np.squeeze(temps))[::-1] def _compute_gravity(self, density, gravity_bottom): """ Calculate the gravity of the planet, based on a density profile. This integrates Poisson's equation in radius, under the assumption that the planet is laterally homogeneous. Used to update the gravity profile in _evaluate_eos() """ start_gravity = gravity_bottom grav = [] for layer in self.layers: grav.extend(layer._compute_gravity( density[layer.n_start: layer.n_end], start_gravity)) start_gravity = grav[-1] return np.array(grav) def _compute_pressure(self, density, gravity, pressure_top): """ Calculate the pressure profile based on density and gravity. This integrates the equation for hydrostatic equilibrium P = rho g z. Used to update the pressure profile in _evaluate_eos() """ start_pressure = pressure_top press = [] for layer in self.layers[::-1]: press.extend(layer._compute_pressure(density[layer.n_start:layer.n_end], gravity[layer.n_start: layer.n_end], start_pressure)[::-1]) start_pressure = press[-1] return np.array(press)[::-1] @property def mass(self): """ calculates the mass of the entire planet [kg] """ return np.sum([layer.mass for layer in self.layers]) @property def average_density(self): """ calculates the average density of the entire planet [kg/m^3] """ return self.mass/self.volume @property def moment_of_inertia(self): """ #Returns the moment of inertia of the planet [kg m^2] """ return np.sum([layer.moment_of_inertia for layer in self.layers]) @property def moment_of_inertia_factor(self): """ #Returns the moment of inertia of the planet [kg m^2] """ moment_factor = self.moment_of_inertia / self.mass / \ self.radius_planet / self.radius_planet return moment_factor @property def depth(self): """ Returns depth of the layer [m] """ return self.evaluate(['depth']) @property def gravity(self): """ Returns gravity of the layer [m s^(-2)] """ return self.evaluate(['gravity']) @property def bullen(self): """ Returns the Bullen parameter """ return self.evaluate(['bullen']) @property def brunt_vasala(self): return self.evaluate(['brunt_vasala']) @property def pressure(self): """ Returns current pressure that was set with :func:`~burnman.Material.set_state`. Aliased with :func:`~burnman.Material.P`. Returns ------- pressure : array of floats Pressure in [Pa]. """ return self.pressures @property def temperature(self): """ Returns current temperature that was set with :func:`~burnman.Material.set_state`. Aliased with :func:`~burnman.Material.T`. Returns ------- temperature : array of floats Temperature in [K]. """ return self.temperatures @material_property def molar_internal_energy(self): """ Returns the molar internal energy of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.energy`. Returns ------- molar_internal_energy : array of floats The internal energy in [J/mol]. """ return self.evaluate(['molar_internal_energy']) @material_property def molar_gibbs(self): """ Returns the molar Gibbs free energy of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.gibbs`. Returns ------- molar_gibbs : array of floats Gibbs free energy in [J/mol]. """ return self.evaluate(['molar_gibbs']) @material_property def molar_helmholtz(self): """ Returns the molar Helmholtz free energy of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.helmholtz`. Returns ------- molar_helmholtz : array of floats Helmholtz free energy in [J/mol]. """ return self.evaluate(['molar_helmholtz']) @material_property def molar_mass(self): """ Returns molar mass of the planet. Needs to be implemented in derived classes. Returns ------- molar_mass : array of floats Molar mass in [kg/mol]. """ return self.evaluate(['molar_mass']) @material_property def molar_volume(self): """ Returns molar volume of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.V`. Returns ------- molar_volume : array of floats Molar volume in [m^3/mol]. """ return self.evaluate(['molar_volume']) @material_property def density(self): """ Returns the density of this planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.rho`. Returns ------- density : array of floats The density of this material in [kg/m^3]. """ return self.evaluate(['density']) @material_property def molar_entropy(self): """ Returns molar entropy of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.S`. Returns ------- molar_entropy : array of floats Entropy in [J/mol]. """ return self.evaluate(['molar_entropy']) @material_property def molar_enthalpy(self): """ Returns molar enthalpy of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.H`. Returns ------- molar_enthalpy : array of floats Enthalpy in [J/mol]. """ return self.evaluate(['molar_enthalpy']) @material_property def isothermal_bulk_modulus(self): """ Returns isothermal bulk modulus of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.K_T`. Returns ------- isothermal_bulk_modulus : array of floats Bulk modulus in [Pa]. """ return self.evaluate(['isothermal_bulk_modulus']) @material_property def adiabatic_bulk_modulus(self): """ Returns the adiabatic bulk modulus of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.K_S`. Returns ------- adiabatic_bulk_modulus : array of floats Adiabatic bulk modulus in [Pa]. """ return self.evaluate(['adiabatic_bulk_modulus']) @material_property def isothermal_compressibility(self): """ Returns isothermal compressibility of the planet (or inverse isothermal bulk modulus). Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.beta_T`. Returns ------- (K_T)^-1 : array of floats Compressibility in [1/Pa]. """ return self.evaluate(['istothermal_compressibility']) @material_property def adiabatic_compressibility(self): """ Returns adiabatic compressibility of the planet (or inverse adiabatic bulk modulus). Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.beta_S`. Returns ------- adiabatic_compressibility : array of floats adiabatic compressibility in [1/Pa]. """ return self.evaluate(['adiabatic_compressibility']) @material_property def shear_modulus(self): """ Returns shear modulus of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.beta_G`. Returns ------- shear_modulus : array of floats Shear modulus in [Pa]. """ return self.evaluate(['shear_modulus']) @material_property def p_wave_velocity(self): """ Returns P wave speed of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.v_p`. Returns ------- p_wave_velocity : array of floats P wave speed in [m/s]. """ return self.evaluate(['p_wave_velocity']) @material_property def bulk_sound_velocity(self): """ Returns bulk sound speed of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.v_phi`. Returns ------- bulk sound velocity: array of floats Sound velocity in [m/s]. """ return self.evaluate(['bulk_sound_velocity']) @material_property def shear_wave_velocity(self): """ Returns shear wave speed of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.v_s`. Returns ------- shear_wave_velocity : array of floats Wave speed in [m/s]. """ return self.evaluate(['shear_wave_velocity']) @material_property def grueneisen_parameter(self): """ Returns the grueneisen parameter of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.gr`. Returns ------- gr : array of floats Grueneisen parameters [unitless]. """ return self.evaluate(['grueneisen_parameter']) @material_property def thermal_expansivity(self): """ Returns thermal expansion coefficient of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.alpha`. Returns ------- alpha : array of floats Thermal expansivity in [1/K]. """ return self.evaluate(['thermal_expansivity']) @material_property def molar_heat_capacity_v(self): """ Returns molar heat capacity at constant volume of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.C_v`. Returns ------- molar_heat_capacity_v : array of floats Heat capacity in [J/K/mol]. """ return self.evaluate(['molar_heat_capacity_v']) @material_property def molar_heat_capacity_p(self): """ Returns molar heat capacity at constant pressure of the planet. Needs to be implemented in derived classes. Aliased with :func:`~burnman.Material.C_p`. Returns ------- molar_heat_capacity_p : array of floats Heat capacity in [J/K/mol]. """ return self.evaluate(['molar_heat_capacity_p']) # # Aliased properties @property def P(self): """Alias for :func:`~burnman.Material.pressure`""" return self.pressure @property def T(self): """Alias for :func:`~burnman.Material.temperature`""" return self.temperature @property def energy(self): """Alias for :func:`~burnman.Material.molar_internal_energy`""" return self.molar_internal_energy @property def helmholtz(self): """Alias for :func:`~burnman.Material.molar_helmholtz`""" return self.molar_helmholtz @property def gibbs(self): """Alias for :func:`~burnman.Material.molar_gibbs`""" return self.molar_gibbs @property def V(self): """Alias for :func:`~burnman.Material.molar_volume`""" return self.molar_volume @property def rho(self): """Alias for :func:`~burnman.Material.density`""" return self.density @property def S(self): """Alias for :func:`~burnman.Material.molar_entropy`""" return self.molar_entropy @property def H(self): """Alias for :func:`~burnman.Material.molar_enthalpy`""" return self.molar_enthalpy @property def K_T(self): """Alias for :func:`~burnman.Material.isothermal_bulk_modulus`""" return self.isothermal_bulk_modulus @property def K_S(self): """Alias for :func:`~burnman.Material.adiabatic_bulk_modulus`""" return self.adiabatic_bulk_modulus @property def beta_T(self): """Alias for :func:`~burnman.Material.isothermal_compressibility`""" return self.isothermal_compressibility @property def beta_S(self): """Alias for :func:`~burnman.Material.adiabatic_compressibility`""" return self.adiabatic_compressibility @property def G(self): """Alias for :func:`~burnman.Material.shear_modulus`""" return self.shear_modulus @property def v_p(self): """Alias for :func:`~burnman.Material.p_wave_velocity`""" return self.p_wave_velocity @property def v_phi(self): """Alias for :func:`~burnman.Material.bulk_sound_velocity`""" return self.bulk_sound_velocity @property def v_s(self): """Alias for :func:`~burnman.Material.shear_wave_velocity`""" return self.shear_wave_velocity @property def gr(self): """Alias for :func:`~burnman.Material.grueneisen_parameter`""" return self.grueneisen_parameter @property def alpha(self): """Alias for :func:`~burnman.Material.thermal_expansivity`""" return self.thermal_expansivity @property def C_v(self): """Alias for :func:`~burnman.Material.molar_heat_capacity_v`""" return self.molar_heat_capacity_v @property def C_p(self): """Alias for :func:`~burnman.Material.molar_heat_capacity_p`""" return self.molar_heat_capacity_p