From abaafc84adba4abbb31fb685ed70e9e39047cc73 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 11 Sep 2024 17:24:17 -0400 Subject: [PATCH 1/2] Add gcpy/examples/geos-chem/check_photolysis.py example script gcpy/examples/geos-chem/check_photolysis.py - Example script that checks the consistency of the photolysis reactions in GEOS-Chem by looking at the FJX_j2j.dat, species_database.yml, and fullchem.eqn files. gcpy/examples/__init__.py CHANGELOG.md - Updated accordingly NOTE: Needs a further commit to add argument parsing. Signed-off-by: Bob Yantosca --- CHANGELOG.md | 8 +- gcpy/examples/__init__.py | 1 + gcpy/examples/geos-chem/check_photolysis.py | 159 ++++++++++++++++++++ 3 files changed, 165 insertions(+), 3 deletions(-) create mode 100644 gcpy/examples/geos-chem/check_photolysis.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e10005a..8675c7e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ## [Unreleased] - TBD ### Added -- Example script `gcpy/examples/hemco/make_hemco_sa_spec.py` (creates the HEMCO standalone configuration file `HEMCO_sa_Spec.rc`) -- Module `benchmark_gcclassic_stats.py` for scraping statistics from GEOS-Chem Classic cloud benchmarks -- Dry deposition velocity comparison plots in 1-month cloud benchmarks +- Added example script `gcpy/examples/hemco/make_hemco_sa_spec.py` (creates the HEMCO standalone configuration file `HEMCO_sa_Spec.rc`) +- Added module `benchmark_gcclassic_stats.py` for scraping statistics from GEOS-Chem Classic cloud benchmarks +- Added dry deposition velocity comparison plots in 1-month cloud benchmarks +- Added `gcpy.examples.geos-chem.check_photolysis.py` script to check consistency of photolysis reactions ### Changed - Changed format of `% diff` column from `12.3e` to `12.3f` in benchmark timing tables @@ -26,6 +27,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - In environment files `read_the_docs_environment.yml` and `read_the_docs_requirements.txt` - Update `jinja` to 3.1.4 (fixes a security issue) - Update `gcpy/setup.py` with the new Python package version numbers +- Updated `gcpy/examples/__init__.py` with to import from `geos-chem` folder ### Fixed - Fixed formatting error in `.github/workflows/stale.yml` that caused the Mark Stale Issues action not to run diff --git a/gcpy/examples/__init__.py b/gcpy/examples/__init__.py index 1d23ef3d..b2aef2c4 100644 --- a/gcpy/examples/__init__.py +++ b/gcpy/examples/__init__.py @@ -5,6 +5,7 @@ #from .bpch_to_nc import * from .diagnostics import * from .dry_run import * +from .geos-chem import * from .plotting import * from .timeseries import * from .working_with_files import * diff --git a/gcpy/examples/geos-chem/check_photolysis.py b/gcpy/examples/geos-chem/check_photolysis.py new file mode 100644 index 00000000..a87683f1 --- /dev/null +++ b/gcpy/examples/geos-chem/check_photolysis.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python3 +""" +Convenience script to verify that photolysis reactions are specified +properly in GEOS-Chem. +""" + +from gcpy.util import read_config_file, verify_variable_type + +def get_photo_species_from_species_db(path): + """ + Returns a list of species that have the "Is_Photolysis: true" + setting in the species_database.yml file. + + Args + path : str : Absolute path to the species_database.yml file + + Returns + species : dict : Names of species having "Is_Photolysis: true" + """ + verify_variable_type(path, str) + + species_database = read_config_file(path) + species = {} + for spc in species_database: + if "Is_Photolysis" in species_database[spc]: + species[spc] = bool(species_database[spc]["Is_Photolysis"]) + + return species + + +def get_photo_species_from_fjx_j2j(path): + """ + Returns a list of photolysis species in the FJX_j2j.dat file. + + Args + path : str : Absolute path to the FJX_j2j.dat file + + Returns + species : list : Tuples containing species name and index + """ + verify_variable_type(path, str) + + species = [] + indices = [] + with open(path, "r", encoding="utf-8") as ifile: + for line in list(ifile): + line = line.strip("\n") + if "PHOTON" in line: + sub_string = line.split() + indices.append(int(sub_string[0])) + species.append(sub_string[1]) + + return list(zip(species, indices)) + + +def get_photo_species_from_kpp_eqn(path): + """ + Returns a list of species having photolysis reactions in + the KPP fullchem.eqn dat file. + + Args + path : str : Absolute path to the fullchem.eqn file + + Returns + species : dict : Species name plus index of the PHOTOL array + """ + verify_variable_type(path, str) + + species = [] + photol = [] + with open(path, "r", encoding="utf-8") as ifile: + for line in list(ifile): + line = line.strip("\n") + name = "" + if "+ hv" in line: + species.append(line.split()[0]) + if "PHOTOL(" in line: + photol.append(int( + line.split(":")[1].split(";")[0]\ + .strip().strip("PHOTOL(").strip(")") + )) + + return list(zip(species, photol)) + + +def get_max_index(multi_list): + """ + Computes the maximum index value in a "multi-index" list + (i.e. a list of tuples). + + Args + multi_list : list : Tuples containing species name and index + + Returns + max_index : int : Maximum index value + """ + max_index = 0 + for (species, index) in multi_list: + if index > max_index: + max_index = index + + return index + + +def print_results(spc_db, fjx_j2j, kpp_eqn): + """ + Prints results + + Args + spc_db : dict : Photolyzing species listed in species_database.yml + fjx_j2j : dict : Photolyzing species listed in FJX_j2j.dat + kpp_eqn : dict : Photolyzing species listed in fullchem.eqn + """ + + # Maximum number of entries in FJX_j2j.dat + max_j2j = get_max_index(fjx_j2j) + max_kpp = get_max_index(kpp_eqn) + print(f"Number of entries in FJX_j2j.dat : {max_j2j:<3}") + print(f"Number of photo rxns in fullchem.eqn : {max_kpp:<3}") + if max_kpp != max_j2j: + msg = "ERROR: Mismatch between FJX_j2j.dat and fullchem.eqn!" + raise ValueError(msg) + + # Print photolyzing species and the corresponding index + # of the KPP PHOTOL array that it uses. Ths should be the + # same order as in the FJX_j2j.dat file. + for (species, index) in kpp_eqn: + if index > 0 and index <= max_j2j: + msg = f"{species:<10} uses PHOTOL({index:<3}) " + msg += f"and has Is_Photolysis={spc_db[species]}" + print(msg) + else: + msg = "{species:<10} has an invalid PHOTOL index " + msg += "({index}) > {max_j2j}!" + print(msg) + + +def main(): + """ + """ + + # Get all species with "Is_Photolysis: true" in the species database + path = "/n/holyscratch01/jacob_lab/ryantosca/tests/14.5.0/gc2457/GCClassic/src/GEOS-Chem/run/shared/species_database.yml" + spc_db = get_photo_species_from_species_db(path) + + # Get all species listed in FJX_j2j.dat + path = "/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/CHEM_INPUTS/CLOUD_J/v2024-09/FJX_j2j.dat" + fjx_j2j = get_photo_species_from_fjx_j2j(path) + + # Get all species with a photo reaction in the f + path = "/n/holyscratch01/jacob_lab/ryantosca/tests/14.5.0/gc2457/GCClassic/src/GEOS-Chem/KPP/fullchem/fullchem.eqn" + kpp_eqn = get_photo_species_from_kpp_eqn(path) + + # Print results + print_results(spc_db, fjx_j2j, kpp_eqn) + + +if __name__ == '__main__': + main() From d4ed2ea51fdde3576413009c64e5fc690989a455 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 12 Sep 2024 12:10:10 -0400 Subject: [PATCH 2/2] Moved check_photolysis.py to examples/geoschem; added arg parsing gcpy/examples/__init__.py - Change geos-chem folder name to geoschem to avoid error gcpy/examples/geoschem/check_photolysis.py - Moved gcpy/examples/geos-chem/check_photolysis.py here - Added argument parsing via argparse library - Updated comments and trimmed trailing whitespace Signed-off-by: Bob Yantosca --- gcpy/examples/__init__.py | 2 +- .../check_photolysis.py | 78 +++++++++++++++---- 2 files changed, 62 insertions(+), 18 deletions(-) rename gcpy/examples/{geos-chem => geoschem}/check_photolysis.py (71%) diff --git a/gcpy/examples/__init__.py b/gcpy/examples/__init__.py index b2aef2c4..16dacf13 100644 --- a/gcpy/examples/__init__.py +++ b/gcpy/examples/__init__.py @@ -5,7 +5,7 @@ #from .bpch_to_nc import * from .diagnostics import * from .dry_run import * -from .geos-chem import * +from .geoschem import * from .plotting import * from .timeseries import * from .working_with_files import * diff --git a/gcpy/examples/geos-chem/check_photolysis.py b/gcpy/examples/geoschem/check_photolysis.py similarity index 71% rename from gcpy/examples/geos-chem/check_photolysis.py rename to gcpy/examples/geoschem/check_photolysis.py index a87683f1..3f4c5d80 100644 --- a/gcpy/examples/geos-chem/check_photolysis.py +++ b/gcpy/examples/geoschem/check_photolysis.py @@ -2,15 +2,21 @@ """ Convenience script to verify that photolysis reactions are specified properly in GEOS-Chem. -""" +Calling sequence: +$ python -m gcpy.examples.geoschem.check_photolysis \ + --spc-db-path /path/to/species_database.yml \ + --fjx-j2j-path /path/to/FJX_j2j.dat \ + --kpp-eqn-path /path/to/fullchem.eqn (or custom.eqn) +""" +import argparse from gcpy.util import read_config_file, verify_variable_type def get_photo_species_from_species_db(path): """ Returns a list of species that have the "Is_Photolysis: true" setting in the species_database.yml file. - + Args path : str : Absolute path to the species_database.yml file @@ -65,7 +71,7 @@ def get_photo_species_from_kpp_eqn(path): species : dict : Species name plus index of the PHOTOL array """ verify_variable_type(path, str) - + species = [] photol = [] with open(path, "r", encoding="utf-8") as ifile: @@ -79,7 +85,7 @@ def get_photo_species_from_kpp_eqn(path): line.split(":")[1].split(";")[0]\ .strip().strip("PHOTOL(").strip(")") )) - + return list(zip(species, photol)) @@ -98,13 +104,13 @@ def get_max_index(multi_list): for (species, index) in multi_list: if index > max_index: max_index = index - + return index def print_results(spc_db, fjx_j2j, kpp_eqn): """ - Prints results + Prints results Args spc_db : dict : Photolyzing species listed in species_database.yml @@ -120,7 +126,7 @@ def print_results(spc_db, fjx_j2j, kpp_eqn): if max_kpp != max_j2j: msg = "ERROR: Mismatch between FJX_j2j.dat and fullchem.eqn!" raise ValueError(msg) - + # Print photolyzing species and the corresponding index # of the KPP PHOTOL array that it uses. Ths should be the # same order as in the FJX_j2j.dat file. @@ -135,25 +141,63 @@ def print_results(spc_db, fjx_j2j, kpp_eqn): print(msg) -def main(): +def check_geos_chem_photolysis(spc_db_path, fjx_j2j_path, kpp_eqn_path): """ """ - + # Get all species with "Is_Photolysis: true" in the species database - path = "/n/holyscratch01/jacob_lab/ryantosca/tests/14.5.0/gc2457/GCClassic/src/GEOS-Chem/run/shared/species_database.yml" - spc_db = get_photo_species_from_species_db(path) + spc_db = get_photo_species_from_species_db(spc_db_path) # Get all species listed in FJX_j2j.dat - path = "/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/data/ExtData/CHEM_INPUTS/CLOUD_J/v2024-09/FJX_j2j.dat" - fjx_j2j = get_photo_species_from_fjx_j2j(path) + fjx_j2j = get_photo_species_from_fjx_j2j(fjx_j2j_path) - # Get all species with a photo reaction in the f - path = "/n/holyscratch01/jacob_lab/ryantosca/tests/14.5.0/gc2457/GCClassic/src/GEOS-Chem/KPP/fullchem/fullchem.eqn" - kpp_eqn = get_photo_species_from_kpp_eqn(path) + # Get all species with a listed photo reaction + # in the the KPP fullchem.eqn ir custom.eqn file + kpp_eqn = get_photo_species_from_kpp_eqn(kpp_eqn_path) # Print results print_results(spc_db, fjx_j2j, kpp_eqn) - + + +def main(): + """ + Parses command-line arguments and calls check_geos_chem_photolysis. + """ + # Tell the parser which arguments to look for + parser = argparse.ArgumentParser( + description="Single-panel plotting example program" + ) + parser.add_argument( + "--spc-db-path", + metavar="spc_db_path", + type=str, + required=True, + help="Path to the species_database.yml file" + ) + parser.add_argument( + "--fjx-j2j-path", + metavar="fjx_j2j_path", + type=str, + required=True, + help="Path to the FJX_j2j.dat file" + ) + parser.add_argument( + "--kpp-eqn-path", + metavar="kpp_eqn_path", + type=str, + required=True, + help="Path to the KPP fullchem.eqn or custom.eqn file" + ) + + # Parse command-line arguments + args = parser.parse_args() + + # Call the plot_single_panel routine + check_geos_chem_photolysis( + args.spc_db_path, + args.fjx_j2j_path, + args.kpp_eqn_path, + ) if __name__ == '__main__': main()