diff --git a/examples/default/init_coupler.toml b/examples/default/init_coupler.toml index 7685801a..58fcd458 100644 --- a/examples/default/init_coupler.toml +++ b/examples/default/init_coupler.toml @@ -90,7 +90,7 @@ author = "Harrison Nicholls, Tim Lichtenberg" mass = 1.0 # M_sun radius = 1.0 # R_sun Teff = 5772.0 # K - Lbol = 1.0 # L_sun + lum_now = 1.0 # L_sun omega = 50.0 # rotation percentile age_now = 4.567 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age diff --git a/examples/hd63433d/init_coupler.toml b/examples/hd63433d/init_coupler.toml index 382e42cc..f76ae534 100644 --- a/examples/hd63433d/init_coupler.toml +++ b/examples/hd63433d/init_coupler.toml @@ -90,7 +90,7 @@ author = "Harrison Nicholls, Tim Lichtenberg" mass = 0.99 # M_sun radius = 0.912 # R_sun Teff = 5640.0 # K - Lbol = 0.753 # L_sun + lum_now = 0.753 # L_sun omega = 50.0 # rotation percentile age_now = 4.140 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age diff --git a/input/default.toml b/input/default.toml index 7685801a..84a46248 100644 --- a/input/default.toml +++ b/input/default.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 1.0 # M_sun radius = 1.0 # R_sun Teff = 5772.0 # K - Lbol = 1.0 # L_sun + mass = 1.0 # M_sun + lum_now = 1.0 # L_sun, Luminosity as observed today. omega = 50.0 # rotation percentile age_now = 4.567 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "spada" # evolution tracks: spada | baraffe diff --git a/input/dummy.toml b/input/dummy.toml index fc43b310..0cab8814 100644 --- a/input/dummy.toml +++ b/input/dummy.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 1.0 # M_sun radius = 1.0 # R_sun Teff = 5772.0 # K - Lbol = 1.0 # L_sun + mass = 1.0 # M_sun + lum_now = 1.0 # L_sun, Luminosity as observed today. omega = 50.0 # rotation percentile age_now = 4.567 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "baraffe" # evolution tracks: spada | baraffe diff --git a/input/dummy_star.toml b/input/dummy_star.toml new file mode 100644 index 00000000..f39c6b62 --- /dev/null +++ b/input/dummy_star.toml @@ -0,0 +1,233 @@ +# PROTEUS configuration file (version 2.0) + +# Root tables should be physical, with the exception of "params" +# Software related options should go within the appropriate physical table + +# The general structure is: +# [root] metadata +# [params] parameters for code execution, output files, time-stepping, convergence +# [star] stellar parameters, model selection +# [orbit] planetary orbital parameters +# [struct] planetary structure (mass, radius) +# [atmos] atmosphere parameters, model selection +# [escape] escape parameters, model selection +# [interior] magma ocean model selection and parameters +# [outgas] outgassing parameters (fO2) and included volatiles +# [delivery] initial volatile inventory, and delivery model selection + +# ---------------------------------------------------- +# Metadata +version = "2.0" +author = "Harrison Nicholls, Tim Lichtenberg" + +# ---------------------------------------------------- +# Parameters +[params] + # output files + [params.out] + path = "dummy_star" + logging = "DEBUG" + plot_mod = 0 # Plotting frequency, 0: wait until completion | n: every n iterations + plot_fmt = "png" # Plotting image file format, "png" or "pdf" recommended + + # time-stepping + [params.dt] + minimum = 1e2 # yr, minimum time-step + maximum = 3e7 # yr, maximum time-step + initial = 1e2 # yr, inital step size + starspec = 1e7 # yr, interval to re-calculate the stellar spectrum + starinst = 1e1 # yr, interval to re-calculate the instellation + method = "adaptive" # proportional | adaptive | maximum + + [params.dt.proportional] + propconst = 52.0 # Proportionality constant + + [params.dt.adaptive] + atol = 0.02 # Step size atol + rtol = 0.10 # Step size rtol + + # termination criteria + [params.stop] + + # required number of iterations + [params.stop.iters] + enabled = true + minimum = 3 + maximum = 9000 + + # required time constraints + [params.stop.time] + enabled = true + minimum = 1.0e3 # yr, model will certainly run to t > minimum + maximum = 4.567e+9 # yr, model will terminate when t > maximum + + # solidification + [params.stop.solid] + enabled = false + phi_crit = 0.005 # non-dim., model will terminate when global melt fraction < phi_crit + + # radiative equilibrium + [params.stop.radeqm] + enabled = false + F_crit = 0.2 # W m-2, model will terminate when |F_atm| < F_crit + + # steady state + [params.stop.steady] + enabled = false + F_crit = 0.8 # Maximum absolute value of F_atm allowed for convergence + dprel = 1.0e-9 # Percentage change in melt fraction over time (dp/p)/dt*100 + + [params.stop.escape] + enabled = true + mass_frac = 3e-4 # Stop when atm_mass < this frac of initial mass + + +# ---------------------------------------------------- +# Star +[star] + + # Physical parameters + radius = 1.0 # R_sun + Teff = 5772.0 # K + mass = 1.0 # M_sun + lum_now = 1.0 # L_sun, Luminosity as observed today. + omega = 50.0 # rotation percentile + age_now = 4.567 # Gyr, current age of star used for scaling + age_ini = 0.100 # Gyr, model initialisation/start age + + module = "dummy" + [star.mors] + tracks = "baraffe" # evolution tracks: spada | baraffe + spec = "stellar_spectra/Named/sun.txt" # stellar spectrum + +# Orbital system +[orbit] + semimajoraxis = 1.0 # AU + eccentricity = 0.0 # dimensionless + zenith_angle = 48.19 # degrees + s0_factor = 0.375 # dimensionless + + # No module specifically for tides / orbital dynamics + module = "none" + + +# Planetary structure - physics table +[struct] + mass = 1.0 # M_earth + radius = 1.0 # R_earth + corefrac = 0.55 # non-dim., radius fraction + + # No module for specifically for internal structure + module = "none" + +# Atmosphere - physics table +[atmos_clim] + prevent_warming = true # do not allow the planet to heat up + surface_d = 0.01 # m, conductive skin thickness + surface_k = 2.0 # W m-1 K-1, conductive skin thermal conductivity + cloud_enabled = false # enable water cloud radiative effects + cloud_alpha = 0.0 # condensate retention fraction (1 -> fully retained) + surf_state = "fixed" # surface scheme: "mixed_layer" | "fixed" | "skin" + surf_albedo = 0.1 # path to file ("string") or grey quantity (float) + albedo_pl = 0.1 # Bond albedo (scattering) + rayleigh = true # enable rayleigh scattering + tmp_minimum = 0.5 # temperature floor on solver + tmp_maximum = 5000.0 # temperature ceiling on solver + + module = "dummy" # Which atmosphere module to use + + [atmos_clim.agni] + p_top = 1.0e-5 # bar, top of atmosphere grid pressure + spectral_group = "Frostflow" # which gas opacities to include + spectral_bands = "48" # how many spectral bands? + num_levels = 40 # Number of atmospheric grid levels + chemistry = "none" # "none" | "eq" | "kin" + + [atmos_clim.janus] + p_top = 1.0e-6 # bar, top of atmosphere grid pressure + spectral_group = "Frostflow" # which gas opacities to include + spectral_bands = "48" # how many spectral bands? + F_atm_bc = 0 # measure outgoing flux at: (0) TOA | (1) Surface + num_levels = 40 # Number of atmospheric grid levels + tropopause = "none" # none | skin | dynamic + + [atmos_clim.dummy] + gamma = 0.7 # atmosphere opacity between 0 and 1 + +# Volatile escape - physics table +[escape] + + module = "dummy" # Which escape module to use + + [escape.zephyrus] + some_parameter = "some_value" + + [escape.dummy] + rate = 2.0e4 # Bulk unfractionated escape rate [kg s-1] + +# Interior - physics table +[interior] + grain_size = 0.1 # crystal settling grain size [m] + F_initial = 1e6 # Initial heat flux guess [W m-2] + + module = "dummy" # Which interior module to use + + [interior.spider] + num_levels = 220 # Number of SPIDER grid levels + mixing_length = 2 # Mixing length parameterization + tolerance = 1.0e-10 # solver tolerance + tsurf_atol = 30.0 # tsurf_poststep_change + tsurf_rtol = 0.02 # tsurf_poststep_change_frac + ini_entropy = 2700.0 # Surface entropy conditions [J K-1 kg-1] + ini_dsdr = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1] + + [interior.aragog] + some_parameter = "some_value" + +# Outgassing - physics table +[outgas] + fO2_shift_IW = 2 # log10(ΔIW), atmosphere/interior boundary oxidation state + + module = "calliope" # Which outgassing module to use + + [outgas.calliope] + include_H2O = true # Include H2O compound + include_CO2 = true # Include CO2 compound + include_N2 = true # Include N2 compound + include_S2 = true # Include S2 compound + include_SO2 = true # Include SO2 compound + include_H2 = true # Include H2 compound + include_CH4 = true # Include CH4 compound + include_CO = true # Include CO compound + + [outgas.atmodeller] + some_parameter = "some_value" + +# Volatile delivery - physics table +[delivery] + + # [Settings here for accretion rate, etc.] + + # Which initial inventory to use? + initial = 'elements' # "elements" | "volatiles" + + # No module for accretion as of yet + module = "none" + + # Set initial volatile inventory by planetary element abundances + [delivery.elements] + CH_ratio = 1.0 # C/H ratio in mantle/atmosphere system + H_oceans = 6.0 # Hydrogen inventory in units of equivalent Earth oceans, by mass + N_ppmw = 2.0 # Nitrogen inventory in ppmw relative to mantle mass, by mass + S_ppmw = 200.0 # Sulfur inventory in ppmw relative to mass of melt + + # Set initial volatile inventory by partial pressures in atmosphere + [delivery.volatiles] + H2O = 0.0 # partial pressure of H2O + CO2 = 0.0 # partial pressure of CO2 + N2 = 0.0 # etc + S2 = 0.0 + SO2 = 0.0 + H2 = 0.0 + CH4 = 0.0 + CO = 0.0 diff --git a/input/escape.toml b/input/escape.toml index 0b4d1474..88be585e 100644 --- a/input/escape.toml +++ b/input/escape.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 1.0 # M_sun radius = 1.0 # R_sun Teff = 5772.0 # K - Lbol = 1.0 # L_sun + mass = 1.0 # M_sun + lum_now = 1.0 # L_sun omega = 50.0 # rotation percentile age_now = 4.567 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "spada" # evolution tracks: spada | baraffe diff --git a/input/hd63433d.toml b/input/hd63433d.toml index 382e42cc..624d0844 100644 --- a/input/hd63433d.toml +++ b/input/hd63433d.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 0.99 # M_sun radius = 0.912 # R_sun Teff = 5640.0 # K - Lbol = 0.753 # L_sun + mass = 0.99 # M_sun + lum_now = 0.753 # L_sun omega = 50.0 # rotation percentile age_now = 4.140 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "baraffe" # evolution tracks: spada | baraffe diff --git a/input/jgr_grid.toml b/input/jgr_grid.toml index 09fddbac..57ee7f91 100644 --- a/input/jgr_grid.toml +++ b/input/jgr_grid.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 1.0 # M_sun radius = 1.0 # R_sun Teff = 5772.0 # K - Lbol = 1.0 # L_sun + mass = 1.0 # M_sun + lum_now = 1.0 # L_sun omega = 60.0 # rotation percentile age_now = 4.567 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "spada" # evolution tracks: spada | baraffe diff --git a/input/ltt3780b.toml b/input/ltt3780b.toml index f808db48..49a103e5 100644 --- a/input/ltt3780b.toml +++ b/input/ltt3780b.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 0.401 # M_sun radius = 0.374 # R_sun Teff = 3331.0 # K - Lbol = 0.167 # L_sun + mass = 0.401 # M_sun + lum_now = 0.167 # L_sun omega = 20.0 # rotation percentile age_now = 3.100 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "spada" # evolution tracks: spada | baraffe diff --git a/input/t1c.toml b/input/t1c.toml index 8487eeb8..95de61cb 100644 --- a/input/t1c.toml +++ b/input/t1c.toml @@ -87,15 +87,14 @@ author = "Harrison Nicholls, Tim Lichtenberg" [star] # Physical parameters - mass = 0.1 # M_sun radius = 0.1192 # R_sun Teff = 2566.0 # K - Lbol = 0.000553 # L_sun + mass = 0.1 # M_sun + lum_now = 0.000553 # L_sun omega = 50.0 # rotation percentile age_now = 7.6 # Gyr, current age of star used for scaling age_ini = 0.100 # Gyr, model initialisation/start age - module = "mors" [star.mors] tracks = "spada" # evolution tracks: spada | baraffe diff --git a/src/proteus/config/_compatibility.py b/src/proteus/config/_compatibility.py index 50608c60..c9cba16c 100644 --- a/src/proteus/config/_compatibility.py +++ b/src/proteus/config/_compatibility.py @@ -50,7 +50,7 @@ def _phi_global(config: Config): 'star_mass': ('star', 'mass'), 'star_radius_modern': ('star', 'radius'), 'star_temperature_modern': ('star', 'Teff'), - 'star_luminosity_modern': ('star', 'Lbol'), + 'star_luminosity_modern': ('star', 'lum_now'), 'star_age_modern': ('star', 'age_now'), 'star_rot_pctle': ('star', 'omega'), 'star_spectrum': ('star', 'mors', 'spec'), diff --git a/src/proteus/config/_star.py b/src/proteus/config/_star.py index 4c847628..16147335 100644 --- a/src/proteus/config/_star.py +++ b/src/proteus/config/_star.py @@ -14,16 +14,16 @@ class Mors: @define class Star: """Stellar parameters, model selection""" - mass: float radius: float + mass: float Teff: float - Lbol: float omega: float + lum_now: float age_now: float age_ini: float module: str | None = field( - validator=validators.in_((None, 'mors')), + validator=validators.in_((None, 'mors', 'dummy')), converter=none_if_none, ) diff --git a/src/proteus/outgas/calliope.py b/src/proteus/outgas/calliope.py index 007da116..9da1b20e 100644 --- a/src/proteus/outgas/calliope.py +++ b/src/proteus/outgas/calliope.py @@ -81,13 +81,8 @@ def calc_surface_pressures(dirs:dict, config:Config, hf_row:dict): if e == "O": continue - # do not allow small values as they will prevent calliope from finding a solution - e_mass = hf_row[e + "_kg_total"] - if e_mass < 1.0e4: - e_mass = 0.0 - # save to dict - solvevol_target[e] = e_mass + solvevol_target[e] = hf_row[e + "_kg_total"] # get atmospheric compositison solvevol_result = equilibrium_atmosphere(solvevol_target, solvevol_inp) diff --git a/src/proteus/outgas/wrapper.py b/src/proteus/outgas/wrapper.py index 1d6e5c2f..7d9cc79c 100644 --- a/src/proteus/outgas/wrapper.py +++ b/src/proteus/outgas/wrapper.py @@ -5,7 +5,7 @@ from typing import TYPE_CHECKING from proteus.outgas.calliope import calc_surface_pressures, calc_target_masses -from proteus.utils.constants import volatile_species +from proteus.utils.constants import element_list, volatile_species if TYPE_CHECKING: from proteus.config import Config @@ -28,6 +28,14 @@ def run_outgassing(dirs:dict, config:Config, hf_row:dict): Run outgassing model to get new volatile surface pressures ''' + # Floating point errors can be problematic here. + # Ensure that zero mass values stay at zero by setting all element mass inventories + # which are less than 10 kg equal to zero. + for e in element_list: + if hf_row[e + "_kg_total"] < 10.0: + hf_row[e + "_kg_total"] = 0.0 + + # Run outgassing calculation if config.outgas.module == 'calliope': calc_surface_pressures(dirs, config, hf_row) else: diff --git a/src/proteus/plot/cpl_sflux.py b/src/proteus/plot/cpl_sflux.py index 5de2a57a..2759ff33 100644 --- a/src/proteus/plot/cpl_sflux.py +++ b/src/proteus/plot/cpl_sflux.py @@ -120,7 +120,7 @@ def plot_sflux(output_dir: str, wl_max: float = 6000.0, plot_format: str="pdf"): log.warning("Only one spectrum was found") ax.set_yscale("log") - ax.set_ylabel("Flux [erg s-1 cm-2 nm-1]") + ax.set_ylabel(r"Flux [erg / (s cm$^2$ nm)]") ax.set_xscale("log") ax.set_xlabel("Wavelength [nm]") diff --git a/src/proteus/plot/cpl_sflux_cross.py b/src/proteus/plot/cpl_sflux_cross.py index 95f8887f..c354c73f 100644 --- a/src/proteus/plot/cpl_sflux_cross.py +++ b/src/proteus/plot/cpl_sflux_cross.py @@ -102,7 +102,7 @@ def plot_sflux_cross( fig,ax = plt.subplots(1,1) ax.set_yscale("log") - ax.set_ylabel("Flux [erg s$^{-1}$ cm$^{-2}$ nm$^{-1}$]") + ax.set_ylabel("Flux [erg / (s cm$^2$ nm)]") ax.set_xscale("log") ax.set_xlabel("Time [yr]") diff --git a/src/proteus/proteus.py b/src/proteus/proteus.py index ee27ed6a..19376a74 100644 --- a/src/proteus/proteus.py +++ b/src/proteus/proteus.py @@ -2,7 +2,6 @@ import logging import os -import shutil import sys from datetime import datetime from pathlib import Path @@ -18,9 +17,17 @@ from proteus.escape.wrapper import RunEscape from proteus.interior import run_interior from proteus.outgas.wrapper import calc_target_elemental_inventories, run_outgassing +from proteus.star.wrapper import ( + get_new_spectrum, + init_star, + scale_spectrum_to_toa, + update_equilibrium_temperature, + update_instellation, + update_stellar_radius, + write_spectrum, +) from proteus.utils.constants import ( AU, - L_sun, M_earth, R_earth, const_G, @@ -28,7 +35,6 @@ volatile_species, ) from proteus.utils.coupler import ( - CalculateEqmTemperature, CreateHelpfileFromDict, CreateLockFile, ExtendHelpfile, @@ -57,11 +63,42 @@ class Proteus: def __init__(self, *, config_path: Path | str) -> None: + + # Read and parse configuration file self.config_path = config_path self.config = read_config_object(config_path) + # Setup directories dictionary self.init_directories() + # Helpfile variables for the current iteration + self.hf_row = None + + # Helpfile variables from all previous iterations + self.hf_all = None + + # Loop counters + self.loops = None + + # Model has finished? + self.finished = False + + # Default values for mors.spada cases + self.star_props = None + self.star_struct = None + + # Default values for mors.baraffe cases + self.baraffe_track = None + self.star_modern_fl = None + self.star_modern_wl = None + + # Stellar spectrum (wavelengths, fluxes) + self.star_wl = None + self.star_fl = None + + # Time at which star was last updated + self.sspec_prev = -np.inf # spectrum + self.sinst_prev = -np.inf # instellation and radius def init_directories(self): """Initialize directories dictionary""" @@ -79,7 +116,6 @@ def start(self, *, resume: bool = False): resume : bool If True, continue from previous simulation """ - import mors UpdateStatusfile(self.directories, 0) @@ -114,7 +150,7 @@ def start(self, *, resume: bool = False): log.info(" ") # Count iterations - loop_counter = { + self.loops = { "total": 0, # Total number of iters performed "total_min": 5, # Minimum number of total loops "total_loops": self.config.params.stop.iters.maximum, # Maximum number of total loops @@ -125,9 +161,6 @@ def start(self, *, resume: bool = False): "steady_check": 15, # Number of iterations to look backwards when checking steady state } - # Model has completed? - finished = False - # Config file paths config_path_backup = os.path.join(self.directories["output"], "init_coupler.toml") @@ -146,29 +179,24 @@ def start(self, *, resume: bool = False): with open(config_path_backup, 'w') as hdl: hdl.writelines(config_raw) - # No previous iterations to be stored. It is therefore important that the - # submodules do not try to read data from the 'past' iterations, since they do - # not yet exist. - hf_all = None - # Create an empty initial row for helpfile - hf_row = ZeroHelpfileRow() + self.hf_row = ZeroHelpfileRow() # Initial time - hf_row["Time"] = 0.0 - hf_row["age_star"] = self.config.star.age_ini * 1e9 + self.hf_row["Time"] = 0.0 + self.hf_row["age_star"] = self.config.star.age_ini * 1e9 # Initial guess for flux - hf_row["F_atm"] = self.config.interior.F_initial - hf_row["F_int"] = hf_row["F_atm"] - hf_row["T_eqm"] = 2000.0 + self.hf_row["F_atm"] = self.config.interior.F_initial + self.hf_row["F_int"] = self.hf_row["F_atm"] + self.hf_row["T_eqm"] = 2000.0 # Planet size conversion, and calculate mantle mass (= liquid + solid) - hf_row["M_int"] = self.config.struct.mass * M_earth - hf_row["R_int"] = self.config.struct.radius * R_earth - hf_row["gravity"] = const_G * hf_row["M_int"] / (hf_row["R_int"] ** 2.0) - hf_row["M_mantle"] = calculate_mantle_mass( - hf_row["R_int"], hf_row["M_int"], self.config.struct.corefrac + self.hf_row["M_int"] = self.config.struct.mass * M_earth + self.hf_row["R_int"] = self.config.struct.radius * R_earth + self.hf_row["gravity"] = const_G * self.hf_row["M_int"] / (self.hf_row["R_int"] ** 2.0) + self.hf_row["M_mantle"] = calculate_mantle_mass( + self.hf_row["R_int"], self.hf_row["M_int"], self.config.struct.corefrac ) # Store partial pressures and list of included volatiles @@ -185,9 +213,9 @@ def start(self, *, resume: bool = False): if include: inc_vols.append(s) - hf_row[s + "_bar"] = max(1.0e-30, float(pp_val)) + self.hf_row[s + "_bar"] = max(1.0e-30, float(pp_val)) else: - hf_row[s + "_bar"] = 0.0 + self.hf_row[s + "_bar"] = 0.0 log.info("Included volatiles: " + str(inc_vols)) else: @@ -202,74 +230,43 @@ def start(self, *, resume: bool = False): IC_INTERIOR = 2 # Read helpfile from disk - hf_all = ReadHelpfileFromCSV(self.directories["output"]) + self.hf_all = ReadHelpfileFromCSV(self.directories["output"]) # Check length - if len(hf_all) <= loop_counter["init_loops"] + 1: + if len(self.hf_all) <= self.loops["init_loops"] + 1: raise Exception("Simulation is too short to be resumed") # Get last row - hf_row = hf_all.iloc[-1].to_dict() + self.hf_row = self.hf_all.iloc[-1].to_dict() # Set loop counters - loop_counter["total"] = len(hf_all) - loop_counter["init"] = loop_counter["init_loops"] + 1 - - # Set instellation - S_0 = hf_row["F_ins"] - F_inst_prev = S_0 + self.loops["total"] = len(self.hf_all) + self.loops["init"] = self.loops["init_loops"] + 1 # Set volatile mass targets f solvevol_target = {} for e in element_list: if e == "O": continue - solvevol_target[e] = hf_row[e + "_kg_total"] + solvevol_target[e] = self.hf_row[e + "_kg_total"] # Download basic data download_sufficient_data(self.config) - # Handle stellar spectrum... - - # Store copy of modern spectrum in memory (1 AU) - sspec_prev = -np.inf - sinst_prev = -np.inf - star_modern_path = os.path.join(self.directories["fwl"], self.config.star.mors.spec) - shutil.copyfile(star_modern_path, os.path.join(self.directories["output"], "-1.sflux")) - - # Prepare stellar models - match self.config.star.mors.tracks: - case 'spada': - # load modern spectrum - star_struct_modern = mors.spec.Spectrum() - star_struct_modern.LoadTSV(star_modern_path) - star_struct_modern.CalcBandFluxes() - - # modern properties - star_props_modern = mors.synthesis.GetProperties( - self.config.star.mass, - self.config.star.omega, - self.config.star.age_now * 1000, - ) - - case 'baraffe': - modern_wl, modern_fl = mors.ModernSpectrumLoad( - star_modern_path, self.directories["output"] + "/-1.sflux" - ) - - baraffe = mors.BaraffeTrack(self.config.star.mass) + # Prepare star stuff + init_star(self) # Create lockfile keepalive_file = CreateLockFile(self.directories["output"]) # Main loop UpdateStatusfile(self.directories, 1) - while not finished: + while not self.finished: # New rows - if loop_counter["total"] > 0: + if self.loops["total"] > 0: # Create new row to hold the updated variables. This will be # overwritten by the routines below. - hf_row = hf_all.iloc[-1].to_dict() + self.hf_row = self.hf_all.iloc[-1].to_dict() log.info(" ") PrintSeparator() @@ -278,12 +275,12 @@ def start(self, *, resume: bool = False): log.info( "%1d/%1d %04d/%04d %1d/%1d" % ( - loop_counter["init"], - loop_counter["init_loops"], - loop_counter["total"], - loop_counter["total_loops"], - loop_counter["steady"], - loop_counter["steady_loops"], + self.loops["init"], + self.loops["init_loops"], + self.loops["total"], + self.loops["total_loops"], + self.loops["steady"], + self.loops["steady_loops"], ) ) @@ -291,7 +288,7 @@ def start(self, *, resume: bool = False): # Calculate time-averaged orbital separation (and convert from AU to metres) # https://physics.stackexchange.com/a/715749 - hf_row["separation"] = self.config.orbit.semimajoraxis * AU * \ + self.hf_row["separation"] = self.config.orbit.semimajoraxis * AU * \ (1 + 0.5 * self.config.orbit.eccentricity**2.0) @@ -301,11 +298,11 @@ def start(self, *, resume: bool = False): # Run interior model dt = run_interior(self.directories, self.config, - loop_counter, IC_INTERIOR, hf_all, hf_row) + self.loops, IC_INTERIOR, self.hf_all, self.hf_row) # Advance current time in main loop according to interior step - hf_row["Time"] += dt # in years - hf_row["age_star"] += dt # in years + self.hf_row["Time"] += dt # in years + self.hf_row["age_star"] += dt # in years ############### / INTERIOR @@ -315,98 +312,57 @@ def start(self, *, resume: bool = False): update_stellar_spectrum = False # Calculate new instellation and radius - if (abs(hf_row["Time"] - sinst_prev) > self.config.params.dt.starinst) or ( - loop_counter["total"] == 0 + if (abs(self.hf_row["Time"] - self.sinst_prev) > self.config.params.dt.starinst) or ( + self.loops["total"] == 0 ): - sinst_prev = hf_row["Time"] + self.sinst_prev = self.hf_row["Time"] - # Get previous instellation - if loop_counter["total"] > 0: - F_inst_prev = S_0 - else: - F_inst_prev = 0.0 - - log.info("Updating instellation and radius") - - match self.config.star.mors.tracks: - case 'spada': - hf_row["R_star"] = ( - mors.Value( - self.config.star.mass, - hf_row["age_star"] / 1e6, - "Rstar" - ) - * mors.const.Rsun - * 1.0e-2 - ) - S_0 = ( - mors.Value( - self.config.star.mass, - hf_row["age_star"] / 1e6, - "Lbol" - ) * L_sun / (4.0 * np.pi * hf_row["separation"]**2.0 ) - ) - - case 'baraffe': - hf_row["R_star"] = ( - baraffe.BaraffeStellarRadius(hf_row["age_star"]) - * mors.const.Rsun - * 1.0e-2 - ) - S_0 = baraffe.BaraffeSolarConstant( - hf_row["age_star"], hf_row["separation"]/AU - ) + # Update value for star's radius + log.info("Update stellar radius") + update_stellar_radius(self.hf_row, self.config, baraffe_track=self.baraffe_track) - # Calculate new eqm temperature - T_eqm_new = CalculateEqmTemperature( - S_0, self.config.orbit.s0_factor, self.config.atmos_clim.albedo_pl - ) + # Update value for instellation flux + log.info("Update instellation") + update_instellation(self.hf_row, self.config, baraffe_track=self.baraffe_track) - hf_row["F_ins"] = S_0 # instellation - hf_row["T_eqm"] = T_eqm_new - hf_row["T_skin"] = ( - T_eqm_new * (0.5**0.25) - ) # Assuming a grey stratosphere in radiative eqm (https://doi.org/10.5194/esd-7-697-2016) + # Calculate new eqm temperature + log.info("Update equilibrium temperature") + update_equilibrium_temperature(self.hf_row, self.config) - log.debug("Instellation change: %+.4e W m-2 (to 4dp)" % abs(S_0 - F_inst_prev)) + # Calculate new skin temperature + # Assuming a grey stratosphere in radiative eqm (https://doi.org/10.5194/esd-7-697-2016) + self.hf_row["T_skin"] = self.hf_row["T_eqm"] * (0.5**0.25) # Calculate a new (historical) stellar spectrum - if (abs(hf_row["Time"] - sspec_prev) > self.config.params.dt.starspec) or ( - loop_counter["total"] == 0 + if (abs(self.hf_row["Time"] - self.sspec_prev) > self.config.params.dt.starspec) or ( + self.loops["total"] == 0 ): - sspec_prev = hf_row["Time"] + self.sspec_prev = self.hf_row["Time"] update_stellar_spectrum = True + # Get the new spectrum using the appropriate module log.info("Updating stellar spectrum") - match self.config.star.mors.tracks: - case 'spada': - synthetic = mors.synthesis.CalcScaledSpectrumFromProps( - star_struct_modern, star_props_modern, hf_row["age_star"] / 1e6 - ) - fl = synthetic.fl # at 1 AU - wl = synthetic.wl - case 'baraffe': - fl = baraffe.BaraffeSpectrumCalc( - hf_row["age_star"], self.config["star_luminosity_modern"], modern_fl - ) - wl = modern_wl + self.star_wl, self.star_fl = get_new_spectrum( + + # Required variables + self.hf_row["age_star"], self.hf_row["R_star"], self.config, + + # Variables needed for mors.spada + star_struct_modern=self.star_struct, + star_props_modern=self.star_props, + + # Variables needed for mors.baraffe + baraffe_track=self.baraffe_track, + modern_wl=self.star_modern_wl, + modern_fl=self.star_modern_fl, + ) # Scale fluxes from 1 AU to TOA - fl *= (AU / hf_row["separation"]) ** 2.0 + self.star_fl = scale_spectrum_to_toa(self.star_fl, self.hf_row["separation"]) # Save spectrum to file - header = ( - "# WL(nm)\t Flux(ergs/cm**2/s/nm) Stellar flux at t_star = %.2e yr" - % hf_row["age_star"] - ) - np.savetxt( - os.path.join(self.directories["output"], "data", "%d.sflux" % hf_row["Time"]), - np.array([wl, fl]).T, - header=header, - comments="", - fmt="%.8e", - delimiter="\t", - ) + write_spectrum(self.star_wl, self.star_fl, + self.hf_row, self.directories["output"]) else: log.info("New spectrum not required at this time") @@ -414,8 +370,8 @@ def start(self, *, resume: bool = False): ############### / STELLAR FLUX MANAGEMENT ############### ESCAPE - if (loop_counter["total"] >= loop_counter["init_loops"]): - RunEscape(self.config, hf_row, dt) + if (self.loops["total"] >= self.loops["init_loops"]): + RunEscape(self.config, self.hf_row, dt) ############### / ESCAPE ############### OUTGASSING @@ -423,19 +379,21 @@ def start(self, *, resume: bool = False): # recalculate mass targets during init phase, since these will be adjusted # depending on the true melt fraction and T_magma found by SPIDER at runtime. - if loop_counter["init"] < loop_counter["init_loops"]: - calc_target_elemental_inventories(self.directories, self.config, hf_row) + if self.loops["init"] < self.loops["init_loops"]: + calc_target_elemental_inventories(self.directories, self.config, self.hf_row) # solve for atmosphere composition - run_outgassing(self.directories, self.config, hf_row) + run_outgassing(self.directories, self.config, self.hf_row) # Add atmosphere mass to interior mass, to get total planet mass - hf_row["M_planet"] = hf_row["M_int"] + hf_row["M_atm"] + self.hf_row["M_planet"] = self.hf_row["M_int"] + self.hf_row["M_atm"] ############### / OUTGASSING - ############### ATMOSPHERE SUB-LOOP - RunAtmosphere(self.config, self.directories, loop_counter, wl, fl, update_stellar_spectrum, hf_all, hf_row) + ############### ATMOSPHERE CLIMATE + RunAtmosphere(self.config, self.directories, self.loops, + self.star_wl, self.star_fl, update_stellar_spectrum, + self.hf_all, self.hf_row) ############### HOUSEKEEPING AND CONVERGENCE CHECK @@ -443,66 +401,66 @@ def start(self, *, resume: bool = False): # Update init loop counter # Next init iter - if loop_counter["init"] < loop_counter["init_loops"]: - loop_counter["init"] += 1 - hf_row["Time"] = 0.0 + if self.loops["init"] < self.loops["init_loops"]: + self.loops["init"] += 1 + self.hf_row["Time"] = 0.0 # Reset restart flag once SPIDER has correct heat flux - if loop_counter["total"] >= loop_counter["init_loops"]: + if self.loops["total"] >= self.loops["init_loops"]: IC_INTERIOR = 2 # Adjust total iteration counters - loop_counter["total"] += 1 + self.loops["total"] += 1 # Update full helpfile - if loop_counter["total"] > 1: + if self.loops["total"] > 1: # append row - hf_all = ExtendHelpfile(hf_all, hf_row) + self.hf_all = ExtendHelpfile(self.hf_all, self.hf_row) else: # first iter => generate new HF from dict - hf_all = CreateHelpfileFromDict(hf_row) + self.hf_all = CreateHelpfileFromDict(self.hf_row) # Write helpfile to disk - WriteHelpfileToCSV(self.directories["output"], hf_all) + WriteHelpfileToCSV(self.directories["output"], self.hf_all) # Print info to terminal and log file - PrintCurrentState(hf_row) + PrintCurrentState(self.hf_row) log.info("Check convergence criteria") # Stop simulation when planet is completely solidified - if self.config.params.stop.solid.enabled and (hf_row["Phi_global"] <= self.config.params.stop.solid.phi_crit): + if self.config.params.stop.solid.enabled and (self.hf_row["Phi_global"] <= self.config.params.stop.solid.phi_crit): UpdateStatusfile(self.directories, 10) log.info("") log.info("===> Planet solidified! <===") log.info("") - finished = True + self.finished = True # Stop simulation when flux is small if ( self.config.params.stop.radeqm.enabled - and (loop_counter["total"] > loop_counter["init_loops"] + 1) - and (abs(hf_row["F_atm"]) <= self.config.params.stop.radeqm.F_crit) + and (self.loops["total"] > self.loops["init_loops"] + 1) + and (abs(self.hf_row["F_atm"]) <= self.config.params.stop.radeqm.F_crit) ): UpdateStatusfile(self.directories, 14) log.info("") log.info("===> Planet no longer cooling! <===") log.info("") - finished = True + self.finished = True # Determine when the simulation enters a steady state if ( self.config.params.stop.steady.enabled - and (loop_counter["total"] > loop_counter["steady_check"] * 2 + 5) - and (loop_counter["steady"] < 0) + and (self.loops["total"] > self.loops["steady_check"] * 2 + 5) + and (self.loops["steady"] < 0) ): # How many iterations to look backwards - lb1 = -int(loop_counter["steady_check"]) + lb1 = -int(self.loops["steady_check"]) lb2 = -1 # Get data - arr_t = np.array(hf_all["Time"]) - arr_f = np.array(hf_all["F_atm"]) - arr_p = np.array(hf_all["Phi_global"]) + arr_t = np.array(self.hf_all["Time"]) + arr_f = np.array(self.hf_all["F_atm"]) + arr_p = np.array(self.hf_all["Phi_global"]) # Time samples t1 = arr_t[lb1] @@ -517,52 +475,53 @@ def start(self, *, resume: bool = False): phi_r = abs(phi_2 - phi_1) / (t2 - t1) # Stop when flux is small and melt fraction is unchanging - if (flx_m < self.config["steady_flux"]) and (phi_r < self.config["steady_dprel"]): + if (flx_m < self.config.params.stop.steady.F_crit) \ + and (phi_r < self.config.params.stop.steady.dprel): log.debug("Steady state declared") - loop_counter["steady"] = 0 + self.loops["steady"] = 0 # Steady-state handling - if loop_counter["steady"] >= 0: + if self.loops["steady"] >= 0: # Force the model to do N=steady_loops more iterations before stopping # This is to make absolutely sure that it has reached a steady state. - loop_counter["steady"] += 1 + self.loops["steady"] += 1 # Stop now? - if loop_counter["steady"] >= loop_counter["steady_loops"]: + if self.loops["steady"] >= self.loops["steady_loops"]: UpdateStatusfile(self.directories, 11) log.info("") log.info("===> Planet entered a steady state! <===") log.info("") - finished = True + self.finished = True # Atmosphere has escaped - if hf_row["M_atm"] <= self.config.params.stop.escape.mass_frac * hf_all.iloc[0]["M_atm"]: + if self.hf_row["M_atm"] <= self.config.params.stop.escape.mass_frac * self.hf_all.iloc[0]["M_atm"]: UpdateStatusfile(self.directories, 15) log.info("") log.info("===> Atmosphere has escaped! <===") log.info("") - finished = True + self.finished = True # Maximum time reached - if hf_row["Time"] >= self.config.params.stop.time.maximum * 1e9: + if self.hf_row["Time"] >= self.config.params.stop.time.maximum: UpdateStatusfile(self.directories, 13) log.info("") log.info("===> Target time reached! <===") log.info("") - finished = True + self.finished = True # Maximum loops reached - if loop_counter["total"] > loop_counter["total_loops"]: + if self.loops["total"] > self.loops["total_loops"]: UpdateStatusfile(self.directories, 12) log.info("") log.info("===> Maximum number of iterations reached! <===") log.info("") - finished = True + self.finished = True # Check if the minimum number of loops have been performed - if finished and (loop_counter["total"] < loop_counter["total_min"]): + if self.finished and (self.loops["total"] < self.loops["total_min"]): log.info("Minimum number of iterations not yet attained; continuing...") - finished = False + self.finished = False UpdateStatusfile(self.directories, 1) # Check if keepalive file has been removed - this means that the model should exit ASAP @@ -571,17 +530,17 @@ def start(self, *, resume: bool = False): log.info("") log.info("===> Model exit was requested! <===") log.info("") - finished = True + self.finished = True # Make plots if required and go to next iteration if ( (self.config.params.out.plot_mod > 0) - and (loop_counter["total"] % self.config.params.out.plot_mod == 0) - and (not finished) + and (self.loops["total"] % self.config.params.out.plot_mod == 0) + and (not self.finished) ): PrintHalfSeparator() log.info("Making plots") - UpdatePlots(hf_all, self.directories, self.config) + UpdatePlots(self.hf_all, self.directories, self.config) ############### / HOUSEKEEPING AND CONVERGENCE CHECK @@ -592,7 +551,7 @@ def start(self, *, resume: bool = False): safe_rm(keepalive_file) # Plot conditions at the end - UpdatePlots(hf_all, self.directories, self.config, end=True) + UpdatePlots(self.hf_all, self.directories, self.config, end=True) end_time = datetime.now() log.info("Simulation stopped at: " + end_time.strftime("%Y-%m-%d_%H:%M:%S")) diff --git a/src/proteus/star/__init__.py b/src/proteus/star/__init__.py new file mode 100644 index 00000000..4d21ee85 --- /dev/null +++ b/src/proteus/star/__init__.py @@ -0,0 +1,3 @@ +from __future__ import annotations + +__all__ = [] diff --git a/src/proteus/star/dummy.py b/src/proteus/star/dummy.py new file mode 100644 index 00000000..6ea1da6e --- /dev/null +++ b/src/proteus/star/dummy.py @@ -0,0 +1,107 @@ +# Dummy star module +from __future__ import annotations + +import logging + +import numpy as np + +from proteus.utils.constants import AU, const_sigma +from proteus.utils.phys import planck_wav + +log = logging.getLogger("fwl."+__name__) + +PLANCK_MIN_TEMPERATURE:float = 0.1 + +def generate_spectrum(tmp:float, R_star:float): + ''' + Get stellar spectrum at 1 AU, assuming that the star emits like a blackbody. + + Parameters + ----------- + tmp : float + Temperature [K] + R_star : float + Stellar radius [m] + + Returns + ----------- + wl_arr : list + Wavelength bin centres [nm] + fl_arr : list + Stellar spectral flux density at 1 AU from star [erg s-1 cm-2 nm-1] + ''' + + log.debug("Generating stellar spectrum at Teff=%.0f K"%tmp) + + # Allocate wavelength array + wl_min = 1e-10 # 0.1 nm + wl_max = 1e-1 # 10 cm + wl_pts = 600 # number of points to sample + wl_arr = np.logspace(np.log10(wl_min), np.log10(wl_max), wl_pts) + + # Allocate flux array with zeros + fl_arr = np.zeros(np.shape(wl_arr)) + + # If Teff=0, keep fluxes at zero. This is equivalent to turning the star 'off'. + + # Else, for non-zero effective temperatures... + if tmp > PLANCK_MIN_TEMPERATURE: + + # Evaluate planck function in each bin + for i,wav in enumerate(wl_arr): + fl_arr[i] = planck_wav(tmp, wav) # W m-2 m-1 at stellar surface + + # Scale from stellar surface to 1 AU + fl_arr *= (R_star/AU)**2 + + # Convert m-1 to nm-1 + fl_arr *= 1e-9 + + # Convert W m-2 to erg s-1 cm-2 + fl_arr *= 1e3 + + # Convert wavelengths to nm + wl_arr *= 1e9 + + # Return as lists + return list(wl_arr), list(fl_arr) + +def calc_star_luminosity(tmp:float, R_star:float): + ''' + Calculate star's bolometric luminosity. + + Assumes that the star emits like a blackbody. + ''' + + # Evaluate stefan-boltzmann law at Teff + lum = 0.0 + if tmp > PLANCK_MIN_TEMPERATURE: + lum = 4 * np.pi * R_star * R_star * const_sigma * (tmp**4) + + return lum + + +def calc_instellation(tmp:float, R_star:float, sep:float): + ''' + Calculate planet's instellation based on the star's luminosity. + + Parameters + ----------- + tmp : float + Star's effective temperature + R_star : float + Star's radius [m] + sep : float + Planet-star separation [m] + + Returns + ----------- + S_0 : float + Instellation [W m-2] + ''' + + # Get luminosity + L_bol = calc_star_luminosity(tmp, R_star) + + # Return flux at planet-star distance + return L_bol / (4 * np.pi * sep * sep) diff --git a/src/proteus/star/wrapper.py b/src/proteus/star/wrapper.py new file mode 100644 index 00000000..fc6302e3 --- /dev/null +++ b/src/proteus/star/wrapper.py @@ -0,0 +1,185 @@ +# Generic stellar evolution wrapper +from __future__ import annotations + +import logging +import os +import shutil +from typing import TYPE_CHECKING + +import numpy as np + +from proteus.utils.constants import AU, L_sun, R_sun, const_sigma + +log = logging.getLogger("fwl."+__name__) + +if TYPE_CHECKING: + from proteus import Proteus + from proteus.config import Config + +def init_star(handler:Proteus): + ''' + Star-related things to be done when the simulation begins. + ''' + + # Path to the modern spectrum + # i.e. that observed by a telescope, or derived from various observations. + # This is what we download from OSF. + star_modern_path = os.path.join(handler.directories["fwl"], + handler.config.star.mors.spec) + + # Copy modern spectrum to output folder, for posterity. + star_backup_path = os.path.join(handler.directories["output"], "-1.sflux") + shutil.copyfile(star_modern_path, star_backup_path) + + # Dummy star modules does not require preparation + + # Prepare MORS + if handler.config.star.module == 'mors': + import mors + + match handler.config.star.mors.tracks: + + case 'spada': + # load modern spectrum + # calculate band-integrated fluxes + handler.star_struct = mors.spec.Spectrum() + handler.star_struct.LoadTSV(star_modern_path) + handler.star_struct.CalcBandFluxes() + + # calculate other properties from modern spectrum + handler.star_props = mors.synthesis.GetProperties( + handler.config.star.mass, + handler.config.star.omega, + handler.config.star.age_now * 1000, # convert Gyr to Myr + ) + + case 'baraffe': + handler.star_modern_wl, handler.star_modern_fl = mors.ModernSpectrumLoad( + star_modern_path, star_backup_path + ) + + handler.baraffe_track = mors.BaraffeTrack(handler.config.star.mass) + +def get_new_spectrum(t_star:float, R_star:float, config:Config, + star_struct_modern=None, star_props_modern=None, + baraffe_track=None, modern_wl=None, modern_fl=None): + ''' + Get new stellar spectrum at 1 AU. + ''' + + log.debug("Get new stellar spectrum (star age = %g Gyr)"%(t_star/1e9)) + + # Dummy case + if config.star.module == 'dummy': + from proteus.star.dummy import generate_spectrum + wl, fl = generate_spectrum(config.star.Teff, R_star) + + # Mors cases + elif config.star.module == 'mors': + + import mors + + match config.star.mors.tracks: + case 'spada': + synthetic = mors.synthesis.CalcScaledSpectrumFromProps( + star_struct_modern, star_props_modern, t_star / 1e6) + fl = synthetic.fl + wl = synthetic.wl + case 'baraffe': + fl = baraffe_track.BaraffeSpectrumCalc( + t_star, config.star.lum_now, modern_fl) + wl = modern_wl + + return wl, fl + +def scale_spectrum_to_toa(fl_arr, sep:float): + ''' + Scale stellar fluxes from 1 AU to top of atmosphere + ''' + return np.array(fl_arr) * ( (AU / sep)**2 ) + +def write_spectrum(wl_arr, fl_arr, hf_row:dict, output_dir:str): + ''' + Write stellar spectrum to file. + ''' + + log.debug("Writing stellar spectrum to file") + + # Header information + header = ( + "# WL(nm)\t Flux(ergs/cm**2/s/nm) Stellar flux at t_star = %.2e yr" + % hf_row["age_star"] + ) + + # Write to TSV file + np.savetxt( + os.path.join(output_dir, "data", "%d.sflux" % hf_row["Time"]), + np.array([wl_arr, fl_arr]).T, + header=header, + comments="", + fmt="%.8e", + delimiter="\t", + ) + +def update_stellar_radius(hf_row:dict, config:Config, baraffe_track=None): + ''' + Update stellar radius in hf_row, stored in SI units. + ''' + + # Dummy case + if config.star.module == 'dummy': + R_star = config.star.radius + + # Mors cases + elif config.star.module == 'mors': + + import mors + + # which track? + match config.star.mors.tracks: + case 'spada': + R_star = mors.Value(config.star.mass, hf_row["age_star"] / 1e6, "Rstar") + case 'baraffe': + R_star = baraffe_track.BaraffeStellarRadius(hf_row["age_star"]) + + # Dimensionalise and store in dictionary + hf_row["R_star"] = R_star * R_sun + +def update_instellation(hf_row:dict, config:Config, baraffe_track=None): + ''' + Update hf_row value of bolometric stellar flux impinging upon the planet. + ''' + + # Dummy case + if config.star.module == 'dummy': + from proteus.star.dummy import calc_instellation + S_0 = calc_instellation(config.star.Teff, hf_row["R_star"], hf_row["separation"]) + + # Mors cases + elif config.star.module == 'mors': + + import mors + + # which track? + match config.star.mors.tracks: + case 'spada': + S_0 = mors.Value(config.star.mass, hf_row["age_star"] / 1e6, "Lbol") \ + * L_sun / (4.0 * np.pi * hf_row["separation"]**2.0 ) + + case 'baraffe': + S_0 = baraffe_track.BaraffeSolarConstant(hf_row["age_star"], + hf_row["separation"]/AU) + + # Update hf_row dictionary + hf_row["F_ins"] = S_0 + +def update_equilibrium_temperature(hf_row:dict, config:Config): + ''' + Calculate planetary equilibrium temperature. + ''' + + # Absorbed stellar flux + F_asf = hf_row["F_ins"] * config.orbit.s0_factor * (1-config.atmos_clim.albedo_pl) + + # Planetary equilibrium temperature + hf_row["T_eqm"] = (F_asf / const_sigma)**(1.0/4.0) diff --git a/src/proteus/utils/coupler.py b/src/proteus/utils/coupler.py index cc821c02..8b6312fd 100644 --- a/src/proteus/utils/coupler.py +++ b/src/proteus/utils/coupler.py @@ -34,7 +34,6 @@ from proteus.plot.cpl_sflux_cross import plot_sflux_cross from proteus.plot.cpl_stacked import plot_stacked from proteus.utils.constants import ( - const_sigma, element_list, volatile_species, ) @@ -63,12 +62,6 @@ def GitRevision(dir:str) -> str: return hash -def CalculateEqmTemperature(I_0, ASF_sf, A_B): - ''' - Calculate planetary equilibrium temperature. - Params: Stellar flux, ASF scale factor, and bond albedo. - ''' - return (I_0 * ASF_sf * (1.0 - A_B) / const_sigma)**(1.0/4.0) def PrintCurrentState(hf_row:dict): ''' diff --git a/src/proteus/utils/data.py b/src/proteus/utils/data.py index 56c296d8..8c683f10 100644 --- a/src/proteus/utils/data.py +++ b/src/proteus/utils/data.py @@ -155,7 +155,7 @@ def download_sufficient_data(config:Config): """ # Star stuff - if config.star.mors.tracks in ('spada', 'baraffe'): + if config.star.module == "mors": download_stellar_spectra() if config.star.mors.tracks == 'spada': download_evolution_tracks("Spada") diff --git a/src/proteus/utils/helper.py b/src/proteus/utils/helper.py index 9bc124e9..28e597ed 100644 --- a/src/proteus/utils/helper.py +++ b/src/proteus/utils/helper.py @@ -194,21 +194,6 @@ def find_nearest(array, target): close = array[idx] return close, idx -def mol_to_ele(mol:str): - ''' - Return the number of each element within a given molecule, as a dictionary - ''' - decomp = re.findall(r'([A-Z][a-z]?)(\d*)', mol) # https://codereview.stackexchange.com/a/232664 - elems = {} - for ev in decomp: - if ev[1] == '': - val = 1 - else: - val = int(ev[1]) - elems[str(ev[0])] = val - return elems - - def recursive_get(d, keys): ''' Function to access nested dictionaries diff --git a/src/proteus/utils/phys.py b/src/proteus/utils/phys.py new file mode 100644 index 00000000..5519074c --- /dev/null +++ b/src/proteus/utils/phys.py @@ -0,0 +1,66 @@ +# Small physics functions that can be used universally +# This file should not depend on too many other files, as this can cause circular import issues + +from __future__ import annotations + +import logging +import re + +import numpy as np + +from proteus.utils.constants import const_c, const_h, const_k + +log = logging.getLogger("fwl."+__name__) + +def mol_to_ele(mol:str): + ''' + Return the number of each element within a given molecule, as a dictionary + ''' + decomp = re.findall(r'([A-Z][a-z]?)(\d*)', mol) # https://codereview.stackexchange.com/a/232664 + elems = {} + for ev in decomp: + if ev[1] == '': + val = 1 + else: + val = int(ev[1]) + elems[str(ev[0])] = val + return elems + +def planck_wav(tmp:float, wav:float): + ''' + Evaluate the planck function at a given temperature and wavelength. + All in SI units. + + Parameters + ----------- + tmp : float + Temperature [K] + wav : float + Wavelength [m] + + Returns + ----------- + flx : float + Spectral flux density [W m-2 m-1] + ''' + + # Optimisation variables + wav5 = wav*wav*wav*wav*wav + hc = const_h * const_c + + # Calculate planck function value [W m-2 sr-1 m-1] + # http://spiff.rit.edu/classes/phys317/lectures/planck.html + with np.errstate(all='raise'): + # Catch overflow error, which can occur when evaluating at short wavelengths + try: + flx = 2.0 * hc * (const_c / wav5) / ( np.exp(hc / (wav * const_k * tmp)) - 1.0) + except FloatingPointError: + flx = 0.0 + + # Integrate solid angle (hemisphere) + flx = flx * np.pi + + # Do not allow zero or near-zero fluxes + flx = max(flx, 1e-40) + + return flx diff --git a/start_proteus.py b/start_proteus.py deleted file mode 100755 index 0130d204..00000000 --- a/start_proteus.py +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/env python3 - -""" -PROTEUS run script -""" -from __future__ import annotations - -from proteus import Proteus - -runner = Proteus(config_path='input/default.toml') -runner.start( - resume=False, -)