From 5dfdca97c33017db8891a1fa5581029c99f3d5de Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Tue, 3 Dec 2024 16:41:44 -0700 Subject: [PATCH 1/2] EAMxx: add a grids manager from pre-built grids This GM is useful for unit tests, but more generally for situations where only grids are avaialble, but interfaces require a grids manager. --- .../dynamics/homme/homme_grids_manager.cpp | 4 +- .../eamxx/src/share/grid/grids_manager.cpp | 30 ++++++++++-- .../eamxx/src/share/grid/grids_manager.hpp | 3 +- .../src/share/grid/library_grids_manager.hpp | 49 +++++++++++++++++++ .../share/grid/mesh_free_grids_manager.cpp | 4 +- 5 files changed, 81 insertions(+), 9 deletions(-) create mode 100644 components/eamxx/src/share/grid/library_grids_manager.hpp diff --git a/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp b/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp index df5de6827f69..7376e9af9dde 100644 --- a/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp +++ b/components/eamxx/src/dynamics/homme/homme_grids_manager.cpp @@ -188,7 +188,7 @@ void HommeGridsManager::build_dynamics_grid () { initialize_vertical_coordinates(dyn_grid); dyn_grid->m_short_name = "dyn"; - add_grid(dyn_grid); + add_nonconst_grid(dyn_grid); } void HommeGridsManager:: @@ -307,7 +307,7 @@ build_physics_grid (const ci_string& type, const ci_string& rebalance) { } phys_grid->m_short_name = type; - add_grid(phys_grid); + add_nonconst_grid(phys_grid); } void HommeGridsManager:: diff --git a/components/eamxx/src/share/grid/grids_manager.cpp b/components/eamxx/src/share/grid/grids_manager.cpp index 77300aa51c90..45494947acad 100644 --- a/components/eamxx/src/share/grid/grids_manager.cpp +++ b/components/eamxx/src/share/grid/grids_manager.cpp @@ -7,7 +7,24 @@ auto GridsManager:: get_grid(const std::string& name) const -> grid_ptr_type { - auto g = get_grid_nonconst(name); + EKAT_REQUIRE_MSG (has_grid(name), + "Error! Grids manager '" + this->name() + "' does not provide grid '" + name + "'.\n" + " Avaialble grids are: " + print_available_grids() + "\n"); + + grid_ptr_type g; + for (const auto& it : m_grids) { + if (it.second->name()==name or + ekat::contains(it.second->aliases(),name)) { + g = it.second; + break; + } + } + + EKAT_REQUIRE_MSG (g!=nullptr, + "Something went wrong while looking up a grid.\n" + " - grids manager: " + this->name() + "\n" + " - grid name : " + name + "\n"); + return g; } @@ -51,7 +68,14 @@ create_remapper (const grid_ptr_type& from_grid, } void GridsManager:: -add_grid (nonconstgrid_ptr_type grid) +add_nonconst_grid (nonconstgrid_ptr_type grid) +{ + add_grid(grid); + m_nonconst_grids[grid->name()] = grid; +} + +void GridsManager:: +add_grid (grid_ptr_type grid) { const auto& name = grid->name(); EKAT_REQUIRE_MSG (not has_grid(name), @@ -59,7 +83,6 @@ add_grid (nonconstgrid_ptr_type grid) " - grids manager: " + this->name() + "\n" " - grid name : " + name + "\n"); - m_nonconst_grids[name] = grid; m_grids[name] = grid; } @@ -86,7 +109,6 @@ get_grid_nonconst (const std::string& name) const " - grid name : " + name + "\n"); return g; - } void GridsManager:: diff --git a/components/eamxx/src/share/grid/grids_manager.hpp b/components/eamxx/src/share/grid/grids_manager.hpp index 4a0e25bc233c..e70bb9678f01 100644 --- a/components/eamxx/src/share/grid/grids_manager.hpp +++ b/components/eamxx/src/share/grid/grids_manager.hpp @@ -56,7 +56,8 @@ class GridsManager protected: - void add_grid (nonconstgrid_ptr_type grid); + void add_nonconst_grid (nonconstgrid_ptr_type grid); + void add_grid (grid_ptr_type grid); void alias_grid (const std::string& grid_name, const std::string& grid_alias); virtual remapper_ptr_type diff --git a/components/eamxx/src/share/grid/library_grids_manager.hpp b/components/eamxx/src/share/grid/library_grids_manager.hpp new file mode 100644 index 000000000000..fd5d3caaef22 --- /dev/null +++ b/components/eamxx/src/share/grid/library_grids_manager.hpp @@ -0,0 +1,49 @@ +#ifndef EAMXX_LIBRARY_GRIDS_MANAGER_HPP +#define EAMXX_LIBRARY_GRIDS_MANAGER_HPP + +#include "share/grid/grids_manager.hpp" + +namespace scream { + +// This class is meant to be used within scopes that need a grids manager +// object, but they only have pre-built grids. The user can then simply +// create a LibraryGridsManager, and add the pre-existing grids. Afterwards, +// it's business as usual with GridsManager's interfaces. + +class LibraryGridsManager : public GridsManager +{ +public: + template + explicit LibraryGridsManager(Pointers&&... ptrs) { + add_grids(std::forward(ptrs)...); + } + + virtual ~LibraryGridsManager () = default; + + std::string name () const { return "Library grids_manager"; } + + void build_grids () override {} + + void add_grids () {} + + template + void add_grids (grid_ptr_type p, Pointers&&... ptrs) { + add_grid(p); + add_grids(std::forward(ptrs)...); + } + +protected: + remapper_ptr_type + do_create_remapper (const grid_ptr_type from_grid, + const grid_ptr_type to_grid) const + { + EKAT_ERROR_MSG ( + "Error! LibraryGridsManager is not capable of creating remappers.\n" + " - from_grid: " + from_grid->name() + "\n" + " - to_grid: " + to_grid->name() + "\n"); + } +}; + +} // namespace scream + +#endif // EAMXX_LIBRARY_GRIDS_MANAGER_HPP diff --git a/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp b/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp index b083dfcfa48b..4d505803e66f 100644 --- a/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp +++ b/components/eamxx/src/share/grid/mesh_free_grids_manager.cpp @@ -106,7 +106,7 @@ build_se_grid (const std::string& name, ekat::ParameterList& params) se_grid->m_short_name = "se"; add_geo_data(se_grid); - add_grid(se_grid); + add_nonconst_grid(se_grid); } void MeshFreeGridsManager:: @@ -132,7 +132,7 @@ build_point_grid (const std::string& name, ekat::ParameterList& params) add_geo_data(pt_grid); pt_grid->m_short_name = "pt"; - add_grid(pt_grid); + add_nonconst_grid(pt_grid); } void MeshFreeGridsManager:: From a851099f111f6cabd7cf225ce1704e2bda65a597 Mon Sep 17 00:00:00 2001 From: Luca Bertagna Date: Tue, 3 Dec 2024 16:44:12 -0700 Subject: [PATCH 2/2] EAMxx: factor diagnostics creation out of IO class * Add free function to build a diagnostic obj from a diag field name * Adapt existing diags so that their name is "static", rather than matching the diag field name * Adapt AtmosphereOutput to use this utility. Enhances encapsulation, and avoids hiding the diag naming convention deep in the IO class --- .../eamxx/src/diagnostics/CMakeLists.txt | 10 +- .../eamxx/src/diagnostics/aerocom_cld.cpp | 5 +- .../eamxx/src/diagnostics/aerocom_cld.hpp | 2 +- .../eamxx/src/diagnostics/field_at_height.cpp | 23 +-- .../eamxx/src/diagnostics/field_at_height.hpp | 2 +- .../eamxx/src/diagnostics/field_at_level.hpp | 2 +- .../diagnostics/field_at_pressure_level.cpp | 32 ++-- .../diagnostics/field_at_pressure_level.hpp | 2 +- .../eamxx/src/diagnostics/number_path.cpp | 4 +- .../eamxx/src/diagnostics/number_path.hpp | 2 +- .../src/diagnostics/potential_temperature.cpp | 7 +- .../src/diagnostics/potential_temperature.hpp | 2 +- .../src/diagnostics/precip_surf_mass_flux.cpp | 2 +- .../src/diagnostics/precip_surf_mass_flux.hpp | 2 +- .../src/diagnostics/register_diagnostics.hpp | 9 +- .../src/diagnostics/tests/CMakeLists.txt | 94 +++++----- .../tests/field_at_height_tests.cpp | 28 ++- .../tests/field_at_pressure_level_tests.cpp | 3 +- .../tests/vertical_layer_tests.cpp | 14 +- .../eamxx/src/diagnostics/vapor_flux.cpp | 8 +- .../eamxx/src/diagnostics/vapor_flux.hpp | 2 +- .../eamxx/src/diagnostics/vertical_layer.cpp | 28 +-- .../eamxx/src/diagnostics/vertical_layer.hpp | 2 +- .../eamxx/src/diagnostics/water_path.cpp | 7 +- .../eamxx/src/diagnostics/water_path.hpp | 2 +- components/eamxx/src/share/io/CMakeLists.txt | 2 +- .../eamxx/src/share/io/scorpio_output.cpp | 148 +++------------- .../eamxx/src/share/io/scream_io_utils.cpp | 87 +++++++++ .../eamxx/src/share/io/scream_io_utils.hpp | 9 + .../eamxx/src/share/io/tests/CMakeLists.txt | 6 + .../eamxx/src/share/io/tests/create_diag.cpp | 166 ++++++++++++++++++ 31 files changed, 418 insertions(+), 294 deletions(-) create mode 100644 components/eamxx/src/share/io/tests/create_diag.cpp diff --git a/components/eamxx/src/diagnostics/CMakeLists.txt b/components/eamxx/src/diagnostics/CMakeLists.txt index be51f4346155..8a8dfca560da 100644 --- a/components/eamxx/src/diagnostics/CMakeLists.txt +++ b/components/eamxx/src/diagnostics/CMakeLists.txt @@ -1,4 +1,7 @@ set(DIAGNOSTIC_SRCS + aerocom_cld.cpp + aodvis.cpp + atm_backtend.cpp atm_density.cpp dry_static_energy.cpp exner.cpp @@ -6,6 +9,7 @@ set(DIAGNOSTIC_SRCS field_at_level.cpp field_at_pressure_level.cpp longwave_cloud_forcing.cpp + number_path.cpp potential_temperature.cpp precip_surf_mass_flux.cpp relative_humidity.cpp @@ -17,15 +21,11 @@ set(DIAGNOSTIC_SRCS virtual_temperature.cpp water_path.cpp wind_speed.cpp - aodvis.cpp - number_path.cpp - aerocom_cld.cpp - atm_backtend.cpp ) add_library(diagnostics ${DIAGNOSTIC_SRCS}) target_link_libraries(diagnostics PUBLIC scream_share) -if (NOT SCREAM_LIB_ONLY) +if (NOT SCREAM_LIB_ONLY AND NOT SCREAM_ONLY_GENERATE_BASELINES) add_subdirectory(tests) endif() diff --git a/components/eamxx/src/diagnostics/aerocom_cld.cpp b/components/eamxx/src/diagnostics/aerocom_cld.cpp index 8606e065b9e3..cca55f6c19f1 100644 --- a/components/eamxx/src/diagnostics/aerocom_cld.cpp +++ b/components/eamxx/src/diagnostics/aerocom_cld.cpp @@ -22,8 +22,6 @@ AeroComCld::AeroComCld(const ekat::Comm &comm, "to be 'Bot' or 'Top' in its input parameters.\n"); } -std::string AeroComCld::name() const { return "AeroComCld" + m_topbot; } - void AeroComCld::set_grids( const std::shared_ptr grids_manager) { using namespace ekat::units; @@ -76,7 +74,8 @@ void AeroComCld::set_grids( m_dz.allocate_view(); // Construct and allocate the output field - FieldIdentifier fid(name(), vector1d_layout, nondim, grid_name); + + FieldIdentifier fid("AeroComCld"+m_topbot, vector1d_layout, nondim, grid_name); m_diagnostic_output = Field(fid); m_diagnostic_output.allocate_view(); diff --git a/components/eamxx/src/diagnostics/aerocom_cld.hpp b/components/eamxx/src/diagnostics/aerocom_cld.hpp index 694eed76f097..53a9a7f8e1ba 100644 --- a/components/eamxx/src/diagnostics/aerocom_cld.hpp +++ b/components/eamxx/src/diagnostics/aerocom_cld.hpp @@ -15,7 +15,7 @@ class AeroComCld : public AtmosphereDiagnostic { AeroComCld(const ekat::Comm &comm, const ekat::ParameterList ¶ms); // The name of the diagnostic - std::string name() const override; + std::string name() const override { return "AeroComCld"; } // Set the grid void set_grids( diff --git a/components/eamxx/src/diagnostics/field_at_height.cpp b/components/eamxx/src/diagnostics/field_at_height.cpp index f61cd3a76c1f..317bbffaf620 100644 --- a/components/eamxx/src/diagnostics/field_at_height.cpp +++ b/components/eamxx/src/diagnostics/field_at_height.cpp @@ -45,21 +45,16 @@ FieldAtHeight (const ekat::Comm& comm, const ekat::ParameterList& params) " - surface reference: " + surf_ref + "\n" " - valid options: sealevel, surface\n"); m_z_name = (surf_ref == "sealevel") ? "z" : "height"; - const auto& location = m_params.get("vertical_location"); - auto chars_start = location.find_first_not_of("0123456789."); - EKAT_REQUIRE_MSG (chars_start!=0 && chars_start!=std::string::npos, - "Error! Invalid string for height value for FieldAtHeight.\n" - " - input string : " + location + "\n" - " - expected format: Nm, with N integer\n"); - const auto z_str = location.substr(0,chars_start); - m_z = std::stod(z_str); - - const auto units = location.substr(chars_start); + + const auto units = m_params.get("height_units"); EKAT_REQUIRE_MSG (units=="m", - "Error! Invalid string for height value for FieldAtHeight.\n" - " - input string : " + location + "\n" - " - expected format: Nm, with N integer\n"); - m_diag_name = m_field_name + "_at_" + m_params.get("vertical_location") + "_above_" + surf_ref; + "Error! Invalid units for FieldAtHeight.\n" + " - input units: " + units + "\n" + " - valid units: m\n"); + + auto z_val = m_params.get("height_value"); + m_z = std::stod(z_val); + m_diag_name = m_field_name + "_at_" + z_val + units + "_above_" + surf_ref; } void FieldAtHeight:: diff --git a/components/eamxx/src/diagnostics/field_at_height.hpp b/components/eamxx/src/diagnostics/field_at_height.hpp index e6198153f940..91f6ae3eb1b5 100644 --- a/components/eamxx/src/diagnostics/field_at_height.hpp +++ b/components/eamxx/src/diagnostics/field_at_height.hpp @@ -18,7 +18,7 @@ class FieldAtHeight : public AtmosphereDiagnostic FieldAtHeight (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const { return m_diag_name; } + std::string name () const { return "FieldAtHeight"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/field_at_level.hpp b/components/eamxx/src/diagnostics/field_at_level.hpp index b63cda0a1a38..3ab9bee15574 100644 --- a/components/eamxx/src/diagnostics/field_at_level.hpp +++ b/components/eamxx/src/diagnostics/field_at_level.hpp @@ -21,7 +21,7 @@ class FieldAtLevel : public AtmosphereDiagnostic FieldAtLevel (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const { return m_diag_name; } + std::string name () const { return "FieldAtLevel"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/field_at_pressure_level.cpp b/components/eamxx/src/diagnostics/field_at_pressure_level.cpp index 21c1ac78dd90..716c8f563af4 100644 --- a/components/eamxx/src/diagnostics/field_at_pressure_level.cpp +++ b/components/eamxx/src/diagnostics/field_at_pressure_level.cpp @@ -15,30 +15,22 @@ FieldAtPressureLevel (const ekat::Comm& comm, const ekat::ParameterList& params) { m_field_name = m_params.get("field_name"); - // Figure out the pressure value - const auto& location = m_params.get("vertical_location"); - auto chars_start = location.find_first_not_of("0123456789."); - EKAT_REQUIRE_MSG (chars_start!=0 && chars_start!=std::string::npos, - "Error! Invalid string for pressure value for FieldAtPressureLevel.\n" - " - input string : " + location + "\n" - " - expected format: Nxyz, with N integer, and xyz='mb', 'hPa', or 'Pa'\n"); - const auto press_str = location.substr(0,chars_start); - m_pressure_level = std::stod(press_str); - - const auto units = location.substr(chars_start); + const auto units = m_params.get("pressure_units"); EKAT_REQUIRE_MSG (units=="mb" or units=="hPa" or units=="Pa", - "Error! Invalid string for pressure value for FieldAtPressureLevel.\n" - " - input string : " + location + "\n" - " - expected format: Nxyz, with N integer, and xyz='mb', 'hPa', or 'Pa'\n"); + "Error! Invalid units for FieldAtPressureLevel.\n" + " - input units: " + units + "\n" + " - valid units: 'mb', 'hPa', 'Pa'\n"); + + // Figure out the pressure value, and convert to Pa if needed + auto p_value = m_params.get("pressure_value"); - // Convert pressure level to Pa, the units of pressure in the simulation if (units=="mb" || units=="hPa") { - m_pressure_level *= 100; + m_pressure_level = std::stod(p_value)*100; + } else { + m_pressure_level = std::stod(p_value); } - m_mask_val = m_params.get("mask_value",Real(constants::DefaultFillValue::value)); - - m_diag_name = m_field_name + "_at_" + location; + m_diag_name = m_field_name + "_at_" + p_value + units; } void FieldAtPressureLevel:: @@ -91,6 +83,8 @@ initialize_impl (const RunType /*run_type*/) // Add a field representing the mask as extra data to the diagnostic field. auto nondim = ekat::units::Units::nondimensional(); const auto& gname = fid.get_grid_name(); + m_mask_val = m_params.get("mask_value",Real(constants::DefaultFillValue::value)); + std::string mask_name = name() + " mask"; FieldLayout mask_layout( {COL}, {num_cols}); diff --git a/components/eamxx/src/diagnostics/field_at_pressure_level.hpp b/components/eamxx/src/diagnostics/field_at_pressure_level.hpp index 950c0c5e2ee9..58e476ec83b9 100644 --- a/components/eamxx/src/diagnostics/field_at_pressure_level.hpp +++ b/components/eamxx/src/diagnostics/field_at_pressure_level.hpp @@ -20,7 +20,7 @@ class FieldAtPressureLevel : public AtmosphereDiagnostic FieldAtPressureLevel (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const { return m_diag_name; } + std::string name () const { return "FieldAtPressureLevel"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/number_path.cpp b/components/eamxx/src/diagnostics/number_path.cpp index d4df2f1bc22f..70e18c6dee30 100644 --- a/components/eamxx/src/diagnostics/number_path.cpp +++ b/components/eamxx/src/diagnostics/number_path.cpp @@ -33,8 +33,6 @@ NumberPathDiagnostic::NumberPathDiagnostic(const ekat::Comm &comm, } } -std::string NumberPathDiagnostic::name() const { return m_kind + "NumberPath"; } - void NumberPathDiagnostic::set_grids( const std::shared_ptr grids_manager) { using namespace ekat::units; @@ -55,7 +53,7 @@ void NumberPathDiagnostic::set_grids( add_field(m_nname, scalar3d, 1 / kg, grid_name); // Construct and allocate the diagnostic field - FieldIdentifier fid(name(), scalar2d, kg/(kg*m2), grid_name); + FieldIdentifier fid(m_kind + "NumberPath", scalar2d, kg/(kg*m2), grid_name); m_diagnostic_output = Field(fid); m_diagnostic_output.allocate_view(); } diff --git a/components/eamxx/src/diagnostics/number_path.hpp b/components/eamxx/src/diagnostics/number_path.hpp index 4888d3601f44..30b383b9452f 100644 --- a/components/eamxx/src/diagnostics/number_path.hpp +++ b/components/eamxx/src/diagnostics/number_path.hpp @@ -16,7 +16,7 @@ class NumberPathDiagnostic : public AtmosphereDiagnostic { const ekat::ParameterList ¶ms); // The name of the diagnostic - std::string name() const; + std::string name() const override { return "NumberPath"; } // Set the grid void set_grids(const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/potential_temperature.cpp b/components/eamxx/src/diagnostics/potential_temperature.cpp index 67260e647f6f..8cc45078de76 100644 --- a/components/eamxx/src/diagnostics/potential_temperature.cpp +++ b/components/eamxx/src/diagnostics/potential_temperature.cpp @@ -24,11 +24,6 @@ PotentialTemperatureDiagnostic::PotentialTemperatureDiagnostic (const ekat::Comm } } -std::string PotentialTemperatureDiagnostic::name() const -{ - return m_ptype; -} - // ========================================================================================= void PotentialTemperatureDiagnostic::set_grids(const std::shared_ptr grids_manager) { @@ -51,7 +46,7 @@ void PotentialTemperatureDiagnostic::set_grids(const std::shared_ptr("qc", scalar3d_layout_mid, kg/kg, grid_name, ps); // Construct and allocate the diagnostic field - FieldIdentifier fid (name(), scalar3d_layout_mid, K, grid_name); + FieldIdentifier fid (m_ptype, scalar3d_layout_mid, K, grid_name); m_diagnostic_output = Field(fid); auto& C_ap = m_diagnostic_output.get_header().get_alloc_properties(); C_ap.request_allocation(ps); diff --git a/components/eamxx/src/diagnostics/potential_temperature.hpp b/components/eamxx/src/diagnostics/potential_temperature.hpp index 37fd1a30806c..0ac1a1201d8f 100644 --- a/components/eamxx/src/diagnostics/potential_temperature.hpp +++ b/components/eamxx/src/diagnostics/potential_temperature.hpp @@ -22,7 +22,7 @@ class PotentialTemperatureDiagnostic : public AtmosphereDiagnostic PotentialTemperatureDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const; + std::string name () const override { return "PotentialTemperature"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/precip_surf_mass_flux.cpp b/components/eamxx/src/diagnostics/precip_surf_mass_flux.cpp index 068456f0522d..329e83b6e9cc 100644 --- a/components/eamxx/src/diagnostics/precip_surf_mass_flux.cpp +++ b/components/eamxx/src/diagnostics/precip_surf_mass_flux.cpp @@ -47,7 +47,7 @@ set_grids(const std::shared_ptr grids_manager) } // Construct and allocate the diagnostic field - FieldIdentifier fid(name(), scalar2d_layout_mid, m/s, grid_name); + FieldIdentifier fid(m_name, scalar2d_layout_mid, m/s, grid_name); m_diagnostic_output = Field(fid); m_diagnostic_output.get_header().get_alloc_properties().request_allocation(); m_diagnostic_output.allocate_view(); diff --git a/components/eamxx/src/diagnostics/precip_surf_mass_flux.hpp b/components/eamxx/src/diagnostics/precip_surf_mass_flux.hpp index 68dd1251646e..6ff3458398da 100644 --- a/components/eamxx/src/diagnostics/precip_surf_mass_flux.hpp +++ b/components/eamxx/src/diagnostics/precip_surf_mass_flux.hpp @@ -17,7 +17,7 @@ class PrecipSurfMassFlux : public AtmosphereDiagnostic PrecipSurfMassFlux (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const { return m_name; } + std::string name () const { return "PrecipSurfMassFlux"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/register_diagnostics.hpp b/components/eamxx/src/diagnostics/register_diagnostics.hpp index 45a04eebeb1f..b4830c39e92c 100644 --- a/components/eamxx/src/diagnostics/register_diagnostics.hpp +++ b/components/eamxx/src/diagnostics/register_diagnostics.hpp @@ -36,13 +36,6 @@ inline void register_diagnostics () { diag_factory.register_product("AtmosphereDensity",&create_atmosphere_diagnostic); diag_factory.register_product("Exner",&create_atmosphere_diagnostic); diag_factory.register_product("VirtualTemperature",&create_atmosphere_diagnostic); - diag_factory.register_product("z_int",&create_atmosphere_diagnostic); - diag_factory.register_product("z_mid",&create_atmosphere_diagnostic); - diag_factory.register_product("geopotential_int",&create_atmosphere_diagnostic); - diag_factory.register_product("geopotential_mid",&create_atmosphere_diagnostic); - diag_factory.register_product("height_int",&create_atmosphere_diagnostic); - diag_factory.register_product("height_mid",&create_atmosphere_diagnostic); - diag_factory.register_product("dz",&create_atmosphere_diagnostic); diag_factory.register_product("DryStaticEnergy",&create_atmosphere_diagnostic); diag_factory.register_product("SeaLevelPressure",&create_atmosphere_diagnostic); diag_factory.register_product("WaterPath",&create_atmosphere_diagnostic); @@ -50,6 +43,7 @@ inline void register_diagnostics () { diag_factory.register_product("LongwaveCloudForcing",&create_atmosphere_diagnostic); diag_factory.register_product("RelativeHumidity",&create_atmosphere_diagnostic); diag_factory.register_product("VaporFlux",&create_atmosphere_diagnostic); + diag_factory.register_product("VerticalLayer",&create_atmosphere_diagnostic); diag_factory.register_product("precip_surf_mass_flux",&create_atmosphere_diagnostic); diag_factory.register_product("surface_upward_latent_heat_flux",&create_atmosphere_diagnostic); diag_factory.register_product("wind_speed",&create_atmosphere_diagnostic); @@ -60,4 +54,5 @@ inline void register_diagnostics () { } } // namespace scream + #endif // SCREAM_REGISTER_DIAGNOSTICS_HPP diff --git a/components/eamxx/src/diagnostics/tests/CMakeLists.txt b/components/eamxx/src/diagnostics/tests/CMakeLists.txt index a684aa248fc0..5b318a049224 100644 --- a/components/eamxx/src/diagnostics/tests/CMakeLists.txt +++ b/components/eamxx/src/diagnostics/tests/CMakeLists.txt @@ -1,4 +1,4 @@ -# NOTE: tests inside this if statement won't be built in a baselines-only build +include(ScreamUtils) function (createDiagTest test_name test_srcs) CreateUnitTest(${test_name} "${test_srcs}" @@ -6,72 +6,68 @@ function (createDiagTest test_name test_srcs) LABELS diagnostics) endfunction () -if (NOT SCREAM_ONLY_GENERATE_BASELINES) - include(ScreamUtils) +# Test extracting a single level of a field +CreateDiagTest(field_at_level "field_at_level_tests.cpp") - # Test extracting a single level of a field - CreateDiagTest(field_at_level "field_at_level_tests.cpp") +# Test interpolating a field onto a single pressure level +CreateDiagTest(field_at_pressure_level "field_at_pressure_level_tests.cpp") - # Test interpolating a field onto a single pressure level - CreateDiagTest(field_at_pressure_level "field_at_pressure_level_tests.cpp") - # Test interpolating a field at a specific height - CreateDiagTest(field_at_height "field_at_height_tests.cpp") +# Test interpolating a field at a specific height +CreateDiagTest(field_at_height "field_at_height_tests.cpp") - # Test potential temperature diagnostic - CreateDiagTest(potential_temperature "potential_temperature_test.cpp") +# Test potential temperature diagnostic +CreateDiagTest(potential_temperature "potential_temperature_test.cpp") - # Test exner diagnostic - CreateDiagTest(exner_function "exner_test.cpp") +# Test exner diagnostic +CreateDiagTest(exner_function "exner_test.cpp") - # Test virtual temperature - CreateDiagTest(virtual_temperature "virtual_temperature_test.cpp") +# Test virtual temperature +CreateDiagTest(virtual_temperature "virtual_temperature_test.cpp") - # Test atmosphere density - CreateDiagTest(atmosphere_density "atm_density_test.cpp") +# Test atmosphere density +CreateDiagTest(atmosphere_density "atm_density_test.cpp") - # Test vertical layer (dz, z_int, z_mid) - CreateDiagTest(vertical_layer "vertical_layer_tests.cpp") +# Test vertical layer (dz, z_int, z_mid) +CreateDiagTest(vertical_layer "vertical_layer_tests.cpp") - # Test dry static energy - CreateDiagTest(dry_static_energy "dry_static_energy_test.cpp") +# Test dry static energy +CreateDiagTest(dry_static_energy "dry_static_energy_test.cpp") - # Test sea level pressure - CreateDiagTest(sea_level_pressure "sea_level_pressure_test.cpp") +# Test sea level pressure +CreateDiagTest(sea_level_pressure "sea_level_pressure_test.cpp") - # Test total water path - CreateDiagTest(water_path "water_path_tests.cpp") +# Test total water path +CreateDiagTest(water_path "water_path_tests.cpp") - # Test shortwave cloud forcing - CreateDiagTest(shortwave_cloud_forcing "shortwave_cloud_forcing_tests.cpp") +# Test shortwave cloud forcing +CreateDiagTest(shortwave_cloud_forcing "shortwave_cloud_forcing_tests.cpp") - # Test longwave cloud forcing - CreateDiagTest(longwave_cloud_forcing "longwave_cloud_forcing_tests.cpp") +# Test longwave cloud forcing +CreateDiagTest(longwave_cloud_forcing "longwave_cloud_forcing_tests.cpp") - # Test Relative Humidity - CreateDiagTest(relative_humidity "relative_humidity_tests.cpp") +# Test Relative Humidity +CreateDiagTest(relative_humidity "relative_humidity_tests.cpp") - # Test Vapor Flux - CreateDiagTest(vapor_flux "vapor_flux_tests.cpp") +# Test Vapor Flux +CreateDiagTest(vapor_flux "vapor_flux_tests.cpp") - # Test precipitation mass surface flux - CreateDiagTest(precip_surf_mass_flux "precip_surf_mass_flux_tests.cpp") +# Test precipitation mass surface flux +CreateDiagTest(precip_surf_mass_flux "precip_surf_mass_flux_tests.cpp") - # Test surface latent heat flux - CreateDiagTest(surface_upward_latent_heat_flux "surf_upward_latent_heat_flux_tests.cpp") +# Test surface latent heat flux +CreateDiagTest(surface_upward_latent_heat_flux "surf_upward_latent_heat_flux_tests.cpp") - # Test wind speed diagnostic - CreateDiagTest(wind_speed "wind_speed_tests.cpp") +# Test wind speed diagnostic +CreateDiagTest(wind_speed "wind_speed_tests.cpp") - # Test AODVIS - CreateDiagTest(aodvis "aodvis_test.cpp") +# Test AODVIS +CreateDiagTest(aodvis "aodvis_test.cpp") - # Test "number" paths - CreateDiagTest(number_paths "number_paths_tests.cpp") +# Test "number" paths +CreateDiagTest(number_paths "number_paths_tests.cpp") - # Test AEROCOM_CLD - CreateDiagTest(aerocom_cld "aerocom_cld_test.cpp") +# Test AEROCOM_CLD +CreateDiagTest(aerocom_cld "aerocom_cld_test.cpp") - # Test atm_tend - CreateDiagTest(atm_backtend "atm_backtend_test.cpp") - -endif() +# Test atm_tend +CreateDiagTest(atm_backtend "atm_backtend_test.cpp") diff --git a/components/eamxx/src/diagnostics/tests/field_at_height_tests.cpp b/components/eamxx/src/diagnostics/tests/field_at_height_tests.cpp index e26d34dd1a55..6e071e4b2013 100644 --- a/components/eamxx/src/diagnostics/tests/field_at_height_tests.cpp +++ b/components/eamxx/src/diagnostics/tests/field_at_height_tests.cpp @@ -104,12 +104,13 @@ TEST_CASE("field_at_height") // Lambda to create and run a diag, and return output auto run_diag = [&](const Field& f, const Field& z, - const std::string& loc, const std::string& surf_ref) { + const double h, const std::string& surf_ref) { util::TimeStamp t0 ({2022,1,1},{0,0,0}); auto& factory = AtmosphereDiagnosticFactory::instance(); ekat::ParameterList pl; pl.set("surface_reference",surf_ref); - pl.set("vertical_location",loc); + pl.set("height_value",std::to_string(h)); + pl.set("height_units",std::string("m")); pl.set("field_name",f.name()); pl.set("grid_name",grid->name()); auto diag = factory.create("FieldAtheight",comm,pl); @@ -173,13 +174,12 @@ TEST_CASE("field_at_height") // Make sure that an unsupported reference height throws an error. print(" -> Testing throws error with unsupported reference height...\n"); { - REQUIRE_THROWS(run_diag (s_mid,h_mid,"1m","foobar")); + REQUIRE_THROWS(run_diag (s_mid,h_mid,1.0,"foobar")); } print(" -> Testing throws error with unsupported reference height... OK\n"); // Run many times int z_tgt; - std::string loc; for (std::string surf_ref : {"sealevel","surface"}) { printf(" -> Testing for a reference height above %s...\n",surf_ref.c_str()); const auto mid_src = surf_ref == "sealevel" ? z_mid : h_mid; @@ -197,32 +197,31 @@ TEST_CASE("field_at_height") // Set target z-slice for testing to a random value. z_tgt = pdf_levs(engine)+max_surf_4test; - loc = std::to_string(z_tgt) + "m"; - printf(" -> test at height of %s.............\n",loc.c_str()); + printf(" -> test at height of %dm............\n",z_tgt); { print(" -> scalar midpoint field...............\n"); - auto d = run_diag(s_mid,mid_src,loc,surf_ref); + auto d = run_diag(s_mid,mid_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,mid_src,s_tgt); REQUIRE (views_are_approx_equal(d,s_tgt,tol)); print(" -> scalar midpoint field............... OK!\n"); } { print(" -> scalar interface field...............\n"); - auto d = run_diag (s_int,int_src,loc,surf_ref); + auto d = run_diag (s_int,int_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,int_src,s_tgt); REQUIRE (views_are_approx_equal(d,s_tgt,tol)); print(" -> scalar interface field............... OK!\n"); } { print(" -> vector midpoint field...............\n"); - auto d = run_diag (v_mid,mid_src,loc,surf_ref); + auto d = run_diag (v_mid,mid_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,mid_src,v_tgt); REQUIRE (views_are_approx_equal(d,v_tgt,tol)); print(" -> vector midpoint field............... OK!\n"); } { print(" -> vector interface field...............\n"); - auto d = run_diag (v_int,int_src,loc,surf_ref); + auto d = run_diag (v_int,int_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,int_src,v_tgt); REQUIRE (views_are_approx_equal(d,v_tgt,tol)); print(" -> vector interface field............... OK!\n"); @@ -230,8 +229,7 @@ TEST_CASE("field_at_height") { print(" -> Forced fail, give incorrect location...............\n"); const int z_tgt_adj = (z_tgt+max_surf_4test)/2; - std::string loc_err = std::to_string(z_tgt_adj) + "m"; - auto d = run_diag(s_int,int_src,loc_err,surf_ref); + auto d = run_diag(s_int,int_src,z_tgt_adj,surf_ref); f_z_tgt(inter,slope,z_tgt,int_src,s_tgt); REQUIRE (!views_are_approx_equal(d,s_tgt,tol,false)); print(" -> Forced fail, give incorrect location............... OK!\n"); @@ -243,15 +241,13 @@ TEST_CASE("field_at_height") auto inter = pdf_y0(engine); f_z_src(inter, slope, int_src, s_int); z_tgt = 2*z_top; - std::string loc = std::to_string(z_tgt) + "m"; - auto dtop = run_diag(s_int,int_src,loc,surf_ref); + auto dtop = run_diag(s_int,int_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,int_src,s_tgt); REQUIRE (views_are_approx_equal(dtop,s_tgt,tol)); print(" -> Forced extrapolation at top............... OK!\n"); print(" -> Forced extrapolation at bot...............\n"); z_tgt = 0; - loc = std::to_string(z_tgt) + "m"; - auto dbot = run_diag(s_int,int_src,loc,surf_ref); + auto dbot = run_diag(s_int,int_src,z_tgt,surf_ref); f_z_tgt(inter,slope,z_tgt,int_src,s_tgt); REQUIRE (views_are_approx_equal(dbot,s_tgt,tol)); print(" -> Forced extrapolation at bot............... OK!\n"); diff --git a/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp b/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp index 4e0deab1dfca..ba733980cade 100644 --- a/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp +++ b/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp @@ -224,7 +224,8 @@ get_test_diag(const ekat::Comm& comm, std::shared_ptr fm, st ekat::ParameterList params; params.set("field_name",field.name()); params.set("grid_name",fm->get_grid()->name()); - params.set("vertical_location",std::to_string(plevel) + "Pa"); + params.set("pressure_value",std::to_string(plevel)); + params.set("pressure_units",std::string("Pa")); auto diag = std::make_shared(comm,params); diag->set_grids(gm); for (const auto& req : diag->get_required_field_requests()) { diff --git a/components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp b/components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp index fe75114611d5..631cb5acd8a8 100644 --- a/components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp +++ b/components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp @@ -50,14 +50,10 @@ void run (const std::string& diag_name, const std::string& location) // Construct the Diagnostic ekat::ParameterList params; - std::string name = diag_name; - if (location=="midpoints") { - name += "_mid"; - } else if (location=="interfaces") { - name += "_int"; - } - params.set("diag_name", name); - auto diag = diag_factory.create(name,comm,params); + + params.set("diag_name", diag_name); + params.set("vert_location",location); + auto diag = diag_factory.create("VerticalLayer",comm,params); diag->set_grids(gm); const bool needs_phis = diag_name=="z" or diag_name=="geopotential"; @@ -180,7 +176,7 @@ TEST_CASE("vertical_layer_test", "vertical_layer_test]"){ std::string msg = " -> Testing diag=dz "; std::string dots (50-msg.size(),'.'); root_print (msg + dots + "\n"); - run("dz", "UNUSED"); + run("dz", "midpoints"); root_print (msg + dots + " PASS!\n"); }; diff --git a/components/eamxx/src/diagnostics/vapor_flux.cpp b/components/eamxx/src/diagnostics/vapor_flux.cpp index ce88c1443015..b6577bf1205a 100644 --- a/components/eamxx/src/diagnostics/vapor_flux.cpp +++ b/components/eamxx/src/diagnostics/vapor_flux.cpp @@ -26,11 +26,6 @@ VaporFluxDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params) } } -std::string VaporFluxDiagnostic::name() const -{ - return m_component==0 ? "ZonalVapFlux" : "MeridionalVapFlux"; -} - void VaporFluxDiagnostic::set_grids(const std::shared_ptr grids_manager) { using namespace ekat::units; @@ -51,7 +46,8 @@ void VaporFluxDiagnostic::set_grids(const std::shared_ptr gr add_field("horiz_winds", vector3d, m/s, grid_name); // Construct and allocate the diagnostic field - FieldIdentifier fid (name(), scalar2d, kg/m/s, grid_name); + std::string dname = m_component==0 ? "ZonalVapFlux" : "MeridionalVapFlux"; + FieldIdentifier fid (dname, scalar2d, kg/m/s, grid_name); m_diagnostic_output = Field(fid); m_diagnostic_output.allocate_view(); } diff --git a/components/eamxx/src/diagnostics/vapor_flux.hpp b/components/eamxx/src/diagnostics/vapor_flux.hpp index 3d82fd882f8c..5bacd78a9b71 100644 --- a/components/eamxx/src/diagnostics/vapor_flux.hpp +++ b/components/eamxx/src/diagnostics/vapor_flux.hpp @@ -17,7 +17,7 @@ class VaporFluxDiagnostic : public AtmosphereDiagnostic VaporFluxDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const; + std::string name () const override { return "VaporFlux"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/vertical_layer.cpp b/components/eamxx/src/diagnostics/vertical_layer.cpp index 32d870da03cc..f9ff526dfcd8 100644 --- a/components/eamxx/src/diagnostics/vertical_layer.cpp +++ b/components/eamxx/src/diagnostics/vertical_layer.cpp @@ -13,25 +13,27 @@ VerticalLayerDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& para : AtmosphereDiagnostic(comm,params) { m_diag_name = params.get("diag_name"); - std::vector supported = { - "z_int", - "z_mid", - "geopotential_int", - "geopotential_mid", - "height_int", - "height_mid", - "dz" - }; + std::vector supported = {"z","geopotential","height","dz"}; EKAT_REQUIRE_MSG(ekat::contains(supported,m_diag_name), "[VerticalLayerDiagnostic] Error! Invalid diag_name.\n" " - diag_name : " + m_diag_name + "\n" " - valid names: " + ekat::join(supported,", ") + "\n"); - m_is_interface_layout = m_diag_name.find("_int") != std::string::npos; + auto vert_pos = params.get("vert_location"); + EKAT_REQUIRE_MSG (vert_pos=="mid" || vert_pos=="int" || + vert_pos=="midpoints" || vert_pos=="interfaces", + "[VerticalLayerDiagnostic] Error! Invalid 'vert_location'.\n" + " - input value: " + vert_pos + "\n" + " - valid names: mid, midpoints, int, interfaces\n"); + m_is_interface_layout = vert_pos=="int" || vert_pos=="interfaces"; + + m_geopotential = m_diag_name=="geopotential"; + m_from_sea_level = m_diag_name=="z" or m_geopotential; - m_geopotential = m_diag_name.substr(0,12)=="geopotential"; - m_from_sea_level = m_diag_name[0]=='z' or m_geopotential; + if (m_diag_name!="dz") { + m_diag_name += m_is_interface_layout ? "_int" : "_mid"; + } } // ======================================================================================== void VerticalLayerDiagnostic:: @@ -88,7 +90,7 @@ initialize_impl (const RunType /*run_type*/) const auto VLEV = m_is_interface_layout ? ILEV : LEV; const auto nlevs = m_is_interface_layout ? m_num_levs+1 : m_num_levs; FieldLayout diag_layout ({COL,VLEV},{m_num_cols,nlevs}); - FieldIdentifier fid (name(), diag_layout, m_geopotential ? m2/s2 : m, grid_name); + FieldIdentifier fid (m_diag_name, diag_layout, m_geopotential ? m2/s2 : m, grid_name); m_diagnostic_output = Field(fid); auto& diag_fap = m_diagnostic_output.get_header().get_alloc_properties(); diff --git a/components/eamxx/src/diagnostics/vertical_layer.hpp b/components/eamxx/src/diagnostics/vertical_layer.hpp index 805fd70028f6..a440bd6a8ee9 100644 --- a/components/eamxx/src/diagnostics/vertical_layer.hpp +++ b/components/eamxx/src/diagnostics/vertical_layer.hpp @@ -24,7 +24,7 @@ class VerticalLayerDiagnostic : public AtmosphereDiagnostic VerticalLayerDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic. - std::string name () const { return m_diag_name; } + std::string name () const { return "VerticalLayer"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/diagnostics/water_path.cpp b/components/eamxx/src/diagnostics/water_path.cpp index 15fa5cdef38c..bca771e6dab6 100644 --- a/components/eamxx/src/diagnostics/water_path.cpp +++ b/components/eamxx/src/diagnostics/water_path.cpp @@ -32,11 +32,6 @@ WaterPathDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params) } } -std::string WaterPathDiagnostic::name() const -{ - return m_kind + "WaterPath"; -} - void WaterPathDiagnostic:: set_grids(const std::shared_ptr grids_manager) { @@ -57,7 +52,7 @@ set_grids(const std::shared_ptr grids_manager) add_field(m_qname, scalar3d, kg/kg, grid_name); // Construct and allocate the diagnostic field - FieldIdentifier fid (name(), scalar2d, kg/m2, grid_name); + FieldIdentifier fid (m_kind + "WaterPath", scalar2d, kg/m2, grid_name); m_diagnostic_output = Field(fid); m_diagnostic_output.allocate_view(); } diff --git a/components/eamxx/src/diagnostics/water_path.hpp b/components/eamxx/src/diagnostics/water_path.hpp index 0b9515e0b421..722a51a31331 100644 --- a/components/eamxx/src/diagnostics/water_path.hpp +++ b/components/eamxx/src/diagnostics/water_path.hpp @@ -17,7 +17,7 @@ class WaterPathDiagnostic : public AtmosphereDiagnostic WaterPathDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params); // The name of the diagnostic - std::string name () const; + std::string name () const override { return "WaterPath"; } // Set the grid void set_grids (const std::shared_ptr grids_manager); diff --git a/components/eamxx/src/share/io/CMakeLists.txt b/components/eamxx/src/share/io/CMakeLists.txt index 6c24d3bc51da..4908e4ae0c95 100644 --- a/components/eamxx/src/share/io/CMakeLists.txt +++ b/components/eamxx/src/share/io/CMakeLists.txt @@ -59,7 +59,7 @@ add_library(scream_io scream_io_utils.cpp ) -target_link_libraries(scream_io PUBLIC scream_share scream_scorpio_interface) +target_link_libraries(scream_io PUBLIC scream_share scream_scorpio_interface diagnostics) if (NOT SCREAM_LIB_ONLY) add_subdirectory(tests) diff --git a/components/eamxx/src/share/io/scorpio_output.cpp b/components/eamxx/src/share/io/scorpio_output.cpp index 47a5616aec55..4e9d5add8004 100644 --- a/components/eamxx/src/share/io/scorpio_output.cpp +++ b/components/eamxx/src/share/io/scorpio_output.cpp @@ -6,6 +6,8 @@ #include "share/util/scream_timing.hpp" #include "share/field/field_utils.hpp" +#include "diagnostics/register_diagnostics.hpp" + #include "ekat/util/ekat_units.hpp" #include "ekat/util/ekat_string_utils.hpp" #include "ekat/std_meta/ekat_std_utils.hpp" @@ -1334,153 +1336,49 @@ void AtmosphereOutput::set_diagnostics() /* ---------------------------------------------------------- */ std::shared_ptr -AtmosphereOutput::create_diagnostic (const std::string& diag_field_name) { - auto& diag_factory = AtmosphereDiagnosticFactory::instance(); +AtmosphereOutput::create_diagnostic (const std::string& diag_field_name) +{ + // We need scream scope resolution, since this->create_diagnostic is hiding it + auto diag = scream::create_diagnostic(diag_field_name,get_field_manager("sim")->get_grid()); - // Construct a diagnostic by this name - ekat::ParameterList params; - std::string diag_name; + // Some diags need some extra setup or trigger extra behaviors std::string diag_avg_cnt_name = ""; - - if (diag_field_name.find("_at_")!=std::string::npos) { - // The diagnostic must be one of - // - ${field_name}_at_lev_${N} <- interface fields still use "_lev_" - // - ${field_name}_at_model_bot - // - ${field_name}_at_model_top - // - ${field_name}_at_${M}X - // where M/N are numbers (N integer), X=Pa, hPa, mb, or m - auto tokens = ekat::split(diag_field_name,"_at_"); - EKAT_REQUIRE_MSG (tokens.size()==2, - "Error! Unexpected diagnostic name: " + diag_field_name + "\n"); - - const auto& fname = tokens.front(); - params.set("field_name",fname); - params.set("grid_name",get_field_manager("sim")->get_grid()->name()); - - params.set("vertical_location", tokens[1]); + auto& params = diag->get_params(); + if (diag->name()=="FieldAtPressureLevel") { params.set("mask_value",m_fill_value); - - // Conventions on notation (N=any integer): - // FieldAtLevel : var_at_lev_N, var_at_model_top, var_at_model_bot - // FieldAtPressureLevel: var_at_Nx, with x=mb,Pa,hPa - // FieldAtHeight : var_at_Nm_above_Y (Y=sealevel or surface) - if (tokens[1].find_first_of("0123456789.")==0) { - auto units_start = tokens[1].find_first_not_of("0123456789."); - auto units = tokens[1].substr(units_start); - if (units.find("_above_") != std::string::npos) { - // The field is at a height above a specific reference. - // Currently we only support FieldAtHeight above "sealevel" or "surface" - auto subtokens = ekat::split(units,"_above_"); - params.set("surface_reference",subtokens[1]); - units = subtokens[0]; - // Need to reset the vertical location to strip the "_above_" part of the string. - params.set("vertical_location", tokens[1].substr(0,units_start)+subtokens[0]); - // If the slice is "above_sealevel" then we need to track the avg cnt uniquely. - // Note, "above_surface" is expected to never have masking and can thus use - // the typical 2d layout avg cnt. - if (subtokens[1]=="sealevel") { - diag_avg_cnt_name = "_" + tokens[1]; // Set avg_cnt tracking for this specific slice - // If we have 2D slices we need to be tracking the average count, - // if m_avg_type is not Instant - m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; - } - } - if (units=="m") { - diag_name = "FieldAtHeight"; - EKAT_REQUIRE_MSG(params.isParameter("surface_reference"),"Error! Output field request for " + diag_field_name + " is missing a surface reference." - " Please add either '_above_sealevel' or '_above_surface' to the field name"); - } else if (units=="mb" or units=="Pa" or units=="hPa") { - diag_name = "FieldAtPressureLevel"; - diag_avg_cnt_name = "_" + tokens[1]; // Set avg_cnt tracking for this specific slice - // If we have 2D slices we need to be tracking the average count, - // if m_avg_type is not Instant - m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; - } else { - EKAT_ERROR_MSG ("Error! Invalid units x for 'field_at_Nx' diagnostic.\n"); - } - } else { - diag_name = "FieldAtLevel"; + diag_avg_cnt_name = "_" + + params.get("pressure_value") + + params.get("pressure_units"); + m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; + } else if (diag->name()=="FieldAtHeight") { + if (params.get("surface_reference")=="sealevel") { + diag_avg_cnt_name = "_" + + params.get("height_value") + + params.get("height_units") + "_above_sealevel"; + m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; } - } else if (diag_field_name=="precip_liq_surf_mass_flux" or - diag_field_name=="precip_ice_surf_mass_flux" or - diag_field_name=="precip_total_surf_mass_flux") { - diag_name = "precip_surf_mass_flux"; - // split will return [X, ''], with X being whatever is before '_surf_mass_flux' - auto type = ekat::split(diag_field_name.substr(7),"_surf_mass_flux").front(); - params.set("precip_type",type); - } else if (diag_field_name=="IceWaterPath" or - diag_field_name=="LiqWaterPath" or - diag_field_name=="RainWaterPath" or - diag_field_name=="RimeWaterPath" or - diag_field_name=="VapWaterPath") { - diag_name = "WaterPath"; - // split will return the list [X, ''], with X being whatever is before 'WaterPath' - params.set("Water Kind",ekat::split(diag_field_name,"WaterPath").front()); - } else if (diag_field_name=="IceNumberPath" or - diag_field_name=="LiqNumberPath" or - diag_field_name=="RainNumberPath") { - diag_name = "NumberPath"; - // split will return the list [X, ''], with X being whatever is before 'NumberPath' - params.set("Number Kind",ekat::split(diag_field_name,"NumberPath").front()); - } else if (diag_field_name=="AeroComCldTop" or - diag_field_name=="AeroComCldBot") { - diag_name = "AeroComCld"; - // split will return the list ['', X], with X being whatever is after 'AeroComCld' - params.set("AeroComCld Kind",ekat::split(diag_field_name,"AeroComCld").back()); - } else if (diag_field_name=="MeridionalVapFlux" or - diag_field_name=="ZonalVapFlux") { - diag_name = "VaporFlux"; - // split will return the list [X, ''], with X being whatever is before 'VapFlux' - params.set("Wind Component",ekat::split(diag_field_name,"VapFlux").front()); - } else if (diag_field_name.find("_atm_backtend")!=std::string::npos) { - diag_name = "AtmBackTendDiag"; - // Set the grid_name - params.set("grid_name",get_field_manager("sim")->get_grid()->name()); - // split will return [X, ''], with X being whatever is before '_atm_tend' - params.set("Tendency Name",ekat::split(diag_field_name,"_atm_backtend").front()); - } else if (diag_field_name=="PotentialTemperature" or - diag_field_name=="LiqPotentialTemperature") { - diag_name = "PotentialTemperature"; - if (diag_field_name == "LiqPotentialTemperature") { - params.set("Temperature Kind", "Liq"); - } else { - params.set("Temperature Kind", "Tot"); - } - } else { - diag_name = diag_field_name; - } - - // These fields are special case of VerticalLayer diagnostic. - // The diagnostics requires the name to be given as param value. - if (diag_name == "z_int" or diag_name == "z_mid" or - diag_name == "geopotential_int" or diag_name == "geopotential_mid" or - diag_name == "height_int" or diag_name == "height_mid" or - diag_name == "dz") { - params.set("diag_name", diag_name); } - // Create the diagnostic - auto diag = diag_factory.create(diag_name,m_comm,params); - diag->set_grids(m_grids_manager); - // Ensure there's an entry in the map for this diag, so .at(diag_name) always works - auto& deps = m_diag_depends_on_diags[diag->name()]; + auto& deps = m_diag_depends_on_diags[diag_field_name]; // Initialize the diagnostic const auto sim_field_mgr = get_field_manager("sim"); for (const auto& freq : diag->get_required_field_requests()) { const auto& fname = freq.fid.name(); if (!sim_field_mgr->has_field(fname)) { + std::cout << diag_field_name << " depends on the diag " << fname << "\n"; // This diag depends on another diag. Create and init the dependency if (m_diagnostics.count(fname)==0) { m_diagnostics[fname] = create_diagnostic(fname); } - auto dep = m_diagnostics.at(fname); deps.push_back(fname); } diag->set_required_field (get_field(fname,"sim")); } + diag->initialize(util::TimeStamp(),RunType::Initial); + // If specified, set avg_cnt tracking for this diagnostic. if (m_track_avg_cnt) { const auto diag_field = diag->get_diagnostic(); diff --git a/components/eamxx/src/share/io/scream_io_utils.cpp b/components/eamxx/src/share/io/scream_io_utils.cpp index 9318728d657a..40b8a97de4b1 100644 --- a/components/eamxx/src/share/io/scream_io_utils.cpp +++ b/components/eamxx/src/share/io/scream_io_utils.cpp @@ -1,6 +1,7 @@ #include "share/io/scream_io_utils.hpp" #include "share/io/scream_scorpio_interface.hpp" +#include "share/grid/library_grids_manager.hpp" #include "share/util/scream_utils.hpp" #include "share/scream_config.hpp" @@ -117,4 +118,90 @@ util::TimeStamp read_timestamp (const std::string& filename, return ts; } +std::shared_ptr +create_diagnostic (const std::string& diag_field_name, + const std::shared_ptr& grid) +{ + // Note: use grouping (the (..) syntax), so you can later query the content + // of each group in the matches output var! + // Note: use raw string syntax R"()" to avoid having to escape the \ character + // Note: the number for field_at_p/h can match positive integer/floating-point numbers + std::regex field_at_l (R"(([A-Za-z0-9_]+)_at_(lev_(\d+)|model_(top|bot))$)"); + std::regex field_at_p (R"(([A-Za-z0-9_]+)_at_(\d+(\.\d+)?)(hPa|mb|Pa)$)"); + std::regex field_at_h (R"(([A-Za-z0-9_]+)_at_(\d+(\.\d+)?)(m)_above_(sealevel|surface)$)"); + std::regex surf_mass_flux ("precip_(liq|ice|total)_surf_mass_flux$"); + std::regex water_path ("(Ice|Liq|Rain|Rime|Vap)WaterPath$"); + std::regex number_path ("(Ice|Liq|Rain)NumberPath$"); + std::regex aerocom_cld ("AeroComCld(Top|Bot)$"); + std::regex vap_flux ("(Meridional|Zonal)VapFlux$"); + std::regex backtend ("([A-Za-z0-9_]+)_atm_backtend$"); + std::regex pot_temp ("(Liq)?PotentialTemperature$"); + std::regex vert_layer ("(z|geopotential|height)_(mid|int)$"); + + std::string diag_name; + std::smatch matches; + ekat::ParameterList params(diag_field_name); + + if (std::regex_search(diag_field_name,matches,field_at_l)) { + params.set("field_name",matches[1].str()); + params.set("grid_name",grid->name()); + params.set("vertical_location", matches[2].str()); + diag_name = "FieldAtLevel"; + } else if (std::regex_search(diag_field_name,matches,field_at_p)) { + params.set("field_name",matches[1].str()); + params.set("grid_name",grid->name()); + params.set("pressure_value",matches[2].str()); + params.set("pressure_units", matches[4].str()); + diag_name = "FieldAtPressureLevel"; + } else if (std::regex_search(diag_field_name,matches,field_at_h)) { + params.set("field_name",matches[1].str()); + params.set("grid_name",grid->name()); + params.set("height_value",matches[2].str()); + params.set("height_units",matches[4].str()); + params.set("surface_reference", matches[5].str()); + diag_name = "FieldAtHeight"; + } else if (std::regex_search(diag_field_name,matches,surf_mass_flux)) { + diag_name = "precip_surf_mass_flux"; + params.set("precip_type",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,water_path)) { + diag_name = "WaterPath"; + params.set("Water Kind",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,number_path)) { + diag_name = "NumberPath"; + params.set("Number Kind",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,aerocom_cld)) { + diag_name = "AeroComCld"; + params.set("AeroComCld Kind",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,vap_flux)) { + diag_name = "VaporFlux"; + params.set("Wind Component",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,backtend)) { + diag_name = "AtmBackTendDiag"; + // Set the grid_name + params.set("grid_name",grid->name()); + params.set("Tendency Name",matches[1].str()); + } else if (std::regex_search(diag_field_name,matches,pot_temp)) { + diag_name = "PotentialTemperature"; + params.set("Temperature Kind", matches[1].str()!="" ? matches[1].str() : std::string("Tot")); + } else if (std::regex_search(diag_field_name,matches,vert_layer)) { + diag_name = "VerticalLayer"; + params.set("diag_name",matches[1].str()); + params.set("vert_location",matches[2].str()); + } else if (diag_field_name=="dz") { + diag_name = "VerticalLayer"; + params.set("diag_name","dz"); + params.set("vert_location","mid"); + } else { + // No existing special regex matches, so we assume that the diag field name IS the diag name. + diag_name = diag_field_name; + } + + auto comm = grid->get_comm(); + auto diag = AtmosphereDiagnosticFactory::instance().create(diag_name,comm,params); + auto gm = std::make_shared(grid); + diag->set_grids(gm); + + return diag; +} + } // namespace scream diff --git a/components/eamxx/src/share/io/scream_io_utils.hpp b/components/eamxx/src/share/io/scream_io_utils.hpp index 9eda82e2bbd3..f725b91b94c8 100644 --- a/components/eamxx/src/share/io/scream_io_utils.hpp +++ b/components/eamxx/src/share/io/scream_io_utils.hpp @@ -3,11 +3,14 @@ #include "scream_io_control.hpp" #include "share/util/scream_time_stamp.hpp" +#include "share/atm_process/atmosphere_diagnostic.hpp" +#include "share/grid/abstract_grid.hpp" #include #include #include +#include namespace scream { @@ -184,5 +187,11 @@ util::TimeStamp read_timestamp (const std::string& filename, const std::string& ts_name, const bool read_nsteps = false); +// Create a diagnostic from a string representation of it. +// E.g., create the diag to compute fieldX_at_500hPa. +std::shared_ptr +create_diagnostic (const std::string& diag_name, + const std::shared_ptr& grid); + } // namespace scream #endif // SCREAM_IO_UTILS_HPP diff --git a/components/eamxx/src/share/io/tests/CMakeLists.txt b/components/eamxx/src/share/io/tests/CMakeLists.txt index 8c819f7896cc..89d92075723d 100644 --- a/components/eamxx/src/share/io/tests/CMakeLists.txt +++ b/components/eamxx/src/share/io/tests/CMakeLists.txt @@ -17,6 +17,12 @@ CreateUnitTest(io_utils "io_utils.cpp" PROPERTIES RESOURCE_LOCK rpointer_file ) +# Test creation of diagnostic from diag_field_name +CreateUnitTest(create_diag "create_diag.cpp" + LIBS diagnostics scream_io + LABELS io diagnostics +) + ## Test basic output (no packs, no diags, all avg types, all freq units) CreateUnitTest(io_basic "io_basic.cpp" LIBS scream_io LABELS io diff --git a/components/eamxx/src/share/io/tests/create_diag.cpp b/components/eamxx/src/share/io/tests/create_diag.cpp new file mode 100644 index 000000000000..bc5091759719 --- /dev/null +++ b/components/eamxx/src/share/io/tests/create_diag.cpp @@ -0,0 +1,166 @@ +#include "catch2/catch.hpp" + +#include "diagnostics/register_diagnostics.hpp" + +#include "share/io/scream_io_utils.hpp" +#include "share/grid/point_grid.hpp" + +namespace scream { + +TEST_CASE("create_diag") +{ + ekat::Comm comm(MPI_COMM_WORLD); + + register_diagnostics(); + + // Create a grid + const int ncols = 3*comm.size(); + const int nlevs = 10; + auto grid = create_point_grid("Physics",ncols,nlevs,comm); + + SECTION ("field_at") { + // FieldAtLevel + auto d1 = create_diagnostic("BlaH_123_at_model_top",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + auto d2 = create_diagnostic("BlaH_123_at_model_bot",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + auto d3 = create_diagnostic("BlaH_123_at_lev_10",grid); + REQUIRE (std::dynamic_pointer_cast(d3)!=nullptr); + + REQUIRE_THROWS(create_diagnostic("BlaH_123_at_modeltop",grid)); // misspelled + + // FieldAtPressureLevel + auto d4 = create_diagnostic("BlaH_123_at_10mb",grid); + REQUIRE (std::dynamic_pointer_cast(d4)!=nullptr); + auto d5 = create_diagnostic("BlaH_123_at_10hPa",grid); + REQUIRE (std::dynamic_pointer_cast(d5)!=nullptr); + auto d6 = create_diagnostic("BlaH_123_at_10Pa",grid); + REQUIRE (std::dynamic_pointer_cast(d6)!=nullptr); + + REQUIRE_THROWS(create_diagnostic("BlaH_123_at_400KPa",grid)); // invalid units + + // FieldAtHeight + auto d7 = create_diagnostic("BlaH_123_at_10m_above_sealevel",grid); + REQUIRE (std::dynamic_pointer_cast(d7)!=nullptr); + auto d8 = create_diagnostic("BlaH_123_at_10m_above_surface",grid); + REQUIRE (std::dynamic_pointer_cast(d8)!=nullptr); + + REQUIRE_THROWS(create_diagnostic("BlaH_123_at_10.5m",grid)); // missing _above_X + REQUIRE_THROWS(create_diagnostic("BlaH_123_at_1km_above_sealevel",grid)); // invalid units + REQUIRE_THROWS(create_diagnostic("BlaH_123_at_1m_above_the_surface",grid)); // invalid reference + } + + SECTION ("precip_mass_flux") { + auto d1 = create_diagnostic("precip_liq_surf_mass_flux",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("precip_type")=="liq"); + + auto d2 = create_diagnostic("precip_ice_surf_mass_flux",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("precip_type")=="ice"); + + auto d3 = create_diagnostic("precip_total_surf_mass_flux",grid); + REQUIRE (std::dynamic_pointer_cast(d3)!=nullptr); + REQUIRE (d3->get_params().get("precip_type")=="total"); + } + + SECTION ("water_and_number_path") { + auto d1 = create_diagnostic("LiqWaterPath",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("Water Kind")=="Liq"); + + auto d2 = create_diagnostic("IceWaterPath",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("Water Kind")=="Ice"); + + auto d3 = create_diagnostic("RainWaterPath",grid); + REQUIRE (std::dynamic_pointer_cast(d3)!=nullptr); + REQUIRE (d3->get_params().get("Water Kind")=="Rain"); + + auto d4 = create_diagnostic("RimeWaterPath",grid); + REQUIRE (std::dynamic_pointer_cast(d4)!=nullptr); + REQUIRE (d4->get_params().get("Water Kind")=="Rime"); + + auto d5 = create_diagnostic("VapWaterPath",grid); + REQUIRE (std::dynamic_pointer_cast(d5)!=nullptr); + REQUIRE (d5->get_params().get("Water Kind")=="Vap"); + + auto d6 = create_diagnostic("LiqNumberPath",grid); + REQUIRE (std::dynamic_pointer_cast(d6)!=nullptr); + REQUIRE (d6->get_params().get("Number Kind")=="Liq"); + + auto d7 = create_diagnostic("IceNumberPath",grid); + REQUIRE (std::dynamic_pointer_cast(d7)!=nullptr); + REQUIRE (d7->get_params().get("Number Kind")=="Ice"); + + auto d8 = create_diagnostic("RainNumberPath",grid); + REQUIRE (std::dynamic_pointer_cast(d8)!=nullptr); + REQUIRE (d8->get_params().get("Number Kind")=="Rain"); + } + + SECTION ("aerocom_cld") { + auto d1 = create_diagnostic("AeroComCldTop",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("AeroComCld Kind")=="Top"); + + auto d2 = create_diagnostic("AeroComCldBot",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("AeroComCld Kind")=="Bot"); + } + + SECTION ("vapor_flux") { + auto d1 = create_diagnostic("MeridionalVapFlux",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("Wind Component")=="Meridional"); + + auto d2 = create_diagnostic("ZonalVapFlux",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("Wind Component")=="Zonal"); + } + + SECTION ("atm_tend") { + auto d1 = create_diagnostic("BlaH_123_atm_backtend",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("Tendency Name")=="BlaH_123"); + } + + SECTION ("pot_temp") { + auto d1 = create_diagnostic("LiqPotentialTemperature",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("Temperature Kind")=="Liq"); + + auto d2 = create_diagnostic("PotentialTemperature",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("Temperature Kind")=="Tot"); + } + + SECTION ("vert_layer") { + auto d1 = create_diagnostic("z_mid",grid); + REQUIRE (std::dynamic_pointer_cast(d1)!=nullptr); + REQUIRE (d1->get_params().get("vert_location")=="mid"); + auto d2 = create_diagnostic("z_int",grid); + REQUIRE (std::dynamic_pointer_cast(d2)!=nullptr); + REQUIRE (d2->get_params().get("vert_location")=="int"); + + auto d3 = create_diagnostic("height_mid",grid); + REQUIRE (std::dynamic_pointer_cast(d3)!=nullptr); + REQUIRE (d3->get_params().get("vert_location")=="mid"); + auto d4 = create_diagnostic("height_int",grid); + REQUIRE (std::dynamic_pointer_cast(d4)!=nullptr); + REQUIRE (d4->get_params().get("vert_location")=="int"); + + auto d5 = create_diagnostic("geopotential_mid",grid); + REQUIRE (std::dynamic_pointer_cast(d5)!=nullptr); + REQUIRE (d5->get_params().get("vert_location")=="mid"); + auto d6 = create_diagnostic("geopotential_int",grid); + REQUIRE (std::dynamic_pointer_cast(d6)!=nullptr); + REQUIRE (d6->get_params().get("vert_location")=="int"); + + auto d7 = create_diagnostic("dz",grid); + REQUIRE (std::dynamic_pointer_cast(d7)!=nullptr); + REQUIRE (d7->get_params().get("vert_location")=="mid"); + } + +} + +} // namespace scream