Conversation
…r species can be set by the user
…c species in the atmosphere and added a new plotting routine which allows comparison of TPS at different time steps between two models
…erent runs are compared at similar timesteps
… on from lavatmos
…0K temperature cap for lavatmos
… allows coupling with LavAtmos,
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #634 +/- ##
==========================================
- Coverage 63.51% 63.49% -0.03%
==========================================
Files 98 98
Lines 9368 9368
Branches 1261 1261
==========================================
- Hits 5950 5948 -2
- Misses 3418 3420 +2
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull request overview
This PR introduces an optional coupling between PROTEUS outgassing (CALLIOPE) and LavAtmos, expands supported gas/species handling, and adds plotting utilities for comparing atmospheric outputs across runs/time steps.
Changes:
- Add LavAtmos ↔ CALLIOPE iterative coupling and make gas lists configurable via
outgas.silicates+outgas.vaplist. - Extend helpfile schema and plotting utilities to include additional species/elements and comparison plots.
- Update environment/config examples for albedo and outgassing test runs.
Reviewed changes
Copilot reviewed 21 out of 21 changed files in this pull request and generated 18 comments.
Show a summary per file
| File | Description |
|---|---|
| tools/get_petsc.sh | Extends OSTYPE handling for Linux PETSc arch selection. |
| src/proteus/utils/plot.py | Adds/adjusts species colour mappings (incl. refractory species). |
| src/proteus/utils/coupler.py | Makes helpfile keys depend on config + adds new output keys. |
| src/proteus/utils/constants.py | Updates supported species/elements lists and related constants. |
| src/proteus/proteus.py | Switches initialization/outgassing flow to use configurable gas lists and LavAtmos loop; updates helpfile read/write call sites. |
| src/proteus/post_processing/compar_atmosphere_models.py | New script to compare atmospheric profiles between two output dirs. |
| src/proteus/plot/cpl_global.py | Uses config-dependent gas list in global plotting. |
| src/proteus/plot/cpl_escape.py | Adds additional logging during escape plotting. |
| src/proteus/plot/cpl_chem_atmosphere.py | Expands default plotted species set to include refractory gases. |
| src/proteus/plot/cpl_atmosphere.py | Tweaks plotting bounds for temperature–pressure/height panels. |
| src/proteus/plot/compar_atmosphere_models.py | New plotting script for comparing runs and selected times (currently includes hardcoded local paths). |
| src/proteus/outgas/wrapper.py | Adds get_gaslist + LavAtmos/CALLIOPE convergence loop. |
| src/proteus/outgas/lavatmos.py | New LavAtmos integration module (currently contains hardcoded paths and several unit/key bugs). |
| src/proteus/outgas/common.py | Makes expected_keys depend on config gas list. |
| src/proteus/outgas/calliope.py | Allows CALLIOPE to use hf_row['fO2_shift'] if present. |
| src/proteus/config/_outgas.py | Adds LavAtmos-related config fields (currently defaults to coupling enabled + absolute paths). |
| src/proteus/cli.py | Minor formatting cleanup. |
| src/proteus/atmos_clim/agni.py | Uses config-dependent gas list and sets AGNI albedo field; adds debug prints. |
| input/run_tests.toml | New comprehensive run config (includes absolute local paths for LavAtmos/FastChem). |
| input/albedo_test.toml | New albedo lookup test config. |
| environment.yml | Pins Python to 3.12. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
src/proteus/proteus.py
Outdated
|
|
||
| # Read helpfile from disk | ||
| self.hf_all = ReadHelpfileFromCSV(self.directories['output']) | ||
| self.hf_all = ReadHelpfileFromCSV(self.directories["output"],self.config) |
There was a problem hiding this comment.
ReadHelpfileFromCSV in proteus.utils.coupler currently only accepts (output_dir: str), but this call passes self.config as a second argument. This will raise a TypeError when resuming, observing, or running offline chemistry. Either update ReadHelpfileFromCSV to accept config (and use it for key validation), or remove the extra argument at call sites.
| self.hf_all = ReadHelpfileFromCSV(self.directories["output"],self.config) | |
| self.hf_all = ReadHelpfileFromCSV(self.directories["output"]) |
src/proteus/outgas/lavatmos.py
Outdated
| hf_row[vol + "_kg_tot"] = hf_row[vol + "_kg_atm"] + hf_row[vol + "_kg_solid"] + hf_row[vol + "_kg_liquid"] | ||
|
|
||
| #elements are not considered as atomic species but just as inventory | ||
| for e in element_list: | ||
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0] * M_atmo_new * species_lib[e].weight/mu_outgassed | ||
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0]+ hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] | ||
| hf_row[e + "_kg_tot"] = hf_row[e + "_kg_atm"] + hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] |
There was a problem hiding this comment.
The LavAtmos output updates write totals to keys named *_kg_tot, but the rest of the codebase/helpfile schema uses *_kg_total (see expected_keys and GetHelpfileKeys). This will leave the canonical *_kg_total fields stale and can trigger KeyErrors/incorrect accounting downstream. Update these assignments to use the existing *_kg_total keys consistently.
| hf_row[vol + "_kg_tot"] = hf_row[vol + "_kg_atm"] + hf_row[vol + "_kg_solid"] + hf_row[vol + "_kg_liquid"] | |
| #elements are not considered as atomic species but just as inventory | |
| for e in element_list: | |
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0] * M_atmo_new * species_lib[e].weight/mu_outgassed | |
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0]+ hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] | |
| hf_row[e + "_kg_tot"] = hf_row[e + "_kg_atm"] + hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] | |
| hf_row[vol + "_kg_total"] = hf_row[vol + "_kg_atm"] + hf_row[vol + "_kg_solid"] + hf_row[vol + "_kg_liquid"] | |
| #elements are not considered as atomic species but just as inventory | |
| for e in element_list: | |
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0] * M_atmo_new * species_lib[e].weight/mu_outgassed | |
| hf_row[e + "_kg_atm"]=new_atmos_abundances[e][0]+ hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] | |
| hf_row[e + "_kg_total"] = hf_row[e + "_kg_atm"] + hf_row[e + "_kg_solid"] + hf_row[e + "_kg_liquid"] |
src/proteus/config/_outgas.py
Outdated
| silicates: bool = field(default=True) | ||
| fastchempath: str = field(default="/data3/leoni/LavAtmos/FastChem/fastchem3/output/") | ||
| vaplist: list = field(default=['SiO','SiO2','MgO','FeO']) |
There was a problem hiding this comment.
These defaults make the LavAtmos coupling effectively “on by default” and embed an absolute, user-specific fastchempath. This is likely to break out-of-the-box runs and CI, and it contradicts the PR description that LavAtmos must be installed separately. Consider defaulting silicates to False, and make fastchempath optional/empty with validation only when coupling is enabled.
| silicates: bool = field(default=True) | |
| fastchempath: str = field(default="/data3/leoni/LavAtmos/FastChem/fastchem3/output/") | |
| vaplist: list = field(default=['SiO','SiO2','MgO','FeO']) | |
| # LavAtmos / silicate coupling is opt-in: default to disabled. | |
| silicates: bool = field(default=False) | |
| # Path to FastChem output directory used by LavAtmos; must be set when silicate coupling is enabled. | |
| fastchempath: str = field(default='') | |
| # List of vapour species for LavAtmos; must be set when silicate coupling is enabled. | |
| vaplist: list = field(factory=list) | |
| def __attrs_post_init__(self) -> None: | |
| """ | |
| Perform cross-field validation after initialization. | |
| When silicate (LavAtmos) coupling is enabled via `silicates=True`, | |
| ensure that `fastchempath` and `vaplist` are explicitly configured. | |
| """ | |
| if self.silicates: | |
| if not self.fastchempath: | |
| raise ValueError( | |
| "Outgas.fastchempath must be set when silicate (LavAtmos) coupling is enabled." | |
| ) | |
| if not self.vaplist: | |
| raise ValueError( | |
| "Outgas.vaplist must be set when silicate (LavAtmos) coupling is enabled." | |
| ) |
| # Supported gases | ||
| vol_list = ['H2O', 'CO2', 'O2', 'H2', 'CH4', 'CO', 'N2', 'NH3', 'S2', 'SO2', 'H2S'] | ||
| vap_list = ['SiO', 'SiO2', 'MgO', 'FeO2'] | ||
| gas_list = vol_list + vap_list | ||
| vap_list = ["SiO", "SiO2", "MgO", "FeO2"] | ||
| vol_list = ["H2O", "CO2", "O2", "H2", "CH4", "CO", "N2", "NH3", "S2", "SO2", "H2S"] | ||
| #vap_list = ["SiO", "SiO2", "MgO", "FeO2"] | ||
| #vap_list = ["SiO", "SiO2", "MgO", "FeO", "Na"] | ||
| #gas_list = vol_list + vap_list | ||
|
|
There was a problem hiding this comment.
gas_list is no longer defined (it’s commented out), but multiple modules/tests still import proteus.utils.constants.gas_list (e.g., atmos_clim/janus.py, inference/*, tests/outgas). This will raise ImportError at runtime and break the test suite. Either restore gas_list = vol_list + vap_list for backward compatibility, or update all imports/usages to use the new get_gaslist(config) pattern consistently (including tests).
There was a problem hiding this comment.
This is a good point, @leojola. Either of the changes suggested by copilot would be fine for me. The second option would be more future-proof.
src/proteus/utils/coupler.py
Outdated
| "M_star", "R_star", "age_star", # [kg], [m], [yr] | ||
| "T_star", # [K] | ||
|
|
||
|
|
||
| # Observational (from infinity) | ||
| "p_obs", # observered radius [bar] | ||
| "R_obs", # observed radius [m] | ||
| "rho_obs", # observed bulk density [kg m-3] | ||
| "transit_depth", "eclipse_depth", # [1], [1] | ||
|
|
||
| "albedo_pl", # input bond albedo [1] | ||
| "bond_albedo", # output bond albedo [1] | ||
|
|
||
| # Imaginary part of k2 Love Number | ||
| "Imk2", # [1] | ||
|
|
||
| # Escape | ||
| "esc_rate_total", "p_xuv", "R_xuv", # [kg s-1], [bar], [m] | ||
|
|
||
| # Atmospheric composition from outgassing | ||
| "M_atm", "P_surf", "atm_kg_per_mol", # [kg], [bar], [kg mol-1] |
There was a problem hiding this comment.
The keys list contains duplicate entries for many fields (e.g., M_star, R_star, age_star, p_obs, M_atm, P_surf, etc.) that are already included earlier in the same list. Duplicate column names in the helpfile DataFrame can lead to ambiguous selection/overwrites and hard-to-debug downstream behavior. Remove the duplicated block and keep a single canonical definition for each key.
| "M_star", "R_star", "age_star", # [kg], [m], [yr] | |
| "T_star", # [K] | |
| # Observational (from infinity) | |
| "p_obs", # observered radius [bar] | |
| "R_obs", # observed radius [m] | |
| "rho_obs", # observed bulk density [kg m-3] | |
| "transit_depth", "eclipse_depth", # [1], [1] | |
| "albedo_pl", # input bond albedo [1] | |
| "bond_albedo", # output bond albedo [1] | |
| # Imaginary part of k2 Love Number | |
| "Imk2", # [1] | |
| # Escape | |
| "esc_rate_total", "p_xuv", "R_xuv", # [kg s-1], [bar], [m] | |
| # Atmospheric composition from outgassing | |
| "M_atm", "P_surf", "atm_kg_per_mol", # [kg], [bar], [kg mol-1] |
| key = gas + '_vmr' | ||
| if key in atm_profile.keys(): | ||
| xarr = list(atm_profile[key]) | ||
| xarr = [xarr[0]] + xarr | ||
| if np.amax(xarr) >= xmin: | ||
| if (np.amax(xarr) >= xmin or key in list(REFRACTORY_GASES) ): | ||
| vmr = float(xarr[-1]) | ||
| ax.plot(xarr, parr, ls='dashed', color=col, lw=_lw, alpha=al) | ||
| ax.plot(xarr, parr, ls = 'dashed', color=col, lw=_lw, alpha=al) | ||
| else: | ||
| print(np.amax(xarr)) | ||
|
|
There was a problem hiding this comment.
The refractory-species inclusion condition is incorrect: key is like 'SiO_vmr', but REFRACTORY_GASES contains bare species names (e.g., 'SiO'). As written, the or key in REFRACTORY_GASES branch will never trigger, and refractory gases below xmin won’t be plotted as intended. Compare gas (or strip the _vmr suffix) instead, and remove the leftover print() debug output.
There was a problem hiding this comment.
This is also correct. Please address this suggestion too.
| output_dir='/data3/leoni/PROTEUS/output/calliope_lavatmos_coupling_xerr001' | ||
| plottimes=[124027,134027,142186,152186] | ||
|
|
||
| plot_atmosphere_selected_times(output_dir, plottimes) | ||
| plot_chemistry_selected_times(output_dir, plottimes) |
There was a problem hiding this comment.
The __main__ section hard-codes a local absolute output directory (/data3/leoni/...) and specific plot times. This makes the script unusable for other users and will fail if someone runs it accidentally. Replace this with sys.argv parsing (as commented out above) or remove the hard-coded block / move it into documentation/examples.
src/proteus/outgas/wrapper.py
Outdated
| log.info("silicates are outgassed") | ||
| log.info("error threshold on fO2 shift : %.6f"%xerr) | ||
| log.info("current error : %.6f"%err) | ||
| while err > xerr: | ||
| old_fO2shift = hf_row['fO2_shift'] | ||
| run_lavatmos(config, hf_row) #in run_lavatmos add a criterion for temperature and melt fraction ? | ||
| log.info("new fO2 shift : %.6f"%hf_row["fO2_shift"]) | ||
| run_outgassing(dirs, config, hf_row) | ||
| err=abs(old_fO2shift - hf_row['fO2_shift']) | ||
| log.info("fO2 shift after runnning lavatmos: %.6f"%hf_row['fO2_shift']) | ||
| log.info("change in fO2 between the last iterations: %.6f"%err) |
There was a problem hiding this comment.
lavatmos_calliope_loop uses an unbounded while err > xerr loop with no maximum iteration count or divergence detection. If LavAtmos/Calliope fail to converge (or oscillate), the simulation can hang indefinitely. Add a max-iterations safeguard (configurable), log the iteration count, and raise a clear error / fall back when convergence isn’t achieved.
| log.info("silicates are outgassed") | |
| log.info("error threshold on fO2 shift : %.6f"%xerr) | |
| log.info("current error : %.6f"%err) | |
| while err > xerr: | |
| old_fO2shift = hf_row['fO2_shift'] | |
| run_lavatmos(config, hf_row) #in run_lavatmos add a criterion for temperature and melt fraction ? | |
| log.info("new fO2 shift : %.6f"%hf_row["fO2_shift"]) | |
| run_outgassing(dirs, config, hf_row) | |
| err=abs(old_fO2shift - hf_row['fO2_shift']) | |
| log.info("fO2 shift after runnning lavatmos: %.6f"%hf_row['fO2_shift']) | |
| log.info("change in fO2 between the last iterations: %.6f"%err) | |
| # Maximum number of LavAtmos/Calliope iterations; allow configuration with a sensible default | |
| max_iterations = getattr(getattr(config, "outgas", config), "max_lavatmos_iterations", 100) | |
| iteration = 0 | |
| log.info("silicates are outgassed") | |
| log.info("error threshold on fO2 shift : %.6f"%xerr) | |
| log.info("initial error : %.6f"%err) | |
| log.info("maximum LavAtmos/Calliope iterations : %d"%max_iterations) | |
| while err > xerr and iteration < max_iterations: | |
| iteration += 1 | |
| log.info("LavAtmos/Calliope iteration %d"%iteration) | |
| old_fO2shift = hf_row['fO2_shift'] | |
| run_lavatmos(config, hf_row) #in run_lavatmos add a criterion for temperature and melt fraction ? | |
| log.info("new fO2 shift : %.6f"%hf_row["fO2_shift"]) | |
| run_outgassing(dirs, config, hf_row) | |
| err=abs(old_fO2shift - hf_row['fO2_shift']) | |
| log.info("fO2 shift after running lavatmos: %.6f"%hf_row['fO2_shift']) | |
| log.info("change in fO2 between the last iterations: %.6f"%err) | |
| if err > xerr: | |
| log.error( | |
| "LavAtmos/Calliope did not converge within %d iterations " | |
| "(final |ΔfO2_shift| = %.6f, threshold = %.6f)", | |
| max_iterations, | |
| err, | |
| xerr, | |
| ) | |
| raise RuntimeError( | |
| "LavAtmos/Calliope convergence failed: " | |
| "maximum number of iterations (%d) exceeded without reaching " | |
| "the fO2_shift error threshold (final error = %.6f, threshold = %.6f)" | |
| % (max_iterations, err, xerr) | |
| ) |
| #This routine needs to be called via python3 src/proteus/plot/compar_atmosphere_models.py. | ||
| #then it needs the folders for outputdir1 and outputdir 2 in the terminal as | ||
| #sys.argv[1] and sys.argv[2]. |
There was a problem hiding this comment.
The module header comment references src/proteus/plot/compar_atmosphere_models.py, but this file lives under src/proteus/post_processing/. This is misleading for users and documentation. Update the invocation instructions to match the actual path and consider adding basic sys.argv length checking with a helpful usage message.
input/run_tests.toml
Outdated
| fO2_shift_IW = 4 # atmosphere/interior boundary oxidation state [log10(ΔIW)] | ||
|
|
||
| silicates = true | ||
| fastchempath = '/data3/leoni/LavAtmos/FastChem/fastchem3/output/' |
There was a problem hiding this comment.
This example config hard-codes fastchempath to an absolute, user-specific location (/data3/leoni/...). If this file is intended for tests or as a shareable example, it won’t work on other machines/CI. Use a relative path under the repo (e.g., a fixture in tests/data) or leave it unset and document how to provide it when silicates=true.
| fastchempath = '/data3/leoni/LavAtmos/FastChem/fastchem3/output/' | |
| # Path to FastChem output directory, relative to the repository root. | |
| fastchempath = "tests/data/fastchem/output" |
… only for element budget. updated the file path for the elemnt abundance file output from lavatmos, such that the correct volatile abundances arae read in each iteration
There was a problem hiding this comment.
Thanks for all your efforts in getting this working, @leojola. Apologies for my week-long delay in addressing your request for a review.
I have a couple of suggestions and requirements before this can be merged.
My main confusion is that the fO2 loop seems semi-implemented still? Additionally, there's inconsistency in how the get_gaslist is used compared to the gas_list global variable. Also, the way LavAtmos and FastChem are accessed is based on paths specific to your user directory - these should be generalised in the same manner as FastChem and SOCRATES are.
There was a problem hiding this comment.
This config file seems quite specific for testing with LavAtmos. Can you rename it to something less generic than run_tests.toml please?
| silicates: bool = field(default=False) | ||
| converge_fO2: bool = field(default=False) | ||
| # Path to FastChem output directory used by LavAtmos; must be set when silicate coupling is enabled. | ||
| fastchempath: str = field(default='') |
There was a problem hiding this comment.
AGNI already requires that a path to Fastchem is set, which is done via the FC_DIR environment variable. This config option is duplicating that functionality, so please remove this fastchempath option and use the FC_DIR variable instead.
| # Path to FastChem output directory used by LavAtmos; must be set when silicate coupling is enabled. | ||
| fastchempath: str = field(default='') | ||
| # Path to FastChem input directory used by LavAtmos; must be set when silicate coupling is enabled. | ||
| elementfile: str = field(default='') |
There was a problem hiding this comment.
Why provide this as a config option? Could it not just be created into the PROTEUS output folder?
|
|
||
| # LavAtmos / silicate coupling is opt-in: default to disabled. | ||
| silicates: bool = field(default=False) | ||
| converge_fO2: bool = field(default=False) |
There was a problem hiding this comment.
Does this converge_fO2 option need to be removed now? Since the coupling is no longer a loop.
| # Supported gases | ||
| vol_list = ['H2O', 'CO2', 'O2', 'H2', 'CH4', 'CO', 'N2', 'NH3', 'S2', 'SO2', 'H2S'] | ||
| vap_list = ['SiO', 'SiO2', 'MgO', 'FeO2'] | ||
| gas_list = vol_list + vap_list | ||
| vap_list = ["SiO", "SiO2", "MgO", "FeO2"] | ||
| vol_list = ["H2O", "CO2", "O2", "H2", "CH4", "CO", "N2", "NH3", "S2", "SO2", "H2S"] | ||
| #vap_list = ["SiO", "SiO2", "MgO", "FeO2"] | ||
| #vap_list = ["SiO", "SiO2", "MgO", "FeO", "Na"] | ||
| #gas_list = vol_list + vap_list | ||
|
|
There was a problem hiding this comment.
This is a good point, @leojola. Either of the changes suggested by copilot would be fine for me. The second option would be more future-proof.
| vmr = float(xarr[-1]) | ||
| ax.plot(xarr, parr, ls='dashed', color=col, lw=_lw, alpha=al) | ||
| else: | ||
| print(np.amax(xarr)) |
There was a problem hiding this comment.
Remove the print() added for debugging here.
| keys.append('runtime') # Simulation wall-clock runtime [s] | ||
|
|
||
| # fmt: on | ||
| keys.append("fO2_shift") #relative to IW buffer |
There was a problem hiding this comment.
See comment above. If this is the fO2 from LavAtmos's vapourisation calculation, please give it a more specific name, since it differs from that in the volatile solubility.
| inc_gases.append(s) | ||
| self.hf_row[s + '_bar'] = 0.0 | ||
| for s in gas_list: | ||
| if s not in vol_list: |
There was a problem hiding this comment.
Why? Doesn't this break the non-vapour cases?
| else: | ||
| run_outgassing(self.directories, self.config, self.hf_row) | ||
| # run_outgassing(self.directories, self.config, self.hf_row) | ||
| if self.config.outgas.converge_fO2 : |
There was a problem hiding this comment.
I thought the loop was removed?
| dependencies: | ||
| # Core Python | ||
| - python=3.13 | ||
| - python=3.12 |
2 participants
Description
A short summary of the changes made in this pull-request.
If silicates=true is set tin the input file LavAtmos and Calliope are run in a loop until they converge I=on the same oxygen fugacity
Checklist
Relevant people
Tag people who should know about this PR here