API Reference

pytesimal.analysis module

Analyse the results of the conductive cooling model for a planetesimal.

This module contains functions to calculate the cooling rates of meteorites based on the empirical relations suggested by Yang et al. (2010); see full references in Murphy Quinlan et al. (2021). It also contains functions to analyse the temperature arrays produced by the pytesimal.numerical_methods module, allowing estimation of the depth of genesis of pallasite meteorites, the relative timing of paleomagnetic recording in meteorites and core dynamo action, and calculation of cooling rates in the mantle of the planetesimal through time.

pytesimal.analysis.cooling_rate(temperature_array, timestep)

Calculate an array of cooling rates from temperature array.

pytesimal.analysis.cooling_rate_cloudyzone_diameter(d)

Cooling rate calculated using cloudy zone particle diameter in nm.

Constants from Yang et al., 2010; obtained by comparing cz particles and tetrataenite bandwidth to modelled Ni diffusion in kamacite and taenite.

Parameters

d (float) – Cloudy zone particle size in nm.

Returns

cz_rate – The cooling rate in K/Myr.

Return type

float

pytesimal.analysis.cooling_rate_tetra_width(tw)

Cooling rate calculated using tetrataenite bandwidth in nm.

Constants from Yang et al., 2010; obtained by comparing cz particles and tetrataenite bandwidth to modelled Ni diffusion in kamacite and taenite.

Parameters

tw (float) – Tetrataenite bandwidth in nm.

Returns

t_rate – The cooling rate in K/Myr.

Return type

float

pytesimal.analysis.cooling_rate_to_seconds(cooling_rate)

Convert cooling rates to seconds.

Parameters

cooling_rate (float) – The cooling rate in K/Myr.

Returns

new_cooling_rate – The cooling rate in K/s.

Return type

float

pytesimal.analysis.core_freezing(coretemp, max_time, times, latent, temp_core_melting, timestep=100000000000.0)

Calculate when the core starts and finishes solidifying.

Takes core temperature and returns boolean array of when the core is below the freezing/melting temperature.

Parameters
  • coretemp (numpy.ndarray) – Array of temperatures in the core, in Kelvin.

  • max_time (float) – Length of time the model runs for, in seconds.

  • times (numpy.ndarray) – Array from 0 to the max time +0.5* the timestep, with a spacing equal to the timestep.

  • latent (list) – List of total latent heat extracted since core began freezing, at each timestep.

  • temp_core_melting (float) – Melting point of core material, in Kelvin.

  • timestep (float, default 1e11) – Discretisation timestep in seconds.

Returns

  • core_frozen (boolean array where temperature <= 1200 K)

  • times_frozen (array of indices of times where the temp <= 1200 K)

  • time_core_frozen (when the core starts to freeze, in seconds)

  • fully_frozen (when the core finished freezing, in seconds)

pytesimal.analysis.meteorite_depth_and_timing(CR, temperatures, dT_by_dt, radii, r_planet, core_size_factor, time_core_frozen, fully_frozen, dr=1000.0, dt=100000000000.0)

Find depth of genesis given the cooling rate.

Function finds the depth, given the cooling rate, and checks if the 593K contour crosses this depth during core solidification, implying whether or not the meteorite is expected to record core dynamo action.

Parameters
  • CR (float) – cooling rate of meteorite, in K/s.

  • temperatures (numpy.ndarray) – Array of mantle temperatures, in Kelvin.

  • dT_by_dt (numpy.ndarray) – Array of mantle cooling rates, in K/dt.

  • radii (numpy.ndarray) – Mantle radii spaced by dr, in m.

  • r_planet (float) – Planetesimal radius, in m.

  • core_size_factor (float, <1) – Radius of the core, expressed as a fraction of r_planet.

  • time_core_frozen (float) – The time the core begins to freeze, in dt.

  • fully_frozen (float) – The time the core is fully frozen, in dt.

  • dr (float, default 1000.0) – Radial step for numerical discretisation, in m.

Returns

  • depth (float) – Depth of genesis of meteorite, in km.

  • string (string) – Relative timing of tetrataenite formation and core crystallisation, in a string format

  • time_core_frozen (float) – The time the core begins to freeze, in dt.

  • Time_of_Crossing (float) – When the meteorite cools through tetrataenite formation temperature, in dt.

  • Critical_Radius (float) – Depth of meteorite genesis given as radius value, in m.

pytesimal.core_function module

Create and track the temperature in the planetesimal core.

This module allows the user to track the temperature of a simple isothermal eutectic core and calculate the change in temperature in the core over a timestep based on the heat extracted across the core-mantle boundary.

The class IsothermalEutecticCore allows a core object to be instantiated. This core object keeps track of its temperature as it cools, and this temperature history can be called at any time in the form of a 1D timeseries or cast across the core radius (as the core is isothermal).

Classes:

IsothermalEutecticCore

Notes

The core object can be replaced with a more complex core that interacts with the mantle in the same way (by extracting energy in Watts across the CMB over a timestep and providing a resulting boundary temperature).

class pytesimal.core_function.IsothermalEutecticCore(initial_temperature, melting_temperature, outer_r, inner_r, rho, cp, core_latent_heat, lat=0)

Bases: object

Core class to represent and manipulate core temperature and latent heat.

initial_temperature

Initial uniform temperature of the core, in K.

Type

float

melting_temperature

Temperature at which core crystallisation initiates, in K.

Type

float

outer_r

Outer core radius, in m.

Type

float

inner_r

Inner core radius, not used by this simple implementation of the core (set to zero), but included so that more complex core models can be coupled to the mantle discretisation function.

Type

float

rho

Core density, kg m^-3.

Type

float

cp

Core heat capacity, J kg^-1 K^-1.

Type

float

core_latent_heat

Latent heat of crystallisation of the core, used to calculate time for core to solidify fully, J kg^-1.

Type

float

lat

Tracks latent heat of core, always initially zero in current implementation but included for forward compatibility with a coupled model where core has already cooled by some degree, in J kg^-1.

Type

float, optional

__init__(initial_temperature, melting_temperature, outer_r, inner_r, rho, cp, core_latent_heat, lat=0)

Create a new core with temperature and latent heat.

extract_heat(power, timestep)

Extract heat (in W) across the core-mantle boundary

Given power (W) and timestep (s), update the boundary_temperature and temperature to reflect the associated cooling. If temperature is equal or less than the melting temperature of the core material, then the core begins to freeze and temperature does not change. Instead, latent heat is tracked (latent). Once latent is greater or equal to the maximum latent heat of the core, the core is fully frozen and begins to cool again as before.

Parameters
  • power (float) – Heat extracted across the CMB in Watts.

  • timestep (float) – The time over which the heat is extracted (in s).

temperature_array_1D()

Return a time-series of core boundary temperatures

Returns

temp_array – Time series of boundary_temperature, in K.

Return type

numpy.ndarray

temperature_array_2D(coretemp_array)

Cast the core boundary temperatures to an array of radii in time

Parameters

coretemp_array (numpy.ndarray) – Array of zeros to be filled wth core temperature history.

Returns

coretemp_array – Array of core temperature history, in K.

Return type

numpy.ndarray

pytesimal.load_plot_save module

Load Plot and Save Module

This module contains functions to load parameters and results from files, to plot results either following a model run or from file, and to save parameters and results to file following a model run.

Example

A directory and parameter file can be quickly generated:

folder = 'path/to/folder'
filename = 'example_param_file.txt'
filepath=f'{folder}/{filename}'

check_folder_exists(folder)
make_default_param_file(filepath=filepath)

This parameters file in json format can then be opened, edited, renamed or moved, and loaded in to set parameter values for a model run.

Notes

Depending on usage, some functions take a folder and filename argument and create an absolute path with these to save or load a file, while some take a full filepath. Please check which argument is required.

pytesimal.load_plot_save.check_folder_exists(folder)

Check directory exists and make directory if not.

pytesimal.load_plot_save.get_million_years_formatters(timestep, maxtime)

Return a matplotlib formatter.

Creates two matplotlib formatters, one to go from timesteps to myrs and one to go from cooling rate per timestep to cooling rate per million years.

Parameters
  • timestep (float) – Numerical timestep, in s.

  • maxtime (float) – Total time for model run, in s.

pytesimal.load_plot_save.load_params_from_file(filepath='example_params.txt')

Load parameters from a json file and return variable values.

Parameters

filepath (str) – Absolute or relative path to load the json parameter file from.

Returns

  • run_ID (str) – Identifier for the specific model run.

  • folder (str) – Absolute path to directory where file is to be saved. Existence of the directory can be checked with the check_folder_exists() function.

  • timestep (float) – The timestep used in numerical method, in s.

  • r_planet (float) – The radius of the planet in m.

  • core_size_factor (float) – The core size as a fraction of the total planet radius.

  • reg_fraction (float) – The regolith thickness as a fraction of the total planet radius.

  • max_time (float) – The total run-time of the model, in millions of years (Myr).

  • temp_core_melting (float) – The melting temperature of the core, in K.

  • mantle_heat_cap_value (float) – The heat capacity of mantle material, in J kg^-1 K^-1.

  • mantle_density_value (float) – The density of mantle material, in kg m^-3.

  • mantle_conductivity_value (float) – The conductivity of the mantle, in W m^-1 K^-1.

  • core_cp (float) – The heat capacity of the core, in J kg^-1 K^-1.

  • core_density (float) – The density of the core, in kg m^-3.

  • temp_init (float, list, numpy array) – The initial temperature of the body, in K.

  • temp_surface (float) – The surface temperature of the planet, in K.

  • core_temp_init (float) – The initial temperature of the core, in K.

  • core_latent_heat (float) – The latent heat of crystallisation of the core, in J kg^-1.

  • kappa_reg (float) – The regolith constant diffusivity, m^2 s^-1.

  • dr (float) – The radial step used in the numerical model, in m.

  • cond_constant (str) – Flag of y or n to specify if mantle conductivity is constant.

  • density_constant (str) – Flag of y or n to specify if mantle density is constant.

  • heat_cap_constant (str) – Flag of y or n to specify if mantle heat capacity is constant.

pytesimal.load_plot_save.make_default_param_file(filepath='example_params.txt')

Save an example parameter json file with default parameters.

Save a dictionary of default parameters as a human-readable json txt file. This example file can then be used as an input file as-is, or can be edited in order to modify the model set up. References for the example values can be found in Murphy Quinlan et al. (2021).

Murphy Quinlan, M., Walker, A. M., Davies, C. J., Mound, J. E., Müller, T., & Harvey, J. (2021). The conductive cooling of planetesimals with temperature-dependent properties. Journal of Geophysical Research: Planets, 126, e2020JE006726. https://doi.org/10.1029/2020JE006726

Parameters

filepath (str, optional) – Absolute or relative path to save the json parameter file to.

pytesimal.load_plot_save.plot_coolingrate_history(dT_by_dt, dT_by_dt_core, timestep, maxtime, ax=None, fig=None, savefile=None, fig_w=8, fig_h=6, show=True)

Generate a heat map of cooling rate vs time.

Generate a heat map of cooling rate vs time, with the colormap showing variation in cooling rate.

Input cooling rate in a n-steps by n-depths array dT_by_dt for the mantle and n-steps by n-depths array dT_by_dt_core for the core. timestep is the length of a timestep and maxtime is the maximum time (both in seconds).

Optional arguments fig and ax can be set to plot on existing matplotlib figure and axis objects. Passing a string via outfile causes the figure to be saved as an image in a file.

pytesimal.load_plot_save.plot_temperature_history(temperatures, coretemp, timestep, maxtime, ax=None, fig=None, savefile=None, fig_w=8, fig_h=6, show=True)

Generate a heat map of depth vs time; colormap shows variation in temp.

Input temperature in a n-steps by n-depths array temperatures for the mantle and n-steps by n-depths array coretemp for the core. timestep is the length of a timestep and maxtime is the maximum time (both in seconds).

Optional arguments fig and ax can be set to plot on existing matplotlib figure and axis objects. Passing a string via outfile causes the figure to be saved as an image in a file.

pytesimal.load_plot_save.read_datafile(filepath)

Read the contents of a model run into numpy arrays.

Reads the content of the numpy ‘npz’ data file representing a model run and returns arrays of the mantle temperature, core temperature, and cooling rates of the mantle and core.

Parameters

filepath (str) – Location of .npz data file, including file name and npz suffix.

Returns

  • temperatures (numpy.ndarray) – Array filled with mantle temperatures, in K.

  • coretemp (numpy.ndarray) – Array filled with core temperatures, in K.

  • dT_by_dt (numpy.ndarray) – Array filled with mantle cooling rates, in K/dt.

  • dT_by_dt_core (numpy.ndarray) – Array filled with core cooling rates, in K/dt.

pytesimal.load_plot_save.read_datafile_with_latent(filepath)

Read contents of a model run into numpy arrays, including latent heat.

Reads the content of the numpy ‘npz’ data file representing a model run and returns arrays of the mantle temperature, core temperature, cooling rates of the mantle and core, and the number of timesteps the core was crystallising for.

Parameters

filepath (str) – Location of .npz data file, including file name and npz suffix.

Returns

  • temperatures (numpy.ndarray) – Array filled with mantle temperatures, in K.

  • coretemp (numpy.ndarray) – Array filled with core temperatures, in K.

  • dT_by_dt (numpy.ndarray) – Array filled with mantle cooling rates, in K/dt.

  • dT_by_dt_core (numpy.ndarray) – Array filled with core cooling rates, in K/dt.

  • latent_array (numpy.ndarray) – Array filled with latent heat of the core, as it crystallises, J kg^-1.

pytesimal.load_plot_save.save_params_and_results(result_filename, run_ID, folder, timestep, r_planet, core_size_factor, reg_fraction, max_time, temp_core_melting, mantle_heat_cap_value, mantle_density_value, mantle_conductivity_value, core_cp, core_density, temp_init, temp_surface, core_temp_init, core_latent_heat, kappa_reg, dr, cond_constant, density_constant, heat_cap_constant, time_core_frozen, fully_frozen, meteorite_results='None given', latent_list_len=0)

Save parameters and results from model run to a json file.

Save the parameters used, core crystallisation timing, and meteorite results to a json file. The resulting file can also be read in as a parameter file to reproduce the same modelling run.

Meteorite results can be excluded, or can be passed in as a string, value, or dictionary of results.

Parameters
  • result_filename (str) – Filename without file suffix; .txt will be appended when the file is saved.

  • run_ID (str) – Identifier for the specific model run.

  • folder (str) – Absolute path to directory where file is to be saved. Existence of the directory can be checked with the check_folder_exists() function.

  • timestep (float) – The timestep used in numerical method, in s.

  • r_planet (float) – The radius of the planet in m.

  • core_size_factor (float) – The core size as a fraction of the total planet radius.

  • reg_fraction (float) – The regolith thickness as a fraction of the total planet radius.

  • max_time (float) – The total run-time of the model, in millions of years (Myr).

  • temp_core_melting (float) – The melting temperature of the core, in K.

  • mantle_heat_cap_value (float) – The heat capacity of mantle material, in J kg^-1 K^-1.

  • mantle_density_value (float) – The density of mantle material, in kg m^-3.

  • mantle_conductivity_value (float) – The conductivity of the mantle, in W m^-1 K^-1.

  • core_cp (float) – The heat capacity of the core, in J kg^-1 K^-1.

  • core_density (float) – The density of the core, in kg m^-3.

  • temp_init (float, list, numpy array) – The initial temperature of the body, in K.

  • temp_surface (float) – The surface temperature of the planet, in K.

  • core_temp_init (float) – The initial temperature of the core, in K.

  • core_latent_heat (float) – The latent heat of crystallisation of the core, in J kg^-1.

  • kappa_reg (float) – The regolith constant diffusivity, m^2 s^-1.

  • dr (float) – The radial step used in the numerical model, in m.

  • cond_constant (str) – Flag of y or n to specify if mantle conductivity is constant.

  • density_constant (str) – Flag of y or n to specify if mantle density is constant.

  • heat_cap_constant (str) – Flag of y or n to specify if mantle heat capacity is constant.

  • time_core_frozen (float) – Time when the core begins to freeze in Myr.

  • fully_frozen (float) – Time when the core finishes freezing, in Myr.

  • meteorite_results (dict, str, float, list, optional) – Depth and timing results of meteorites, can be passed as a dictionary of results for different samples, as a list of results, or as a single result in the form of a float or string.

  • latent_list_len (float, optional) – The length of the latent heat list, needed for further analysis of core crystallisation duration at a later point.

Returns

Return type

File is saved to folder/result_filename.txt in the json format.

pytesimal.load_plot_save.save_result_arrays(result_filename, folder, mantle_temperature_array, core_temperature_array, mantle_cooling_rates, core_cooling_rates, latent=[])

Save results as a compressed Numpy array (npz).

Result arrays of temperatures and cooling rates for both the mantle and the core (numpy arrays) are saved to a specified file.

Parameters
  • result_filename (str) – Filename without file suffix; .npz will be appended when the file is saved.

  • folder (str) – Absolute path to directory where file is to be saved. Existence of the directory can be checked with the check_folder_exists() function.

  • mantle_temperature_array (numpy.ndarray) – Temperatures in the mantle for all radii through time, in K.

  • core_temperature_array (numpy.ndarray) – Temperatures in the core through time, in K.

  • mantle_cooling_rates (numpy.ndarray) – Cooling rates in the mantle for all radii through time, in K/dt.

  • core_cooling_rates (numpy.ndarray) – Cooling rates in the core through time, in K/dt.

  • latent (list, optional) – List of latent heat values for the core; needed to calculate timing of core crystallisation, in J kg^-1.

Returns

  • File is saved to folder/result_filename.npz in the compressed Numpy

  • array format.

pytesimal.load_plot_save.two_in_one(fig_w, fig_h, temperatures, coretemp, dT_by_dt, dT_by_dt_core, savefile=None, timestep=100000000000.0)

Return a heat map of depth vs time; colormap shows variation in temp.

Change save=”n” to save=”y” when function is called to produce a png image named after the data filename.

pytesimal.mantle_properties module

Define mantle properties as constant or temperature-dependent.

Set mantle conductivity, density and heat capacity as constant values or define them as functions of temperature. The MantleProperties class has methods which define constant values for mantle properties, which can then individually be overridden by the VariableConductivity, VariableDensity and VariableHeatCapacity subclasses. The VariableConductivity subclass also contains a method to calculate the derivative of the conductivity with respect to time.

The functions used for variable properties are based on experimental results and mineral physics theory, discussed in Murphy Quinlan et al. (2021) and the references therein.

class pytesimal.mantle_properties.MantleProperties(rho=3341.0, cp=819.0, k=3.0)

Bases: object

Mantle properties class to define thermal properties of mantle material.

The value of conductivity, heat capacity or density can be called using the get methods, and can be changed using the set methods. They can be overridden with temperature-dependent functions using subclasses for each individual property.

Temperature and pressure are optional arguments for the get methods; these are not used when the values are temperature-independent but allow for easy insertion of temperature or pressure dependent functions into pre-existing code with minimal changes.

rho

The density of the mantle material (constant), in kg m^-3.

Type

float, default 3341.0

cp

The heat capacity of mantle material (constant), in J kg^-1 K^-1.

Type

float, default 819.0

k

The conductivity of mantle material (constant), in W m^-1 K^-1.

Type

float, default 3.0

__init__(rho=3341.0, cp=819.0, k=3.0)

Initialise mantle properties.

getcp(T=295, P=0.1)

Get heat capacity.

getdkdT(T=295, P=0.1)

Get gradient of conductivity.

getk(T=295, P=0.1)

Get conductivity.

getkappa()

Get diffusivity.

getrho(T=295, P=0.1)

Get density.

setcp(value)

Set heat capacity.

setk(value)

Set conductivity.

setrho(value)

Set density.

class pytesimal.mantle_properties.VariableConductivity(rho=3341.0, cp=819.0, k=3.0)

Bases: pytesimal.mantle_properties.MantleProperties

Make conductivity T-dependent.

getdkdT(T=295)

Get derivative of conductivity with respect to temperature.

getk(T=295, P=0.1)

Get conductivity.

class pytesimal.mantle_properties.VariableDensity(rho=3341.0, cp=819.0, k=3.0)

Bases: pytesimal.mantle_properties.MantleProperties

Make density T-dependent.

getrho(T=295.0)

Get density.

class pytesimal.mantle_properties.VariableHeatCapacity(rho=3341.0, cp=819.0, k=3.0)

Bases: pytesimal.mantle_properties.MantleProperties

Make heat capacity T-dependent.

getcp(T=295)

Get heat capacity.

pytesimal.mantle_properties.set_up_mantle_properties(cond_constant='y', density_constant='y', heat_cap_constant='y', mantle_density=3341.0, mantle_heat_capacity=819.0, mantle_conductivity=3.0)

Define mantle properties quickly

A quick set-up function that can read parameters and flags from loaded parameter files to set up mantle properties.

Parameters
  • cond_constant (str, default 'y') – Flag to define if conductivity is constant in temperature or not. Default ‘y’ results in constant conductivity, while any other string produces variable conductivity.

  • density_constant (str, default 'y') – Flag to define if density is constant in temperature or not. Default ‘y’ results in constant density, while any other string produces variable density.

  • heat_cap_constant (str, default 'y') – Flag to define if heat capacity is constant in temperature or not. Default ‘y’ results in constant heat capacity, while any other string produces variable heat capacity.

  • mantle_density (float, default 3341.0) – Constant value for mantle density, in kg m^-3.

  • mantle_heat_capacity (float, default 819.0) – Constant value for mantle heat capacity, in J kg^-1 K^-1.

  • mantle_conductivity (float) – Constant value for mantle conductivity, in W m^-1 K^-1.

Returns

  • conductivity (object) – Conductivity object, with constant or temperature dependent value

  • heat_capacity (object) – Heat capacity object, with constant or temperature dependent value

  • density (object) – Density object, with constant or temperature dependent value

pytesimal.numerical_methods module

Forward-Time Central-Space (FTCS) discretisation for the mantle.

Set up boundary conditions, calculate the heat extracted across the core-mantle boundary in a timestep, and numerically discretise the conductive cooling of the mantle of a planetesimal. This module also contains functions to calculate the thermal diffusivity of a material form the thermal conductivity, heat capacity and the density, and to check whether the diffusivity, timestep and radial discretisation meet Von Neumann stability criteria.

class pytesimal.numerical_methods.EnergyExtractedAcrossCMB(outer_r, timestep, radius_step)

Bases: object

Class to calculate the energy extracted across the cmb in one timestep.

outer_r

Core radius, in m.

Type

float

timestep

Time over which heat is extracted, in s.

Type

float

radius_step

Radial step for numerical discretisation, in m.

Type

float

__init__(outer_r, timestep, radius_step)

Initialize self. See help(type(self)) for accurate signature.

energy_extracted(mantle_temperatures, i, k)

Calculate energy extracted in one timestep

power(mantle_temperatures, i, k)

Calculate heat (power) extracted in one timestep

pytesimal.numerical_methods.calculate_diffusivity(conductivity, heat_capacity, density)

Calculate diffusivity from conductivity, heat capacity and density.

Returns a value for diffusivity at a certain temperature, given float values for conductivity, heat capacity and density.

Parameters
  • conductivity (float) – Thermal conductivity of the material, in W m^-1 K^-1.

  • heat_capacity (float) – Heat capacity of the material, in J kg^-1 K^-1.

  • density (float) – Density of the material, in kg m^-3.

Returns

diffusivity – The calculated diffusivity, in m^2 s^-1.

Return type

float

pytesimal.numerical_methods.check_stability(max_diffusivity, timestep, dr)

Check adherence to Von Neumann stability criteria.

Use the maximum diffusivity of the system to return the most restrictive criteria. Diffusivity can be calculated using the calculate_diffusivity function, with max diffusivity where conductivity is maximised and density and heat capacity are minimised.

Parameters
  • max_diffusivity (float) – The highest diffusivity of the system to impose the most restrictive conditions, , in m^2 s^-1.

  • timestep (float) – The timestep used for the numerical scheme, in s.

  • dr (float) – The radial step used for the numerical scheme, in m.

Returns

result – Boolean, True if parameters pass stability criteria and false if they fail.

Return type

bool

pytesimal.numerical_methods.cmb_dirichlet_bc(temperatures, core_boundary_temperature, i)

Set a fixed temperature boundary condition at the base of the mantle.

Parameters
  • temperatures (numpy.ndarray) – Numpy array of mantle temperatures to apply condition to, in K.

  • core_boundary_temperature (float) – The temperature at the core mantle boundary, in K.

  • i (int) – Index along time axis where condition is to be set.

Returns

temperatures – Temperature array with condition applied, in K.

Return type

numpy.ndarray

pytesimal.numerical_methods.cmb_neumann_bc(temperatures, core_boundary_temperature, i)

Set a zero flux boundary condition at the base of the mantle

Note that core radius must be set to zero for this to approximate the analytical solution of a conductively cooling sphere or to model an undifferentiated meteorite parent body. Sets zero flux at base of the mantle by approximating dT/dr using forward differences and finding the necessary temperature. See eq. 6.31 of http://folk.ntnu.no/leifh/teaching/tkt4140/._main056.html

Parameters
  • temperatures (numpy.ndarray) – Numpy array of mantle temperatures to apply condition to, in K.

  • core_boundary_temperature (float) – The temperature at the core mantle boundary; this is not used by this boundary condition but inclusion allows functions to be easily swapped. In K.

  • i (int) – Index along time axis where condition is to be set.

Returns

temperatures – Temperature array with condition applied, in K.

Return type

numpy.ndarray

pytesimal.numerical_methods.discretisation(core_values, latent, temp_init, core_temp_init, top_mantle_bc, bottom_mantle_bc, temp_surface, temperatures, dr, coretemp_array, timestep, r_core, radii, times, where_regolith, kappa_reg, cond, heatcap, dens, non_lin_term='y')

Finite difference solver with variable k.

Uses variable heat capacity, conductivity, density as required.

Uses diffusivity for regolith layer.

Parameters
  • core_values (core object) – An object that represents the state of the layer inside the current layer (normally a metallic core). The object must provide one method and two attributes. The method extract_heat(power, timestep) is called on each timestep and represents the amount of heat that is lost from the the inner layer to the present layer (power, in W) over an amount of time (timestep, in s). The attribute temperature gives the temperature at the top of the inner layer and this is used (after calling extract_heat) as input to the basal boundary condition calculation. The attribute latent reports any latent heat released by freezing and this is not explicitly used in the evaluation of mantle temperatures.

  • latent (list) – Empty list (unless coupling two models) of latent heat extracted from the core.

  • temp_init (float, numpy.ndarray) – The initial temperature (in K) of the mantle with float implying initial homogeneous temperature distribution.

  • core_temp_init (float, numpy.ndarray) – Initial temperature of the core; current core object is isothermal so only accepts float but more complex core models could track the temperature distribution in the core. In K.

  • top_mantle_bc (callable) – Calleable function that defines the boundary condition at the planetesimal surface. The calling signature is top_mantle_bc(temperatures, surface_temperature, timestep_index) where temperatures is the temperatures array to be updated with the boundary condition, core_boundary_temperature is the temperature (that may be involved in the calculation) and timestep_index is the column index of the current timestep. The function must return an updated temperatures array. See surface_dirichlet_bc for an example.

  • bottom_mantle_bc (callable) – Calleable function that defines the boundary condition at the base of the planetesimal mantle. The calling signature is bottom_mantle_bc(temperatures, core_boundary_temperature, timestep_index) where temperatures is the temperatures array to be updated with the boundary condition, core_boundary_temperature is the temperature (that may be involved in the calculation) and timestep_index is the column index of the current timestep. The function must return an updated temperatures array. See cmb_neumann_bc for an example.

  • temp_surface (float) – Temperature at the surface of the planetesimal, in K.

  • temperatures (numpy.ndarray) – Numpy array to fill with mantle temperatures in K.

  • dr (float) – Radial step for numerical discretisation, in m.

  • coretemp_array (numpy.ndarray) – Numpy array to fill with core temperatures

  • timestep (float) – Timestep for numerical discretisation, in s.

  • r_core (float) – Radius of the core in m.

  • radii (numpy.ndarray) – Numpy array of radii values in the mantle, with spacing defined by dr, in m.

  • times (numpy.ndarray) – Numpy array of time values in s, up to the maximum time, with spacing controlled by timestep, in s.

  • where_regolith (numpy.ndarray) – Boolean array recording presence of regolith.

  • kappa_reg (float) – Constant diffusivity of the regolith, in m^2 s^-1.

  • cond (callable) – Callable function or method that defines the mantle conductivity. The calling signature is cond.getk(temperatures[radial_index, timestep_index]), where temperatures is the temperatures array, radial_index is the row index of the radius, and timestep_index is the column index of the timestep, that define the value in temperatures at which conductivity should be evaluated. The function must return a value for conductivity in in W m^-1 K^-1.

  • heatcap (callable) – Callable function or method that defines the mantle heat capacity. The calling signature is heatcap.getcp(temperatures[radial_index, timestep_index]), where temperatures is the temperatures array, radial_index is the row index of the radius, and timestep_index is the column index of the timestep, that define the value in temperatures at which heat capacity should be evaluated. The function must return a value for heat capacity in J kg^-1 K^-1.

  • dens (callable) – Callable function or method that defines the mantle density. The calling signature is dens.getrho(temperatures[radial_index, timestep_index]), where temperatures is the temperatures array, radial_index is the row index of the radius, and timestep_index is the column index of the timestep, that define the value in temperatures at which heat capacity should be evaluated. The function must return a value for density in kg m^-3.

  • non_lin_term (str, default ‘y’) – Flag to switch off the non-linear term when temperature-dependent conductivity is being used.

Returns

  • temperatures (numpy.ndarray) – Array filled with mantle temperatures, in K.

  • coretemp (numpy.ndarray) – Array filled with core temperatures, in K.

  • latent (list) – List of latent heat values during core crystallisation, in J kg^-1.

pytesimal.numerical_methods.surface_dirichlet_bc(temperatures, temp_surface, i)

Set a fixed temperature boundary condition at the planetesimal’s surface.

Parameters
  • temperatures (numpy.ndarray) – Numpy array of mantle temperatures to apply condition to, in K.

  • temp_surface (float) – The temperature at the surface boundary, in K.

  • i (int) – Index along time axis where condition is to be set.

Returns

temperatures – Temperature array with condition applied, in K.

Return type

numpy.ndarray

pytesimal.quick_workflow module

Set up complete runs with a single function call and a parameter file.

This module provides a simple workflow function which follows a basic default workflow, with parameters set by an input parameter file, and results saved as a json file and as compressed numpy arrays. These results files can then be loaded and analysed, with meteorite results calculated and the results plotted.

The function takes two arguments, filename which is the name of a parameters file to be loaded, without the extension (the extension of the file must be .txt). The location of this input file is given by folder_path which should be the relative or absolute path to the location of the parameters file. Within this parameters file, the “folder” field defines the path of the directory where results

pytesimal.quick_workflow.workflow(filename, folder_path)

Run model in full with parameters set by an input file.

Saves a results file (json format, .txt) and a results array file (.npz) to the folder specified in the parameter file. If you want the results to save to the same folder that the parameter file is in, ensure that the field “folder” in the json parameter file is the same as folder_path.

Results files will be saved under the filename with _results and the appropriate file extension appended.

The json (.txt) file contains a list of parameter names and values, while the array file (.npz) contains four different arrays: mantle_temperature_array - an array of mantle temperatures in K, core_temperature_array - an array of core temperatures in K, mantle_cooling_rates - an array of mantle cooling rates in K/dt, core_cooling_rates - an array of core cooling rates in K/dt.

Parameters
  • filename (str) – The filename of the parameters file to read in, without file extension. File must be in json format with a .txt extension. It’s recommended to generate an example parameters file using the pytesimal.load_plot_save.make_default_param_file(filepath) function and edit the default settings in this file

  • folder_path (str) – The absolute path to the directory that holds the parameters file. If the “folder” field of the parameters file == folder_path, the results file will be saved alongside the parameters file.

pytesimal.setup_functions module

Define the geometry of planetesimal and set up required empty arrays.

This module allows the user to set up a basic geometry based on parameters instead of manually defining ‘numpy.ndarrays’.

pytesimal.setup_functions.set_up(timestep=100000000000.0, r_planet=250000.0, core_size_factor=0.5, reg_fraction=0.032, max_time=400.0, dr=1000.0)

Define the geometry and set up corresponding arrays.

Parameters
  • timestep (float, default 1e11) – A timestep for the numerical discretisation, in s

  • r_planet (float, default 250000.0) – The radius of the planetesimal, in m

  • core_size_factor (float, < 1.0, default 0.5) – The core radius expressed as a fraction of r_planet

  • reg_fraction (float, <1.0, default 0.032) – The core thickness expressed as a fraction of r_planet

  • max_time (float, default 400.0) – Total time for model to run, in millions of years (Myr)

  • dr (float, default 1000.0) – Radial step for the numerical discretisation, in m

Returns

  • r_core (float,) – Radius of the core, in m

  • radii (numpy.ndarray) – Numpy array of radius values in m for the mantle, with spacing defined by dr

  • core_radii (numpy.ndarray) – Numpy array of radius values in m for the core, with spacing defined by dr

  • reg_thickness (float) – Regolith thickness in m

  • where_regolith (numpy.ndarray) – Boolean array with location of regolith

  • times (numpy.ndarray) – Numpy array starting at 0 and going to 400 Myr, with timestep controlling the spacing

  • mantle_temperature_array (numpy.ndarray) – Numpy array of zeros to be filled with mantle temperatures in K

  • core_temperature_array (numpy.ndarray) – Numpy array of zeros to be filled with core temperatures in K