diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py index 1f45d02da..17483e71f 100644 --- a/examples/all-sky/make_problem_size_loop.py +++ b/examples/all-sky/make_problem_size_loop.py @@ -3,15 +3,14 @@ # This script... # import argparse -import csv -import os +import csv -# template: exe ncol nlay nloops output_file k_dist clouds aeorsols +# template: exe ncol nlay nloops output_file k_dist clouds aeorsols -# specify: kdist, optional clouds, aerosols. specify nloops -# Toggle clouds and aerosols? -# Loop over sets of ncol, nlay, -# output name +# specify: kdist, optional clouds, aerosols. specify nloops +# Toggle clouds and aerosols? +# Loop over sets of ncol, nlay, +# output name if __name__ == '__main__': @@ -19,44 +18,47 @@ description="Description here ") # Argument parseing described at # https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse - parser.add_argument("-x", "--executable", - type=str, + parser.add_argument("-x", "--executable", + type=str, default="./rrtmgp_allsky", help="Path to exectuable") - parser.add_argument("-c", "--ncol", - type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + parser.add_argument("-c", "--ncol", + type=lambda items: [int(i) for i in + list(csv.reader([items]))[0]], default="2,4,8,16", help="Number of columns (multiple e.g. 2,4,8,16)") - parser.add_argument("-l", "--nlay", - type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + parser.add_argument("-l", "--nlay", + type=lambda items: [int(i) for i in + list(csv.reader([items]))[0]], default="32, 64, 96", help="Number of layers (multiple e.g. 32,64.96)") - parser.add_argument("-i", "--nloops", + parser.add_argument("-i", "--nloops", type=int, default=1, help="Number of loops (same for all ncol)") - parser.add_argument("-o", "--output_file", - type=str, + parser.add_argument("-o", "--output_file", + type=str, default="rrtmgp-allsky-output.nc", help="Path to output file") - parser.add_argument("-k", "--k_distribution", - type=str, - required = True, + parser.add_argument("-k", "--k_distribution", + type=str, + required=True, help="Path to gas optics file [required]") - parser.add_argument("-cl", "--cloud_optics", + parser.add_argument("-cl", "--cloud_optics", type=str, default="", help="Path to cloud optics file") - parser.add_argument("-a", "--aerosol_optics", + parser.add_argument("-a", "--aerosol_optics", type=str, default="", help="Path to aerosol optics file") args = parser.parse_args() - # Can't supply aerosols without clouds - if(args.cloud_optics == "" and args.aerosol_optics != ""): - raise AssertionError("Need to supply cloud optics if providing aerosol optics") + # Can't supply aerosols without clouds + if args.cloud_optics == "" and args.aerosol_optics != "": + raise AssertionError( + "Need to supply cloud optics if providing aerosol optics") - # Every combo of ncol, nlay - for l in args.nlay: - for i in args.ncol: - print(f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + \ - f"{args.output_file} {args.k_distribution}" + \ - f"{args.cloud_optics} {args.aerosol_optics}") + # Every combo of ncol, nlay + for l in args.nlay: + for i in args.ncol: + print(f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + + f"{args.output_file} {args.k_distribution}" + + f"{args.cloud_optics} {args.aerosol_optics}") diff --git a/examples/all-sky/run-allsky-example.py b/examples/all-sky/run-allsky-example.py index 2b7f604a4..dd9b95a25 100644 --- a/examples/all-sky/run-allsky-example.py +++ b/examples/all-sky/run-allsky-example.py @@ -3,9 +3,7 @@ # This script runs runs the RTE+RRTMGP all-sky examples # import argparse -import glob import os -import shutil import subprocess rte_rrtmgp_dir = os.environ["RRTMGP_DATA"] @@ -38,8 +36,8 @@ "cloudy columns") args = parser.parse_args() - ncol_str = '{0:5d}'.format(args.ncol) - nlay_str = '{0:5d}'.format(args.nlay) + ncol_str = '{0:5d}'.format(args.ncol) + nlay_str = '{0:5d}'.format(args.nlay) nloops_str = '{0:5d}'.format(args.nloops) if args.run_command: print("using the run command") @@ -47,8 +45,12 @@ os.chdir(all_sky_dir) # Remove cloudy-sky fluxes from the file containing the atmospheric profiles subprocess.run( - [all_sky_exe_name, ncol_str, nlay_str, nloops_str, "rrtmgp-allsky-lw-no-aerosols.nc", lw_gas_coeffs_file, lw_clouds_coeff_file]) + [all_sky_exe_name, ncol_str, nlay_str, nloops_str, + "rrtmgp-allsky-lw-no-aerosols.nc", lw_gas_coeffs_file, + lw_clouds_coeff_file]) subprocess.run( - [all_sky_exe_name, ncol_str, nlay_str, nloops_str, "rrtmgp-allsky-sw-no-aerosols.nc", sw_gas_coeffs_file, sw_clouds_coeff_file]) + [all_sky_exe_name, ncol_str, nlay_str, nloops_str, + "rrtmgp-allsky-sw-no-aerosols.nc", sw_gas_coeffs_file, + sw_clouds_coeff_file]) # end main diff --git a/examples/compare-to-reference.py b/examples/compare-to-reference.py index 702efc4bc..bb7d847a3 100644 --- a/examples/compare-to-reference.py +++ b/examples/compare-to-reference.py @@ -9,43 +9,45 @@ # Thresholds come from environement variables when set? # # -import argparse -import os import sys -import warnings +import argparse import numpy as np +import os +import warnings import xarray as xr -REPORT = os.getenv("REPORTING_THRESHOLD") if os.getenv("REPORTING_THRESHOLD") is not None else 0. -FAILURE = os.getenv("FAILURE_THRESHOLD") if os.getenv("FAILURE_THRESHOLD") is not None else 1.e-5 # # Comparing reference and test results # if __name__ == '__main__': warnings.simplefilter("ignore", xr.SerializationWarning) - + parser = argparse.ArgumentParser( description="Compares example output to file in reference directory") - parser.add_argument("--ref_dir", type=str, + parser.add_argument("--ref_dir", type=str, help="Directory containing reference files") - parser.add_argument("--tst_dir", type=str, default = ".", + parser.add_argument("--tst_dir", type=str, default=".", help="Directory contining test values") - parser.add_argument("--file_names", type=str, nargs='+', default = [], + parser.add_argument("--file_names", type=str, nargs='+', default=[], help="Name[s] of files to compare") - parser.add_argument("--variables", type=str, nargs='+', default = None, + parser.add_argument("--variables", type=str, nargs='+', default=None, help="Name[s] of files to compare") - parser.add_argument("--report_threshold", type=float, - default=REPORT, + parser.add_argument("--report_threshold", type=float, + default=os.getenv("REPORTING_THRESHOLD", 0.), help="Threshold for reporting differences") - parser.add_argument("--failure_threshold", type=float, - default=FAILURE, + parser.add_argument("--failure_threshold", type=float, + default=os.getenv("FAILURE_THRESHOLD", 1.e-5), help="Threshold at which differences cause failure " "(for continuous integration)") args = parser.parse_args() - tst = xr.open_mfdataset([os.path.join(args.tst_dir, f) for f in args.file_names],combine='by_coords') - ref = xr.open_mfdataset([os.path.join(args.ref_dir, f) for f in args.file_names],combine='by_coords') + tst = xr.open_mfdataset( + [os.path.join(args.tst_dir, f) for f in args.file_names], + combine='by_coords') + ref = xr.open_mfdataset( + [os.path.join(args.ref_dir, f) for f in args.file_names], + combine='by_coords') variables = args.variables if args.variables is not None else tst.variables failed = False @@ -62,23 +64,25 @@ # # Reporting # - if not np.allclose(tst[v], ref[v], atol=args.report_threshold, rtol=0): - diff = abs( (tst - ref)[v].values) + if not np.allclose(tst[v], ref[v], atol=args.report_threshold, rtol=0): + diff = abs((tst - ref)[v].values) avg = 0.5 * (tst + ref)[v].values - # Division raises a runtime warning when we divide by zero even if the - # values in those locations will be ignored. + # Division raises a runtime warning when we divide by zero even if + # the values in those locations will be ignored. with np.errstate(divide='ignore', invalid='ignore'): frac_diff = np.where( (avg > 2. * np.finfo(float).eps), diff / avg, 0) print('Variable %s differs (max abs difference: %e; ' - 'max percent difference: %e%%)' % (v, diff.max(), 100.0 * frac_diff.max())) + 'max percent difference: %e%%)' % ( + v, diff.max(), 100.0 * frac_diff.max())) else: print('Variable %s: No diffs' % v) # # Failure # - if not np.allclose(tst[v], ref[v], atol=args.failure_threshold, rtol=0): failed = True + if not np.allclose(tst[v], ref[v], atol=args.failure_threshold, + rtol=0): failed = True - if failed: print ("Tests failed") + if failed: + print("Tests failed") sys.exit(1) if failed else sys.exit(0) - diff --git a/examples/rfmip-clear-sky/run-rfmip-examples.py b/examples/rfmip-clear-sky/run-rfmip-examples.py index 3ac480b5f..c747cefd1 100755 --- a/examples/rfmip-clear-sky/run-rfmip-examples.py +++ b/examples/rfmip-clear-sky/run-rfmip-examples.py @@ -3,7 +3,6 @@ # This script runs runs RRTMGP on the RFMIP off-line test cases # import argparse -import glob import os import subprocess @@ -15,11 +14,13 @@ rfmip_dir = "." # Code should be run in the rfmip_dir directory -conds_file = os.path.join(os.environ["RRTMGP_DATA"], - "examples", "rfmip-clear-sky", "inputs", - "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc") -lw_gas_coeffs_file = os.path.join(os.environ["RRTMGP_DATA"], "rrtmgp-gas-lw-g256.nc") -sw_gas_coeffs_file = os.path.join(os.environ["RRTMGP_DATA"], "rrtmgp-gas-sw-g224.nc") +conds_file = os.path.join( + os.environ["RRTMGP_DATA"], "examples", "rfmip-clear-sky", "inputs", + "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc") +lw_gas_coeffs_file = os.path.join( + os.environ["RRTMGP_DATA"], "rrtmgp-gas-lw-g256.nc") +sw_gas_coeffs_file = os.path.join( + os.environ["RRTMGP_DATA"], "rrtmgp-gas-sw-g224.nc") rfmip_lw_exe_name = "./rrtmgp_rfmip_lw" rfmip_sw_exe_name = "./rrtmgp_rfmip_sw" diff --git a/examples/rfmip-clear-sky/stage_files.py b/examples/rfmip-clear-sky/stage_files.py index e9efc24d8..b23a6986e 100755 --- a/examples/rfmip-clear-sky/stage_files.py +++ b/examples/rfmip-clear-sky/stage_files.py @@ -3,8 +3,9 @@ # This script downloads and/or creates files needed for the RFMIP off-line test # cases # -import subprocess import sys + +import subprocess import urllib.request from pathlib import Path @@ -39,13 +40,6 @@ print("Downloading RFMIP input files") urllib.request.urlretrieve(conds_url, conds_file) -generate_templates = False -if generate_templates: - print("Downloading scripts for generating output templates") - urllib.request.urlretrieve(templ_scr_url, templ_scr) - subprocess.run( - [sys.executable, templ_scr, "--source_id", "RTE-RRTMGP-181204"]) -else: - print("Downloading output templates") - for f in output_files: - urllib.request.urlretrieve(conds_url.replace(conds_file, f), f) +print("Downloading output templates") +for f in output_files: + urllib.request.urlretrieve(conds_url.replace(conds_file, f), f) diff --git a/tests/validation-plots.py b/tests/validation-plots.py index ea121efdf..4dcc9434b 100644 --- a/tests/validation-plots.py +++ b/tests/validation-plots.py @@ -1,10 +1,9 @@ -import warnings -import urllib.request import colorcet as cc import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np -import seaborn as sb +import urllib.request +import warnings import xarray as xr from matplotlib.backends.backend_pdf import PdfPages @@ -23,35 +22,37 @@ def rms(diff, col_dim): return np.sqrt(np.square(diff).mean(dim=col_dim)) -def make_comparison_plot(variants, labels, reference, vscale, col_dim="site", - lay_dim="layer"): +def make_comparison_plot(variants, labels, reference, vscale, colors, + col_dim="site"): # # Make a plot comparing differences with respect to reference # if type(variants) is not list: - make_comparison_plot([variants], [labels], reference, vscale) + make_comparison_plot([variants], [labels], reference, vscale, colors, + col_dim) else: for i in np.arange(0, len(variants)): delta = variants[i] - reference plt.plot(mae(delta, col_dim), vscale, '-', - color=cols[i], label=labels[i]) + color=colors[i], label=labels[i]) plt.plot(rms(delta, col_dim), vscale, '--', - color=cols[i]) + color=colors[i]) # Reverse vertical ordering plt.legend() # Reverse vertical ordering plt.ylim(vscale.max(), vscale.min()) + def construct_lbl_esgf_root(var, esgf_node="llnl"): # # For a given variable name, provide the https URL for the LBLRM # line-by-line results # - model="LBLRTM-12-8" - prefix = ("http://esgf3.dkrz.de/thredds/fileServer/cmip6/RFMIP/AER/" + model + - "/rad-irf/r1i1p1f1/Efx/") + model = "LBLRTM-12-8" + prefix = ("http://esgf3.dkrz.de/thredds/fileServer/cmip6/RFMIP/AER/" + + model + "/rad-irf/r1i1p1f1/Efx/") if esgf_node == "llnl": prefix = ("http://esgf-data1.llnl.gov/thredds/fileServer/css03_data/" "CMIP6/RFMIP/AER/" + model + "/rad-irf/r1i1p1f1/Efx/") @@ -60,21 +61,25 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): ######################################################################## -if __name__ == '__main__': +def main(): warnings.simplefilter("ignore", xr.SerializationWarning) # - # Reference values from LBLRTM - download locally, since OpenDAP access is so inconsistent + # Reference values from LBLRTM - download locally, since OpenDAP access is + # so inconsistent # fluxes = ["rsd", "rsu", "rld", "rlu"] lbl_suffix = "_Efx_LBLRTM-12-8_rad-irf_r1i1p1f1_gn.nc" for v in fluxes: - try: + try: try: - urllib.request.urlretrieve(construct_lbl_esgf_root(v), v+lbl_suffix) + urllib.request.urlretrieve(construct_lbl_esgf_root(v), + v + lbl_suffix) except: - urllib.request.urlretrieve(construct_lbl_esgf_root(v), v+lbl_suffix, node="dkrz") + urllib.request.urlretrieve( + construct_lbl_esgf_root(v, esgf_node="dkrz"), + v + lbl_suffix) except: - raise Exception("Failed to download {0}".format(v+lbl_suffix)) + raise Exception("Failed to download {0}".format(v + lbl_suffix)) lbl = xr.open_mfdataset([v + lbl_suffix for v in fluxes], combine="by_coords").sel(expt=0) @@ -145,7 +150,8 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): "fewer points + optimal-angle", "3-angle", "2-stream", "rescaled"], reference=r, - vscale=plev / 100.) + vscale=plev / 100., + colors=cols) plt.ylabel("Pressure (Pa)") plt.xlabel("Error (W/m$^2$), solid=mean, dash=RMS") plt.title(t) @@ -167,7 +173,8 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): make_comparison_plot(v, labels=["no-tlev", "2str"], reference=r, - vscale=plev / 100.) + vscale=plev / 100., + colors=cols) plt.ylabel("Pressure (Pa)") plt.xlabel("Difference (W/m$^2$), solid=mean, dash=RMS") plt.title(t) @@ -189,9 +196,14 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): make_comparison_plot(v, labels=["default", "fewer-g-points"], reference=r, - vscale=plev / 100.) + vscale=plev / 100., + colors=cols) plt.ylabel("Pressure (Pa)") plt.xlabel("Error (W/m$^2$), solid=mean, dash=RMS") plt.title(t) pdf.savefig() plt.close() + + +if __name__ == '__main__': + main()