Skip to content

Commit

Permalink
EAMxx: cleanup the VerticalLayer diagnostic
Browse files Browse the repository at this point in the history 
* Rename altitude with height
* Clean up includes
* Clean up unit tests
  • Loading branch information
@bartgol
bartgol committed May 20, 2024
1 parent 0da3b22 commit 2f4be63
Show file tree
Hide file tree
Showing 5 changed files with 215 additions and 166 deletions.
4 changes: 3 additions & 1 deletion components/eamxx/src/diagnostics/register_diagnostics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,9 +36,11 @@ inline void register_diagnostics () {
diag_factory.register_product("Exner",&create_atmosphere_diagnostic<ExnerDiagnostic>);
diag_factory.register_product("VirtualTemperature",&create_atmosphere_diagnostic<VirtualTemperatureDiagnostic>);
diag_factory.register_product("z_int",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("geopotential_int",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("z_mid",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("geopotential_int",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("geopotential_mid",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("height_int",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("height_mid",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("dz",&create_atmosphere_diagnostic<VerticalLayerDiagnostic>);
diag_factory.register_product("DryStaticEnergy",&create_atmosphere_diagnostic<DryStaticEnergyDiagnostic>);
diag_factory.register_product("SeaLevelPressure",&create_atmosphere_diagnostic<SeaLevelPressureDiagnostic>);
Expand Down
176 changes: 88 additions & 88 deletions components/eamxx/src/diagnostics/tests/vertical_layer_tests.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,8 @@
#include "catch2/catch.hpp"

#include "share/grid/mesh_free_grids_manager.hpp"
#include "diagnostics/vertical_layer.hpp"
#include "diagnostics/register_diagnostics.hpp"

#include "physics/share/physics_constants.hpp"

#include "share/util/scream_setup_random_test.hpp"
#include "share/util/scream_common_physics_functions.hpp"
#include "share/field/field_utils.hpp"

#include "ekat/ekat_pack.hpp"
#include "ekat/kokkos/ekat_kokkos_utils.hpp"
#include "ekat/util/ekat_test_utils.hpp"

#include <iomanip>
#include "share/grid/mesh_free_grids_manager.hpp"

namespace scream {

Expand All @@ -39,24 +27,13 @@ create_gm (const ekat::Comm& comm, const int ncols, const int nlevs) {
}

//-----------------------------------------------------------------------------------------------//
template<typename DeviceT, typename Engine>
void run( Engine& engine,
std::string diag_name,
const std::string& location)
template<typename DeviceT, int N>
void run (const std::string& diag_name, const std::string& location)
{
using PC = scream::physics::Constants<Real>;
using Pack = ekat::Pack<Real,SCREAM_PACK_SIZE>;
using KT = ekat::KokkosTypes<DeviceT>;
using ExecSpace = typename KT::ExeSpace;
using MemberType = typename KT::MemberType;
using view_1d = typename KT::template view_1d<Pack>;
using rview_1d = typename KT::template view_1d<Real>;
using view_2d = typename KT::template view_2d<Pack>;

const int packsize = SCREAM_PACK_SIZE;
using PC = scream::physics::Constants<Real>;

const int packsize = N;
constexpr int num_levs = packsize*2 + 1; // Number of levels to use for tests, make sure the last pack can also have some empty slots (packsize>1).
const int num_mid_packs = ekat::npack<Pack>(num_levs);
const int num_mid_packs_p1 = ekat::npack<Pack>(num_levs+1);

// A world comm
ekat::Comm comm(MPI_COMM_WORLD);
Expand All @@ -68,31 +45,19 @@ void run( Engine& engine,
// A time stamp
util::TimeStamp t0 ({2022,1,1},{0,0,0});

// Kokkos Policy
auto policy = ekat::ExeSpaceUtils<ExecSpace>::get_default_team_policy(ncols, num_mid_packs);

// Input (randomized) views
view_1d temperature("temperature",num_mid_packs),
pseudodensity("pseudo_density",num_mid_packs),
pressure("pressure",num_mid_packs),
watervapor("watervapor",num_mid_packs);

auto dview_as_real = [&] (const view_1d& v) -> rview_1d {
return rview_1d(reinterpret_cast<Real*>(v.data()),v.size()*packsize);
};

register_diagnostics();
auto& diag_factory = AtmosphereDiagnosticFactory::instance();

// Construct the Diagnostic
ekat::ParameterList params;
std::string name = diag_name;
if (location=="midpoints") {
diag_name += "_mid";
name += "_mid";
} else if (location=="interfaces") {
diag_name += "_int";
name += "_int";
}
params.set<std::string>("diag_name", diag_name);
auto diag = diag_factory.create(diag_name,comm,params);
params.set<std::string>("diag_name", name);
auto diag = diag_factory.create(name,comm,params);
diag->set_grids(gm);

const bool needs_phis = diag_name=="z" or diag_name=="geopotential";
Expand All @@ -107,21 +72,31 @@ void run( Engine& engine,
const auto name = f.name();
f.get_header().get_tracking().update_time_stamp(t0);
diag->set_required_field(f.get_const());
REQUIRE_THROWS(diag->set_computed_field(f));
input_fields.emplace(name,f);
}

// Can't set computed fields in the diag
REQUIRE_THROWS(diag->set_computed_field(input_fields.begin()->second));

// Note: we are not testing the calculate_dz utility. We are testing
// the diag class, so use some inputs that make checking results easier
// With these inputs, T_virt=T, and dz=8*rd/g
const Real dz = PC::RD/PC::gravit;
const Real phis = 3;
input_fields["T_mid"].deep_copy(Real(4));
input_fields["p_mid"].deep_copy(Real(2));
input_fields["pseudo_density"].deep_copy(Real(4));
input_fields["qv"].deep_copy(Real(0));
const Real g = PC::gravit;
const Real rho_val = 4;
const Real qv_val = 0;
const Real p_val = 2;
const Real T_val = 4;
const Real c1 = -PC::ONE + PC::ONE / PC::ep_2;
const Real Tvirt_val = T_val*(PC::ONE + c1*qv_val);
const Real dz_val = (PC::RD/g) * rho_val*Tvirt_val / p_val;
const Real phis_val = 3;

input_fields["T_mid"].deep_copy(T_val);
input_fields["p_mid"].deep_copy(p_val);
input_fields["pseudo_density"].deep_copy(rho_val);
input_fields["qv"].deep_copy(qv_val);
if (needs_phis) {
input_fields["phis"].deep_copy(phis);
input_fields["phis"].deep_copy(phis_val);
}

// Initialize and run the diagnostic
Expand All @@ -132,28 +107,43 @@ void run( Engine& engine,
auto d_h = diag_out.get_view<Real**,Host>();

// Compare against expecte value
const auto last_lev = location=="interfaces" ? num_levs : num_levs-1;
const auto last_int = num_levs;
const auto last_mid = last_int-1;

// Precompute surface value and increment depending on the diag type
Real delta, surf_val;
if (diag_name=="altitude") {
surf_val = 0;
delta = dz_val;
} else if (diag_name=="z") {
surf_val = phis_val/g;
delta = dz_val;
} else {
surf_val = phis_val;
delta = dz_val*g;
}

for (int icol=0; icol<ncols; ++icol) {
for (int ilev=last_lev; ilev>=0; --ilev) {
int num_mid_levs_below = ilev-last_lev;
Real tgt_mid, tgt_int;
if (diag_name=="dz") {
tgt_mid = dz;
} else if (diag_name=="altitude") {
tgt_int = phis/PC::gravit + num_mid_levs_below*dz;
tgt_mid = tgt_int + dz/2;
} else if (diag_name=="geopotential") {
tgt_int = phis + num_mid_levs_below*dz*PC::gravit;
tgt_mid = tgt_int + PC::gravit*dz/2;
} else {
tgt_int = num_mid_levs_below*dz;
tgt_mid = tgt_int + dz/2;
}
Real prev_int_val = surf_val;

if (location=="interfaces") {
REQUIRE (d_h(icol,ilev)==tgt_int);
if (location=="interfaces") {
// Check surface value
REQUIRE (d_h(icol,num_levs)==prev_int_val);
}

for (int ilev=last_mid; ilev>=0; --ilev) {
if (diag_name=="dz") {
REQUIRE (d_h(icol,ilev)==dz_val);
} else {
REQUIRE (d_h(icol,ilev)==tgt_mid);
// If interface, check value, otherwise perform int->mid averaging and check value
auto int_val = prev_int_val + delta;
if (location=="interfaces") {
REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(int_val,1e-5));
} else {
auto mid_val = (int_val + prev_int_val) / 2;
REQUIRE_THAT(d_h(icol,ilev), Catch::Matchers::WithinRel(mid_val,1e-5));
}
prev_int_val = int_val;
}
}
}
Expand All @@ -168,26 +158,36 @@ TEST_CASE("vertical_layer_test", "vertical_layer_test]"){
using scream::Real;
using Device = scream::DefaultDevice;

constexpr int num_runs = 5;

auto engine = scream::setup_random_test();

printf("Test specs\n");
printf(" - number of randomized runs: %d\n",num_runs);
printf(" - scalar type: Pack<Real,%d>\n\n", SCREAM_PACK_SIZE);

for (int irun=0; irun<num_runs; ++irun) {
ekat::Comm comm(MPI_COMM_WORLD);
auto root_print = [&](const std::string& msg) {
if (comm.am_i_root()) {
printf("%s",msg.c_str());
}
};
auto do_run = [&] (auto int_const) {
constexpr int N = decltype(int_const)::value;
root_print("\n");
root_print(" -> Testing diagnostic for pack_size=" + std::to_string(N) + "\n");
for (std::string loc : {"midpoints","interfaces"}) {
for (std::string diag : {"geopotential","altitude","z"}) {
printf(" -> Testing diag=%s at %s ...\n",diag.c_str(),loc.c_str());
run<Device>(engine, loc, diag);
printf(" -> Testing diag=%s at %s ... PASS!\n",diag.c_str(),loc.c_str());
std::string msg = " -> Testing diag=" + diag + " at " + loc + " ";
std::string dots (50-msg.size(),'.');
root_print (msg + dots + "\n");
run<Device,N>(diag, loc);
root_print (msg + dots + " PASS!\n");
}
}
printf(" -> Testing diag=dz ...\n");
run<Device>(engine, "", "dz");
printf(" -> Testing diag=dz ... PASS!\n");
std::string msg = " -> Testing diag=dz ";
std::string dots (50-msg.size(),'.');
root_print (msg + dots + "\n");
run<Device, N>("dz", "UNUSED");
root_print (msg + dots + " PASS!\n");
};

if (SCREAM_PACK_SIZE!=1) {
do_run(std::integral_constant<int,1>());
}
do_run(std::integral_constant<int,SCREAM_PACK_SIZE>());

} // TEST_CASE

Expand Down
Loading

0 comments on commit 2f4be63

Please sign in to comment.