Validate escape config

.github/workflows/tests.yaml

@@ -122,25 +122,15 @@ jobs:
       - name: Test with pytest
         run: coverage run -m pytest
 
-      # Upload output run images if tests fail
-      - name: Upload PNG plots as artifacts
+      # Upload result if tests fail
+      - name: Upload result as artifact
         if: failure()
         uses: actions/upload-artifact@v3
         with:
-          name: png-plots
+          name: integration-physical
           path: |
-            output/physical/plot_atmosphere.png
-            output/physical/plot_escape.png
-            output/physical/plot_global_log.png
-            output/physical/plot_observables.png
-            output/physical/plot_stacked.png
-            output/physical/plot_elements.png
-            output/physical/plot_fluxes_atmosphere.png
-            output/physical/plot_interior_cmesh.png
-            output/physical/plot_sflux_cross.png
-            output/physical/plot_emission.png
-            output/physical/plot_interior.png
-            output/physical/plot_sflux.png
+            output/physical/plot_*.png
+            output/physical/runtime_helpfile.csv
 
       - name: Report coverage
         run: |

pyproject.toml

@@ -49,7 +49,7 @@ dependencies = [
     "fwl-mors>=24.11.18",
     "fwl-calliope>=24.9.10",
     "fwl-zephyrus>=24.10.15",
-    "aragog>=0.1.6a0",
+    "aragog==0.1.6a0",
     "cmcrameri",
     "juliacall",
     "matplotlib",

src/proteus/atmos_clim/agni.py

@@ -10,9 +10,9 @@
 from juliacall import Main as jl
 from scipy.interpolate import PchipInterpolator
 
-from proteus.atmos_clim.common import get_spfile_path
+from proteus.atmos_clim.common import get_height_from_pressure, get_spfile_path
 from proteus.utils.constants import gas_list
-from proteus.utils.helper import UpdateStatusfile, create_tmp_folder, find_nearest, safe_rm
+from proteus.utils.helper import UpdateStatusfile, create_tmp_folder, safe_rm
 from proteus.utils.logs import GetCurrentLogfileIndex, GetLogfilePath
 
 if TYPE_CHECKING:
@@ -469,8 +469,13 @@ def run_agni(atmos, loops_total:int, dirs:dict, config:Config, hf_row:dict):
                                         (net_flux[-1], net_flux[0] ,LW_flux_up[0]))
 
     # XUV height in atm
-    p_xuv, idx_xuv =  find_nearest(atmos.p, config.escape.zephyrus.Pxuv*1e5)
-    z_xuv = atmos.z[idx_xuv]
+    if config.escape.module == 'zephyrus':
+        # escape level set by zephyrus config
+        p_xuv = config.escape.zephyrus.Pxuv # [bar]
+    else:
+        # escape level set to surface
+        p_xuv = hf_row["P_surf"] # [bar]
+    p_xuv, z_xuv = get_height_from_pressure(atmos.p, atmos.z, p_xuv*1e5) # [Pa], [m]
 
     # final things to store
     output = {}

src/proteus/atmos_clim/common.py

@@ -8,6 +8,8 @@
 import netCDF4 as nc
 import numpy as np
 
+from proteus.utils.helper import find_nearest
+
 if TYPE_CHECKING:
     from proteus.config import Config
 
@@ -114,3 +116,29 @@ def get_spfile_path(fwl_dir:str, config:Config):
 
     # Construct file path
     return os.path.join(fwl_dir,"spectral_files",group,bands,group)+".sf"
+
+def get_height_from_pressure(p_arr, z_arr, p_tgt):
+    """
+    Get the geometric height [m] corresponding to a given pressure.
+
+    Parameters:
+    ----------------
+        p_arr: list
+            Pressure array
+        z_arr: list
+            Height array
+        p_tgt: float
+            Target pressure
+
+    Returns:
+    ----------------
+        p_close: float
+            Closest pressure in the array
+        z_close: float
+            Closest height in the array
+    """
+
+    p_close, idx = find_nearest(p_arr, p_tgt)
+    z_close = z_arr[idx]
+
+    return float(p_close), z_close

src/proteus/atmos_clim/dummy.py

@@ -102,6 +102,12 @@ def _resid(x):
     log.info("    F_olr  =  %.3e  W m-2" % fluxes["fl_U_LW"])
     log.info("    F_sct  =  %.3e  W m-2" % fluxes["fl_U_SW"])
 
+    # Escape level always at surface for dummy atmosphere
+    if config.escape.module == 'zephyrus':
+        log.warning("Setting escape level to surface because dummy atmosphere is used")
+    p_xuv = P_surf # already in bar
+    z_xuv = 0.0
+
     output = {}
     output["T_surf"]  = T_surf_atm
     output["F_atm"]   = F_atm_lim             # Net flux at TOA
@@ -110,7 +116,7 @@ def _resid(x):
     output["z_obs"]   = 0.0
     output["rho_obs"] = 3 * M_int / (4*np.pi*R_int**3)
     output["albedo"]  = fluxes["fl_U_SW"]/fluxes["fl_D_SW"]
-    output["p_xuv"]   = P_surf
-    output["z_xuv"]   = 0.0
+    output["p_xuv"]   = p_xuv
+    output["z_xuv"]   = z_xuv
 
     return output

src/proteus/atmos_clim/janus.py

@@ -11,8 +11,9 @@
 import numpy as np
 import pandas as pd
 
+from proteus.atmos_clim.common import get_height_from_pressure
 from proteus.utils.constants import vap_list, vol_list
-from proteus.utils.helper import UpdateStatusfile, create_tmp_folder, find_nearest
+from proteus.utils.helper import UpdateStatusfile, create_tmp_folder
 
 if TYPE_CHECKING:
     from proteus.config import Config
@@ -269,18 +270,25 @@ def RunJANUS(atm, dirs:dict, config:Config, hf_row:dict, hf_all:pd.DataFrame,
         log.warning("    %g  ->  %g" % (F_atm_new , F_atm_lim))
 
     # observables
-    z_obs = -1.0
+    z_obs = 0.0
     rho_obs = -1.0
-    if not atm.height_error:
-        # find 1 mbar level
-        idx = find_nearest(atm.p, 1e2)[1]
-        z_obs = atm.z[idx]
-        # calc observed density
+    if atm.height_error:
+        log.error("Hydrostatic integration failed in JANUS!")
+    else:
+        # find observed level [m] at 1 mbar
+        _, z_obs = get_height_from_pressure(atm.p, atm.z, 1e2)
+
+        # calc observed density [kg m-3]
         rho_obs = calc_observed_rho(atm)
 
     # XUV height in atm
-    p_xuv, idx_xuv =  find_nearest(atm.p, config.escape.zephyrus.Pxuv*1e5)
-    z_xuv = atm.z[idx_xuv]
+    if config.escape.module == 'zephyrus':
+        # escape level set by zephyrus config
+        p_xuv = config.escape.zephyrus.Pxuv # [bar]
+    else:
+        # escape level set to surface
+        p_xuv = hf_row["P_surf"] # [bar]
+    p_xuv, z_xuv = get_height_from_pressure(atm.p, atm.z, p_xuv*1e5) # [Pa], [m]
 
     # final things to store
     output={}

src/proteus/atmos_clim/wrapper.py

@@ -152,12 +152,12 @@ def RunAtmosphere(config:Config, dirs:dict, loop_counter:dict,
     hf_row["T_surf"] = atm_output["T_surf"]
     hf_row["F_net"]  = hf_row["F_int"] - hf_row["F_atm"]
     hf_row["bond_albedo"]= atm_output["albedo"]
-    hf_row["p_xuv"]  = atm_output["p_xuv"]                  # Closest pressure from Pxuv    [Pa]
+    hf_row["p_xuv"]  = atm_output["p_xuv"]                  # Closest pressure from Pxuv    [bar]
     hf_row["z_xuv"]  = atm_output["z_xuv"]                  # Height at Pxuv                [m]
-    hf_row["R_xuv"]  = hf_row["R_int"] + hf_row["z_xuv"]    # Computed XUV planetary radius
+    hf_row["R_xuv"]  = hf_row["R_int"] + hf_row["z_xuv"]    # Radius at z_xuv [m]
 
     # Calculate observables (measured at infinite distance)
-    R_obs = hf_row["z_obs"] + hf_row["R_int"] # observed radius
+    R_obs = hf_row["z_obs"] + hf_row["R_int"] # observed radius [m]
     hf_row["transit_depth"] =  (R_obs / hf_row["R_star"])**2.0
     hf_row["contrast_ratio"] = ((hf_row["F_olr"]+hf_row["F_sct"])/hf_row["F_ins"]) * \
                                  (R_obs / hf_row["separation"])**2.0

src/proteus/config/_config.py

@@ -18,6 +18,13 @@
 
 log = logging.getLogger('fwl.' + __name__)
 
+def spada_zephyrus(instance, attribute, value):
+    # using zephyrus
+    #     zephyrus requires MORS + Spada
+    if (instance.escape.module == 'zephyrus') and \
+        not ( (instance.star.module == 'mors') and (instance.star.mors.tracks == 'spada')):
+        raise ValueError('ZEPHYRUS must be used with MORS and the Spada evolution tracks')
+
 @define
 class Config:
     """Root config parameters.
@@ -56,7 +63,7 @@ class Config:
     orbit: Orbit
     struct: Struct
     atmos_clim: AtmosClim
-    escape: Escape
+    escape: Escape = field(validator=(spada_zephyrus,))
     interior: Interior
     outgas: Outgas
     delivery: Delivery

src/proteus/orbit/wrapper.py

@@ -47,8 +47,15 @@ def update_period(hf_row:dict, sma:float):
             Semimajor axis [AU]
     '''
 
+    # Total mass of system, kg
+    M_total = hf_row["M_star"] + hf_row["M_tot"]
+
+    # Sanity check
+    if M_total < 1e3:
+       log.error("Unreasonable star+planet mass: %.5e kg"%M_total)
+
     # Standard gravitational parameter (planet mass + star mass)
-    mu = const_G * (hf_row["M_star"] + hf_row["M_tot"])
+    mu = const_G * M_total
 
     # Semimajor axis in SI units
     a = sma * AU

src/proteus/proteus.py

@@ -17,6 +17,7 @@
     get_new_spectrum,
     init_star,
     scale_spectrum_to_toa,
+    update_stellar_mass,
     update_stellar_quantities,
     write_spectrum,
 )
@@ -171,6 +172,9 @@ def start(self, *, resume: bool = False):
             # Create an empty initial row for helpfile
             self.hf_row = ZeroHelpfileRow()
 
+            # Stellar mass
+            update_stellar_mass(self.hf_row, self.config)
+
             # Initial time
             self.hf_row["Time"] = 0.0
             self.hf_row["age_star"] = self.config.star.age_ini * 1e9